1 /* mpfr_get_ld, mpfr_get_ld_2exp -- convert a multiple precision floating-point 2 number to a machine long double 3 4 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 5 Contributed by the Arenaire 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 #include "mpfr-impl.h" 27 28 #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE 29 30 long double 31 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) 32 { 33 34 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 35 return (long double) mpfr_get_d (x, rnd_mode); 36 else /* now x is a normal non-zero number */ 37 { 38 long double r; /* result */ 39 long double m; 40 double s; /* part of result */ 41 mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */ 42 mpfr_t y, z; 43 int sign; 44 45 /* first round x to the target long double precision, so that 46 all subsequent operations are exact (this avoids double rounding 47 problems) */ 48 mpfr_init2 (y, MPFR_LDBL_MANT_DIG); 49 mpfr_init2 (z, MPFR_LDBL_MANT_DIG); 50 /* Note about the precision of z: even though IEEE_DBL_MANT_DIG is 51 sufficient, z has been set to the same precision as y so that 52 the mpfr_sub below calls mpfr_sub1sp, which is faster than the 53 generic subtraction, even in this particular case (from tests 54 done by Patrick Pelissier on a 64-bit Core2 Duo against r7285). 55 But here there is an important cancellation in the subtraction. 56 TODO: get more information about what has been tested. */ 57 58 mpfr_set (y, x, rnd_mode); 59 sh = MPFR_GET_EXP (y); 60 sign = MPFR_SIGN (y); 61 MPFR_SET_EXP (y, 0); 62 MPFR_SET_POS (y); 63 64 r = 0.0; 65 do { 66 s = mpfr_get_d (y, MPFR_RNDN); /* high part of y */ 67 r += (long double) s; 68 mpfr_set_d (z, s, MPFR_RNDN); /* exact */ 69 mpfr_sub (y, y, z, MPFR_RNDN); /* exact */ 70 } while (!MPFR_IS_ZERO (y)); 71 72 mpfr_clear (z); 73 mpfr_clear (y); 74 75 /* we now have to multiply back by 2^sh */ 76 MPFR_ASSERTD (r > 0); 77 if (sh != 0) 78 { 79 /* An overflow may occurs (example: 0.5*2^1024) */ 80 while (r < 1.0) 81 { 82 r += r; 83 sh--; 84 } 85 86 if (sh > 0) 87 m = 2.0; 88 else 89 { 90 m = 0.5; 91 sh = -sh; 92 } 93 94 for (;;) 95 { 96 if (sh % 2) 97 r = r * m; 98 sh >>= 1; 99 if (sh == 0) 100 break; 101 m = m * m; 102 } 103 } 104 if (sign < 0) 105 r = -r; 106 return r; 107 } 108 } 109 110 #else 111 112 long double 113 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) 114 { 115 mpfr_long_double_t ld; 116 mpfr_t tmp; 117 int inex; 118 MPFR_SAVE_EXPO_DECL (expo); 119 120 MPFR_SAVE_EXPO_MARK (expo); 121 122 mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG); 123 inex = mpfr_set (tmp, x, rnd_mode); 124 125 mpfr_set_emin (-16382-63); 126 mpfr_set_emax (16384); 127 mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode); 128 mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */ 129 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp))) 130 ld.ld = (long double) mpfr_get_d (tmp, rnd_mode); 131 else 132 { 133 mp_limb_t *tmpmant; 134 mpfr_exp_t e, denorm; 135 136 tmpmant = MPFR_MANT (tmp); 137 e = MPFR_GET_EXP (tmp); 138 /* the smallest normal number is 2^(-16382), which is 0.5*2^(-16381) 139 in MPFR, thus any exponent <= -16382 corresponds to a subnormal 140 number */ 141 denorm = MPFR_UNLIKELY (e <= -16382) ? - e - 16382 + 1 : 0; 142 #if GMP_NUMB_BITS >= 64 143 ld.s.manl = (tmpmant[0] >> denorm); 144 ld.s.manh = (tmpmant[0] >> denorm) >> 32; 145 #elif GMP_NUMB_BITS == 32 146 if (MPFR_LIKELY (denorm == 0)) 147 { 148 ld.s.manl = tmpmant[0]; 149 ld.s.manh = tmpmant[1]; 150 } 151 else if (denorm < 32) 152 { 153 ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm)); 154 ld.s.manh = tmpmant[1] >> denorm; 155 } 156 else /* 32 <= denorm <= 64 */ 157 { 158 ld.s.manl = tmpmant[1] >> (denorm - 32); 159 ld.s.manh = 0; 160 } 161 #else 162 # error "GMP_NUMB_BITS must be 32 or >= 64" 163 /* Other values have never been supported anyway. */ 164 #endif 165 if (MPFR_LIKELY (denorm == 0)) 166 { 167 ld.s.exph = (e + 0x3FFE) >> 8; 168 ld.s.expl = (e + 0x3FFE); 169 } 170 else 171 ld.s.exph = ld.s.expl = 0; 172 ld.s.sign = MPFR_IS_NEG (x); 173 } 174 175 mpfr_clear (tmp); 176 MPFR_SAVE_EXPO_FREE (expo); 177 return ld.ld; 178 } 179 180 #endif 181 182 /* contributed by Damien Stehle */ 183 long double 184 mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode) 185 { 186 long double ret; 187 mpfr_exp_t exp; 188 mpfr_t tmp; 189 190 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src))) 191 return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode); 192 193 tmp[0] = *src; /* Hack copy mpfr_t */ 194 MPFR_SET_EXP (tmp, 0); 195 ret = mpfr_get_ld (tmp, rnd_mode); 196 197 if (MPFR_IS_PURE_FP(src)) 198 { 199 exp = MPFR_GET_EXP (src); 200 201 /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */ 202 if (ret == 1.0) 203 { 204 ret = 0.5; 205 exp ++; 206 } 207 else if (ret == -1.0) 208 { 209 ret = -0.5; 210 exp ++; 211 } 212 213 MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0) 214 || (ret <= -0.5 && ret > -1.0)); 215 MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX); 216 } 217 else 218 exp = 0; 219 220 *expptr = exp; 221 return ret; 222 } 223