1 /* mpfr_set_ld -- convert a machine long double to 2 a multiple precision floating-point number 3 4 Copyright 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 /* Various i386 systems have been seen with <float.h> LDBL constants equal 30 to the DBL ones, whereas they ought to be bigger, reflecting the 10-byte 31 IEEE extended format on that processor. gcc 3.2.1 on FreeBSD and Solaris 32 has been seen with the problem, and gcc 2.95.4 on FreeBSD 4.7. */ 33 34 #if HAVE_LDOUBLE_IEEE_EXT_LITTLE 35 static const union { 36 char bytes[10]; 37 long double d; 38 } ldbl_max_struct = { 39 { '\377','\377','\377','\377', 40 '\377','\377','\377','\377', 41 '\376','\177' } 42 }; 43 #define MPFR_LDBL_MAX (ldbl_max_struct.d) 44 #else 45 #define MPFR_LDBL_MAX LDBL_MAX 46 #endif 47 48 #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE 49 50 /* Generic code */ 51 int 52 mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode) 53 { 54 mpfr_t t, u; 55 int inexact, shift_exp; 56 long double x; 57 MPFR_SAVE_EXPO_DECL (expo); 58 59 /* Check for NAN */ 60 LONGDOUBLE_NAN_ACTION (d, goto nan); 61 62 /* Check for INF */ 63 if (d > MPFR_LDBL_MAX) 64 { 65 mpfr_set_inf (r, 1); 66 return 0; 67 } 68 else if (d < -MPFR_LDBL_MAX) 69 { 70 mpfr_set_inf (r, -1); 71 return 0; 72 } 73 /* Check for ZERO */ 74 else if (d == 0.0) 75 return mpfr_set_d (r, (double) d, rnd_mode); 76 77 mpfr_init2 (t, MPFR_LDBL_MANT_DIG); 78 mpfr_init2 (u, IEEE_DBL_MANT_DIG); 79 80 MPFR_SAVE_EXPO_MARK (expo); 81 82 convert: 83 x = d; 84 MPFR_SET_ZERO (t); /* The sign doesn't matter. */ 85 shift_exp = 0; /* invariant: remainder to deal with is d*2^shift_exp */ 86 while (x != (long double) 0.0) 87 { 88 /* Check overflow of double */ 89 if (x > (long double) DBL_MAX || (-x) > (long double) DBL_MAX) 90 { 91 long double div9, div10, div11, div12, div13; 92 93 #define TWO_64 18446744073709551616.0 /* 2^64 */ 94 #define TWO_128 (TWO_64 * TWO_64) 95 #define TWO_256 (TWO_128 * TWO_128) 96 div9 = (long double) (double) (TWO_256 * TWO_256); /* 2^(2^9) */ 97 div10 = div9 * div9; 98 div11 = div10 * div10; /* 2^(2^11) */ 99 div12 = div11 * div11; /* 2^(2^12) */ 100 div13 = div12 * div12; /* 2^(2^13) */ 101 if (ABS (x) >= div13) 102 { 103 x /= div13; /* exact */ 104 shift_exp += 8192; 105 mpfr_div_2si (t, t, 8192, MPFR_RNDZ); 106 } 107 if (ABS (x) >= div12) 108 { 109 x /= div12; /* exact */ 110 shift_exp += 4096; 111 mpfr_div_2si (t, t, 4096, MPFR_RNDZ); 112 } 113 if (ABS (x) >= div11) 114 { 115 x /= div11; /* exact */ 116 shift_exp += 2048; 117 mpfr_div_2si (t, t, 2048, MPFR_RNDZ); 118 } 119 if (ABS (x) >= div10) 120 { 121 x /= div10; /* exact */ 122 shift_exp += 1024; 123 mpfr_div_2si (t, t, 1024, MPFR_RNDZ); 124 } 125 /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024, 126 therefore we have one extra exponent reduction step */ 127 if (ABS (x) >= div9) 128 { 129 x /= div9; /* exact */ 130 shift_exp += 512; 131 mpfr_div_2si (t, t, 512, MPFR_RNDZ); 132 } 133 } /* Check overflow of double */ 134 else /* no overflow on double */ 135 { 136 long double div9, div10, div11; 137 138 div9 = (long double) (double) 7.4583407312002067432909653e-155; 139 /* div9 = 2^(-2^9) */ 140 div10 = div9 * div9; /* 2^(-2^10) */ 141 div11 = div10 * div10; /* 2^(-2^11) if extended precision */ 142 /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not 143 overflow here */ 144 if (ABS(x) < div10 && 145 div11 != (long double) 0.0 && 146 div11 / div10 == div10) /* possible underflow */ 147 { 148 long double div12, div13; 149 /* After the divisions, any bit of x must be >= div10, 150 hence the possible division by div9. */ 151 div12 = div11 * div11; /* 2^(-2^12) */ 152 div13 = div12 * div12; /* 2^(-2^13) */ 153 if (ABS (x) <= div13) 154 { 155 x /= div13; /* exact */ 156 shift_exp -= 8192; 157 mpfr_mul_2si (t, t, 8192, MPFR_RNDZ); 158 } 159 if (ABS (x) <= div12) 160 { 161 x /= div12; /* exact */ 162 shift_exp -= 4096; 163 mpfr_mul_2si (t, t, 4096, MPFR_RNDZ); 164 } 165 if (ABS (x) <= div11) 166 { 167 x /= div11; /* exact */ 168 shift_exp -= 2048; 169 mpfr_mul_2si (t, t, 2048, MPFR_RNDZ); 170 } 171 if (ABS (x) <= div10) 172 { 173 x /= div10; /* exact */ 174 shift_exp -= 1024; 175 mpfr_mul_2si (t, t, 1024, MPFR_RNDZ); 176 } 177 if (ABS(x) <= div9) 178 { 179 x /= div9; /* exact */ 180 shift_exp -= 512; 181 mpfr_mul_2si (t, t, 512, MPFR_RNDZ); 182 } 183 } 184 else /* no underflow */ 185 { 186 inexact = mpfr_set_d (u, (double) x, MPFR_RNDZ); 187 MPFR_ASSERTD (inexact == 0); 188 if (mpfr_add (t, t, u, MPFR_RNDZ) != 0) 189 { 190 if (!mpfr_number_p (t)) 191 break; 192 /* Inexact. This cannot happen unless the C implementation 193 "lies" on the precision or when long doubles are 194 implemented with FP expansions like under Mac OS X. */ 195 if (MPFR_PREC (t) != MPFR_PREC (r) + 1) 196 { 197 /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX. 198 The precision MPFR_PREC (r) + 1 allows us to 199 deduce the rounding bit and the sticky bit. */ 200 mpfr_set_prec (t, MPFR_PREC (r) + 1); 201 goto convert; 202 } 203 else 204 { 205 mp_limb_t *tp; 206 int rb_mask; 207 208 /* Since mpfr_add was inexact, the sticky bit is 1. */ 209 tp = MPFR_MANT (t); 210 rb_mask = MPFR_LIMB_ONE << 211 (GMP_NUMB_BITS - 1 - 212 (MPFR_PREC (r) & (GMP_NUMB_BITS - 1))); 213 if (rnd_mode == MPFR_RNDN) 214 rnd_mode = (*tp & rb_mask) ^ MPFR_IS_NEG (t) ? 215 MPFR_RNDU : MPFR_RNDD; 216 *tp |= rb_mask; 217 break; 218 } 219 } 220 x -= (long double) mpfr_get_d1 (u); /* exact */ 221 } 222 } 223 } 224 inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode); 225 mpfr_clear (t); 226 mpfr_clear (u); 227 228 MPFR_SAVE_EXPO_FREE (expo); 229 return mpfr_check_range (r, inexact, rnd_mode); 230 231 nan: 232 MPFR_SET_NAN(r); 233 MPFR_RET_NAN; 234 } 235 236 #else /* IEEE Extended Little Endian Code */ 237 238 int 239 mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode) 240 { 241 int inexact, i, k, cnt; 242 mpfr_t tmp; 243 mp_limb_t tmpmant[MPFR_LIMBS_PER_LONG_DOUBLE]; 244 mpfr_long_double_t x; 245 mpfr_exp_t exp; 246 int signd; 247 MPFR_SAVE_EXPO_DECL (expo); 248 249 /* Check for NAN */ 250 if (MPFR_UNLIKELY (d != d)) 251 { 252 MPFR_SET_NAN (r); 253 MPFR_RET_NAN; 254 } 255 /* Check for INF */ 256 else if (MPFR_UNLIKELY (d > MPFR_LDBL_MAX)) 257 { 258 MPFR_SET_INF (r); 259 MPFR_SET_POS (r); 260 return 0; 261 } 262 else if (MPFR_UNLIKELY (d < -MPFR_LDBL_MAX)) 263 { 264 MPFR_SET_INF (r); 265 MPFR_SET_NEG (r); 266 return 0; 267 } 268 /* Check for ZERO */ 269 else if (MPFR_UNLIKELY (d == 0.0)) 270 { 271 x.ld = d; 272 MPFR_SET_ZERO (r); 273 if (x.s.sign == 1) 274 MPFR_SET_NEG(r); 275 else 276 MPFR_SET_POS(r); 277 return 0; 278 } 279 280 /* now d is neither 0, nor NaN nor Inf */ 281 MPFR_SAVE_EXPO_MARK (expo); 282 283 MPFR_MANT (tmp) = tmpmant; 284 MPFR_PREC (tmp) = 64; 285 286 /* Extract sign */ 287 x.ld = d; 288 signd = MPFR_SIGN_POS; 289 if (x.ld < 0.0) 290 { 291 signd = MPFR_SIGN_NEG; 292 x.ld = -x.ld; 293 } 294 295 /* Extract mantissa */ 296 #if GMP_NUMB_BITS >= 64 297 tmpmant[0] = ((mp_limb_t) x.s.manh << 32) | ((mp_limb_t) x.s.manl); 298 #else 299 tmpmant[0] = (mp_limb_t) x.s.manl; 300 tmpmant[1] = (mp_limb_t) x.s.manh; 301 #endif 302 303 /* Normalize mantissa */ 304 i = MPFR_LIMBS_PER_LONG_DOUBLE; 305 MPN_NORMALIZE_NOT_ZERO (tmpmant, i); 306 k = MPFR_LIMBS_PER_LONG_DOUBLE - i; 307 count_leading_zeros (cnt, tmpmant[i - 1]); 308 if (MPFR_LIKELY (cnt != 0)) 309 mpn_lshift (tmpmant + k, tmpmant, i, cnt); 310 else if (k != 0) 311 MPN_COPY (tmpmant + k, tmpmant, i); 312 if (MPFR_UNLIKELY (k != 0)) 313 MPN_ZERO (tmpmant, k); 314 315 /* Set exponent */ 316 exp = (mpfr_exp_t) ((x.s.exph << 8) + x.s.expl); /* 15-bit unsigned int */ 317 if (MPFR_UNLIKELY (exp == 0)) 318 exp -= 0x3FFD; 319 else 320 exp -= 0x3FFE; 321 322 MPFR_SET_EXP (tmp, exp - cnt - k * GMP_NUMB_BITS); 323 324 /* tmp is exact */ 325 inexact = mpfr_set4 (r, tmp, rnd_mode, signd); 326 327 MPFR_SAVE_EXPO_FREE (expo); 328 return mpfr_check_range (r, inexact, rnd_mode); 329 } 330 331 #endif 332