1 /* mpfr_set_d -- convert a machine double precision float to 2 a multiple precision floating-point number 3 4 Copyright 1999-2004, 2006-2018 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 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 finite and <> 0. 31 Returns the corresponding exponent such that |d| = {rp, n} * 2^exp, 32 with the value of {rp, n} in [1/2, 1). 33 The int type should be sufficient for exp. 34 */ 35 static int 36 extract_double (mpfr_limb_ptr rp, double d) 37 { 38 int exp; 39 mp_limb_t manl; 40 #if GMP_NUMB_BITS == 32 41 mp_limb_t manh; 42 #endif 43 44 /* FIXME: Generalize to handle all GMP_NUMB_BITS. */ 45 46 MPFR_ASSERTD(!DOUBLE_ISNAN(d)); 47 MPFR_ASSERTD(!DOUBLE_ISINF(d)); 48 MPFR_ASSERTD(d != 0.0); 49 50 #if _MPFR_IEEE_FLOATS 51 52 { 53 union mpfr_ieee_double_extract x; 54 x.d = d; 55 56 exp = x.s.exp; 57 if (exp) 58 { 59 #if GMP_NUMB_BITS >= 64 60 manl = ((MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1)) | 61 ((mp_limb_t) x.s.manh << (GMP_NUMB_BITS - 21)) | 62 ((mp_limb_t) x.s.manl << (GMP_NUMB_BITS - 53))); 63 #else 64 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 32); 65 manh = (MPFR_LIMB_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21); 66 manl = x.s.manl << 11; 67 #endif 68 exp -= 1022; 69 } 70 else /* subnormal number */ 71 { 72 int cnt; 73 exp = -1021; 74 #if GMP_NUMB_BITS >= 64 75 manl = (((mp_limb_t) x.s.manh << (GMP_NUMB_BITS - 21)) | 76 ((mp_limb_t) x.s.manl << (GMP_NUMB_BITS - 53))); 77 count_leading_zeros (cnt, manl); 78 #else 79 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 32); 80 manh = (x.s.manh << 11) /* high 21 bits */ 81 | (x.s.manl >> 21); /* middle 11 bits */ 82 manl = x.s.manl << 11; /* low 21 bits */ 83 if (manh == 0) 84 { 85 manh = manl; 86 manl = 0; 87 exp -= GMP_NUMB_BITS; 88 } 89 count_leading_zeros (cnt, manh); 90 manh = (manh << cnt) | 91 (cnt != 0 ? manl >> (GMP_NUMB_BITS - cnt) : 0); 92 #endif 93 manl <<= cnt; 94 exp -= cnt; 95 } 96 } 97 98 #else /* _MPFR_IEEE_FLOATS */ 99 100 { 101 /* Unknown (or known to be non-IEEE) double format. */ 102 exp = 0; 103 d = ABS (d); 104 if (d >= 1.0) 105 { 106 while (d >= 32768.0) 107 { 108 d *= (1.0 / 65536.0); 109 exp += 16; 110 } 111 while (d >= 1.0) 112 { 113 d *= 0.5; 114 exp += 1; 115 } 116 } 117 else if (d < 0.5) 118 { 119 while (d < (1.0 / 65536.0)) 120 { 121 d *= 65536.0; 122 exp -= 16; 123 } 124 while (d < 0.5) 125 { 126 d *= 2.0; 127 exp -= 1; 128 } 129 } 130 131 d *= MP_BASE_AS_DOUBLE; 132 #if GMP_NUMB_BITS >= 64 133 #ifndef __clang__ 134 manl = d; 135 #else 136 /* clang produces an invalid exception when d >= 2^63, 137 see https://bugs.llvm.org//show_bug.cgi?id=17686. 138 Since this is always the case, here, we use the following patch. */ 139 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 64); 140 manl = 0x8000000000000000 + (mp_limb_t) (d - 0x8000000000000000); 141 #endif /* __clang__ */ 142 #else 143 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 32); 144 manh = (mp_limb_t) d; 145 manl = (mp_limb_t) ((d - manh) * MP_BASE_AS_DOUBLE); 146 #endif 147 } 148 149 #endif /* _MPFR_IEEE_FLOATS */ 150 151 rp[0] = manl; 152 #if GMP_NUMB_BITS == 32 153 rp[1] = manh; 154 #endif 155 156 MPFR_ASSERTD((rp[MPFR_LIMBS_PER_DOUBLE - 1] & MPFR_LIMB_HIGHBIT) != 0); 157 158 return exp; 159 } 160 161 /* End of part included from gmp-2.0.2 */ 162 163 int 164 mpfr_set_d (mpfr_ptr r, double d, mpfr_rnd_t rnd_mode) 165 { 166 int inexact; 167 mpfr_t tmp; 168 mp_limb_t tmpmant[MPFR_LIMBS_PER_DOUBLE]; 169 MPFR_SAVE_EXPO_DECL (expo); 170 171 if (MPFR_UNLIKELY(DOUBLE_ISNAN(d))) 172 { 173 MPFR_SET_NAN(r); 174 MPFR_RET_NAN; 175 } 176 else if (MPFR_UNLIKELY(d == 0)) 177 { 178 #if _MPFR_IEEE_FLOATS 179 union mpfr_ieee_double_extract x; 180 181 MPFR_SET_ZERO(r); 182 /* set correct sign */ 183 x.d = d; 184 if (x.s.sig == 1) 185 MPFR_SET_NEG(r); 186 else 187 MPFR_SET_POS(r); 188 #else /* _MPFR_IEEE_FLOATS */ 189 MPFR_SET_ZERO(r); 190 { 191 /* This is to get the sign of zero on non-IEEE hardware 192 Some systems support +0.0, -0.0, and unsigned zero. 193 Some other systems may just have an unsigned zero. 194 We can't use d == +0.0 since it should be always true, 195 so we check that the memory representation of d is the 196 same than +0.0, etc. 197 Note: r is set to -0 only if d is detected as a negative zero 198 *and*, for the double type, -0 has a different representation 199 from +0. If -0.0 has several representations, the code below 200 may not work as expected, but this is hardly fixable in a 201 portable way (without depending on a math library) and only 202 the sign could be incorrect. Such systems should be taken 203 into account on a case-by-case basis. If the code is changed 204 here, set_d64.c code should be updated too. */ 205 double poszero = +0.0, negzero = DBL_NEG_ZERO; 206 if (memcmp(&d, &poszero, sizeof(double)) == 0) 207 MPFR_SET_POS(r); 208 else if (memcmp(&d, &negzero, sizeof(double)) == 0) 209 MPFR_SET_NEG(r); 210 else 211 MPFR_SET_POS(r); 212 } 213 #endif /* _MPFR_IEEE_FLOATS */ 214 return 0; /* 0 is exact */ 215 } 216 else if (MPFR_UNLIKELY(DOUBLE_ISINF(d))) 217 { 218 MPFR_SET_INF(r); 219 if (d > 0) 220 MPFR_SET_POS(r); 221 else 222 MPFR_SET_NEG(r); 223 return 0; /* infinity is exact */ 224 } 225 226 /* now d is neither 0, nor NaN nor Inf */ 227 228 MPFR_SAVE_EXPO_MARK (expo); 229 230 /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE, 231 since PREC(r) may be different from PREC(tmp), and then both variables 232 would have same precision in the mpfr_set4 call below. */ 233 MPFR_MANT(tmp) = tmpmant; 234 MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG; 235 236 /* don't use MPFR_SET_EXP here since the exponent may be out of range */ 237 MPFR_EXP(tmp) = extract_double (tmpmant, d); 238 239 /* tmp is exact since PREC(tmp)=53 */ 240 inexact = mpfr_set4 (r, tmp, rnd_mode, 241 (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS); 242 243 MPFR_SAVE_EXPO_FREE (expo); 244 return mpfr_check_range (r, inexact, rnd_mode); 245 } 246 247 248 249