1/* mpfr_subnormalize -- Subnormalize a floating point number
2   emulating sub-normal numbers.
3
4Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5Contributed by the AriC and Caramel projects, INRIA.
6
7This file is part of the GNU MPFR Library.
8
9The GNU MPFR Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MPFR Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2251 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23
24#include "mpfr-impl.h"
25
26/* For MPFR_RNDN, we can have a problem of double rounding.
27   In such a case, this table helps to conclude what to do (y positive):
28     Rounding Bit |  Sticky Bit | inexact  | Action    | new inexact
29     0            |   ?         |  ?       | Trunc     | sticky
30     1            |   0         |  1       | Trunc     |
31     1            |   0         |  0       | Trunc if even |
32     1            |   0         | -1       | AddOneUlp |
33     1            |   1         |  ?       | AddOneUlp |
34
35   For other rounding mode, there isn't such a problem.
36   Just round it again and merge the ternary values.
37
38   Set the inexact flag if the returned ternary value is non-zero.
39   Set the underflow flag if a second rounding occurred (whether this
40   rounding is exact or not). See
41     https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00000.html
42     https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00008.html
43     https://sympa.inria.fr/sympa/arc/mpfr/2009-06/msg00010.html
44*/
45
46int
47mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
48{
49  int sign;
50
51  /* The subnormal exponent range is [ emin, emin + MPFR_PREC(y) - 2 ] */
52  if (MPFR_LIKELY (MPFR_IS_SINGULAR (y)
53                   || (MPFR_GET_EXP (y) >=
54                       __gmpfr_emin + (mpfr_exp_t) MPFR_PREC (y) - 1)))
55    MPFR_RET (old_inexact);
56
57  mpfr_set_underflow ();
58  sign = MPFR_SIGN (y);
59
60  /* We have to emulate one bit rounding if EXP(y) = emin */
61  if (MPFR_GET_EXP (y) == __gmpfr_emin)
62    {
63      /* If this is a power of 2, we don't need rounding.
64         It handles cases when |y| = 0.1 * 2^emin */
65      if (mpfr_powerof2_raw (y))
66        MPFR_RET (old_inexact);
67
68      /* We keep the same sign for y.
69         Assuming Y is the real value and y the approximation
70         and since y is not a power of 2:  0.5*2^emin < Y < 1*2^emin
71         We also know the direction of the error thanks to ternary value. */
72
73      if (rnd == MPFR_RNDN)
74        {
75          mp_limb_t *mant, rb ,sb;
76          mp_size_t s;
77          /* We need the rounding bit and the sticky bit. Read them
78             and use the previous table to conclude. */
79          s = MPFR_LIMB_SIZE (y) - 1;
80          mant = MPFR_MANT (y) + s;
81          rb = *mant & (MPFR_LIMB_HIGHBIT >> 1);
82          if (rb == 0)
83            goto set_min;
84          sb = *mant & ((MPFR_LIMB_HIGHBIT >> 1) - 1);
85          while (sb == 0 && s-- != 0)
86            sb = *--mant;
87          if (sb != 0)
88            goto set_min_p1;
89          /* Rounding bit is 1 and sticky bit is 0.
90             We need to examine old inexact flag to conclude. */
91          if ((old_inexact > 0 && sign > 0) ||
92              (old_inexact < 0 && sign < 0))
93            goto set_min;
94          /* If inexact != 0, return 0.1*2^(emin+1).
95             Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0
96             So we have 0.1100000000000000000000000*2^emin exactly.
97             We return 0.1*2^(emin+1) according to the even-rounding
98             rule on subnormals. */
99          goto set_min_p1;
100        }
101      else if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y)))
102        {
103        set_min:
104          mpfr_setmin (y, __gmpfr_emin);
105          MPFR_RET (-sign);
106        }
107      else
108        {
109        set_min_p1:
110          /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */
111          mpfr_setmin (y, __gmpfr_emin + 1);
112          MPFR_RET (sign);
113        }
114    }
115  else /* Hard case: It is more or less the same problem than mpfr_cache */
116    {
117      mpfr_t dest;
118      mpfr_prec_t q;
119      int inexact, inex2;
120
121      MPFR_ASSERTD (MPFR_GET_EXP (y) > __gmpfr_emin);
122
123      /* Compute the intermediary precision */
124      q = (mpfr_uexp_t) MPFR_GET_EXP (y) - __gmpfr_emin + 1;
125      MPFR_ASSERTD (q >= MPFR_PREC_MIN && q < MPFR_PREC (y));
126
127      /* TODO: perform the rounding in place. */
128      mpfr_init2 (dest, q);
129      /* Round y in dest */
130      MPFR_SET_EXP (dest, MPFR_GET_EXP (y));
131      MPFR_SET_SIGN (dest, sign);
132      MPFR_RNDRAW_EVEN (inexact, dest,
133                        MPFR_MANT (y), MPFR_PREC (y), rnd, sign,
134                        MPFR_SET_EXP (dest, MPFR_GET_EXP (dest) + 1));
135      if (MPFR_LIKELY (old_inexact != 0))
136        {
137          if (MPFR_UNLIKELY (rnd == MPFR_RNDN &&
138                             (inexact == MPFR_EVEN_INEX ||
139                              inexact == -MPFR_EVEN_INEX)))
140            {
141              /* if both roundings are in the same direction, we have to go
142                 back in the other direction */
143              if (SAME_SIGN (inexact, old_inexact))
144                {
145                  if (SAME_SIGN (inexact, MPFR_INT_SIGN (y)))
146                    mpfr_nexttozero (dest);
147                  else
148                    mpfr_nexttoinf (dest);
149                  inexact = -inexact;
150                }
151            }
152          else if (MPFR_UNLIKELY (inexact == 0))
153            inexact = old_inexact;
154        }
155
156      inex2 = mpfr_set (y, dest, rnd);
157      MPFR_ASSERTN (inex2 == 0);
158      MPFR_ASSERTN (MPFR_IS_PURE_FP (y));
159      mpfr_clear (dest);
160
161      MPFR_RET (inexact);
162    }
163}
164