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 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 2.1 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; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  20 MA 02110-1301, USA. */
  21
  22 #include "gmp.h"
  23 #include "gmp-impl.h"
  24
  25 #ifdef XDEBUG
  26 #undef _GMP_IEEE_FLOATS
  27 #endif
  28
  29 #ifndef _GMP_IEEE_FLOATS
  30 #define _GMP_IEEE_FLOATS 0
  31 #endif
  32
  33 #define BITS_IN_MANTISSA 53
  34
  35 /* Extract a non-negative double in d.  */
  36
  37 int
  38 __gmp_extract_double (mp_ptr rp, double d)
  39 {
  40   long exp;
  41   unsigned sc;
  42 #ifdef _LONG_LONG_LIMB
  43 #define BITS_PER_PART 64        /* somewhat bogus */
  44   unsigned long long int manl;
  45 #else
  46 #define BITS_PER_PART GMP_LIMB_BITS
  47   unsigned long int manh, manl;
  48 #endif
  49
  50   /* BUGS
  51
  52      1. Should handle Inf and NaN in IEEE specific code.
  53      2. Handle Inf and NaN also in default code, to avoid hangs.
  54      3. Generalize to handle all GMP_LIMB_BITS >= 32.
  55      4. This lits is incomplete and misspelled.
  56    */
  57
  58   ASSERT (d >= 0.0);
  59
  60   if (d == 0.0)
  61     {
  62       MPN_ZERO (rp, LIMBS_PER_DOUBLE);
  63       return 0;
  64     }
  65
  66 #if _GMP_IEEE_FLOATS
  67   {
  68 #if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8
  69     /* Work around alpha-specific bug in GCC 2.8.x.  */
  70     volatile
  71 #endif
  72     union ieee_double_extract x;
  73     x.d = d;
  74     exp = x.s.exp;
  75 #if BITS_PER_PART == 64         /* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */
  76     manl = (((mp_limb_t) 1 << 63)
  77             | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
  78     if (exp == 0)
  79       {
  80         /* Denormalized number.  Don't try to be clever about this,
  81            since it is not an important case to make fast.  */
  82         exp = 1;
  83         do
  84           {
  85             manl = manl << 1;
  86             exp--;
  87           }
  88         while ((manl & GMP_LIMB_HIGHBIT) == 0);
  89       }
  90 #endif
  91 #if BITS_PER_PART == 32
  92     manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
  93     manl = x.s.manl << 11;
  94     if (exp == 0)
  95       {
  96         /* Denormalized number.  Don't try to be clever about this,
  97            since it is not an important case to make fast.  */
  98         exp = 1;
  99         do
 100           {
 101             manh = (manh << 1) | (manl >> 31);
 102             manl = manl << 1;
 103             exp--;
 104           }
 105         while ((manh & GMP_LIMB_HIGHBIT) == 0);
 106       }
 107 #endif
 108 #if BITS_PER_PART != 32 && BITS_PER_PART != 64
 109   You need to generalize the code above to handle this.
 110 #endif
 111     exp -= 1022;                /* Remove IEEE bias.  */
 112   }
 113 #else
 114   {
 115     /* Unknown (or known to be non-IEEE) double format.  */
 116     exp = 0;
 117     if (d >= 1.0)
 118       {
 119         ASSERT_ALWAYS (d * 0.5 != d);
 120
 121         while (d >= 32768.0)
 122           {
 123             d *= (1.0 / 65536.0);
 124             exp += 16;
 125           }
 126         while (d >= 1.0)
 127           {
 128             d *= 0.5;
 129             exp += 1;
 130           }
 131       }
 132     else if (d < 0.5)
 133       {
 134         while (d < (1.0 / 65536.0))
 135           {
 136             d *=  65536.0;
 137             exp -= 16;
 138           }
 139         while (d < 0.5)
 140           {
 141             d *= 2.0;
 142             exp -= 1;
 143           }
 144       }
 145
 146     d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
 147 #if BITS_PER_PART == 64
 148     manl = d;
 149 #else
 150     manh = d;
 151     manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
 152 #endif
 153   }
 154 #endif /* IEEE */
 155
 156   /* Up until here, we have ignored the actual limb size.  Remains
 157      to split manh,,manl into an array of LIMBS_PER_DOUBLE limbs.
 158   */
 159
 160   sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS;
 161
 162   /* We add something here to get rounding right.  */
 163   exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1;
 164
 165 #if LIMBS_PER_DOUBLE == 2
 166 #if GMP_NAIL_BITS == 0
 167   if (sc != 0)
 168     {
 169       rp[1] = manl >> (GMP_LIMB_BITS - sc);
 170       rp[0] = manl << sc;
 171     }
 172   else
 173     {
 174       rp[1] = manl;
 175       rp[0] = 0;
 176       exp--;
 177     }
 178 #else
 179   if (sc > GMP_NAIL_BITS)
 180     {
 181       rp[1] = manl >> (GMP_LIMB_BITS - sc);
 182       rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK;
 183     }
 184   else
 185     {
 186       if (sc == 0)
 187         {
 188           rp[1] = manl >> GMP_NAIL_BITS;
 189           rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
 190           exp--;
 191         }
 192       else
 193         {
 194           rp[1] = manl >> (GMP_LIMB_BITS - sc);
 195           rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
 196         }
 197     }
 198 #endif
 199 #endif
 200
 201 #if LIMBS_PER_DOUBLE == 3
 202 #if GMP_NAIL_BITS == 0
 203   if (sc != 0)
 204     {
 205       rp[2] = manh >> (GMP_LIMB_BITS - sc);
 206       rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
 207       rp[0] = manl << sc;
 208     }
 209   else
 210     {
 211       rp[2] = manh;
 212       rp[1] = manl;
 213       rp[0] = 0;
 214       exp--;
 215     }
 216 #else
 217   if (sc > GMP_NAIL_BITS)
 218     {
 219       rp[2] = (manh >> (GMP_LIMB_BITS - sc));
 220       rp[1] = ((manh << (sc - GMP_NAIL_BITS)) |
 221                (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK;
 222       if (sc >= 2 * GMP_NAIL_BITS)
 223         rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
 224       else
 225         rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
 226     }
 227   else
 228     {
 229       if (sc == 0)
 230         {
 231           rp[2] = manh >> GMP_NAIL_BITS;
 232           rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK;
 233           rp[0] = 0;
 234           exp--;
 235         }
 236       else
 237         {
 238           rp[2] = (manh >> (GMP_LIMB_BITS - sc));
 239           rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
 240           rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc))
 241                    | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK;
 242         }
 243     }
 244 #endif
 245 #endif
 246
 247 #if LIMBS_PER_DOUBLE > 3
 248   if (sc == 0)
 249     {
 250       int i;
 251
 252       for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--)
 253         {
 254           rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
 255           manh = ((manh << GMP_NUMB_BITS)
 256                   | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
 257           manl = manl << GMP_NUMB_BITS;
 258         }
 259       exp--;
 260     }
 261   else
 262     {
 263       int i;
 264
 265       rp[LIMBS_PER_DOUBLE - 1] = (manh >> (GMP_LIMB_BITS - sc));
 266       manh = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
 267       manl = (manl << sc);
 268       for (i = LIMBS_PER_DOUBLE - 2; i >= 0; i--)
 269         {
 270           rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
 271           manh = ((manh << GMP_NUMB_BITS)
 272                   | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
 273           manl = manl << GMP_NUMB_BITS;
 274         }
 275   }
 276 #endif
 277
 278   return exp;
 279 }