1 /* mpfr_subnormalize -- Subnormalize a floating point number
2 emulating sub-normal numbers.
3
4 Copyright 2005-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba projects, INRIA.
6
7 This file is part of the GNU MPFR Library.
8
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 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 modes, 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
46 int
47 mpfr_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 || rnd == MPFR_RNDNA)
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. Note the same holds for RNDNA. */
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 as mpfr_cache */
116 {
117 mpfr_t dest;
118 mpfr_prec_t q;
119 mpfr_rnd_t rnd2;
120 int inexact, inex2;
121
122 MPFR_ASSERTD (MPFR_GET_EXP (y) > __gmpfr_emin);
123
124 /* Compute the intermediary precision */
125 q = (mpfr_uexp_t) MPFR_GET_EXP (y) - __gmpfr_emin + 1;
126 MPFR_ASSERTD (q >= MPFR_PREC_MIN && q < MPFR_PREC (y));
127
128 /* TODO: perform the rounding in place. */
129 mpfr_init2 (dest, q);
130 /* Round y in dest */
131 MPFR_SET_EXP (dest, MPFR_GET_EXP (y));
132 MPFR_SET_SIGN (dest, sign);
133 rnd2 = rnd == MPFR_RNDNA ? MPFR_RNDN : rnd;
134 MPFR_RNDRAW_EVEN (inexact, dest,
135 MPFR_MANT (y), MPFR_PREC (y), rnd2, sign,
136 MPFR_SET_EXP (dest, MPFR_GET_EXP (dest) + 1));
137 if (MPFR_LIKELY (old_inexact != 0))
138 {
139 if (MPFR_UNLIKELY (rnd2 == MPFR_RNDN &&
140 (inexact == MPFR_EVEN_INEX ||
141 inexact == -MPFR_EVEN_INEX)))
142 {
143 /* If both roundings are in the same direction,
144 we have to go back in the other direction.
145 For MPFR_RNDNA it is the same, since we are not
146 exactly in the middle case (old_inexact != 0). */
147 if (SAME_SIGN (inexact, old_inexact))
148 {
149 if (SAME_SIGN (inexact, MPFR_INT_SIGN (y)))
150 mpfr_nexttozero (dest);
151 else /* subnormal range, thus no overflow */
152 {
153 mpfr_nexttoinf (dest);
154 MPFR_ASSERTD(!MPFR_IS_INF (dest));
155 }
156 inexact = -inexact;
157 }
158 }
159 else if (MPFR_UNLIKELY (inexact == 0))
160 inexact = old_inexact;
161 }
162 else if (rnd == MPFR_RNDNA &&
163 (inexact == MPFR_EVEN_INEX || inexact == -MPFR_EVEN_INEX))
164 {
165 /* We are in the middle case: since we used RNDN to round, we should
166 round in the opposite direction when inexact has the opposite
167 sign of y. */
168 if (!SAME_SIGN (inexact, MPFR_INT_SIGN (y)))
169 {
170 mpfr_nexttoinf (dest);
171 MPFR_ASSERTD(!MPFR_IS_INF (dest));
172 inexact = -inexact;
173 }
174 }
175
176 inex2 = mpfr_set (y, dest, rnd);
177 MPFR_ASSERTN (inex2 == 0);
178 MPFR_ASSERTN (MPFR_IS_PURE_FP (y));
179 mpfr_clear (dest);
180
181 MPFR_RET (inexact);
182 }
183 }