1 /* mpfr_get_ld, mpfr_get_ld_2exp -- convert a multiple precision floating-point
2 number to a machine long double
3
4 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel 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 http://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 <float.h>
25
26 #include "mpfr-impl.h"
27
28 #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE
29
30 long double
mpfr_get_ld(mpfr_srcptr x,mpfr_rnd_t rnd_mode)31 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
32 {
33
34 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
35 return (long double) mpfr_get_d (x, rnd_mode);
36 else /* now x is a normal non-zero number */
37 {
38 long double r; /* result */
39 long double m;
40 double s; /* part of result */
41 mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */
42 mpfr_t y, z;
43 int sign;
44
45 /* first round x to the target long double precision, so that
46 all subsequent operations are exact (this avoids double rounding
47 problems) */
48 mpfr_init2 (y, MPFR_LDBL_MANT_DIG);
49 mpfr_init2 (z, MPFR_LDBL_MANT_DIG);
50 /* Note about the precision of z: even though IEEE_DBL_MANT_DIG is
51 sufficient, z has been set to the same precision as y so that
52 the mpfr_sub below calls mpfr_sub1sp, which is faster than the
53 generic subtraction, even in this particular case (from tests
54 done by Patrick Pelissier on a 64-bit Core2 Duo against r7285).
55 But here there is an important cancellation in the subtraction.
56 TODO: get more information about what has been tested. */
57
58 mpfr_set (y, x, rnd_mode);
59 sh = MPFR_GET_EXP (y);
60 sign = MPFR_SIGN (y);
61 MPFR_SET_EXP (y, 0);
62 MPFR_SET_POS (y);
63
64 r = 0.0;
65 do {
66 s = mpfr_get_d (y, MPFR_RNDN); /* high part of y */
67 r += (long double) s;
68 mpfr_set_d (z, s, MPFR_RNDN); /* exact */
69 mpfr_sub (y, y, z, MPFR_RNDN); /* exact */
70 } while (!MPFR_IS_ZERO (y));
71
72 mpfr_clear (z);
73 mpfr_clear (y);
74
75 /* we now have to multiply back by 2^sh */
76 MPFR_ASSERTD (r > 0);
77 if (sh != 0)
78 {
79 /* An overflow may occurs (example: 0.5*2^1024) */
80 while (r < 1.0)
81 {
82 r += r;
83 sh--;
84 }
85
86 if (sh > 0)
87 m = 2.0;
88 else
89 {
90 m = 0.5;
91 sh = -sh;
92 }
93
94 for (;;)
95 {
96 if (sh % 2)
97 r = r * m;
98 sh >>= 1;
99 if (sh == 0)
100 break;
101 m = m * m;
102 }
103 }
104 if (sign < 0)
105 r = -r;
106 return r;
107 }
108 }
109
110 #else
111
112 long double
mpfr_get_ld(mpfr_srcptr x,mpfr_rnd_t rnd_mode)113 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
114 {
115 mpfr_long_double_t ld;
116 mpfr_t tmp;
117 int inex;
118 MPFR_SAVE_EXPO_DECL (expo);
119
120 MPFR_SAVE_EXPO_MARK (expo);
121
122 mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
123 inex = mpfr_set (tmp, x, rnd_mode);
124
125 mpfr_set_emin (-16382-63);
126 mpfr_set_emax (16384);
127 mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode);
128 mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */
129 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
130 ld.ld = (long double) mpfr_get_d (tmp, rnd_mode);
131 else
132 {
133 mp_limb_t *tmpmant;
134 mpfr_exp_t e, denorm;
135
136 tmpmant = MPFR_MANT (tmp);
137 e = MPFR_GET_EXP (tmp);
138 /* the smallest normal number is 2^(-16382), which is 0.5*2^(-16381)
139 in MPFR, thus any exponent <= -16382 corresponds to a subnormal
140 number */
141 denorm = MPFR_UNLIKELY (e <= -16382) ? - e - 16382 + 1 : 0;
142 #if GMP_NUMB_BITS >= 64
143 ld.s.manl = (tmpmant[0] >> denorm);
144 ld.s.manh = (tmpmant[0] >> denorm) >> 32;
145 #elif GMP_NUMB_BITS == 32
146 if (MPFR_LIKELY (denorm == 0))
147 {
148 ld.s.manl = tmpmant[0];
149 ld.s.manh = tmpmant[1];
150 }
151 else if (denorm < 32)
152 {
153 ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
154 ld.s.manh = tmpmant[1] >> denorm;
155 }
156 else /* 32 <= denorm <= 64 */
157 {
158 ld.s.manl = tmpmant[1] >> (denorm - 32);
159 ld.s.manh = 0;
160 }
161 #else
162 # error "GMP_NUMB_BITS must be 32 or >= 64"
163 /* Other values have never been supported anyway. */
164 #endif
165 if (MPFR_LIKELY (denorm == 0))
166 {
167 ld.s.exph = (e + 0x3FFE) >> 8;
168 ld.s.expl = (e + 0x3FFE);
169 }
170 else
171 ld.s.exph = ld.s.expl = 0;
172 ld.s.sign = MPFR_IS_NEG (x);
173 }
174
175 mpfr_clear (tmp);
176 MPFR_SAVE_EXPO_FREE (expo);
177 return ld.ld;
178 }
179
180 #endif
181
182 /* contributed by Damien Stehle */
183 long double
mpfr_get_ld_2exp(long * expptr,mpfr_srcptr src,mpfr_rnd_t rnd_mode)184 mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
185 {
186 long double ret;
187 mpfr_exp_t exp;
188 mpfr_t tmp;
189
190 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
191 return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode);
192
193 tmp[0] = *src; /* Hack copy mpfr_t */
194 MPFR_SET_EXP (tmp, 0);
195 ret = mpfr_get_ld (tmp, rnd_mode);
196
197 if (MPFR_IS_PURE_FP(src))
198 {
199 exp = MPFR_GET_EXP (src);
200
201 /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
202 if (ret == 1.0)
203 {
204 ret = 0.5;
205 exp ++;
206 }
207 else if (ret == -1.0)
208 {
209 ret = -0.5;
210 exp ++;
211 }
212
213 MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
214 || (ret <= -0.5 && ret > -1.0));
215 MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
216 }
217 else
218 exp = 0;
219
220 *expptr = exp;
221 return ret;
222 }
223