xref: /netbsd-src/external/lgpl3/gmp/dist/extract-dbl.c (revision fa28c6faa16e0b00edee7acdcaf4899797043def)
1 /* __gmp_extract_double -- convert from double to array of mp_limb_t.
2 
3 Copyright 1996, 1999, 2000, 2001, 2002, 2006, 2012 Free Software Foundation,
4 Inc.
5 
6 This file is part of the GNU MP Library.
7 
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 The GNU MP Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16 License for more details.
17 
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
20 
21 #include "gmp.h"
22 #include "gmp-impl.h"
23 
24 #ifdef XDEBUG
25 #undef _GMP_IEEE_FLOATS
26 #endif
27 
28 #ifndef _GMP_IEEE_FLOATS
29 #define _GMP_IEEE_FLOATS 0
30 #endif
31 
32 /* Extract a non-negative double in d.  */
33 
34 int
35 __gmp_extract_double (mp_ptr rp, double d)
36 {
37   long exp;
38   unsigned sc;
39 #ifdef _LONG_LONG_LIMB
40 #define BITS_PER_PART 64	/* somewhat bogus */
41   unsigned long long int manl;
42 #else
43 #define BITS_PER_PART GMP_LIMB_BITS
44   unsigned long int manh, manl;
45 #endif
46 
47   /* BUGS
48 
49      1. Should handle Inf and NaN in IEEE specific code.
50      2. Handle Inf and NaN also in default code, to avoid hangs.
51      3. Generalize to handle all GMP_LIMB_BITS >= 32.
52      4. This lits is incomplete and misspelled.
53    */
54 
55   ASSERT (d >= 0.0);
56 
57   if (d == 0.0)
58     {
59       MPN_ZERO (rp, LIMBS_PER_DOUBLE);
60       return 0;
61     }
62 
63 #if _GMP_IEEE_FLOATS
64   {
65 #if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8
66     /* Work around alpha-specific bug in GCC 2.8.x.  */
67     volatile
68 #endif
69     union ieee_double_extract x;
70     x.d = d;
71     exp = x.s.exp;
72 #if BITS_PER_PART == 64		/* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */
73     manl = (((mp_limb_t) 1 << 63)
74 	    | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
75     if (exp == 0)
76       {
77 	/* Denormalized number.  Don't try to be clever about this,
78 	   since it is not an important case to make fast.  */
79 	exp = 1;
80 	do
81 	  {
82 	    manl = manl << 1;
83 	    exp--;
84 	  }
85 	while ((manl & GMP_LIMB_HIGHBIT) == 0);
86       }
87 #endif
88 #if BITS_PER_PART == 32
89     manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
90     manl = x.s.manl << 11;
91     if (exp == 0)
92       {
93 	/* Denormalized number.  Don't try to be clever about this,
94 	   since it is not an important case to make fast.  */
95 	exp = 1;
96 	do
97 	  {
98 	    manh = (manh << 1) | (manl >> 31);
99 	    manl = manl << 1;
100 	    exp--;
101 	  }
102 	while ((manh & GMP_LIMB_HIGHBIT) == 0);
103       }
104 #endif
105 #if BITS_PER_PART != 32 && BITS_PER_PART != 64
106   You need to generalize the code above to handle this.
107 #endif
108     exp -= 1022;		/* Remove IEEE bias.  */
109   }
110 #else
111   {
112     /* Unknown (or known to be non-IEEE) double format.  */
113     exp = 0;
114     if (d >= 1.0)
115       {
116 	ASSERT_ALWAYS (d * 0.5 != d);
117 
118 	while (d >= 32768.0)
119 	  {
120 	    d *= (1.0 / 65536.0);
121 	    exp += 16;
122 	  }
123 	while (d >= 1.0)
124 	  {
125 	    d *= 0.5;
126 	    exp += 1;
127 	  }
128       }
129     else if (d < 0.5)
130       {
131 	while (d < (1.0 / 65536.0))
132 	  {
133 	    d *=  65536.0;
134 	    exp -= 16;
135 	  }
136 	while (d < 0.5)
137 	  {
138 	    d *= 2.0;
139 	    exp -= 1;
140 	  }
141       }
142 
143     d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
144 #if BITS_PER_PART == 64
145     manl = d;
146 #endif
147 #if BITS_PER_PART == 32
148     manh = d;
149     manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
150 #endif
151   }
152 #endif /* IEEE */
153 
154   sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS;
155 
156   /* We add something here to get rounding right.  */
157   exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1;
158 
159 #if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 2
160 #if GMP_NAIL_BITS == 0
161   if (sc != 0)
162     {
163       rp[1] = manl >> (GMP_LIMB_BITS - sc);
164       rp[0] = manl << sc;
165     }
166   else
167     {
168       rp[1] = manl;
169       rp[0] = 0;
170       exp--;
171     }
172 #else
173   if (sc > GMP_NAIL_BITS)
174     {
175       rp[1] = manl >> (GMP_LIMB_BITS - sc);
176       rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK;
177     }
178   else
179     {
180       if (sc == 0)
181 	{
182 	  rp[1] = manl >> GMP_NAIL_BITS;
183 	  rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
184 	  exp--;
185 	}
186       else
187 	{
188 	  rp[1] = manl >> (GMP_LIMB_BITS - sc);
189 	  rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
190 	}
191     }
192 #endif
193 #endif
194 
195 #if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 3
196   if (sc > GMP_NAIL_BITS)
197     {
198       rp[2] = manl >> (GMP_LIMB_BITS - sc);
199       rp[1] = (manl << sc - GMP_NAIL_BITS) & GMP_NUMB_MASK;
200       if (sc >= 2 * GMP_NAIL_BITS)
201 	rp[0] = 0;
202       else
203 	rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK;
204     }
205   else
206     {
207       if (sc == 0)
208 	{
209 	  rp[2] = manl >> GMP_NAIL_BITS;
210 	  rp[1] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
211 	  rp[0] = 0;
212 	  exp--;
213 	}
214       else
215 	{
216 	  rp[2] = manl >> (GMP_LIMB_BITS - sc);
217 	  rp[1] = (manl >> GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
218 	  rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK;
219 	}
220     }
221 #endif
222 
223 #if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE == 3
224 #if GMP_NAIL_BITS == 0
225   if (sc != 0)
226     {
227       rp[2] = manh >> (GMP_LIMB_BITS - sc);
228       rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
229       rp[0] = manl << sc;
230     }
231   else
232     {
233       rp[2] = manh;
234       rp[1] = manl;
235       rp[0] = 0;
236       exp--;
237     }
238 #else
239   if (sc > GMP_NAIL_BITS)
240     {
241       rp[2] = (manh >> (GMP_LIMB_BITS - sc));
242       rp[1] = ((manh << (sc - GMP_NAIL_BITS)) |
243 	       (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK;
244       if (sc >= 2 * GMP_NAIL_BITS)
245 	rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
246       else
247 	rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
248     }
249   else
250     {
251       if (sc == 0)
252 	{
253 	  rp[2] = manh >> GMP_NAIL_BITS;
254 	  rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK;
255 	  rp[0] = (manl << GMP_NUMB_BITS - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
256 	  exp--;
257 	}
258       else
259 	{
260 	  rp[2] = (manh >> (GMP_LIMB_BITS - sc));
261 	  rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
262 	  rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc))
263 		   | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK;
264 	}
265     }
266 #endif
267 #endif
268 
269 #if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE > 3
270   if (sc == 0)
271     {
272       int i;
273 
274       for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--)
275 	{
276 	  rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
277 	  manh = ((manh << GMP_NUMB_BITS)
278 		  | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
279 	  manl = manl << GMP_NUMB_BITS;
280 	}
281       exp--;
282     }
283   else
284     {
285       int i;
286 
287       rp[LIMBS_PER_DOUBLE - 1] = (manh >> (GMP_LIMB_BITS - sc));
288       manh = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
289       manl = (manl << sc);
290       for (i = LIMBS_PER_DOUBLE - 2; i >= 0; i--)
291 	{
292 	  rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
293 	  manh = ((manh << GMP_NUMB_BITS)
294 		  | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
295 	  manl = manl << GMP_NUMB_BITS;
296 	}
297   }
298 #endif
299 
300   return exp;
301 }
302