contrib/mpfr/src/set_d.c

   1 /* mpfr_set_d -- convert a machine double precision float to
   2                  a multiple precision floating-point number
   3
   4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   5 Contributed by the AriC 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>  /* For DOUBLE_ISINF and DOUBLE_ISNAN */
  25
  26 #define MPFR_NEED_LONGLONG_H
  27 #include "mpfr-impl.h"
  28
  29 /* extracts the bits of d in rp[0..n-1] where n=ceil(53/GMP_NUMB_BITS).
  30    Assumes d is neither 0 nor NaN nor Inf. */
  31 static long
  32 __gmpfr_extract_double (mpfr_limb_ptr rp, double d)
  33      /* e=0 iff GMP_NUMB_BITS=32 and rp has only one limb */
  34 {
  35   long exp;
  36   mp_limb_t manl;
  37 #if GMP_NUMB_BITS == 32
  38   mp_limb_t manh;
  39 #endif
  40
  41   /* BUGS
  42      1. Should handle Inf and NaN in IEEE specific code.
  43      2. Handle Inf and NaN also in default code, to avoid hangs.
  44      3. Generalize to handle all GMP_NUMB_BITS.
  45      4. This lits is incomplete and misspelled.
  46    */
  47
  48   MPFR_ASSERTD(!DOUBLE_ISNAN(d));
  49   MPFR_ASSERTD(!DOUBLE_ISINF(d));
  50   MPFR_ASSERTD(d != 0.0);
  51
  52 #if _GMP_IEEE_FLOATS
  53
  54   {
  55     union ieee_double_extract x;
  56     x.d = d;
  57
  58     exp = x.s.exp;
  59     if (exp)
  60       {
  61 #if GMP_NUMB_BITS >= 64
  62         manl = ((MPFR_LIMB_ONE << 63)
  63                 | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
  64 #else
  65         manh = (MPFR_LIMB_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
  66         manl = x.s.manl << 11;
  67 #endif
  68       }
  69     else /* denormalized number */
  70       {
  71 #if GMP_NUMB_BITS >= 64
  72         manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11);
  73 #else
  74         manh = (x.s.manh << 11) /* high 21 bits */
  75           | (x.s.manl >> 21); /* middle 11 bits */
  76         manl = x.s.manl << 11; /* low 21 bits */
  77 #endif
  78       }
  79
  80     if (exp)
  81       exp -= 1022;
  82     else
  83       exp = -1021;
  84   }
  85
  86 #else /* _GMP_IEEE_FLOATS */
  87
  88   {
  89     /* Unknown (or known to be non-IEEE) double format.  */
  90     exp = 0;
  91     if (d >= 1.0)
  92       {
  93         MPFR_ASSERTN (d * 0.5 != d);
  94         while (d >= 32768.0)
  95           {
  96             d *= (1.0 / 65536.0);
  97             exp += 16;
  98           }
  99         while (d >= 1.0)
 100           {
 101             d *= 0.5;
 102             exp += 1;
 103           }
 104       }
 105     else if (d < 0.5)
 106       {
 107         while (d < (1.0 / 65536.0))
 108           {
 109             d *=  65536.0;
 110             exp -= 16;
 111           }
 112         while (d < 0.5)
 113           {
 114             d *= 2.0;
 115             exp -= 1;
 116           }
 117       }
 118
 119     d *= MP_BASE_AS_DOUBLE;
 120 #if GMP_NUMB_BITS >= 64
 121     manl = d;
 122 #else
 123     manh = (mp_limb_t) d;
 124     manl = (mp_limb_t) ((d - manh) * MP_BASE_AS_DOUBLE);
 125 #endif
 126   }
 127
 128 #endif /* _GMP_IEEE_FLOATS */
 129
 130 #if GMP_NUMB_BITS >= 64
 131   rp[0] = manl;
 132 #else
 133   rp[1] = manh;
 134   rp[0] = manl;
 135 #endif
 136
 137   return exp;
 138 }
 139
 140 /* End of part included from gmp-2.0.2 */
 141
 142 int
 143 mpfr_set_d (mpfr_ptr r, double d, mpfr_rnd_t rnd_mode)
 144 {
 145   int signd, inexact;
 146   unsigned int cnt;
 147   mp_size_t i, k;
 148   mpfr_t tmp;
 149   mp_limb_t tmpmant[MPFR_LIMBS_PER_DOUBLE];
 150   MPFR_SAVE_EXPO_DECL (expo);
 151
 152   if (MPFR_UNLIKELY(DOUBLE_ISNAN(d)))
 153     {
 154       MPFR_SET_NAN(r);
 155       MPFR_RET_NAN;
 156     }
 157   else if (MPFR_UNLIKELY(d == 0))
 158     {
 159 #if _GMP_IEEE_FLOATS
 160       union ieee_double_extract x;
 161
 162       MPFR_SET_ZERO(r);
 163       /* set correct sign */
 164       x.d = d;
 165       if (x.s.sig == 1)
 166         MPFR_SET_NEG(r);
 167       else
 168         MPFR_SET_POS(r);
 169 #else /* _GMP_IEEE_FLOATS */
 170       MPFR_SET_ZERO(r);
 171       {
 172         /* This is to get the sign of zero on non-IEEE hardware
 173            Some systems support +0.0, -0.0 and unsigned zero.
 174            We can't use d==+0.0 since it should be always true,
 175            so we check that the memory representation of d is the
 176            same than +0.0. etc */
 177         /* FIXME: consider the case where +0.0 or -0.0 may have several
 178            representations. */
 179         double poszero = +0.0, negzero = DBL_NEG_ZERO;
 180         if (memcmp(&d, &poszero, sizeof(double)) == 0)
 181           MPFR_SET_POS(r);
 182         else if (memcmp(&d, &negzero, sizeof(double)) == 0)
 183           MPFR_SET_NEG(r);
 184         else
 185           MPFR_SET_POS(r);
 186       }
 187 #endif
 188       return 0; /* 0 is exact */
 189     }
 190   else if (MPFR_UNLIKELY(DOUBLE_ISINF(d)))
 191     {
 192       MPFR_SET_INF(r);
 193       if (d > 0)
 194         MPFR_SET_POS(r);
 195       else
 196         MPFR_SET_NEG(r);
 197       return 0; /* infinity is exact */
 198     }
 199
 200   /* now d is neither 0, nor NaN nor Inf */
 201
 202   MPFR_SAVE_EXPO_MARK (expo);
 203
 204   /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE,
 205      since PREC(r) may be different from PREC(tmp), and then both variables
 206      would have same precision in the mpfr_set4 call below. */
 207   MPFR_MANT(tmp) = tmpmant;
 208   MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG;
 209
 210   signd = (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS;
 211   d = ABS (d);
 212
 213   /* don't use MPFR_SET_EXP here since the exponent may be out of range */
 214   MPFR_EXP(tmp) = __gmpfr_extract_double (tmpmant, d);
 215
 216 #ifdef WANT_ASSERT
 217   /* Failed assertion if the stored value is 0 (e.g., if the exponent range
 218      has been reduced at the wrong moment and an underflow to 0 occurred).
 219      Probably a bug in the C implementation if this happens. */
 220   i = 0;
 221   while (tmpmant[i] == 0)
 222     {
 223       i++;
 224       MPFR_ASSERTN(i < MPFR_LIMBS_PER_DOUBLE);
 225     }
 226 #endif
 227
 228   /* determine the index i-1 of the most significant non-zero limb
 229      and the number k of zero high limbs */
 230   i = MPFR_LIMBS_PER_DOUBLE;
 231   MPN_NORMALIZE_NOT_ZERO(tmpmant, i);
 232   k = MPFR_LIMBS_PER_DOUBLE - i;
 233
 234   count_leading_zeros (cnt, tmpmant[i - 1]);
 235
 236   if (MPFR_LIKELY(cnt != 0))
 237     mpn_lshift (tmpmant + k, tmpmant, i, cnt);
 238   else if (k != 0)
 239     MPN_COPY (tmpmant + k, tmpmant, i);
 240
 241   if (MPFR_UNLIKELY(k != 0))
 242     MPN_ZERO (tmpmant, k);
 243
 244   /* don't use MPFR_SET_EXP here since the exponent may be out of range */
 245   MPFR_EXP(tmp) -= (mpfr_exp_t) (cnt + k * GMP_NUMB_BITS);
 246
 247   /* tmp is exact since PREC(tmp)=53 */
 248   inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
 249
 250   MPFR_SAVE_EXPO_FREE (expo);
 251   return mpfr_check_range (r, inexact, rnd_mode);
 252 }
 253
 254
 255