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