contrib/mpfr/src/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, 2010, 2011 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Caramel 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 3 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.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. */
  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
  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.  */
  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, mpfr_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               mpfr_div_2si (t, t, 8192, MPFR_RNDZ);
 106             }
 107           if (ABS (x) >= div12)
 108             {
 109               x /= div12; /* exact */
 110               shift_exp += 4096;
 111               mpfr_div_2si (t, t, 4096, MPFR_RNDZ);
 112             }
 113           if (ABS (x) >= div11)
 114             {
 115               x /= div11; /* exact */
 116               shift_exp += 2048;
 117               mpfr_div_2si (t, t, 2048, MPFR_RNDZ);
 118             }
 119           if (ABS (x) >= div10)
 120             {
 121               x /= div10; /* exact */
 122               shift_exp += 1024;
 123               mpfr_div_2si (t, t, 1024, MPFR_RNDZ);
 124             }
 125           /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024,
 126              therefore we have one extra exponent reduction step */
 127           if (ABS (x) >= div9)
 128             {
 129               x /= div9; /* exact */
 130               shift_exp += 512;
 131               mpfr_div_2si (t, t, 512, MPFR_RNDZ);
 132             }
 133         } /* Check overflow of double */
 134       else /* no overflow on double */
 135         {
 136           long double div9, div10, div11;
 137
 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
 143              overflow here */
 144           if (ABS(x) < div10 &&
 145               div11 != (long double) 0.0 &&
 146               div11 / div10 == div10) /* possible underflow */
 147             {
 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)
 154                 {
 155                   x /= div13; /* exact */
 156                   shift_exp -= 8192;
 157                   mpfr_mul_2si (t, t, 8192, MPFR_RNDZ);
 158                 }
 159               if (ABS (x) <= div12)
 160                 {
 161                   x /= div12; /* exact */
 162                   shift_exp -= 4096;
 163                   mpfr_mul_2si (t, t, 4096, MPFR_RNDZ);
 164                 }
 165               if (ABS (x) <= div11)
 166                 {
 167                   x /= div11; /* exact */
 168                   shift_exp -= 2048;
 169                   mpfr_mul_2si (t, t, 2048, MPFR_RNDZ);
 170                 }
 171               if (ABS (x) <= div10)
 172                 {
 173                   x /= div10; /* exact */
 174                   shift_exp -= 1024;
 175                   mpfr_mul_2si (t, t, 1024, MPFR_RNDZ);
 176                 }
 177               if (ABS(x) <= div9)
 178                 {
 179                   x /= div9;  /* exact */
 180                   shift_exp -= 512;
 181                   mpfr_mul_2si (t, t, 512, MPFR_RNDZ);
 182                 }
 183             }
 184           else /* no underflow */
 185             {
 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)
 189                 {
 190                   if (!mpfr_number_p (t))
 191                     break;
 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)
 196                     {
 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);
 201                       goto convert;
 202                     }
 203                   else
 204                     {
 205                       mp_limb_t *tp;
 206                       int rb_mask;
 207
 208                       /* Since mpfr_add was inexact, the sticky bit is 1. */
 209                       tp = MPFR_MANT (t);
 210                       rb_mask = MPFR_LIMB_ONE <<
 211                         (GMP_NUMB_BITS - 1 -
 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;
 216                       *tp |= rb_mask;
 217                       break;
 218                     }
 219                 }
 220               x -= (long double) mpfr_get_d1 (u); /* exact */
 221             }
 222         }
 223     }
 224   inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode);
 225   mpfr_clear (t);
 226   mpfr_clear (u);
 227
 228   MPFR_SAVE_EXPO_FREE (expo);
 229   return mpfr_check_range (r, inexact, rnd_mode);
 230
 231  nan:
 232   MPFR_SET_NAN(r);
 233   MPFR_RET_NAN;
 234 }
 235
 236 #else /* IEEE Extended Little Endian Code */
 237
 238 int
 239 mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode)
 240 {
 241   int inexact, i, k, cnt;
 242   mpfr_t tmp;
 243   mp_limb_t tmpmant[MPFR_LIMBS_PER_LONG_DOUBLE];
 244   mpfr_long_double_t x;
 245   mpfr_exp_t exp;
 246   int signd;
 247   MPFR_SAVE_EXPO_DECL (expo);
 248
 249   /* Check for NAN */
 250   if (MPFR_UNLIKELY (d != d))
 251     {
 252       MPFR_SET_NAN (r);
 253       MPFR_RET_NAN;
 254     }
 255   /* Check for INF */
 256   else if (MPFR_UNLIKELY (d > MPFR_LDBL_MAX))
 257     {
 258       MPFR_SET_INF (r);
 259       MPFR_SET_POS (r);
 260       return 0;
 261     }
 262   else if (MPFR_UNLIKELY (d < -MPFR_LDBL_MAX))
 263     {
 264       MPFR_SET_INF (r);
 265       MPFR_SET_NEG (r);
 266       return 0;
 267     }
 268   /* Check for ZERO */
 269   else if (MPFR_UNLIKELY (d == 0.0))
 270     {
 271       x.ld = d;
 272       MPFR_SET_ZERO (r);
 273       if (x.s.sign == 1)
 274         MPFR_SET_NEG(r);
 275       else
 276         MPFR_SET_POS(r);
 277       return 0;
 278     }
 279
 280   /* now d is neither 0, nor NaN nor Inf */
 281   MPFR_SAVE_EXPO_MARK (expo);
 282
 283   MPFR_MANT (tmp) = tmpmant;
 284   MPFR_PREC (tmp) = 64;
 285
 286   /* Extract sign */
 287   x.ld = d;
 288   signd = MPFR_SIGN_POS;
 289   if (x.ld < 0.0)
 290     {
 291       signd = MPFR_SIGN_NEG;
 292       x.ld = -x.ld;
 293     }
 294
 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);
 298 #else
 299   tmpmant[0] = (mp_limb_t) x.s.manl;
 300   tmpmant[1] = (mp_limb_t) x.s.manh;
 301 #endif
 302
 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);
 310   else if (k != 0)
 311     MPN_COPY (tmpmant + k, tmpmant, i);
 312   if (MPFR_UNLIKELY (k != 0))
 313     MPN_ZERO (tmpmant, k);
 314
 315   /* Set exponent */
 316   exp = (mpfr_exp_t) ((x.s.exph << 8) + x.s.expl);  /* 15-bit unsigned int */
 317   if (MPFR_UNLIKELY (exp == 0))
 318     exp -= 0x3FFD;
 319   else
 320     exp -= 0x3FFE;
 321
 322   MPFR_SET_EXP (tmp, exp - cnt - k * GMP_NUMB_BITS);
 323
 324   /* tmp is exact */
 325   inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
 326
 327   MPFR_SAVE_EXPO_FREE (expo);
 328   return mpfr_check_range (r, inexact, rnd_mode);
 329 }
 330
 331 #endif