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 } 184