1 /* __gmp_extract_double -- convert from double to array of mp_limb_t. 2 3 Copyright 1996, 1999, 2000, 2001, 2002, 2006, 2012 Free Software Foundation, 4 Inc. 5 6 This file is part of the GNU MP Library. 7 8 The GNU MP Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MP Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 20 21 #include "gmp.h" 22 #include "gmp-impl.h" 23 24 #ifdef XDEBUG 25 #undef _GMP_IEEE_FLOATS 26 #endif 27 28 #ifndef _GMP_IEEE_FLOATS 29 #define _GMP_IEEE_FLOATS 0 30 #endif 31 32 /* Extract a non-negative double in d. */ 33 34 int 35 __gmp_extract_double (mp_ptr rp, double d) 36 { 37 long exp; 38 unsigned sc; 39 #ifdef _LONG_LONG_LIMB 40 #define BITS_PER_PART 64 /* somewhat bogus */ 41 unsigned long long int manl; 42 #else 43 #define BITS_PER_PART GMP_LIMB_BITS 44 unsigned long int manh, manl; 45 #endif 46 47 /* BUGS 48 49 1. Should handle Inf and NaN in IEEE specific code. 50 2. Handle Inf and NaN also in default code, to avoid hangs. 51 3. Generalize to handle all GMP_LIMB_BITS >= 32. 52 4. This lits is incomplete and misspelled. 53 */ 54 55 ASSERT (d >= 0.0); 56 57 if (d == 0.0) 58 { 59 MPN_ZERO (rp, LIMBS_PER_DOUBLE); 60 return 0; 61 } 62 63 #if _GMP_IEEE_FLOATS 64 { 65 #if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8 66 /* Work around alpha-specific bug in GCC 2.8.x. */ 67 volatile 68 #endif 69 union ieee_double_extract x; 70 x.d = d; 71 exp = x.s.exp; 72 #if BITS_PER_PART == 64 /* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */ 73 manl = (((mp_limb_t) 1 << 63) 74 | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11)); 75 if (exp == 0) 76 { 77 /* Denormalized number. Don't try to be clever about this, 78 since it is not an important case to make fast. */ 79 exp = 1; 80 do 81 { 82 manl = manl << 1; 83 exp--; 84 } 85 while ((manl & GMP_LIMB_HIGHBIT) == 0); 86 } 87 #endif 88 #if BITS_PER_PART == 32 89 manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21); 90 manl = x.s.manl << 11; 91 if (exp == 0) 92 { 93 /* Denormalized number. Don't try to be clever about this, 94 since it is not an important case to make fast. */ 95 exp = 1; 96 do 97 { 98 manh = (manh << 1) | (manl >> 31); 99 manl = manl << 1; 100 exp--; 101 } 102 while ((manh & GMP_LIMB_HIGHBIT) == 0); 103 } 104 #endif 105 #if BITS_PER_PART != 32 && BITS_PER_PART != 64 106 You need to generalize the code above to handle this. 107 #endif 108 exp -= 1022; /* Remove IEEE bias. */ 109 } 110 #else 111 { 112 /* Unknown (or known to be non-IEEE) double format. */ 113 exp = 0; 114 if (d >= 1.0) 115 { 116 ASSERT_ALWAYS (d * 0.5 != d); 117 118 while (d >= 32768.0) 119 { 120 d *= (1.0 / 65536.0); 121 exp += 16; 122 } 123 while (d >= 1.0) 124 { 125 d *= 0.5; 126 exp += 1; 127 } 128 } 129 else if (d < 0.5) 130 { 131 while (d < (1.0 / 65536.0)) 132 { 133 d *= 65536.0; 134 exp -= 16; 135 } 136 while (d < 0.5) 137 { 138 d *= 2.0; 139 exp -= 1; 140 } 141 } 142 143 d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2))); 144 #if BITS_PER_PART == 64 145 manl = d; 146 #endif 147 #if BITS_PER_PART == 32 148 manh = d; 149 manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2))); 150 #endif 151 } 152 #endif /* IEEE */ 153 154 sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS; 155 156 /* We add something here to get rounding right. */ 157 exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1; 158 159 #if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 2 160 #if GMP_NAIL_BITS == 0 161 if (sc != 0) 162 { 163 rp[1] = manl >> (GMP_LIMB_BITS - sc); 164 rp[0] = manl << sc; 165 } 166 else 167 { 168 rp[1] = manl; 169 rp[0] = 0; 170 exp--; 171 } 172 #else 173 if (sc > GMP_NAIL_BITS) 174 { 175 rp[1] = manl >> (GMP_LIMB_BITS - sc); 176 rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK; 177 } 178 else 179 { 180 if (sc == 0) 181 { 182 rp[1] = manl >> GMP_NAIL_BITS; 183 rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK; 184 exp--; 185 } 186 else 187 { 188 rp[1] = manl >> (GMP_LIMB_BITS - sc); 189 rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK; 190 } 191 } 192 #endif 193 #endif 194 195 #if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 3 196 if (sc > GMP_NAIL_BITS) 197 { 198 rp[2] = manl >> (GMP_LIMB_BITS - sc); 199 rp[1] = (manl << sc - GMP_NAIL_BITS) & GMP_NUMB_MASK; 200 if (sc >= 2 * GMP_NAIL_BITS) 201 rp[0] = 0; 202 else 203 rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK; 204 } 205 else 206 { 207 if (sc == 0) 208 { 209 rp[2] = manl >> GMP_NAIL_BITS; 210 rp[1] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK; 211 rp[0] = 0; 212 exp--; 213 } 214 else 215 { 216 rp[2] = manl >> (GMP_LIMB_BITS - sc); 217 rp[1] = (manl >> GMP_NAIL_BITS - sc) & GMP_NUMB_MASK; 218 rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK; 219 } 220 } 221 #endif 222 223 #if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE == 3 224 #if GMP_NAIL_BITS == 0 225 if (sc != 0) 226 { 227 rp[2] = manh >> (GMP_LIMB_BITS - sc); 228 rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc)); 229 rp[0] = manl << sc; 230 } 231 else 232 { 233 rp[2] = manh; 234 rp[1] = manl; 235 rp[0] = 0; 236 exp--; 237 } 238 #else 239 if (sc > GMP_NAIL_BITS) 240 { 241 rp[2] = (manh >> (GMP_LIMB_BITS - sc)); 242 rp[1] = ((manh << (sc - GMP_NAIL_BITS)) | 243 (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK; 244 if (sc >= 2 * GMP_NAIL_BITS) 245 rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK; 246 else 247 rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK; 248 } 249 else 250 { 251 if (sc == 0) 252 { 253 rp[2] = manh >> GMP_NAIL_BITS; 254 rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK; 255 rp[0] = (manl << GMP_NUMB_BITS - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK; 256 exp--; 257 } 258 else 259 { 260 rp[2] = (manh >> (GMP_LIMB_BITS - sc)); 261 rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK; 262 rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)) 263 | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK; 264 } 265 } 266 #endif 267 #endif 268 269 #if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE > 3 270 if (sc == 0) 271 { 272 int i; 273 274 for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--) 275 { 276 rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS); 277 manh = ((manh << GMP_NUMB_BITS) 278 | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS))); 279 manl = manl << GMP_NUMB_BITS; 280 } 281 exp--; 282 } 283 else 284 { 285 int i; 286 287 rp[LIMBS_PER_DOUBLE - 1] = (manh >> (GMP_LIMB_BITS - sc)); 288 manh = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc)); 289 manl = (manl << sc); 290 for (i = LIMBS_PER_DOUBLE - 2; i >= 0; i--) 291 { 292 rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS); 293 manh = ((manh << GMP_NUMB_BITS) 294 | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS))); 295 manl = manl << GMP_NUMB_BITS; 296 } 297 } 298 #endif 299 300 return exp; 301 } 302