lib/gmp/extract-dbl.c

   1 /* __gmp_extract_double -- convert from double to array of mp_limb_t.
   2
   3 Copyright (C) 1996 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 Library General Public License as published by
   9 the Free Software Foundation; either version 2 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 Library General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Library 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 <config.h>
  23 #include "gmp.h"
  24 #include "gmp-impl.h"
  25
  26 #ifdef XDEBUG
  27 #undef _GMP_IEEE_FLOATS
  28 #endif
  29
  30 #ifndef _GMP_IEEE_FLOATS
  31 #define _GMP_IEEE_FLOATS 0
  32 #endif
  33
  34 #define MP_BASE_AS_DOUBLE (2.0 * ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1)))
  35
  36 /* Extract a non-negative double in d.  */
  37
  38 int
  39 #if __STDC__
  40 __gmp_extract_double (mp_ptr rp, double d)
  41 #else
  42 __gmp_extract_double (rp, d)
  43      mp_ptr rp;
  44      double d;
  45 #endif
  46 {
  47   long exp;
  48   unsigned sc;
  49   mp_limb_t manh, manl;
  50
  51   /* BUGS
  52
  53      1. Should handle Inf and NaN in IEEE specific code.
  54      2. Handle Inf and NaN also in default code, to avoid hangs.
  55      3. Generalize to handle all BITS_PER_MP_LIMB >= 32.
  56      4. This lits is incomplete and misspelled.
  57    */
  58
  59   if (d == 0.0)
  60     {
  61       rp[0] = 0;
  62       rp[1] = 0;
  63 #if BITS_PER_MP_LIMB == 32
  64       rp[2] = 0;
  65 #endif
  66       return 0;
  67     }
  68
  69 #if _GMP_IEEE_FLOATS
  70   {
  71     union ieee_double_extract x;
  72     x.d = d;
  73
  74     exp = x.s.exp;
  75     sc = (unsigned) (exp + 2) % BITS_PER_MP_LIMB;
  76 #if BITS_PER_MP_LIMB == 64
  77     manl = (((mp_limb_t) 1 << 63)
  78             | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
  79 #else
  80     manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
  81     manl = x.s.manl << 11;
  82 #endif
  83   }
  84 #else
  85   {
  86     /* Unknown (or known to be non-IEEE) double format.  */
  87     exp = 0;
  88     if (d >= 1.0)
  89       {
  90         if (d * 0.5 == d)
  91           abort ();
  92
  93         while (d >= 32768.0)
  94           {
  95             d *= (1.0 / 65536.0);
  96             exp += 16;
  97           }
  98         while (d >= 1.0)
  99           {
 100             d *= 0.5;
 101             exp += 1;
 102           }
 103       }
 104     else if (d < 0.5)
 105       {
 106         while (d < (1.0 / 65536.0))
 107           {
 108             d *=  65536.0;
 109             exp -= 16;
 110           }
 111         while (d < 0.5)
 112           {
 113             d *= 2.0;
 114             exp -= 1;
 115           }
 116       }
 117
 118     sc = (unsigned) exp % BITS_PER_MP_LIMB;
 119
 120     d *= MP_BASE_AS_DOUBLE;
 121 #if BITS_PER_MP_LIMB == 64
 122     manl = d;
 123 #else
 124     manh = d;
 125     manl = (d - manh) * MP_BASE_AS_DOUBLE;
 126 #endif
 127
 128     exp += 1022;
 129   }
 130 #endif
 131
 132   exp = (unsigned) (exp + 1) / BITS_PER_MP_LIMB - 1024 / BITS_PER_MP_LIMB + 1;
 133
 134 #if BITS_PER_MP_LIMB == 64
 135   if (sc != 0)
 136     {
 137       rp[1] = manl >> (BITS_PER_MP_LIMB - sc);
 138       rp[0] = manl << sc;
 139     }
 140   else
 141     {
 142       rp[1] = manl;
 143       rp[0] = 0;
 144     }
 145 #else
 146   if (sc != 0)
 147     {
 148       rp[2] = manh >> (BITS_PER_MP_LIMB - sc);
 149       rp[1] = (manl >> (BITS_PER_MP_LIMB - sc)) | (manh << sc);
 150       rp[0] = manl << sc;
 151     }
 152   else
 153     {
 154       rp[2] = manh;
 155       rp[1] = manl;
 156       rp[0] = 0;
 157     }
 158 #endif
 159
 160   return exp;
 161 }