1 /* mpfr_subnormalize -- Subnormalize a floating point number
2 emulating sub-normal numbers.
4 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
5 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to
21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
22 MA 02110-1301, USA. */
24 #include "mpfr-impl.h"
26 /* For GMP_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 inexact flags.
40 mpfr_subnormalize (mpfr_ptr y
, int old_inexact
, mp_rnd_t rnd
)
44 /* The subnormal exponent range are from:
45 mpfr_emin to mpfr_emin + MPFR_PREC(y) - 1 */
46 if (MPFR_LIKELY (MPFR_IS_SINGULAR (y
)
47 || (MPFR_GET_EXP (y
) >=
48 __gmpfr_emin
+ (mp_exp_t
) MPFR_PREC (y
) - 1)))
49 inexact
= old_inexact
;
51 /* We have to emulate one bit rounding if EXP(y) = emin */
52 else if (MPFR_GET_EXP (y
) == __gmpfr_emin
)
54 /* If this is a power of 2, we don't need rounding.
55 It handles cases when rouding away and y=0.1*2^emin */
56 if (mpfr_powerof2_raw (y
))
57 inexact
= old_inexact
;
58 /* We keep the same sign for y.
59 Assuming Y is the real value and y the approximation
60 and since y is not a power of 2: 0.5*2^emin < Y < 1*2^emin
61 We also know the direction of the error thanks to inexact flag */
62 else if (rnd
== GMP_RNDN
)
64 mp_limb_t
*mant
, rb
,sb
;
66 /* We need the rounding bit and the sticky bit. Read them
67 and use the previous table to conclude. */
68 s
= MPFR_LIMB_SIZE (y
) - 1;
69 mant
= MPFR_MANT (y
) + s
;
70 rb
= *mant
& (MPFR_LIMB_HIGHBIT
>> 1);
73 sb
= *mant
& ((MPFR_LIMB_HIGHBIT
>> 1) - 1);
74 while (sb
== 0 && s
-- != 0)
78 /* Rounding bit is 1 and sticky bit is 0.
79 We need to examine old inexact flag to conclude. */
80 if ((old_inexact
> 0 && MPFR_SIGN (y
) > 0) ||
81 (old_inexact
< 0 && MPFR_SIGN (y
) < 0))
83 /* If inexact != 0, return 0.1*2^(emin+1).
84 Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0
85 So we have 0.1100000000000000000000000*2^emin exactly.
86 We return 0.1*2^(emin+1) according to the even-rounding
87 rule on subnormals. */
90 else if (MPFR_IS_LIKE_RNDZ (rnd
, MPFR_IS_NEG (y
)))
93 mpfr_setmin (y
, __gmpfr_emin
);
94 inexact
= -MPFR_SIGN (y
);
99 /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */
100 mpfr_setmin (y
, __gmpfr_emin
+ 1);
101 inexact
= MPFR_SIGN (y
);
105 else /* Hard case: It is more or less the same problem than mpfr_cache */
111 /* Compute the intermediary precision */
112 q
= (mpfr_uexp_t
) MPFR_GET_EXP (y
) - __gmpfr_emin
+ 1;
113 mpfr_init2 (dest
, q
);
114 /* Round y in dest */
115 sign
= MPFR_SIGN (y
);
116 MPFR_SET_EXP (dest
, MPFR_GET_EXP (y
));
117 MPFR_SET_SIGN (dest
, sign
);
118 MPFR_RNDRAW_EVEN (inexact
, dest
,
119 MPFR_MANT (y
), MPFR_PREC (y
), rnd
, sign
,
120 MPFR_SET_EXP (dest
, MPFR_GET_EXP (dest
)+1));
121 if (MPFR_LIKELY (old_inexact
!= 0))
123 if (MPFR_UNLIKELY(rnd
== GMP_RNDN
&& (inexact
== MPFR_EVEN_INEX
124 || inexact
== -MPFR_EVEN_INEX
)))
126 /* if both roundings are in the same direction, we have to go
127 back in the other direction */
128 if (SAME_SIGN (inexact
, old_inexact
))
130 if (SAME_SIGN (inexact
, MPFR_INT_SIGN (y
)))
131 mpfr_nexttozero (dest
);
133 mpfr_nexttoinf (dest
);
137 else if (MPFR_UNLIKELY (inexact
== 0))
138 inexact
= old_inexact
;
140 old_inexact
= mpfr_set (y
, dest
, rnd
);
141 MPFR_ASSERTN (old_inexact
== 0);
142 MPFR_ASSERTN (MPFR_IS_PURE_FP (y
));