Home | History | Annotate | Line # | Download | only in src 
      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-2023 Free Software Foundation, Inc.
      5 Contributed by the AriC and Caramba 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 https://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         {
     46           /* we don't propagate the sign bit */
     47           return MPFR_DBL_NAN;
     48         }
     49 
     50       negative = MPFR_IS_NEG (src);
     51 
     52       if (MPFR_IS_INF (src))
     53         return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
     54 
     55       MPFR_ASSERTD (MPFR_IS_ZERO(src));
     56       return negative ? DBL_NEG_ZERO : 0.0;
     57     }
     58 
     59   e = MPFR_GET_EXP (src);
     60   negative = MPFR_IS_NEG (src);
     61 
     62   if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
     63     rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;
     64 
     65   /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
     66      subnormal is 2^(-1074)=0.1e-1073 */
     67   if (MPFR_UNLIKELY (e < -1073))
     68     {
     69       /* Note: Avoid using a constant expression DBL_MIN * DBL_EPSILON
     70          as this gives 0 instead of the correct result with gcc on some
     71          Alpha machines. */
     72       d = negative ?
     73         (rnd_mode == MPFR_RNDD ||
     74          (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
     75          ? -DBL_MIN : DBL_NEG_ZERO) :
     76         (rnd_mode == MPFR_RNDU ||
     77          (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
     78          ? DBL_MIN : 0.0);
     79       if (d != 0.0) /* we multiply DBL_MIN = 2^(-1022) by DBL_EPSILON = 2^(-52)
     80                        to get +-2^(-1074) */
     81         d *= DBL_EPSILON;
     82     }
     83   /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
     84   else if (MPFR_UNLIKELY (e > 1024))
     85     {
     86       d = negative ?
     87         (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDU ?
     88          -DBL_MAX : MPFR_DBL_INFM) :
     89         (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD ?
     90          DBL_MAX : MPFR_DBL_INFP);
     91     }
     92   else
     93     {
     94       int nbits;
     95       mp_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
     96       int carry;
     97 
     98       nbits = IEEE_DBL_MANT_DIG; /* 53 */
     99       if (MPFR_UNLIKELY (e < -1021))
    100         /*In the subnormal case, compute the exact number of significant bits*/
    101         {
    102           nbits += 1021 + e;
    103           MPFR_ASSERTD (1 <= nbits && nbits < IEEE_DBL_MANT_DIG);
    104         }
    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 #if MPFR_LIMBS_PER_DOUBLE == 1
    112           d = (double) tp[0] / MP_BASE_AS_DOUBLE;
    113 #else
    114           mp_size_t np, i;
    115           MPFR_ASSERTD (nbits <= IEEE_DBL_MANT_DIG);
    116           np = MPFR_PREC2LIMBS (nbits);
    117           MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
    118           /* The following computations are exact thanks to the previous
    119              mpfr_round_raw. */
    120           d = (double) tp[0] / MP_BASE_AS_DOUBLE;
    121           for (i = 1 ; i < np ; i++)
    122             d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
    123           /* d is the mantissa (between 1/2 and 1) of the argument rounded
    124              to 53 bits */
    125 #endif
    126         }
    127       d = mpfr_scale2 (d, e);
    128       if (negative)
    129         d = -d;
    130     }
    131 
    132   return d;
    133 }
    134 
    135 #undef mpfr_get_d1
    136 double
    137 mpfr_get_d1 (mpfr_srcptr src)
    138 {
    139   return mpfr_get_d (src, __gmpfr_default_rounding_mode);
    140 }
    141 
    142 double
    143 mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
    144 {
    145   double ret;
    146   mpfr_exp_t exp;
    147   mpfr_t tmp;
    148 
    149   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
    150     {
    151       int negative;
    152       *expptr = 0;
    153       if (MPFR_IS_NAN (src))
    154         {
    155           /* we don't propagate the sign bit */
    156           return MPFR_DBL_NAN;
    157         }
    158       negative = MPFR_IS_NEG (src);
    159       if (MPFR_IS_INF (src))
    160         return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
    161       MPFR_ASSERTD (MPFR_IS_ZERO(src));
    162       return negative ? DBL_NEG_ZERO : 0.0;
    163     }
    164 
    165   MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0);
    166   ret = mpfr_get_d (tmp, rnd_mode);
    167 
    168   exp = MPFR_GET_EXP (src);
    169 
    170   /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
    171   if (ret == 1.0)
    172     {
    173       ret = 0.5;
    174       exp++;
    175     }
    176   else if (ret == -1.0)
    177     {
    178       ret = -0.5;
    179       exp++;
    180     }
    181 
    182   MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
    183                 || (ret <= -0.5 && ret > -1.0));
    184   MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
    185 
    186   *expptr = exp;
    187   return ret;
    188 }
    189