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