Annotation of OpenXM_contrib/gmp/mpfr/get_d.c, Revision 1.1.1.1
1.1 ohara 1: /* mpfr_get_d -- convert a multiple precision floating-point number
2: to a machine double precision float
3:
4: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
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
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
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
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 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 <float.h>
24: #ifndef NO_MATH_DEFS
25: #include <math.h>
26: #endif
27: #include "gmp.h"
28: #include "gmp-impl.h"
29: #include "longlong.h"
30: #include "mpfr.h"
31: #include "mpfr-impl.h"
32: #include "mpfr-math.h"
33:
34: static double mpfr_scale2 _PROTO ((double, int));
35:
36: #define IEEE_DBL_MANT_DIG 53
37:
38: /* multiplies 1/2 <= d <= 1 by 2^exp */
39: static double
40: mpfr_scale2 (double d, int exp)
41: {
42: #if _GMP_IEEE_FLOATS
43: {
44: union ieee_double_extract x;
45:
46: if (d == 1.0)
47: {
48: d = 0.5;
49: exp ++;
50: }
51:
52: /* now 1/2 <= d < 1 */
53:
54: /* infinities and zeroes have already been checked */
55: MPFR_ASSERTN(-1073 <= exp && exp <= 1025);
56:
57: x.d = d;
58: if (exp < -1021) /* subnormal case */
59: {
60: x.s.exp += exp + 52;
61: x.d *= DBL_EPSILON;
62: }
63: else /* normalized case */
64: {
65: x.s.exp += exp;
66: }
67: return x.d;
68: }
69: #else
70: {
71: double factor;
72:
73: if (exp < 0)
74: {
75: factor = 0.5;
76: exp = -exp;
77: }
78: else
79: {
80: factor = 2.0;
81: }
82: while (exp != 0)
83: {
84: if ((exp & 1) != 0)
85: d *= factor;
86: exp >>= 1;
87: factor *= factor;
88: }
89: return d;
90: }
91: #endif
92: }
93:
94: /* Assumes IEEE-754 double precision; otherwise, only an approximated
95: result will be returned, without any guaranty (and special cases
96: such as NaN must be avoided if not supported). */
97:
98: double
99: mpfr_get_d3 (mpfr_srcptr src, mp_exp_t e, mp_rnd_t rnd_mode)
100: {
101: double d;
102: int negative;
103:
104: if (MPFR_IS_NAN(src))
105: return MPFR_DBL_NAN;
106:
107: negative = MPFR_SIGN(src) < 0;
108:
109: if (MPFR_IS_INF(src))
110: return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
111:
112: if (MPFR_IS_ZERO(src))
113: return negative ? -0.0 : 0.0;
114:
115: /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
116: subnormal is 2^(-1074)=0.1e-1073 */
117: if (e < -1073)
118: {
119: d = negative ?
120: (rnd_mode == GMP_RNDD ||
121: (rnd_mode == GMP_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
122: ? -DBL_MIN : -0.0) :
123: (rnd_mode == GMP_RNDU ||
124: (rnd_mode == GMP_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
125: ? DBL_MIN : 0.0);
126: if (d != 0.0)
127: d *= DBL_EPSILON;
128: }
129: /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
130: else if (e > 1024)
131: {
132: d = negative ?
133: (rnd_mode == GMP_RNDZ || rnd_mode == GMP_RNDU ?
134: -DBL_MAX : MPFR_DBL_INFM) :
135: (rnd_mode == GMP_RNDZ || rnd_mode == GMP_RNDD ?
136: DBL_MAX : MPFR_DBL_INFP);
137: }
138: else
139: {
140: int nbits;
141: mp_size_t np, i;
142: mp_ptr tp;
143: int carry;
144:
145: nbits = IEEE_DBL_MANT_DIG; /* 53 */
146: if (e < -1021) /* in the subnormal case, compute the exact number of
147: significant bits */
148: {
149: nbits += (1021 + e);
150: MPFR_ASSERTN(nbits >= 1);
151: }
152: np = (nbits - 1) / BITS_PER_MP_LIMB;
153: tp = (mp_ptr) (*__gmp_allocate_func) ((np + 1) * BYTES_PER_MP_LIMB);
154: carry = mpfr_round_raw (tp, MPFR_MANT(src), MPFR_PREC(src), negative,
155: nbits, rnd_mode, (int *) 0);
156: if (carry)
157: d = 1.0;
158: else
159: {
160: /* Warning: the rounding may still be incorrect in the rounding
161: to the nearest mode when the result is a subnormal because of
162: a double rounding (-> 53 bits -> final precision). */
163: d = (double) tp[0] / MP_BASE_AS_DOUBLE;
164: for (i = 1; i <= np; i++)
165: d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
166: /* d is the mantissa (between 1/2 and 1) of the argument rounded
167: to 53 bits */
168: }
169: d = mpfr_scale2 (d, e);
170: if (negative)
171: d = -d;
172:
173: (*__gmp_free_func) (tp, (np + 1) * BYTES_PER_MP_LIMB);
174: }
175:
176: return d;
177: }
178:
179: /* Note: do not read the exponent if it has no meaning (avoid possible
180: traps on some implementations). */
181:
182: double
183: mpfr_get_d (mpfr_srcptr src, mp_rnd_t rnd_mode)
184: {
185: return mpfr_get_d3 (src, MPFR_IS_FP(src) && MPFR_NOTZERO(src) ?
186: MPFR_EXP(src) : 0, rnd_mode);
187: }
188:
189: double
190: mpfr_get_d1 (mpfr_srcptr src)
191: {
192: return mpfr_get_d3 (src, MPFR_IS_FP(src) && MPFR_NOTZERO(src) ?
193: MPFR_EXP(src) : 0, __gmp_default_rounding_mode);
194: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to get_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