contrib/mpfr/src/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, 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>
  25
  26 #define MPFR_NEED_LONGLONG_H
  27 #include "mpfr-impl.h"
  28
  29 #include "ieee_floats.h"
  30
  31 /* Assumes IEEE-754 double precision; otherwise, only an approximated
  32    result will be returned, without any guaranty (and special cases
  33    such as NaN must be avoided if not supported). */
  34
  35 double
  36 mpfr_get_d (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
  37 {
  38   double d;
  39   int negative;
  40   mpfr_exp_t e;
  41
  42   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
  43     {
  44       if (MPFR_IS_NAN (src))
  45         return MPFR_DBL_NAN;
  46
  47       negative = MPFR_IS_NEG (src);
  48
  49       if (MPFR_IS_INF (src))
  50         return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
  51
  52       MPFR_ASSERTD (MPFR_IS_ZERO(src));
  53       return negative ? DBL_NEG_ZERO : 0.0;
  54     }
  55
  56   e = MPFR_GET_EXP (src);
  57   negative = MPFR_IS_NEG (src);
  58
  59   if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
  60     rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;
  61
  62   /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
  63      subnormal is 2^(-1074)=0.1e-1073 */
  64   if (MPFR_UNLIKELY (e < -1073))
  65     {
  66       /* Note: Avoid using a constant expression DBL_MIN * DBL_EPSILON
  67          as this gives 0 instead of the correct result with gcc on some
  68          Alpha machines. */
  69       d = negative ?
  70         (rnd_mode == MPFR_RNDD ||
  71          (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
  72          ? -DBL_MIN : DBL_NEG_ZERO) :
  73         (rnd_mode == MPFR_RNDU ||
  74          (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
  75          ? DBL_MIN : 0.0);
  76       if (d != 0.0) /* we multiply DBL_MIN = 2^(-1022) by DBL_EPSILON = 2^(-52)
  77                        to get +-2^(-1074) */
  78         d *= DBL_EPSILON;
  79     }
  80   /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
  81   else if (MPFR_UNLIKELY (e > 1024))
  82     {
  83       d = negative ?
  84         (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDU ?
  85          -DBL_MAX : MPFR_DBL_INFM) :
  86         (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD ?
  87          DBL_MAX : MPFR_DBL_INFP);
  88     }
  89   else
  90     {
  91       int nbits;
  92       mp_size_t np, i;
  93       mp_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
  94       int carry;
  95
  96       nbits = IEEE_DBL_MANT_DIG; /* 53 */
  97       if (MPFR_UNLIKELY (e < -1021))
  98         /*In the subnormal case, compute the exact number of significant bits*/
  99         {
 100           nbits += (1021 + e);
 101           MPFR_ASSERTD (nbits >= 1);
 102         }
 103       np = MPFR_PREC2LIMBS (nbits);
 104       MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
 105       carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative,
 106                                 nbits, rnd_mode);
 107       if (MPFR_UNLIKELY(carry))
 108         d = 1.0;
 109       else
 110         {
 111           /* The following computations are exact thanks to the previous
 112              mpfr_round_raw. */
 113           d = (double) tp[0] / MP_BASE_AS_DOUBLE;
 114           for (i = 1 ; i < np ; i++)
 115             d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
 116           /* d is the mantissa (between 1/2 and 1) of the argument rounded
 117              to 53 bits */
 118         }
 119       d = mpfr_scale2 (d, e);
 120       if (negative)
 121         d = -d;
 122     }
 123
 124   return d;
 125 }
 126
 127 #undef mpfr_get_d1
 128 double
 129 mpfr_get_d1 (mpfr_srcptr src)
 130 {
 131   return mpfr_get_d (src, __gmpfr_default_rounding_mode);
 132 }
 133
 134 double
 135 mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
 136 {
 137   double ret;
 138   mpfr_exp_t exp;
 139   mpfr_t tmp;
 140
 141   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
 142     {
 143       int negative;
 144       *expptr = 0;
 145       if (MPFR_IS_NAN (src))
 146         return MPFR_DBL_NAN;
 147       negative = MPFR_IS_NEG (src);
 148       if (MPFR_IS_INF (src))
 149         return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
 150       MPFR_ASSERTD (MPFR_IS_ZERO(src));
 151       return negative ? DBL_NEG_ZERO : 0.0;
 152     }
 153
 154   tmp[0] = *src;        /* Hack copy mpfr_t */
 155   MPFR_SET_EXP (tmp, 0);
 156   ret = mpfr_get_d (tmp, rnd_mode);
 157
 158   if (MPFR_IS_PURE_FP(src))
 159     {
 160       exp = MPFR_GET_EXP (src);
 161
 162       /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
 163       if (ret == 1.0)
 164         {
 165           ret = 0.5;
 166           exp++;
 167         }
 168       else if (ret == -1.0)
 169         {
 170           ret = -0.5;
 171           exp++;
 172         }
 173
 174       MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
 175                     || (ret <= -0.5 && ret > -1.0));
 176       MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
 177     }
 178   else
 179     exp = 0;
 180
 181   *expptr = exp;
 182   return ret;
 183 }