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-2020 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 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_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ]; 93 int carry; 94 95 nbits = IEEE_DBL_MANT_DIG; /* 53 */ 96 if (MPFR_UNLIKELY (e < -1021)) 97 /*In the subnormal case, compute the exact number of significant bits*/ 98 { 99 nbits += 1021 + e; 100 MPFR_ASSERTD (1 <= nbits && nbits < IEEE_DBL_MANT_DIG); 101 } 102 carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative, 103 nbits, rnd_mode); 104 if (MPFR_UNLIKELY(carry)) 105 d = 1.0; 106 else 107 { 108 #if MPFR_LIMBS_PER_DOUBLE == 1 109 d = (double) tp[0] / MP_BASE_AS_DOUBLE; 110 #else 111 mp_size_t np, i; 112 MPFR_ASSERTD (nbits <= IEEE_DBL_MANT_DIG); 113 np = MPFR_PREC2LIMBS (nbits); 114 MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE ); 115 /* The following computations are exact thanks to the previous 116 mpfr_round_raw. */ 117 d = (double) tp[0] / MP_BASE_AS_DOUBLE; 118 for (i = 1 ; i < np ; i++) 119 d = (d + tp[i]) / MP_BASE_AS_DOUBLE; 120 /* d is the mantissa (between 1/2 and 1) of the argument rounded 121 to 53 bits */ 122 #endif 123 } 124 d = mpfr_scale2 (d, e); 125 if (negative) 126 d = -d; 127 } 128 129 return d; 130 } 131 132 #undef mpfr_get_d1 133 double 134 mpfr_get_d1 (mpfr_srcptr src) 135 { 136 return mpfr_get_d (src, __gmpfr_default_rounding_mode); 137 } 138 139 double 140 mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode) 141 { 142 double ret; 143 mpfr_exp_t exp; 144 mpfr_t tmp; 145 146 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src))) 147 { 148 int negative; 149 *expptr = 0; 150 if (MPFR_IS_NAN (src)) 151 return MPFR_DBL_NAN; 152 negative = MPFR_IS_NEG (src); 153 if (MPFR_IS_INF (src)) 154 return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP; 155 MPFR_ASSERTD (MPFR_IS_ZERO(src)); 156 return negative ? DBL_NEG_ZERO : 0.0; 157 } 158 159 MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0); 160 ret = mpfr_get_d (tmp, rnd_mode); 161 162 exp = MPFR_GET_EXP (src); 163 164 /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */ 165 if (ret == 1.0) 166 { 167 ret = 0.5; 168 exp++; 169 } 170 else if (ret == -1.0) 171 { 172 ret = -0.5; 173 exp++; 174 } 175 176 MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0) 177 || (ret <= -0.5 && ret > -1.0)); 178 MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX); 179 180 *expptr = exp; 181 return ret; 182 } 183