1 /* mpfr_set_ld -- convert a machine long double to
2 a multiple precision floating-point number
4 Copyright 2002-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 /* Various i386 systems have been seen with <float.h> LDBL constants equal
30 to the DBL ones, whereas they ought to be bigger, reflecting the 10-byte
31 IEEE extended format on that processor. gcc 3.2.1 on FreeBSD and Solaris
32 has been seen with the problem, and gcc 2.95.4 on FreeBSD 4.7. */
34 #if HAVE_LDOUBLE_IEEE_EXT_LITTLE
39 { '\377','\377','\377','\377',
40 '\377','\377','\377','\377',
43 #define MPFR_LDBL_MAX (ldbl_max_struct.d)
45 #define MPFR_LDBL_MAX LDBL_MAX
48 #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE
52 mpfr_set_ld (mpfr_ptr r
, long double d
, mpfr_rnd_t rnd_mode
)
55 int inexact
, shift_exp
;
57 MPFR_SAVE_EXPO_DECL (expo
);
60 LONGDOUBLE_NAN_ACTION (d
, goto nan
);
63 if (d
> MPFR_LDBL_MAX
)
68 else if (d
< -MPFR_LDBL_MAX
)
75 return mpfr_set_d (r
, (double) d
, rnd_mode
);
77 mpfr_init2 (t
, MPFR_LDBL_MANT_DIG
);
78 mpfr_init2 (u
, IEEE_DBL_MANT_DIG
);
80 MPFR_SAVE_EXPO_MARK (expo
);
84 MPFR_SET_ZERO (t
); /* The sign doesn't matter. */
85 shift_exp
= 0; /* invariant: remainder to deal with is d*2^shift_exp */
86 while (x
!= (long double) 0.0)
88 /* Check overflow of double */
89 if (x
> (long double) DBL_MAX
|| (-x
) > (long double) DBL_MAX
)
91 long double div9
, div10
, div11
, div12
, div13
;
93 #define TWO_64 18446744073709551616.0 /* 2^64 */
94 #define TWO_128 (TWO_64 * TWO_64)
95 #define TWO_256 (TWO_128 * TWO_128)
96 div9
= (long double) (double) (TWO_256
* TWO_256
); /* 2^(2^9) */
98 div11
= div10
* div10
; /* 2^(2^11) */
99 div12
= div11
* div11
; /* 2^(2^12) */
100 div13
= div12
* div12
; /* 2^(2^13) */
101 if (ABS (x
) >= div13
)
103 x
/= div13
; /* exact */
105 mpfr_div_2si (t
, t
, 8192, MPFR_RNDZ
);
107 if (ABS (x
) >= div12
)
109 x
/= div12
; /* exact */
111 mpfr_div_2si (t
, t
, 4096, MPFR_RNDZ
);
113 if (ABS (x
) >= div11
)
115 x
/= div11
; /* exact */
117 mpfr_div_2si (t
, t
, 2048, MPFR_RNDZ
);
119 if (ABS (x
) >= div10
)
121 x
/= div10
; /* exact */
123 mpfr_div_2si (t
, t
, 1024, MPFR_RNDZ
);
125 /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024,
126 therefore we have one extra exponent reduction step */
129 x
/= div9
; /* exact */
131 mpfr_div_2si (t
, t
, 512, MPFR_RNDZ
);
133 } /* Check overflow of double */
134 else /* no overflow on double */
136 long double div9
, div10
, div11
;
138 div9
= (long double) (double) 7.4583407312002067432909653e-155;
139 /* div9 = 2^(-2^9) */
140 div10
= div9
* div9
; /* 2^(-2^10) */
141 div11
= div10
* div10
; /* 2^(-2^11) if extended precision */
142 /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not
144 if (ABS(x
) < div10
&&
145 div11
!= (long double) 0.0 &&
146 div11
/ div10
== div10
) /* possible underflow */
148 long double div12
, div13
;
149 /* After the divisions, any bit of x must be >= div10,
150 hence the possible division by div9. */
151 div12
= div11
* div11
; /* 2^(-2^12) */
152 div13
= div12
* div12
; /* 2^(-2^13) */
153 if (ABS (x
) <= div13
)
155 x
/= div13
; /* exact */
157 mpfr_mul_2si (t
, t
, 8192, MPFR_RNDZ
);
159 if (ABS (x
) <= div12
)
161 x
/= div12
; /* exact */
163 mpfr_mul_2si (t
, t
, 4096, MPFR_RNDZ
);
165 if (ABS (x
) <= div11
)
167 x
/= div11
; /* exact */
169 mpfr_mul_2si (t
, t
, 2048, MPFR_RNDZ
);
171 if (ABS (x
) <= div10
)
173 x
/= div10
; /* exact */
175 mpfr_mul_2si (t
, t
, 1024, MPFR_RNDZ
);
179 x
/= div9
; /* exact */
181 mpfr_mul_2si (t
, t
, 512, MPFR_RNDZ
);
184 else /* no underflow */
186 inexact
= mpfr_set_d (u
, (double) x
, MPFR_RNDZ
);
187 MPFR_ASSERTD (inexact
== 0);
188 if (mpfr_add (t
, t
, u
, MPFR_RNDZ
) != 0)
190 if (!mpfr_number_p (t
))
192 /* Inexact. This cannot happen unless the C implementation
193 "lies" on the precision or when long doubles are
194 implemented with FP expansions like under Mac OS X. */
195 if (MPFR_PREC (t
) != MPFR_PREC (r
) + 1)
197 /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX.
198 The precision MPFR_PREC (r) + 1 allows us to
199 deduce the rounding bit and the sticky bit. */
200 mpfr_set_prec (t
, MPFR_PREC (r
) + 1);
208 /* Since mpfr_add was inexact, the sticky bit is 1. */
210 rb_mask
= MPFR_LIMB_ONE
<<
212 (MPFR_PREC (r
) & (GMP_NUMB_BITS
- 1)));
213 if (rnd_mode
== MPFR_RNDN
)
214 rnd_mode
= (*tp
& rb_mask
) ^ MPFR_IS_NEG (t
) ?
215 MPFR_RNDU
: MPFR_RNDD
;
220 x
-= (long double) mpfr_get_d1 (u
); /* exact */
224 inexact
= mpfr_mul_2si (r
, t
, shift_exp
, rnd_mode
);
228 MPFR_SAVE_EXPO_FREE (expo
);
229 return mpfr_check_range (r
, inexact
, rnd_mode
);
236 #else /* IEEE Extended Little Endian Code */
239 mpfr_set_ld (mpfr_ptr r
, long double d
, mpfr_rnd_t rnd_mode
)
241 int inexact
, i
, k
, cnt
;
243 mp_limb_t tmpmant
[MPFR_LIMBS_PER_LONG_DOUBLE
];
244 mpfr_long_double_t x
;
247 MPFR_SAVE_EXPO_DECL (expo
);
250 if (MPFR_UNLIKELY (d
!= d
))
256 else if (MPFR_UNLIKELY (d
> MPFR_LDBL_MAX
))
262 else if (MPFR_UNLIKELY (d
< -MPFR_LDBL_MAX
))
269 else if (MPFR_UNLIKELY (d
== 0.0))
280 /* now d is neither 0, nor NaN nor Inf */
281 MPFR_SAVE_EXPO_MARK (expo
);
283 MPFR_MANT (tmp
) = tmpmant
;
284 MPFR_PREC (tmp
) = 64;
288 signd
= MPFR_SIGN_POS
;
291 signd
= MPFR_SIGN_NEG
;
295 /* Extract mantissa */
296 #if GMP_NUMB_BITS >= 64
297 tmpmant
[0] = ((mp_limb_t
) x
.s
.manh
<< 32) | ((mp_limb_t
) x
.s
.manl
);
299 tmpmant
[0] = (mp_limb_t
) x
.s
.manl
;
300 tmpmant
[1] = (mp_limb_t
) x
.s
.manh
;
303 /* Normalize mantissa */
304 i
= MPFR_LIMBS_PER_LONG_DOUBLE
;
305 MPN_NORMALIZE_NOT_ZERO (tmpmant
, i
);
306 k
= MPFR_LIMBS_PER_LONG_DOUBLE
- i
;
307 count_leading_zeros (cnt
, tmpmant
[i
- 1]);
308 if (MPFR_LIKELY (cnt
!= 0))
309 mpn_lshift (tmpmant
+ k
, tmpmant
, i
, cnt
);
311 MPN_COPY (tmpmant
+ k
, tmpmant
, i
);
312 if (MPFR_UNLIKELY (k
!= 0))
313 MPN_ZERO (tmpmant
, k
);
316 exp
= (mpfr_exp_t
) ((x
.s
.exph
<< 8) + x
.s
.expl
); /* 15-bit unsigned int */
317 if (MPFR_UNLIKELY (exp
== 0))
322 MPFR_SET_EXP (tmp
, exp
- cnt
- k
* GMP_NUMB_BITS
);
325 inexact
= mpfr_set4 (r
, tmp
, rnd_mode
, signd
);
327 MPFR_SAVE_EXPO_FREE (expo
);
328 return mpfr_check_range (r
, inexact
, rnd_mode
);