rts/gmp/extract-dbl.c

   1 /* __gmp_extract_double -- convert from double to array of mp_limb_t.
   2
   3 Copyright (C) 1996, 1999, 2000 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., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, 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 /* Extract a non-negative double in d.  */
  34
  35 int
  36 #if __STDC__
  37 __gmp_extract_double (mp_ptr rp, double d)
  38 #else
  39 __gmp_extract_double (rp, d)
  40      mp_ptr rp;
  41      double d;
  42 #endif
  43 {
  44   long exp;
  45   unsigned sc;
  46   mp_limb_t manh, manl;
  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 BITS_PER_MP_LIMB >= 32.
  53      4. This lits is incomplete and misspelled.
  54    */
  55
  56   if (d == 0.0)
  57     {
  58       rp[0] = 0;
  59       rp[1] = 0;
  60 #if BITS_PER_MP_LIMB == 32
  61       rp[2] = 0;
  62 #endif
  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_MP_LIMB == 64
  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 ((mp_limb_signed_t) manl >= 0);
  89       }
  90 #else
  91     manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
  92     manl = x.s.manl << 11;
  93     if (exp == 0)
  94       {
  95         /* Denormalized number.  Don't try to be clever about this,
  96            since it is not an important case to make fast.  */
  97         exp = 1;
  98         do
  99           {
 100             manh = (manh << 1) | (manl >> 31);
 101             manl = manl << 1;
 102             exp--;
 103           }
 104         while ((mp_limb_signed_t) manh >= 0);
 105       }
 106 #endif
 107     exp -= 1022;                /* Remove IEEE bias.  */
 108   }
 109 #else
 110   {
 111     /* Unknown (or known to be non-IEEE) double format.  */
 112     exp = 0;
 113     if (d >= 1.0)
 114       {
 115         if (d * 0.5 == d)
 116           abort ();
 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 *= MP_BASE_AS_DOUBLE;
 144 #if BITS_PER_MP_LIMB == 64
 145     manl = d;
 146 #else
 147     manh = d;
 148     manl = (d - manh) * MP_BASE_AS_DOUBLE;
 149 #endif
 150   }
 151 #endif
 152
 153   sc = (unsigned) exp % BITS_PER_MP_LIMB;
 154
 155   /* We add something here to get rounding right.  */
 156   exp = (exp + 2048) / BITS_PER_MP_LIMB - 2048 / BITS_PER_MP_LIMB + 1;
 157
 158 #if BITS_PER_MP_LIMB == 64
 159   if (sc != 0)
 160     {
 161       rp[1] = manl >> (BITS_PER_MP_LIMB - sc);
 162       rp[0] = manl << sc;
 163     }
 164   else
 165     {
 166       rp[1] = manl;
 167       rp[0] = 0;
 168       exp--;
 169     }
 170 #else
 171   if (sc != 0)
 172     {
 173       rp[2] = manh >> (BITS_PER_MP_LIMB - sc);
 174       rp[1] = (manl >> (BITS_PER_MP_LIMB - sc)) | (manh << sc);
 175       rp[0] = manl << sc;
 176     }
 177   else
 178     {
 179       rp[2] = manh;
 180       rp[1] = manl;
 181       rp[0] = 0;
 182       exp--;
 183     }
 184 #endif
 185
 186   return exp;
 187 }