extract-dbl.c

   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 }