Annotation of OpenXM_contrib/gmp/extract-dbl.c, Revision 1.1.1.2
1.1 maekawa 1: /* __gmp_extract_double -- convert from double to array of mp_limb_t.
2:
1.1.1.2 ! ohara 3: Copyright 1996, 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
1.1 maekawa 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:
1.1.1.2 ! ohara 33: #define BITS_IN_MANTISSA 53
! 34:
1.1 maekawa 35: /* Extract a non-negative double in d. */
36:
37: int
38: __gmp_extract_double (mp_ptr rp, double d)
39: {
40: long exp;
41: unsigned sc;
1.1.1.2 ! ohara 42: #ifdef _LONG_LONG_LIMB
! 43: #define BITS_PER_PART 64 /* somewhat bogus */
! 44: unsigned long long int manl;
! 45: #else
! 46: #define BITS_PER_PART GMP_LIMB_BITS
! 47: unsigned long int manh, manl;
! 48: #endif
1.1 maekawa 49:
50: /* BUGS
51:
52: 1. Should handle Inf and NaN in IEEE specific code.
53: 2. Handle Inf and NaN also in default code, to avoid hangs.
1.1.1.2 ! ohara 54: 3. Generalize to handle all GMP_LIMB_BITS >= 32.
1.1 maekawa 55: 4. This lits is incomplete and misspelled.
56: */
57:
1.1.1.2 ! ohara 58: ASSERT (d >= 0.0);
! 59:
1.1 maekawa 60: if (d == 0.0)
61: {
1.1.1.2 ! ohara 62: MPN_ZERO (rp, LIMBS_PER_DOUBLE);
1.1 maekawa 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;
1.1.1.2 ! ohara 75: #if BITS_PER_PART == 64 /* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */
1.1 maekawa 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: }
1.1.1.2 ! ohara 90: #endif
! 91: #if BITS_PER_PART == 32
1.1 maekawa 92: manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
93: manl = x.s.manl << 11;
94: if (exp == 0)
95: {
96: /* Denormalized number. Don't try to be clever about this,
97: since it is not an important case to make fast. */
98: exp = 1;
99: do
100: {
101: manh = (manh << 1) | (manl >> 31);
102: manl = manl << 1;
103: exp--;
104: }
105: while ((mp_limb_signed_t) manh >= 0);
106: }
107: #endif
1.1.1.2 ! ohara 108: #if BITS_PER_PART != 32 && BITS_PER_PART != 64
! 109: You need to generalize the code above to handle this.
! 110: #endif
1.1 maekawa 111: exp -= 1022; /* Remove IEEE bias. */
112: }
113: #else
114: {
115: /* Unknown (or known to be non-IEEE) double format. */
116: exp = 0;
117: if (d >= 1.0)
118: {
119: if (d * 0.5 == d)
120: abort ();
121:
122: while (d >= 32768.0)
123: {
124: d *= (1.0 / 65536.0);
125: exp += 16;
126: }
127: while (d >= 1.0)
128: {
129: d *= 0.5;
130: exp += 1;
131: }
132: }
133: else if (d < 0.5)
134: {
135: while (d < (1.0 / 65536.0))
136: {
137: d *= 65536.0;
138: exp -= 16;
139: }
140: while (d < 0.5)
141: {
142: d *= 2.0;
143: exp -= 1;
144: }
145: }
146:
1.1.1.2 ! ohara 147: d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
! 148: #if BITS_PER_PART == 64
1.1 maekawa 149: manl = d;
150: #else
151: manh = d;
1.1.1.2 ! ohara 152: manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
1.1 maekawa 153: #endif
154: }
1.1.1.2 ! ohara 155: #endif /* IEEE */
! 156:
! 157: /* Up until here, we have ignored the actual limb size. Remains
! 158: to split manh,,manl into an array of LIMBS_PER_DOUBLE limbs.
! 159: */
1.1 maekawa 160:
1.1.1.2 ! ohara 161: sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS;
1.1 maekawa 162:
163: /* We add something here to get rounding right. */
1.1.1.2 ! ohara 164: exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1;
1.1 maekawa 165:
1.1.1.2 ! ohara 166: #if LIMBS_PER_DOUBLE == 2
! 167: #if GMP_NAIL_BITS == 0
1.1 maekawa 168: if (sc != 0)
169: {
1.1.1.2 ! ohara 170: rp[1] = manl >> (GMP_LIMB_BITS - sc);
1.1 maekawa 171: rp[0] = manl << sc;
172: }
173: else
174: {
175: rp[1] = manl;
176: rp[0] = 0;
177: exp--;
178: }
179: #else
1.1.1.2 ! ohara 180: if (sc > GMP_NAIL_BITS)
! 181: {
! 182: rp[1] = manl >> (GMP_LIMB_BITS - sc);
! 183: rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK;
! 184: }
! 185: else
! 186: {
! 187: if (sc == 0)
! 188: {
! 189: rp[1] = manl >> GMP_NAIL_BITS;
! 190: rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
! 191: exp--;
! 192: }
! 193: else
! 194: {
! 195: rp[1] = manl >> (GMP_LIMB_BITS - sc);
! 196: rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
! 197: }
! 198: }
! 199: #endif
! 200: #endif
! 201:
! 202: #if LIMBS_PER_DOUBLE == 3
! 203: #if GMP_NAIL_BITS == 0
1.1 maekawa 204: if (sc != 0)
205: {
1.1.1.2 ! ohara 206: rp[2] = manh >> (GMP_LIMB_BITS - sc);
! 207: rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
1.1 maekawa 208: rp[0] = manl << sc;
209: }
210: else
211: {
212: rp[2] = manh;
213: rp[1] = manl;
214: rp[0] = 0;
215: exp--;
216: }
1.1.1.2 ! ohara 217: #else
! 218: if (sc > GMP_NAIL_BITS)
! 219: {
! 220: rp[2] = (manh >> (GMP_LIMB_BITS - sc));
! 221: rp[1] = ((manh << (sc - GMP_NAIL_BITS)) |
! 222: (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK;
! 223: if (sc >= 2 * GMP_NAIL_BITS)
! 224: rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
! 225: else
! 226: rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
! 227: }
! 228: else
! 229: {
! 230: if (sc == 0)
! 231: {
! 232: rp[2] = manh >> GMP_NAIL_BITS;
! 233: rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK;
! 234: rp[0] = 0;
! 235: exp--;
! 236: }
! 237: else
! 238: {
! 239: rp[2] = (manh >> (GMP_LIMB_BITS - sc));
! 240: rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
! 241: rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc))
! 242: | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK;
! 243: }
! 244: }
! 245: #endif
! 246: #endif
! 247:
! 248: #if LIMBS_PER_DOUBLE > 3
! 249: /* Insert code for splitting manh,,manl into LIMBS_PER_DOUBLE
! 250: mp_limb_t's at rp. */
! 251: if (sc != 0)
! 252: {
! 253: /* This is not perfect, and would fail for GMP_LIMB_BITS == 16.
! 254: The ASSERT_ALWAYS should catch the problematic cases. */
! 255: ASSERT_ALWAYS ((manl << sc) == 0);
! 256: manl = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
! 257: manh = manh >> (GMP_LIMB_BITS - sc);
! 258: }
! 259: {
! 260: int i;
! 261: for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--)
! 262: {
! 263: rp[i] = manh >> (BITS_PER_LONGINT - GMP_LIMB_BITS);
! 264: manh = ((manh << GMP_LIMB_BITS)
! 265: | (manl >> (BITS_PER_LONGINT - GMP_LIMB_BITS)));
! 266: manl = manl << GMP_LIMB_BITS;
! 267: }
! 268: }
1.1 maekawa 269: #endif
270:
271: return exp;
272: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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