source/libs/gmp/gmp-src/extract-dbl.c

   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 }