contrib/mpfr/set_ld.c

   1 /* mpfr_set_ld -- convert a machine long double to
   2                   a multiple precision floating-point number
   3
   4 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  22 MA 02110-1301, USA. */
  23
  24 #include <float.h>
  25
  26 #define MPFR_NEED_LONGLONG_H
  27 #include "mpfr-impl.h"
  28
  29 /* Various i386 systems have been seen with float.h LDBL constants equal to
  30    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.  */
  33
  34 #if HAVE_LDOUBLE_IEEE_EXT_LITTLE
  35 static const union {
  36   char         bytes[10];
  37   long double  d;
  38 } ldbl_max_struct = {
  39   { '\377','\377','\377','\377',
  40     '\377','\377','\377','\377',
  41     '\376','\177' }
  42 };
  43 #define MPFR_LDBL_MAX   (ldbl_max_struct.d)
  44 #else
  45 #define MPFR_LDBL_MAX   LDBL_MAX
  46 #endif
  47
  48 #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE
  49
  50 /* Generic code */
  51 int
  52 mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode)
  53 {
  54   mpfr_t t, u;
  55   int inexact, shift_exp;
  56   long double x;
  57   MPFR_SAVE_EXPO_DECL (expo);
  58
  59   /* Check for NAN */
  60   LONGDOUBLE_NAN_ACTION (d, goto nan);
  61
  62   /* Check for INF */
  63   if (d > MPFR_LDBL_MAX)
  64     {
  65       mpfr_set_inf (r, 1);
  66       return 0;
  67     }
  68   else if (d < -MPFR_LDBL_MAX)
  69     {
  70       mpfr_set_inf (r, -1);
  71       return 0;
  72     }
  73   /* Check for ZERO */
  74   else if (d == 0.0)
  75     return mpfr_set_d (r, (double) d, rnd_mode);
  76
  77   mpfr_init2 (t, MPFR_LDBL_MANT_DIG);
  78   mpfr_init2 (u, IEEE_DBL_MANT_DIG);
  79
  80   MPFR_SAVE_EXPO_MARK (expo);
  81
  82  convert:
  83   x = d;
  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)
  87     {
  88       /* Check overflow of double */
  89       if (x > (long double) DBL_MAX || (-x) > (long double) DBL_MAX)
  90         {
  91           long double div9, div10, div11, div12, div13;
  92
  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) */
  97           div10 = div9 * div9;
  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)
 102             {
 103               x /= div13; /* exact */
 104               shift_exp += 8192;
 105             }
 106           if (ABS (x) >= div12)
 107             {
 108               x /= div12; /* exact */
 109               shift_exp += 4096;
 110             }
 111           if (ABS (x) >= div11)
 112             {
 113               x /= div11; /* exact */
 114               shift_exp += 2048;
 115             }
 116           if (ABS (x) >= div10)
 117             {
 118               x /= div10; /* exact */
 119               shift_exp += 1024;
 120             }
 121           /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024,
 122              therefore we have one extra exponent reduction step */
 123           if (ABS (x) >= div9)
 124             {
 125               x /= div9; /* exact */
 126               shift_exp += 512;
 127             }
 128         } /* Check overflow of double */
 129       else
 130         {
 131           long double div9, div10, div11;
 132
 133           div9 = (long double) (double) 7.4583407312002067432909653e-155;
 134           /* div9 = 2^(-2^9) */
 135           div10 = div9  * div9;  /* 2^(-2^10) */
 136           div11 = div10 * div10; /* 2^(-2^11) if extended precision */
 137           /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not
 138              overflow here */
 139           if (ABS(x) < div10 &&
 140               div11 != (long double) 0.0 &&
 141               div11 / div10 == div10) /* possible underflow */
 142             {
 143               long double div12, div13;
 144               /* After the divisions, any bit of x must be >= div10,
 145                  hence the possible division by div9. */
 146               div12 = div11 * div11; /* 2^(-2^12) */
 147               div13 = div12 * div12; /* 2^(-2^13) */
 148               if (ABS (x) <= div13)
 149                 {
 150                   x /= div13; /* exact */
 151                   shift_exp -= 8192;
 152                 }
 153               if (ABS (x) <= div12)
 154                 {
 155                   x /= div12; /* exact */
 156                   shift_exp -= 4096;
 157                 }
 158               if (ABS (x) <= div11)
 159                 {
 160                   x /= div11; /* exact */
 161                   shift_exp -= 2048;
 162                 }
 163               if (ABS (x) <= div10)
 164                 {
 165                   x /= div10; /* exact */
 166                   shift_exp -= 1024;
 167                 }
 168               if (ABS(x) <= div9)
 169                 {
 170                   x /= div9;  /* exact */
 171                   shift_exp -= 512;
 172                 }
 173             }
 174           else
 175             {
 176               inexact = mpfr_set_d (u, (double) x, GMP_RNDZ);
 177               MPFR_ASSERTD (inexact == 0);
 178               if (mpfr_add (t, t, u, GMP_RNDZ) != 0)
 179                 {
 180                   if (!mpfr_number_p (t))
 181                     break;
 182                   /* Inexact. This cannot happen unless the C implementation
 183                      "lies" on the precision or when long doubles are
 184                      implemented with FP expansions like under Mac OS X. */
 185                   if (MPFR_PREC (t) != MPFR_PREC (r) + 1)
 186                     {
 187                       /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX.
 188                          The precision MPFR_PREC (r) + 1 allows us to
 189                          deduce the rounding bit and the sticky bit. */
 190                       mpfr_set_prec (t, MPFR_PREC (r) + 1);
 191                       goto convert;
 192                     }
 193                   else
 194                     {
 195                       mp_limb_t *tp;
 196                       int rb_mask;
 197
 198                       /* Since mpfr_add was inexact, the sticky bit is 1. */
 199                       tp = MPFR_MANT (t);
 200                       rb_mask = MPFR_LIMB_ONE <<
 201                         (BITS_PER_MP_LIMB - 1 -
 202                          (MPFR_PREC (r) & (BITS_PER_MP_LIMB - 1)));
 203                       if (rnd_mode == GMP_RNDN)
 204                         rnd_mode = (*tp & rb_mask) ^ MPFR_IS_NEG (t) ?
 205                           GMP_RNDU : GMP_RNDD;
 206                       *tp |= rb_mask;
 207                       break;
 208                     }
 209                 }
 210               x -= (long double) mpfr_get_d1 (u); /* exact */
 211             }
 212         }
 213     }
 214   inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode);
 215   mpfr_clear (t);
 216   mpfr_clear (u);
 217
 218   MPFR_SAVE_EXPO_FREE (expo);
 219   return mpfr_check_range (r, inexact, rnd_mode);
 220
 221  nan:
 222   MPFR_SET_NAN(r);
 223   MPFR_RET_NAN;
 224 }
 225
 226 #else /* IEEE Extended Little Endian Code */
 227
 228 int
 229 mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode)
 230 {
 231   int inexact, i, k, cnt;
 232   mpfr_t tmp;
 233   mp_limb_t tmpmant[MPFR_LIMBS_PER_LONG_DOUBLE];
 234   mpfr_long_double_t x;
 235   mp_exp_t exp;
 236   int signd;
 237   MPFR_SAVE_EXPO_DECL (expo);
 238
 239   /* Check for NAN */
 240   if (MPFR_UNLIKELY (d != d))
 241     {
 242       MPFR_SET_NAN (r);
 243       MPFR_RET_NAN;
 244     }
 245   /* Check for INF */
 246   else if (MPFR_UNLIKELY (d > MPFR_LDBL_MAX))
 247     {
 248       MPFR_SET_INF (r);
 249       MPFR_SET_POS (r);
 250       return 0;
 251     }
 252   else if (MPFR_UNLIKELY (d < -MPFR_LDBL_MAX))
 253     {
 254       MPFR_SET_INF (r);
 255       MPFR_SET_NEG (r);
 256       return 0;
 257     }
 258   /* Check for ZERO */
 259   else if (MPFR_UNLIKELY (d == 0.0))
 260     {
 261       x.ld = d;
 262       MPFR_SET_ZERO (r);
 263       if (x.s.sign == 1)
 264         MPFR_SET_NEG(r);
 265       else
 266         MPFR_SET_POS(r);
 267       return 0;
 268     }
 269
 270   /* now d is neither 0, nor NaN nor Inf */
 271   MPFR_SAVE_EXPO_MARK (expo);
 272
 273   MPFR_MANT (tmp) = tmpmant;
 274   MPFR_PREC (tmp) = 64;
 275
 276   /* Extract sign */
 277   x.ld = d;
 278   signd = MPFR_SIGN_POS;
 279   if (x.ld < 0.0)
 280     {
 281       signd = MPFR_SIGN_NEG;
 282       x.ld = -x.ld;
 283     }
 284
 285   /* Extract mantissa */
 286 #if BITS_PER_MP_LIMB >= 64
 287   tmpmant[0] = ((mp_limb_t) x.s.manh << 32) | ((mp_limb_t) x.s.manl);
 288 #else
 289   tmpmant[0] = (mp_limb_t) x.s.manl;
 290   tmpmant[1] = (mp_limb_t) x.s.manh;
 291 #endif
 292
 293   /* Normalize mantissa */
 294   i = MPFR_LIMBS_PER_LONG_DOUBLE;
 295   MPN_NORMALIZE_NOT_ZERO (tmpmant, i);
 296   k = MPFR_LIMBS_PER_LONG_DOUBLE - i;
 297   count_leading_zeros (cnt, tmpmant[i - 1]);
 298   if (MPFR_LIKELY (cnt != 0))
 299     mpn_lshift (tmpmant + k, tmpmant, i, cnt);
 300   else if (k != 0)
 301     MPN_COPY (tmpmant + k, tmpmant, i);
 302   if (MPFR_UNLIKELY (k != 0))
 303     MPN_ZERO (tmpmant, k);
 304
 305   /* Set exponent */
 306   exp = (mp_exp_t) ((x.s.exph << 8) + x.s.expl);  /* 15-bit unsigned int */
 307   if (MPFR_UNLIKELY (exp == 0))
 308     exp -= 0x3FFD;
 309   else
 310     exp -= 0x3FFE;
 311
 312   MPFR_SET_EXP (tmp, exp - cnt - k * BITS_PER_MP_LIMB);
 313
 314   /* tmp is exact */
 315   inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
 316
 317   MPFR_SAVE_EXPO_FREE (expo);
 318   return mpfr_check_range (r, inexact, rnd_mode);
 319 }
 320
 321 #endif