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