1 /* mpn_get_d -- limbs to double conversion.
2
3 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
4 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
5 FUTURE GNU MP RELEASES.
6
7 Copyright 2003, 2004 Free Software Foundation, Inc.
8
9 This file is part of the GNU MP Library.
10
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 2.1 of the License, or (at your
14 option) any later version.
15
16 The GNU MP Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 License for more details.
20
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
23 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
24 MA 02110-1301, USA. */
25
26 #include "mpir.h"
27 #include "gmp-impl.h"
28 #include "longlong.h"
29
30
31 #define CONST_1024 (1024)
32 #define CONST_NEG_1023 (-1023)
33 #define CONST_NEG_1022_SUB_53 (-1022 - 53)
34
35 /* Return the value {ptr,size}*2^exp, and negative if sign<0.
36 Must have size>=1, and a non-zero high limb ptr[size-1].
37
38 {ptr,size} is truncated towards zero. This is consistent with other gmp
39 conversions, like mpz_set_f or mpz_set_q, and is easy to implement and
40 test.
41
42 In the past conversions had attempted (imperfectly) to let the hardware
43 float rounding mode take effect, but that gets tricky since multiple
44 roundings need to be avoided, or taken into account, and denorms mean the
45 effective precision of the mantissa is not constant. (For reference,
46 mpz_get_d on IEEE systems was ok, except it operated on the absolute
47 value. mpf_get_d and mpq_get_d suffered from multiple roundings and from
48 not always using enough bits to get the rounding right.)
49
50 It's felt that GMP is not primarily concerned with hardware floats, and
51 really isn't enhanced by getting involved with hardware rounding modes
52 (which could even be some weird unknown style), so something unambiguous
53 and straightforward is best.
54
55
56 The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
57 limb and is done with shifts and masks. The 64-bit case in particular
58 should come out nice and compact.
59
60 The generic code works one bit at a time, which will be quite slow, but
61 should support any binary-based "double" and be safe against any rounding
62 mode. Note in particular it works on IEEE systems too.
63
64
65 Traps:
66
67 Hardware traps for overflow to infinity, underflow to zero, or
68 unsupported denorms may or may not be taken. The IEEE code works bitwise
69 and so probably won't trigger them, the generic code works by float
70 operations and so probably will. This difference might be thought less
71 than ideal, but again its felt straightforward code is better than trying
72 to get intimate with hardware exceptions (of perhaps unknown nature).
73
74
75 Not done:
76
77 mpz_get_d in the past handled size==1 with a cast limb->double. This
78 might still be worthwhile there (for up to the mantissa many bits), but
79 for mpn_get_d here, the cost of applying "exp" to the resulting exponent
80 would probably use up any benefit a cast may have over bit twiddling.
81 Also, if the exponent is pushed into denorm range then bit twiddling is
82 the only option, to ensure the desired truncation is obtained.
83
84
85 Other:
86
87 For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
88 to the kernel for values >= 2^63. This makes it slow, and worse the
89 Linux kernel (what versions?) apparently uses untested code in its trap
90 handling routines, and gets the sign wrong. We don't use such a limb to
91 double cast, neither in the IEEE or generic code. */
92
93
94 double
mpn_get_d(mp_srcptr ptr,mp_size_t size,mp_size_t sign,long exp)95 mpn_get_d (mp_srcptr ptr, mp_size_t size, mp_size_t sign, long exp)
96 {
97 ASSERT (size >= 0);
98 ASSERT_MPN (ptr, size);
99 ASSERT (size == 0 || ptr[size-1] != 0);
100
101 if (size == 0)
102 return 0.0;
103
104 /* Adjust exp to a radix point just above {ptr,size}, guarding against
105 overflow. After this exp can of course be reduced to anywhere within
106 the {ptr,size} region without underflow. */
107 if (UNLIKELY ((mpir_ui) (GMP_NUMB_BITS * size)
108 > (mpir_ui) (LONG_MAX - exp)))
109 {
110 goto ieee_infinity;
111
112 /* generic */
113 exp = LONG_MAX;
114 }
115 else
116 {
117 exp += GMP_NUMB_BITS * size;
118 }
119
120 #define ONE_LIMB (GMP_LIMB_BITS == 64 && 2*GMP_NUMB_BITS >= 53)
121 #define TWO_LIMBS (GMP_LIMB_BITS == 32 && 3*GMP_NUMB_BITS >= 53)
122
123 if (ONE_LIMB || TWO_LIMBS)
124 {
125 union ieee_double_extract u;
126 mp_limb_t m0, m1, m2, rmask;
127 int lshift, rshift;
128
129 m0 = ptr[size-1]; /* high limb */
130 m1 = (size >= 2 ? ptr[size-2] : 0); /* second highest limb */
131 count_leading_zeros (lshift, m0);
132
133 /* relative to just under high non-zero bit */
134 exp -= (lshift - GMP_NAIL_BITS) + 1;
135
136 if (ONE_LIMB)
137 {
138 /* lshift to have high of m0 non-zero, and collapse nails */
139 rshift = GMP_LIMB_BITS - lshift;
140 m1 <<= GMP_NAIL_BITS;
141 rmask = GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX;
142 m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
143
144 /* rshift back to have bit 53 of m0 the high non-zero */
145 m0 >>= 11;
146 }
147 else /* TWO_LIMBS */
148 {
149 m2 = (size >= 3 ? ptr[size-3] : 0); /* third highest limb */
150
151 /* collapse nails from m1 and m2 */
152 #if GMP_NAIL_BITS != 0
153 m1 = (m1 << GMP_NAIL_BITS) | (m2 >> (GMP_NUMB_BITS-GMP_NAIL_BITS));
154 m2 <<= 2*GMP_NAIL_BITS;
155 #endif
156
157 /* lshift to have high of m0:m1 non-zero, collapse nails from m0 */
158 rshift = GMP_LIMB_BITS - lshift;
159 rmask = (GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX);
160 m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
161 m1 = (m1 << lshift) | ((m2 >> rshift) & rmask);
162
163 /* rshift back to have bit 53 of m0:m1 the high non-zero */
164 m1 = (m1 >> 11) | (m0 << (GMP_LIMB_BITS-11));
165 m0 >>= 11;
166 }
167
168 if (UNLIKELY (exp >= CONST_1024))
169 {
170 /* overflow, return infinity */
171 ieee_infinity:
172 m0 = 0;
173 m1 = 0;
174 exp = 1024;
175 }
176 else if (UNLIKELY (exp <= CONST_NEG_1023))
177 {
178 if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
179 return 0.0; /* denorm underflows to zero */
180
181 rshift = -1022 - exp;
182 ASSERT (rshift > 0 && rshift < 53);
183 if (ONE_LIMB)
184 {
185 m0 >>= rshift;
186 }
187 else /* TWO_LIMBS */
188 {
189 if (rshift >= 32)
190 {
191 m1 = m0;
192 m0 = 0;
193 rshift -= 32;
194 }
195 lshift = GMP_LIMB_BITS - rshift;
196 m1 = (m1 >> rshift) | (rshift == 0 ? 0 : m0 << lshift);
197 m0 >>= rshift;
198 }
199 exp = -1023;
200 }
201
202 if (ONE_LIMB)
203 {
204 #if GMP_LIMB_BITS > 32 /* avoid compiler warning about big shift */
205 u.s.manh = m0 >> 32;
206 #endif
207 u.s.manl = m0;
208 }
209 else /* TWO_LIMBS */
210 {
211 u.s.manh = m0;
212 u.s.manl = m1;
213 }
214
215 u.s.exp = exp + 1023;
216 u.s.sig = (sign < 0);
217 return u.d;
218 }
219 else
220 {
221 /* Non-IEEE or strange limb size, do something generic. */
222
223 mp_size_t i;
224 mp_limb_t limb, bit;
225 int shift;
226 double base, factor, prev_factor, d, new_d, diff;
227
228 /* "limb" is "ptr[i]" the limb being examined, "bit" is a mask for the
229 bit being examined, initially the highest non-zero bit. */
230 i = size-1;
231 limb = ptr[i];
232 count_leading_zeros (shift, limb);
233 bit = GMP_LIMB_HIGHBIT >> shift;
234
235 /* relative to just under high non-zero bit */
236 exp -= (shift - GMP_NAIL_BITS) + 1;
237
238 /* Power up "factor" to 2^exp, being the value of the "bit" in "limb"
239 being examined. */
240 base = (exp >= 0 ? 2.0 : 0.5);
241 exp = ABS (exp);
242 factor = 1.0;
243 for (;;)
244 {
245 if (exp & 1)
246 {
247 prev_factor = factor;
248 factor *= base;
249 if (factor == 0.0)
250 return 0.0; /* underflow */
251 if (factor == prev_factor)
252 {
253 d = factor; /* overflow, apparent infinity */
254 goto generic_done;
255 }
256 }
257 exp >>= 1;
258 if (exp == 0)
259 break;
260 base *= base;
261 }
262
263 /* Add a "factor" for each non-zero bit, working from high to low.
264 Stop if any rounding occurs, hence implementing a truncation.
265
266 Note no attention is paid to DBL_MANT_DIG, since the effective
267 number of bits in the mantissa isn't constant when in denorm range.
268 We also encountered an ARM system with apparently somewhat doubtful
269 software floats where DBL_MANT_DIG claimed 53 bits but only 32
270 actually worked. */
271
272 d = factor; /* high bit */
273 for (;;)
274 {
275 factor *= 0.5; /* next bit */
276 bit >>= 1;
277 if (bit == 0)
278 {
279 /* next limb, if any */
280 i--;
281 if (i < 0)
282 break;
283 limb = ptr[i];
284 bit = GMP_NUMB_HIGHBIT;
285 }
286
287 if (bit & limb)
288 {
289 new_d = d + factor;
290 diff = new_d - d;
291 if (diff != factor)
292 break; /* rounding occured, stop now */
293 d = new_d;
294 }
295 }
296
297 generic_done:
298 return (sign >= 0 ? d : -d);
299 }
300 }
301