contrib/mpfr/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 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>  /* 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/BITS_PER_MP_LIMB).
  30    Assumes d is neither 0 nor NaN nor Inf. */
  31 static long
  32 __gmpfr_extract_double (mp_ptr rp, double d)
  33      /* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */
  34 {
  35   long exp;
  36   mp_limb_t manl;
  37 #if BITS_PER_MP_LIMB == 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 BITS_PER_MP_LIMB.
  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 BITS_PER_MP_LIMB >= 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 BITS_PER_MP_LIMB >= 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 BITS_PER_MP_LIMB >= 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 BITS_PER_MP_LIMB >= 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, mp_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   MPFR_CLEAR_FLAGS(r);
 153
 154   if (MPFR_UNLIKELY(DOUBLE_ISNAN(d)))
 155     {
 156       MPFR_SET_NAN(r);
 157       MPFR_RET_NAN;
 158     }
 159   else if (MPFR_UNLIKELY(d == 0))
 160     {
 161 #if _GMP_IEEE_FLOATS
 162       union ieee_double_extract x;
 163
 164       MPFR_SET_ZERO(r);
 165       /* set correct sign */
 166       x.d = d;
 167       if (x.s.sig == 1)
 168         MPFR_SET_NEG(r);
 169       else
 170         MPFR_SET_POS(r);
 171 #else /* _GMP_IEEE_FLOATS */
 172       MPFR_SET_ZERO(r);
 173       {
 174         /* This is to get the sign of zero on non-IEEE hardware
 175            Some systems support +0.0, -0.0 and unsigned zero.
 176            We can't use d==+0.0 since it should be always true,
 177            so we check that the memory representation of d is the
 178            same than +0.0. etc */
 179         /* FIXME: consider the case where +0.0 or -0.0 may have several
 180            representations. */
 181         double poszero = +0.0, negzero = DBL_NEG_ZERO;
 182         if (memcmp(&d, &poszero, sizeof(double)) == 0)
 183           MPFR_SET_POS(r);
 184         else if (memcmp(&d, &negzero, sizeof(double)) == 0)
 185           MPFR_SET_NEG(r);
 186         else
 187           MPFR_SET_POS(r);
 188       }
 189 #endif
 190       return 0; /* 0 is exact */
 191     }
 192   else if (MPFR_UNLIKELY(DOUBLE_ISINF(d)))
 193     {
 194       MPFR_SET_INF(r);
 195       if (d > 0)
 196         MPFR_SET_POS(r);
 197       else
 198         MPFR_SET_NEG(r);
 199       return 0; /* infinity is exact */
 200     }
 201
 202   /* now d is neither 0, nor NaN nor Inf */
 203
 204   MPFR_SAVE_EXPO_MARK (expo);
 205
 206   /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE,
 207      since PREC(r) may be different from PREC(tmp), and then both variables
 208      would have same precision in the mpfr_set4 call below. */
 209   MPFR_MANT(tmp) = tmpmant;
 210   MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG;
 211
 212   signd = (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS;
 213   d = ABS (d);
 214
 215   /* don't use MPFR_SET_EXP here since the exponent may be out of range */
 216   MPFR_EXP(tmp) = __gmpfr_extract_double (tmpmant, d);
 217
 218 #ifdef WANT_ASSERT
 219   /* Failed assertion if the stored value is 0 (e.g., if the exponent range
 220      has been reduced at the wrong moment and an underflow to 0 occurred).
 221      Probably a bug in the C implementation if this happens. */
 222   i = 0;
 223   while (tmpmant[i] == 0)
 224     {
 225       i++;
 226       MPFR_ASSERTN(i < MPFR_LIMBS_PER_DOUBLE);
 227     }
 228 #endif
 229
 230   /* determine the index i-1 of the most significant non-zero limb
 231      and the number k of zero high limbs */
 232   i = MPFR_LIMBS_PER_DOUBLE;
 233   MPN_NORMALIZE_NOT_ZERO(tmpmant, i);
 234   k = MPFR_LIMBS_PER_DOUBLE - i;
 235
 236   count_leading_zeros (cnt, tmpmant[i - 1]);
 237
 238   if (MPFR_LIKELY(cnt != 0))
 239     mpn_lshift (tmpmant + k, tmpmant, i, cnt);
 240   else if (k != 0)
 241     MPN_COPY (tmpmant + k, tmpmant, i);
 242
 243   if (MPFR_UNLIKELY(k != 0))
 244     MPN_ZERO (tmpmant, k);
 245
 246   /* don't use MPFR_SET_EXP here since the exponent may be out of range */
 247   MPFR_EXP(tmp) -= (mp_exp_t) (cnt + k * BITS_PER_MP_LIMB);
 248
 249   /* tmp is exact since PREC(tmp)=53 */
 250   inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
 251
 252   MPFR_SAVE_EXPO_FREE (expo);
 253   return mpfr_check_range (r, inexact, rnd_mode);
 254 }
 255
 256
 257