1 /* mpfr_get_d, mpfr_get_d_2exp -- convert a multiple precision floating-point
2 number to a machine double precision float
3
4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel projects, INRIA.
6
7 This file is part of the GNU MPFR Library.
8
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23
24 #include <float.h>
25
26 #define MPFR_NEED_LONGLONG_H
27 #include "mpfr-impl.h"
28
29 #include "ieee_floats.h"
30
31 /* Assumes IEEE-754 double precision; otherwise, only an approximated
32 result will be returned, without any guaranty (and special cases
33 such as NaN must be avoided if not supported). */
34
35 double
mpfr_get_d(mpfr_srcptr src,mpfr_rnd_t rnd_mode)36 mpfr_get_d (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
37 {
38 double d;
39 int negative;
40 mpfr_exp_t e;
41
42 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
43 {
44 if (MPFR_IS_NAN (src))
45 return MPFR_DBL_NAN;
46
47 negative = MPFR_IS_NEG (src);
48
49 if (MPFR_IS_INF (src))
50 return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
51
52 MPFR_ASSERTD (MPFR_IS_ZERO(src));
53 return negative ? DBL_NEG_ZERO : 0.0;
54 }
55
56 e = MPFR_GET_EXP (src);
57 negative = MPFR_IS_NEG (src);
58
59 if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
60 rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;
61
62 /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
63 subnormal is 2^(-1074)=0.1e-1073 */
64 if (MPFR_UNLIKELY (e < -1073))
65 {
66 /* Note: Avoid using a constant expression DBL_MIN * DBL_EPSILON
67 as this gives 0 instead of the correct result with gcc on some
68 Alpha machines. */
69 d = negative ?
70 (rnd_mode == MPFR_RNDD ||
71 (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
72 ? -DBL_MIN : DBL_NEG_ZERO) :
73 (rnd_mode == MPFR_RNDU ||
74 (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
75 ? DBL_MIN : 0.0);
76 if (d != 0.0) /* we multiply DBL_MIN = 2^(-1022) by DBL_EPSILON = 2^(-52)
77 to get +-2^(-1074) */
78 d *= DBL_EPSILON;
79 }
80 /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
81 else if (MPFR_UNLIKELY (e > 1024))
82 {
83 d = negative ?
84 (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDU ?
85 -DBL_MAX : MPFR_DBL_INFM) :
86 (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD ?
87 DBL_MAX : MPFR_DBL_INFP);
88 }
89 else
90 {
91 int nbits;
92 mp_size_t np, i;
93 mp_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
94 int carry;
95
96 nbits = IEEE_DBL_MANT_DIG; /* 53 */
97 if (MPFR_UNLIKELY (e < -1021))
98 /*In the subnormal case, compute the exact number of significant bits*/
99 {
100 nbits += (1021 + e);
101 MPFR_ASSERTD (nbits >= 1);
102 }
103 np = MPFR_PREC2LIMBS (nbits);
104 MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
105 carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative,
106 nbits, rnd_mode);
107 if (MPFR_UNLIKELY(carry))
108 d = 1.0;
109 else
110 {
111 /* The following computations are exact thanks to the previous
112 mpfr_round_raw. */
113 d = (double) tp[0] / MP_BASE_AS_DOUBLE;
114 for (i = 1 ; i < np ; i++)
115 d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
116 /* d is the mantissa (between 1/2 and 1) of the argument rounded
117 to 53 bits */
118 }
119 d = mpfr_scale2 (d, e);
120 if (negative)
121 d = -d;
122 }
123
124 return d;
125 }
126
127 #undef mpfr_get_d1
128 double
mpfr_get_d1(mpfr_srcptr src)129 mpfr_get_d1 (mpfr_srcptr src)
130 {
131 return mpfr_get_d (src, __gmpfr_default_rounding_mode);
132 }
133
134 double
mpfr_get_d_2exp(long * expptr,mpfr_srcptr src,mpfr_rnd_t rnd_mode)135 mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
136 {
137 double ret;
138 mpfr_exp_t exp;
139 mpfr_t tmp;
140
141 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
142 {
143 int negative;
144 *expptr = 0;
145 if (MPFR_IS_NAN (src))
146 return MPFR_DBL_NAN;
147 negative = MPFR_IS_NEG (src);
148 if (MPFR_IS_INF (src))
149 return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
150 MPFR_ASSERTD (MPFR_IS_ZERO(src));
151 return negative ? DBL_NEG_ZERO : 0.0;
152 }
153
154 tmp[0] = *src; /* Hack copy mpfr_t */
155 MPFR_SET_EXP (tmp, 0);
156 ret = mpfr_get_d (tmp, rnd_mode);
157
158 if (MPFR_IS_PURE_FP(src))
159 {
160 exp = MPFR_GET_EXP (src);
161
162 /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
163 if (ret == 1.0)
164 {
165 ret = 0.5;
166 exp++;
167 }
168 else if (ret == -1.0)
169 {
170 ret = -0.5;
171 exp++;
172 }
173
174 MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
175 || (ret <= -0.5 && ret > -1.0));
176 MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
177 }
178 else
179 exp = 0;
180
181 *expptr = exp;
182 return ret;
183 }
184