![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to extract-dbl.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp|
Return to extract-dbl.c CVS log | Up to [local] / OpenXM_contrib / gmp |
1.1 maekawa 1: /* __gmp_extract_double -- convert from double to array of mp_limb_t.
2:
3: Copyright (C) 1996, 1999, 2000 Free Software Foundation, Inc.
4:
5: This file is part of the GNU MP Library.
6:
7: The GNU MP Library is free software; you can redistribute it and/or modify
8: it under the terms of the GNU Lesser General Public License as published by
9: the Free Software Foundation; either version 2.1 of the License, or (at your
10: option) any later version.
11:
12: The GNU MP Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Lesser General Public License
18: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include "gmp.h"
23: #include "gmp-impl.h"
24:
25: #ifdef XDEBUG
26: #undef _GMP_IEEE_FLOATS
27: #endif
28:
29: #ifndef _GMP_IEEE_FLOATS
30: #define _GMP_IEEE_FLOATS 0
31: #endif
32:
33: /* Extract a non-negative double in d. */
34:
35: int
36: #if __STDC__
37: __gmp_extract_double (mp_ptr rp, double d)
38: #else
39: __gmp_extract_double (rp, d)
40: mp_ptr rp;
41: double d;
42: #endif
43: {
44: long exp;
45: unsigned sc;
46: mp_limb_t manh, manl;
47:
48: /* BUGS
49:
50: 1. Should handle Inf and NaN in IEEE specific code.
51: 2. Handle Inf and NaN also in default code, to avoid hangs.
52: 3. Generalize to handle all BITS_PER_MP_LIMB >= 32.
53: 4. This lits is incomplete and misspelled.
54: */
55:
56: if (d == 0.0)
57: {
58: rp[0] = 0;
59: rp[1] = 0;
60: #if BITS_PER_MP_LIMB == 32
61: rp[2] = 0;
62: #endif
63: return 0;
64: }
65:
66: #if _GMP_IEEE_FLOATS
67: {
68: #if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8
69: /* Work around alpha-specific bug in GCC 2.8.x. */
70: volatile
71: #endif
72: union ieee_double_extract x;
73: x.d = d;
74: exp = x.s.exp;
75: #if BITS_PER_MP_LIMB == 64
76: manl = (((mp_limb_t) 1 << 63)
77: | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
78: if (exp == 0)
79: {
80: /* Denormalized number. Don't try to be clever about this,
81: since it is not an important case to make fast. */
82: exp = 1;
83: do
84: {
85: manl = manl << 1;
86: exp--;
87: }
88: while ((mp_limb_signed_t) manl >= 0);
89: }
90: #else
91: manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
92: manl = x.s.manl << 11;
93: if (exp == 0)
94: {
95: /* Denormalized number. Don't try to be clever about this,
96: since it is not an important case to make fast. */
97: exp = 1;
98: do
99: {
100: manh = (manh << 1) | (manl >> 31);
101: manl = manl << 1;
102: exp--;
103: }
104: while ((mp_limb_signed_t) manh >= 0);
105: }
106: #endif
107: exp -= 1022; /* Remove IEEE bias. */
108: }
109: #else
110: {
111: /* Unknown (or known to be non-IEEE) double format. */
112: exp = 0;
113: if (d >= 1.0)
114: {
115: if (d * 0.5 == d)
116: abort ();
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 *= MP_BASE_AS_DOUBLE;
144: #if BITS_PER_MP_LIMB == 64
145: manl = d;
146: #else
147: manh = d;
148: manl = (d - manh) * MP_BASE_AS_DOUBLE;
149: #endif
150: }
151: #endif
152:
153: sc = (unsigned) exp % BITS_PER_MP_LIMB;
154:
155: /* We add something here to get rounding right. */
156: exp = (exp + 2048) / BITS_PER_MP_LIMB - 2048 / BITS_PER_MP_LIMB + 1;
157:
158: #if BITS_PER_MP_LIMB == 64
159: if (sc != 0)
160: {
161: rp[1] = manl >> (BITS_PER_MP_LIMB - sc);
162: rp[0] = manl << sc;
163: }
164: else
165: {
166: rp[1] = manl;
167: rp[0] = 0;
168: exp--;
169: }
170: #else
171: if (sc != 0)
172: {
173: rp[2] = manh >> (BITS_PER_MP_LIMB - sc);
174: rp[1] = (manl >> (BITS_PER_MP_LIMB - sc)) | (manh << sc);
175: rp[0] = manl << sc;
176: }
177: else
178: {
179: rp[2] = manh;
180: rp[1] = manl;
181: rp[0] = 0;
182: exp--;
183: }
184: #endif
185:
186: return exp;
187: }