1/* mpfr_set_ld -- convert a machine long double to 2 a multiple precision floating-point number 3 4Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 5Contributed by the AriC and Caramel projects, INRIA. 6 7This file is part of the GNU MPFR Library. 8 9The GNU MPFR Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MPFR Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2251 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 35static 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 */ 51int 52mpfr_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 238int 239mpfr_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