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