1 /* mpfr_subnormalize -- Subnormalize a floating point number
2 emulating sub-normal numbers.
4 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel projects, INRIA.
7 This file is part of the GNU MPFR Library.
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.
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.
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 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 #include "mpfr-impl.h"
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
31 1 | 0 | 0 | Trunc if even |
32 1 | 0 | -1 | AddOneUlp |
33 1 | 1 | ? | AddOneUlp |
35 For other rounding mode, there isn't such a problem.
36 Just round it again and merge the ternary values.
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
47 mpfr_subnormalize (mpfr_ptr y
, int old_inexact
, mpfr_rnd_t rnd
)
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
);
57 mpfr_set_underflow ();
60 /* We have to emulate one bit rounding if EXP(y) = emin */
61 if (MPFR_GET_EXP (y
) == __gmpfr_emin
)
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
);
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. */
75 mp_limb_t
*mant
, rb
,sb
;
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);
84 sb
= *mant
& ((MPFR_LIMB_HIGHBIT
>> 1) - 1);
85 while (sb
== 0 && s
-- != 0)
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))
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. */
101 else if (MPFR_IS_LIKE_RNDZ (rnd
, MPFR_IS_NEG (y
)))
104 mpfr_setmin (y
, __gmpfr_emin
);
110 /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */
111 mpfr_setmin (y
, __gmpfr_emin
+ 1);
115 else /* Hard case: It is more or less the same problem than mpfr_cache */
121 MPFR_ASSERTD (MPFR_GET_EXP (y
) > __gmpfr_emin
);
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
));
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))
137 if (MPFR_UNLIKELY (rnd
== MPFR_RNDN
&&
138 (inexact
== MPFR_EVEN_INEX
||
139 inexact
== -MPFR_EVEN_INEX
)))
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
))
145 if (SAME_SIGN (inexact
, MPFR_INT_SIGN (y
)))
146 mpfr_nexttozero (dest
);
148 mpfr_nexttoinf (dest
);
152 else if (MPFR_UNLIKELY (inexact
== 0))
153 inexact
= old_inexact
;
156 inex2
= mpfr_set (y
, dest
, rnd
);
157 MPFR_ASSERTN (inex2
== 0);
158 MPFR_ASSERTN (MPFR_IS_PURE_FP (y
));