Annotation of OpenXM_contrib/gmp/mpfr/set_d.c, Revision 1.1.1.2
1.1.1.2 ! ohara 1: /* mpfr_set_d -- convert a machine double precision float to
! 2: a multiple precision floating-point number
1.1 maekawa 3:
1.1.1.2 ! ohara 4: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
1.1 maekawa 5:
6: This file is part of the MPFR Library.
7:
8: The MPFR Library is free software; you can redistribute it and/or modify
1.1.1.2 ! ohara 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 11: option) any later version.
12:
13: The MPFR Library is distributed in the hope that it will be useful, but
14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! ohara 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 16: License for more details.
17:
1.1.1.2 ! ohara 18: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 19: along with the MPFR Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: #include "gmp.h"
24: #include "gmp-impl.h"
25: #include "longlong.h"
26: #include "mpfr.h"
1.1.1.2 ! ohara 27: #include "mpfr-impl.h"
! 28:
! 29: #if (BITS_PER_MP_LIMB==32)
! 30: #define MPFR_LIMBS_PER_DOUBLE 2
! 31: #elif (BITS_PER_MP_LIMB >= 64)
! 32: #define MPFR_LIMBS_PER_DOUBLE 1
! 33: #elif (BITS_PER_MP_LIMB == 16)
! 34: #define MPFR_LIMBS_PER_DOUBLE 4
! 35: #endif
1.1 maekawa 36:
1.1.1.2 ! ohara 37: static int __mpfr_extract_double _PROTO ((mp_ptr, double));
1.1 maekawa 38:
39: /* Included from gmp-2.0.2, patched to support denorms */
40:
41: #ifdef XDEBUG
42: #undef _GMP_IEEE_FLOATS
43: #endif
44:
45: #ifndef _GMP_IEEE_FLOATS
46: #define _GMP_IEEE_FLOATS 0
47: #endif
48:
1.1.1.2 ! ohara 49: static int
! 50: __mpfr_extract_double (mp_ptr rp, double d)
1.1 maekawa 51: /* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */
52: {
53: long exp;
54: mp_limb_t manl;
55: #if BITS_PER_MP_LIMB == 32
56: mp_limb_t manh;
57: #endif
58:
59: /* BUGS
60:
61: 1. Should handle Inf and NaN in IEEE specific code.
62: 2. Handle Inf and NaN also in default code, to avoid hangs.
63: 3. Generalize to handle all BITS_PER_MP_LIMB >= 32.
64: 4. This lits is incomplete and misspelled.
65: */
66:
67: if (d == 0.0)
68: {
69: rp[0] = 0;
70: return 0;
71: }
72:
73: #if _GMP_IEEE_FLOATS
74: {
75: union ieee_double_extract x;
76: x.d = d;
77:
78: exp = x.s.exp;
1.1.1.2 ! ohara 79: if (exp)
1.1 maekawa 80: {
81: #if BITS_PER_MP_LIMB == 64
1.1.1.2 ! ohara 82: manl = ((MP_LIMB_T_ONE << 63)
1.1 maekawa 83: | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
84: #else
1.1.1.2 ! ohara 85: manh = (MP_LIMB_T_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
! 86: manl = x.s.manl << 11;
1.1 maekawa 87: #endif
88: }
1.1.1.2 ! ohara 89: else /* denormalized number */
1.1 maekawa 90: {
91: #if BITS_PER_MP_LIMB == 64
92: manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11);
93: #else
1.1.1.2 ! ohara 94: manh = (x.s.manh << 11) /* high 21 bits */
! 95: | (x.s.manl >> 21); /* middle 11 bits */
! 96: manl = x.s.manl << 11; /* low 21 bits */
1.1 maekawa 97: #endif
98: }
99: }
100: #else
101: {
102: /* Unknown (or known to be non-IEEE) double format. */
103: exp = 0;
104: if (d >= 1.0)
105: {
106: if (d * 0.5 == d)
107: abort ();
108:
109: while (d >= 32768.0)
110: {
111: d *= (1.0 / 65536.0);
112: exp += 16;
113: }
114: while (d >= 1.0)
115: {
116: d *= 0.5;
117: exp += 1;
118: }
119: }
120: else if (d < 0.5)
121: {
122: while (d < (1.0 / 65536.0))
123: {
124: d *= 65536.0;
125: exp -= 16;
126: }
127: while (d < 0.5)
128: {
129: d *= 2.0;
130: exp -= 1;
131: }
132: }
133:
134: d *= MP_BASE_AS_DOUBLE;
135: #if BITS_PER_MP_LIMB == 64
136: manl = d;
137: #else
138: manh = d;
139: manl = (d - manh) * MP_BASE_AS_DOUBLE;
140: #endif
141:
142: exp += 1022;
143: }
144: #endif
145:
1.1.1.2 ! ohara 146: if (exp) exp = (unsigned) exp - 1022; else exp = -1021;
1.1 maekawa 147:
148: #if BITS_PER_MP_LIMB == 64
1.1.1.2 ! ohara 149: rp[0] = manl;
1.1 maekawa 150: #else
1.1.1.2 ! ohara 151: rp[1] = manh;
! 152: rp[0] = manl;
1.1 maekawa 153: #endif
154:
155: return exp;
156: }
157:
158: /* End of part included from gmp-2.0.2 */
159:
1.1.1.2 ! ohara 160: int
! 161: mpfr_set_d (mpfr_ptr r, double d, mp_rnd_t rnd_mode)
! 162: {
! 163: int signd, sizetmp, inexact;
! 164: unsigned int cnt, k;
! 165: mpfr_ptr tmp;
! 166: TMP_DECL(marker);
1.1 maekawa 167:
1.1.1.2 ! ohara 168: TMP_MARK(marker);
! 169: MPFR_CLEAR_FLAGS(r);
1.1 maekawa 170:
1.1.1.2 ! ohara 171: if (d == 0)
! 172: {
! 173: union ieee_double_extract x;
1.1 maekawa 174:
1.1.1.2 ! ohara 175: MPFR_SET_ZERO(r);
! 176: /* set correct sign */
! 177: x.d = d;
! 178: if (((x.s.sig == 1) && (MPFR_SIGN(r) > 0))
! 179: || ((x.s.sig == 0) && (MPFR_SIGN(r) < 0)))
! 180: MPFR_CHANGE_SIGN(r);
! 181: return 0; /* 0 is exact */
! 182: }
! 183: else if (DOUBLE_ISNAN(d))
! 184: {
! 185: MPFR_SET_NAN(r);
! 186: MPFR_RET_NAN;
! 187: }
! 188: else if (DOUBLE_ISINF(d))
! 189: {
! 190: MPFR_SET_INF(r);
! 191: if ((d > 0 && (MPFR_SIGN(r) == -1)) || (d < 0 && (MPFR_SIGN(r) == 1)))
! 192: MPFR_CHANGE_SIGN(r);
! 193: return 0; /* infinity is exact */
! 194: }
1.1 maekawa 195:
1.1.1.2 ! ohara 196: /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE,
! 197: since PREC(r) may be different from PREC(tmp), and then both variables
! 198: would have same precision in the mpfr_set4 call below. */
! 199: tmp = (mpfr_ptr) TMP_ALLOC(sizeof(mpfr_t));
! 200: MPFR_MANT(tmp) = TMP_ALLOC(MPFR_LIMBS_PER_DOUBLE * BYTES_PER_MP_LIMB);
! 201: MPFR_PREC(tmp) = 53;
! 202: MPFR_SIZE(tmp) = MPFR_LIMBS_PER_DOUBLE;
! 203: sizetmp = MPFR_LIMBS_PER_DOUBLE;
1.1 maekawa 204:
205: signd = (d < 0) ? -1 : 1;
206: d = ABS (d);
207:
1.1.1.2 ! ohara 208: MPFR_EXP(tmp) = __mpfr_extract_double (MPFR_MANT(tmp), d);
1.1 maekawa 209:
1.1.1.2 ! ohara 210: /* determine number k of zero high limbs */
! 211: for (k = 0; k < sizetmp && MPFR_MANT(tmp)[sizetmp - 1 - k] == 0; k++);
1.1 maekawa 212:
1.1.1.2 ! ohara 213: count_leading_zeros (cnt, MPFR_MANT(tmp)[sizetmp - 1 - k]);
1.1 maekawa 214:
1.1.1.2 ! ohara 215: if (cnt)
! 216: mpn_lshift (MPFR_MANT(tmp) + k, MPFR_MANT(tmp), sizetmp - k, cnt);
! 217: else if (k)
! 218: MPN_COPY (MPFR_MANT(tmp) + k, MPFR_MANT(tmp), sizetmp - k);
! 219: if (k)
! 220: MPN_ZERO (MPFR_MANT(tmp), k);
1.1 maekawa 221:
1.1.1.2 ! ohara 222: MPFR_EXP(tmp) -= cnt + k * BITS_PER_MP_LIMB;
! 223:
! 224: /* tmp is exact since PREC(tmp)=53 */
! 225: inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
1.1 maekawa 226:
1.1.1.2 ! ohara 227: TMP_FREE(marker);
! 228: return inexact;
! 229: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to set_d.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr