1 /* mpfr_set_d -- convert a machine double precision float to 2 a multiple precision floating-point number 3 4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 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> /* For DOUBLE_ISINF and DOUBLE_ISNAN */ 25 26 #define MPFR_NEED_LONGLONG_H 27 #include "mpfr-impl.h" 28 29 /* extracts the bits of d in rp[0..n-1] where n=ceil(53/GMP_NUMB_BITS). 30 Assumes d is neither 0 nor NaN nor Inf. */ 31 static long 32 __gmpfr_extract_double (mpfr_limb_ptr rp, double d) 33 /* e=0 iff GMP_NUMB_BITS=32 and rp has only one limb */ 34 { 35 long exp; 36 mp_limb_t manl; 37 #if GMP_NUMB_BITS == 32 38 mp_limb_t manh; 39 #endif 40 41 /* BUGS 42 1. Should handle Inf and NaN in IEEE specific code. 43 2. Handle Inf and NaN also in default code, to avoid hangs. 44 3. Generalize to handle all GMP_NUMB_BITS. 45 4. This lits is incomplete and misspelled. 46 */ 47 48 MPFR_ASSERTD(!DOUBLE_ISNAN(d)); 49 MPFR_ASSERTD(!DOUBLE_ISINF(d)); 50 MPFR_ASSERTD(d != 0.0); 51 52 #if _GMP_IEEE_FLOATS 53 54 { 55 union ieee_double_extract x; 56 x.d = d; 57 58 exp = x.s.exp; 59 if (exp) 60 { 61 #if GMP_NUMB_BITS >= 64 62 manl = ((MPFR_LIMB_ONE << 63) 63 | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11)); 64 #else 65 manh = (MPFR_LIMB_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21); 66 manl = x.s.manl << 11; 67 #endif 68 } 69 else /* denormalized number */ 70 { 71 #if GMP_NUMB_BITS >= 64 72 manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11); 73 #else 74 manh = (x.s.manh << 11) /* high 21 bits */ 75 | (x.s.manl >> 21); /* middle 11 bits */ 76 manl = x.s.manl << 11; /* low 21 bits */ 77 #endif 78 } 79 80 if (exp) 81 exp -= 1022; 82 else 83 exp = -1021; 84 } 85 86 #else /* _GMP_IEEE_FLOATS */ 87 88 { 89 /* Unknown (or known to be non-IEEE) double format. */ 90 exp = 0; 91 if (d >= 1.0) 92 { 93 MPFR_ASSERTN (d * 0.5 != d); 94 while (d >= 32768.0) 95 { 96 d *= (1.0 / 65536.0); 97 exp += 16; 98 } 99 while (d >= 1.0) 100 { 101 d *= 0.5; 102 exp += 1; 103 } 104 } 105 else if (d < 0.5) 106 { 107 while (d < (1.0 / 65536.0)) 108 { 109 d *= 65536.0; 110 exp -= 16; 111 } 112 while (d < 0.5) 113 { 114 d *= 2.0; 115 exp -= 1; 116 } 117 } 118 119 d *= MP_BASE_AS_DOUBLE; 120 #if GMP_NUMB_BITS >= 64 121 manl = d; 122 #else 123 manh = (mp_limb_t) d; 124 manl = (mp_limb_t) ((d - manh) * MP_BASE_AS_DOUBLE); 125 #endif 126 } 127 128 #endif /* _GMP_IEEE_FLOATS */ 129 130 #if GMP_NUMB_BITS >= 64 131 rp[0] = manl; 132 #else 133 rp[1] = manh; 134 rp[0] = manl; 135 #endif 136 137 return exp; 138 } 139 140 /* End of part included from gmp-2.0.2 */ 141 142 int 143 mpfr_set_d (mpfr_ptr r, double d, mpfr_rnd_t rnd_mode) 144 { 145 int signd, inexact; 146 unsigned int cnt; 147 mp_size_t i, k; 148 mpfr_t tmp; 149 mp_limb_t tmpmant[MPFR_LIMBS_PER_DOUBLE]; 150 MPFR_SAVE_EXPO_DECL (expo); 151 152 if (MPFR_UNLIKELY(DOUBLE_ISNAN(d))) 153 { 154 MPFR_SET_NAN(r); 155 MPFR_RET_NAN; 156 } 157 else if (MPFR_UNLIKELY(d == 0)) 158 { 159 #if _GMP_IEEE_FLOATS 160 union ieee_double_extract x; 161 162 MPFR_SET_ZERO(r); 163 /* set correct sign */ 164 x.d = d; 165 if (x.s.sig == 1) 166 MPFR_SET_NEG(r); 167 else 168 MPFR_SET_POS(r); 169 #else /* _GMP_IEEE_FLOATS */ 170 MPFR_SET_ZERO(r); 171 { 172 /* This is to get the sign of zero on non-IEEE hardware 173 Some systems support +0.0, -0.0 and unsigned zero. 174 We can't use d==+0.0 since it should be always true, 175 so we check that the memory representation of d is the 176 same than +0.0. etc */ 177 /* FIXME: consider the case where +0.0 or -0.0 may have several 178 representations. */ 179 double poszero = +0.0, negzero = DBL_NEG_ZERO; 180 if (memcmp(&d, &poszero, sizeof(double)) == 0) 181 MPFR_SET_POS(r); 182 else if (memcmp(&d, &negzero, sizeof(double)) == 0) 183 MPFR_SET_NEG(r); 184 else 185 MPFR_SET_POS(r); 186 } 187 #endif 188 return 0; /* 0 is exact */ 189 } 190 else if (MPFR_UNLIKELY(DOUBLE_ISINF(d))) 191 { 192 MPFR_SET_INF(r); 193 if (d > 0) 194 MPFR_SET_POS(r); 195 else 196 MPFR_SET_NEG(r); 197 return 0; /* infinity is exact */ 198 } 199 200 /* now d is neither 0, nor NaN nor Inf */ 201 202 MPFR_SAVE_EXPO_MARK (expo); 203 204 /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE, 205 since PREC(r) may be different from PREC(tmp), and then both variables 206 would have same precision in the mpfr_set4 call below. */ 207 MPFR_MANT(tmp) = tmpmant; 208 MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG; 209 210 signd = (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS; 211 d = ABS (d); 212 213 /* don't use MPFR_SET_EXP here since the exponent may be out of range */ 214 MPFR_EXP(tmp) = __gmpfr_extract_double (tmpmant, d); 215 216 #ifdef WANT_ASSERT 217 /* Failed assertion if the stored value is 0 (e.g., if the exponent range 218 has been reduced at the wrong moment and an underflow to 0 occurred). 219 Probably a bug in the C implementation if this happens. */ 220 i = 0; 221 while (tmpmant[i] == 0) 222 { 223 i++; 224 MPFR_ASSERTN(i < MPFR_LIMBS_PER_DOUBLE); 225 } 226 #endif 227 228 /* determine the index i-1 of the most significant non-zero limb 229 and the number k of zero high limbs */ 230 i = MPFR_LIMBS_PER_DOUBLE; 231 MPN_NORMALIZE_NOT_ZERO(tmpmant, i); 232 k = MPFR_LIMBS_PER_DOUBLE - i; 233 234 count_leading_zeros (cnt, tmpmant[i - 1]); 235 236 if (MPFR_LIKELY(cnt != 0)) 237 mpn_lshift (tmpmant + k, tmpmant, i, cnt); 238 else if (k != 0) 239 MPN_COPY (tmpmant + k, tmpmant, i); 240 241 if (MPFR_UNLIKELY(k != 0)) 242 MPN_ZERO (tmpmant, k); 243 244 /* don't use MPFR_SET_EXP here since the exponent may be out of range */ 245 MPFR_EXP(tmp) -= (mpfr_exp_t) (cnt + k * GMP_NUMB_BITS); 246 247 /* tmp is exact since PREC(tmp)=53 */ 248 inexact = mpfr_set4 (r, tmp, rnd_mode, signd); 249 250 MPFR_SAVE_EXPO_FREE (expo); 251 return mpfr_check_range (r, inexact, rnd_mode); 252 } 253 254 255 256