1 /* mpfr_get_d, mpfr_get_d_2exp -- convert a multiple precision floating-point
2 number to a machine double precision float
4 Copyright 1999-2015 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel projects, INRIA.
7 This file is part of the GNU MPFR Library.
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.
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.
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. */
26 #define MPFR_NEED_LONGLONG_H
27 #include "mpfr-impl.h"
29 #include "ieee_floats.h"
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). */
36 mpfr_get_d (mpfr_srcptr src
, mpfr_rnd_t rnd_mode
)
42 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src
)))
44 if (MPFR_IS_NAN (src
))
47 negative
= MPFR_IS_NEG (src
);
49 if (MPFR_IS_INF (src
))
50 return negative
? MPFR_DBL_INFM
: MPFR_DBL_INFP
;
52 MPFR_ASSERTD (MPFR_IS_ZERO(src
));
53 return negative
? DBL_NEG_ZERO
: 0.0;
56 e
= MPFR_GET_EXP (src
);
57 negative
= MPFR_IS_NEG (src
);
59 if (MPFR_UNLIKELY(rnd_mode
== MPFR_RNDA
))
60 rnd_mode
= negative
? MPFR_RNDD
: MPFR_RNDU
;
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))
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
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)
76 if (d
!= 0.0) /* we multiply DBL_MIN = 2^(-1022) by DBL_EPSILON = 2^(-52)
80 /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
81 else if (MPFR_UNLIKELY (e
> 1024))
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
);
93 mp_limb_t tp
[ MPFR_LIMBS_PER_DOUBLE
];
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*/
101 MPFR_ASSERTD (nbits
>= 1);
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
,
107 if (MPFR_UNLIKELY(carry
))
111 /* The following computations are exact thanks to the previous
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
119 d
= mpfr_scale2 (d
, e
);
129 mpfr_get_d1 (mpfr_srcptr src
)
131 return mpfr_get_d (src
, __gmpfr_default_rounding_mode
);
135 mpfr_get_d_2exp (long *expptr
, mpfr_srcptr src
, mpfr_rnd_t rnd_mode
)
141 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src
)))
145 if (MPFR_IS_NAN (src
))
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;
154 tmp
[0] = *src
; /* Hack copy mpfr_t */
155 MPFR_SET_EXP (tmp
, 0);
156 ret
= mpfr_get_d (tmp
, rnd_mode
);
158 if (MPFR_IS_PURE_FP(src
))
160 exp
= MPFR_GET_EXP (src
);
162 /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
168 else if (ret
== -1.0)
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
);