contrib/mpfr/get_d.c

   1 /* mpfr_get_d, mpfr_get_d_2exp -- convert a multiple precision floating-point
   2                                   number to a machine double precision float
   3
   4 Copyright 1999, 2000, 2001, 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 /* "double" NaN and infinities are written as explicit bytes to be sure of
  30    getting what we want, and to be sure of not depending on libm.
  31
  32    Could use 4-byte "float" values and let the code convert them, but it
  33    seems more direct to give exactly what we want.  Certainly for gcc 3.0.2
  34    on alphaev56-unknown-freebsd4.3 the NaN must be 8-bytes, since that
  35    compiler+system was seen incorrectly converting from a "float" NaN.  */
  36
  37 #if _GMP_IEEE_FLOATS
  38
  39 /* The "d" field guarantees alignment to a suitable boundary for a double.
  40    Could use a union instead, if we checked the compiler supports union
  41    initializers.  */
  42 struct dbl_bytes {
  43   unsigned char b[8];
  44   double d;
  45 };
  46
  47 #define MPFR_DBL_INFP  (* (const double *) dbl_infp.b)
  48 #define MPFR_DBL_INFM  (* (const double *) dbl_infm.b)
  49 #define MPFR_DBL_NAN   (* (const double *) dbl_nan.b)
  50
  51 #if HAVE_DOUBLE_IEEE_LITTLE_ENDIAN
  52 static const struct dbl_bytes dbl_infp =
  53   { { 0, 0, 0, 0, 0, 0, 0xF0, 0x7F }, 0.0 };
  54 static const struct dbl_bytes dbl_infm =
  55   { { 0, 0, 0, 0, 0, 0, 0xF0, 0xFF }, 0.0 };
  56 static const struct dbl_bytes dbl_nan  =
  57   { { 0, 0, 0, 0, 0, 0, 0xF8, 0x7F }, 0.0 };
  58 #endif
  59 #if HAVE_DOUBLE_IEEE_LITTLE_SWAPPED
  60 static const struct dbl_bytes dbl_infp =
  61   { { 0, 0, 0xF0, 0x7F, 0, 0, 0, 0 }, 0.0 };
  62 static const struct dbl_bytes dbl_infm =
  63   { { 0, 0, 0xF0, 0xFF, 0, 0, 0, 0 }, 0.0 };
  64 static const struct dbl_bytes dbl_nan  =
  65   { { 0, 0, 0xF8, 0x7F, 0, 0, 0, 0 }, 0.0 };
  66 #endif
  67 #if HAVE_DOUBLE_IEEE_BIG_ENDIAN
  68 static const struct dbl_bytes dbl_infp =
  69   { { 0x7F, 0xF0, 0, 0, 0, 0, 0, 0 }, 0.0 };
  70 static const struct dbl_bytes dbl_infm =
  71   { { 0xFF, 0xF0, 0, 0, 0, 0, 0, 0 }, 0.0 };
  72 static const struct dbl_bytes dbl_nan  =
  73   { { 0x7F, 0xF8, 0, 0, 0, 0, 0, 0 }, 0.0 };
  74 #endif
  75
  76 #else /* _GMP_IEEE_FLOATS */
  77
  78 #define MPFR_DBL_INFP DBL_POS_INF
  79 #define MPFR_DBL_INFM DBL_NEG_INF
  80 #define MPFR_DBL_NAN DBL_NAN
  81
  82 #endif /* _GMP_IEEE_FLOATS */
  83
  84
  85 /* multiplies 1/2 <= d <= 1 by 2^exp */
  86 static double
  87 mpfr_scale2 (double d, int exp)
  88 {
  89 #if _GMP_IEEE_FLOATS
  90   {
  91     union ieee_double_extract x;
  92
  93     if (MPFR_UNLIKELY (d == 1.0))
  94       {
  95         d = 0.5;
  96         exp ++;
  97       }
  98
  99     /* now 1/2 <= d < 1 */
 100
 101     /* infinities and zeroes have already been checked */
 102     MPFR_ASSERTD (-1073 <= exp && exp <= 1025);
 103
 104     x.d = d;
 105     if (MPFR_UNLIKELY (exp < -1021)) /* subnormal case */
 106       {
 107         x.s.exp += exp + 52;
 108         x.d *= DBL_EPSILON;
 109       }
 110     else /* normalized case */
 111       {
 112         x.s.exp += exp;
 113       }
 114     return x.d;
 115   }
 116 #else /* _GMP_IEEE_FLOATS */
 117   {
 118     double factor;
 119
 120     /* An overflow may occurs (example: 0.5*2^1024) */
 121     if (d < 1.0)
 122       {
 123         d += d;
 124         exp--;
 125       }
 126     /* Now 1.0 <= d < 2.0 */
 127
 128     if (exp < 0)
 129       {
 130         factor = 0.5;
 131         exp = -exp;
 132       }
 133     else
 134       {
 135         factor = 2.0;
 136       }
 137     while (exp != 0)
 138       {
 139         if ((exp & 1) != 0)
 140           d *= factor;
 141         exp >>= 1;
 142         factor *= factor;
 143       }
 144     return d;
 145   }
 146 #endif
 147 }
 148
 149 /* Assumes IEEE-754 double precision; otherwise, only an approximated
 150    result will be returned, without any guaranty (and special cases
 151    such as NaN must be avoided if not supported). */
 152
 153 double
 154 mpfr_get_d (mpfr_srcptr src, mp_rnd_t rnd_mode)
 155 {
 156   double d;
 157   int negative;
 158   mp_exp_t e;
 159
 160   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
 161     {
 162       if (MPFR_IS_NAN (src))
 163         return MPFR_DBL_NAN;
 164
 165       negative = MPFR_IS_NEG (src);
 166
 167       if (MPFR_IS_INF (src))
 168         return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
 169
 170       MPFR_ASSERTD (MPFR_IS_ZERO(src));
 171       return negative ? DBL_NEG_ZERO : 0.0;
 172     }
 173
 174   e = MPFR_GET_EXP (src);
 175   negative = MPFR_IS_NEG (src);
 176
 177   /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
 178      subnormal is 2^(-1074)=0.1e-1073 */
 179   if (MPFR_UNLIKELY (e < -1073))
 180     {
 181       /* Note: Avoid using a constant expression DBL_MIN * DBL_EPSILON
 182          as this gives 0 instead of the correct result with gcc on some
 183          Alpha machines. */
 184       d = negative ?
 185         (rnd_mode == GMP_RNDD ||
 186          (rnd_mode == GMP_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
 187          ? -DBL_MIN : DBL_NEG_ZERO) :
 188         (rnd_mode == GMP_RNDU ||
 189          (rnd_mode == GMP_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
 190          ? DBL_MIN : 0.0);
 191       if (d != 0.0)
 192         d *= DBL_EPSILON;
 193     }
 194   /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
 195   else if (MPFR_UNLIKELY (e > 1024))
 196     {
 197       d = negative ?
 198         (rnd_mode == GMP_RNDZ || rnd_mode == GMP_RNDU ?
 199          -DBL_MAX : MPFR_DBL_INFM) :
 200         (rnd_mode == GMP_RNDZ || rnd_mode == GMP_RNDD ?
 201          DBL_MAX : MPFR_DBL_INFP);
 202     }
 203   else
 204     {
 205       int nbits;
 206       mp_size_t np, i;
 207       mp_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
 208       int carry;
 209
 210       nbits = IEEE_DBL_MANT_DIG; /* 53 */
 211       if (MPFR_UNLIKELY (e < -1021))
 212         /*In the subnormal case, compute the exact number of significant bits*/
 213         {
 214           nbits += (1021 + e);
 215           MPFR_ASSERTD (nbits >= 1);
 216         }
 217       np = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
 218       MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
 219       carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative,
 220                                 nbits, rnd_mode);
 221       if (MPFR_UNLIKELY(carry))
 222         d = 1.0;
 223       else
 224         {
 225           /* The following computations are exact thanks to the previous
 226              mpfr_round_raw. */
 227           d = (double) tp[0] / MP_BASE_AS_DOUBLE;
 228           for (i = 1 ; i < np ; i++)
 229             d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
 230           /* d is the mantissa (between 1/2 and 1) of the argument rounded
 231              to 53 bits */
 232         }
 233       d = mpfr_scale2 (d, e);
 234       if (negative)
 235         d = -d;
 236     }
 237
 238   return d;
 239 }
 240
 241 #undef mpfr_get_d1
 242 double
 243 mpfr_get_d1 (mpfr_srcptr src)
 244 {
 245   return mpfr_get_d (src, __gmpfr_default_rounding_mode);
 246 }
 247
 248 double
 249 mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mp_rnd_t rnd_mode)
 250 {
 251   double ret;
 252   mp_exp_t exp;
 253   mpfr_t tmp;
 254
 255   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
 256     {
 257       int negative;
 258       *expptr = 0;
 259       if (MPFR_IS_NAN (src))
 260         return MPFR_DBL_NAN;
 261       negative = MPFR_IS_NEG (src);
 262       if (MPFR_IS_INF (src))
 263         return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
 264       MPFR_ASSERTD (MPFR_IS_ZERO(src));
 265       return negative ? DBL_NEG_ZERO : 0.0;
 266     }
 267
 268   tmp[0] = *src;        /* Hack copy mpfr_t */
 269   MPFR_SET_EXP (tmp, 0);
 270   ret = mpfr_get_d (tmp, rnd_mode);
 271
 272   if (MPFR_IS_PURE_FP(src))
 273     {
 274       exp = MPFR_GET_EXP (src);
 275
 276       /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
 277       if (ret == 1.0)
 278         {
 279           ret = 0.5;
 280           exp++;
 281         }
 282       else if (ret == -1.0)
 283         {
 284           ret = -0.5;
 285           exp++;
 286         }
 287
 288       MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
 289                     || (ret <= -0.5 && ret > -1.0));
 290       MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
 291     }
 292   else
 293     exp = 0;
 294
 295   *expptr = exp;
 296   return ret;
 297 }