contrib/gmp/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 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 }