--- /dev/null
+/* 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
projects / msp430-gcc.git / blobdiff