| /* mpfr_set_ld -- convert a machine long double to |
| a multiple precision floating-point number |
| |
| Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. |
| Contributed by the Arenaire and Cacao projects, INRIA. |
| |
| This file is part of the GNU MPFR Library. |
| |
| The GNU MPFR Library is free software; you can redistribute it and/or modify |
| it under the terms of the GNU Lesser General Public License as published by |
| the Free Software Foundation; either version 2.1 of the License, or (at your |
| option) any later version. |
| |
| The GNU MPFR Library is distributed in the hope that it will be useful, but |
| WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
| or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
| License for more details. |
| |
| You should have received a copy of the GNU Lesser General Public License |
| along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to |
| the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
| MA 02110-1301, USA. */ |
| |
| #include <float.h> |
| |
| #define MPFR_NEED_LONGLONG_H |
| #include "mpfr-impl.h" |
| |
| /* Various i386 systems have been seen with float.h LDBL constants equal to |
| the DBL ones, whereas they ought to be bigger, reflecting the 10-byte |
| IEEE extended format on that processor. gcc 3.2.1 on FreeBSD and Solaris |
| has been seen with the problem, and gcc 2.95.4 on FreeBSD 4.7. */ |
| |
| #if HAVE_LDOUBLE_IEEE_EXT_LITTLE |
| static const struct { |
| char bytes[10]; |
| long double dummy; /* for memory alignment */ |
| } ldbl_max_struct = { |
| { '\377','\377','\377','\377', |
| '\377','\377','\377','\377', |
| '\376','\177' }, 0.0 |
| }; |
| #define MPFR_LDBL_MAX (* (const long double *) ldbl_max_struct.bytes) |
| #else |
| #define MPFR_LDBL_MAX LDBL_MAX |
| #endif |
| |
| #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE |
| |
| /* Generic code */ |
| int |
| mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode) |
| { |
| mpfr_t t, u; |
| int inexact, shift_exp; |
| long double x; |
| MPFR_SAVE_EXPO_DECL (expo); |
| |
| /* Check for NAN */ |
| LONGDOUBLE_NAN_ACTION (d, goto nan); |
| |
| /* Check for INF */ |
| if (d > MPFR_LDBL_MAX) |
| { |
| mpfr_set_inf (r, 1); |
| return 0; |
| } |
| else if (d < -MPFR_LDBL_MAX) |
| { |
| mpfr_set_inf (r, -1); |
| return 0; |
| } |
| /* Check for ZERO */ |
| else if (d == 0.0) |
| return mpfr_set_d (r, (double) d, rnd_mode); |
| |
| mpfr_init2 (t, MPFR_LDBL_MANT_DIG); |
| mpfr_init2 (u, IEEE_DBL_MANT_DIG); |
| |
| MPFR_SAVE_EXPO_MARK (expo); |
| |
| convert: |
| x = d; |
| MPFR_SET_ZERO (t); /* The sign doesn't matter. */ |
| shift_exp = 0; /* invariant: remainder to deal with is d*2^shift_exp */ |
| while (x != (long double) 0.0) |
| { |
| /* Check overflow of double */ |
| if (x > (long double) DBL_MAX || (-x) > (long double) DBL_MAX) |
| { |
| long double div9, div10, div11, div12, div13; |
| |
| #define TWO_64 18446744073709551616.0 /* 2^64 */ |
| #define TWO_128 (TWO_64 * TWO_64) |
| #define TWO_256 (TWO_128 * TWO_128) |
| div9 = (long double) (double) (TWO_256 * TWO_256); /* 2^(2^9) */ |
| div10 = div9 * div9; |
| div11 = div10 * div10; /* 2^(2^11) */ |
| div12 = div11 * div11; /* 2^(2^12) */ |
| div13 = div12 * div12; /* 2^(2^13) */ |
| if (ABS (x) >= div13) |
| { |
| x /= div13; /* exact */ |
| shift_exp += 8192; |
| } |
| if (ABS (x) >= div12) |
| { |
| x /= div12; /* exact */ |
| shift_exp += 4096; |
| } |
| if (ABS (x) >= div11) |
| { |
| x /= div11; /* exact */ |
| shift_exp += 2048; |
| } |
| if (ABS (x) >= div10) |
| { |
| x /= div10; /* exact */ |
| shift_exp += 1024; |
| } |
| /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024, |
| therefore we have one extra exponent reduction step */ |
| if (ABS (x) >= div9) |
| { |
| x /= div9; /* exact */ |
| shift_exp += 512; |
| } |
| } /* Check overflow of double */ |
| else |
| { |
| long double div9, div10, div11; |
| |
| div9 = (long double) (double) 7.4583407312002067432909653e-155; |
| /* div9 = 2^(-2^9) */ |
| div10 = div9 * div9; /* 2^(-2^10) */ |
| div11 = div10 * div10; /* 2^(-2^11) if extended precision */ |
| /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not |
| overflow here */ |
| if (ABS(x) < div10 && |
| div11 != (long double) 0.0 && |
| div11 / div10 == div10) /* possible underflow */ |
| { |
| long double div12, div13; |
| /* After the divisions, any bit of x must be >= div10, |
| hence the possible division by div9. */ |
| div12 = div11 * div11; /* 2^(-2^12) */ |
| div13 = div12 * div12; /* 2^(-2^13) */ |
| if (ABS (x) <= div13) |
| { |
| x /= div13; /* exact */ |
| shift_exp -= 8192; |
| } |
| if (ABS (x) <= div12) |
| { |
| x /= div12; /* exact */ |
| shift_exp -= 4096; |
| } |
| if (ABS (x) <= div11) |
| { |
| x /= div11; /* exact */ |
| shift_exp -= 2048; |
| } |
| if (ABS (x) <= div10) |
| { |
| x /= div10; /* exact */ |
| shift_exp -= 1024; |
| } |
| if (ABS(x) <= div9) |
| { |
| x /= div9; /* exact */ |
| shift_exp -= 512; |
| } |
| } |
| else |
| { |
| inexact = mpfr_set_d (u, (double) x, GMP_RNDZ); |
| MPFR_ASSERTD (inexact == 0); |
| if (mpfr_add (t, t, u, GMP_RNDZ) != 0) |
| { |
| if (!mpfr_number_p (t)) |
| break; |
| /* Inexact. This cannot happen unless the C implementation |
| "lies" on the precision or when long doubles are |
| implemented with FP expansions like under Mac OS X. */ |
| if (MPFR_PREC (t) != MPFR_PREC (r) + 1) |
| { |
| /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX. |
| The precision MPFR_PREC (r) + 1 allows us to |
| deduce the rounding bit and the sticky bit. */ |
| mpfr_set_prec (t, MPFR_PREC (r) + 1); |
| goto convert; |
| } |
| else |
| { |
| mp_limb_t *tp; |
| int rb_mask; |
| |
| /* Since mpfr_add was inexact, the sticky bit is 1. */ |
| tp = MPFR_MANT (t); |
| rb_mask = MPFR_LIMB_ONE << |
| (BITS_PER_MP_LIMB - 1 - |
| (MPFR_PREC (r) & (BITS_PER_MP_LIMB - 1))); |
| if (rnd_mode == GMP_RNDN) |
| rnd_mode = (*tp & rb_mask) ^ MPFR_IS_NEG (t) ? |
| GMP_RNDU : GMP_RNDD; |
| *tp |= rb_mask; |
| break; |
| } |
| } |
| x -= (long double) mpfr_get_d1 (u); /* exact */ |
| } |
| } |
| } |
| inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode); |
| mpfr_clear (t); |
| mpfr_clear (u); |
| |
| MPFR_SAVE_EXPO_FREE (expo); |
| return mpfr_check_range (r, inexact, rnd_mode); |
| |
| nan: |
| MPFR_SET_NAN(r); |
| MPFR_RET_NAN; |
| } |
| |
| #else /* IEEE Extended Little Endian Code */ |
| |
| int |
| mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode) |
| { |
| int inexact, i, k, cnt; |
| mpfr_t tmp; |
| mp_limb_t tmpmant[MPFR_LIMBS_PER_LONG_DOUBLE]; |
| mpfr_long_double_t x; |
| mp_exp_t exp; |
| int signd; |
| MPFR_SAVE_EXPO_DECL (expo); |
| |
| /* Check for NAN */ |
| if (MPFR_UNLIKELY (d != d)) |
| { |
| MPFR_SET_NAN (r); |
| MPFR_RET_NAN; |
| } |
| /* Check for INF */ |
| else if (MPFR_UNLIKELY (d > MPFR_LDBL_MAX)) |
| { |
| MPFR_SET_INF (r); |
| MPFR_SET_POS (r); |
| return 0; |
| } |
| else if (MPFR_UNLIKELY (d < -MPFR_LDBL_MAX)) |
| { |
| MPFR_SET_INF (r); |
| MPFR_SET_NEG (r); |
| return 0; |
| } |
| /* Check for ZERO */ |
| else if (MPFR_UNLIKELY (d == 0.0)) |
| { |
| x.ld = d; |
| MPFR_SET_ZERO (r); |
| if (x.s.sign == 1) |
| MPFR_SET_NEG(r); |
| else |
| MPFR_SET_POS(r); |
| return 0; |
| } |
| |
| /* now d is neither 0, nor NaN nor Inf */ |
| MPFR_SAVE_EXPO_MARK (expo); |
| |
| MPFR_MANT (tmp) = tmpmant; |
| MPFR_PREC (tmp) = 64; |
| |
| /* Extract sign */ |
| x.ld = d; |
| signd = MPFR_SIGN_POS; |
| if (x.ld < 0.0) |
| { |
| signd = MPFR_SIGN_NEG; |
| x.ld = -x.ld; |
| } |
| |
| /* Extract mantissa */ |
| #if BITS_PER_MP_LIMB >= 64 |
| tmpmant[0] = ((mp_limb_t) x.s.manh << 32) | ((mp_limb_t) x.s.manl); |
| #else |
| tmpmant[0] = (mp_limb_t) x.s.manl; |
| tmpmant[1] = (mp_limb_t) x.s.manh; |
| #endif |
| |
| /* Normalize mantissa */ |
| i = MPFR_LIMBS_PER_LONG_DOUBLE; |
| MPN_NORMALIZE_NOT_ZERO (tmpmant, i); |
| k = MPFR_LIMBS_PER_LONG_DOUBLE - i; |
| count_leading_zeros (cnt, tmpmant[i - 1]); |
| if (MPFR_LIKELY (cnt != 0)) |
| mpn_lshift (tmpmant + k, tmpmant, i, cnt); |
| else if (k != 0) |
| MPN_COPY (tmpmant + k, tmpmant, i); |
| if (MPFR_UNLIKELY (k != 0)) |
| MPN_ZERO (tmpmant, k); |
| |
| /* Set exponent */ |
| exp = (mp_exp_t) ((x.s.exph << 8) + x.s.expl); /* 15-bit unsigned int */ |
| if (MPFR_UNLIKELY (exp == 0)) |
| exp -= 0x3FFD; |
| else |
| exp -= 0x3FFE; |
| |
| MPFR_SET_EXP (tmp, exp - cnt - k * BITS_PER_MP_LIMB); |
| |
| /* tmp is exact */ |
| inexact = mpfr_set4 (r, tmp, rnd_mode, signd); |
| |
| MPFR_SAVE_EXPO_FREE (expo); |
| return mpfr_check_range (r, inexact, rnd_mode); |
| } |
| |
| #endif |
