[BACK]Return to extract-double.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2

Annotation of OpenXM/src/kan96xx/gmp-2.0.2/extract-double.c, Revision 1.1.1.1

1.1       maekawa     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 "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 MP_BASE_AS_DOUBLE (2.0 * ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1)))
                     34:
                     35: /* Extract a non-negative double in d.  */
                     36:
                     37: int
                     38: #if __STDC__
                     39: __gmp_extract_double (mp_ptr rp, double d)
                     40: #else
                     41: __gmp_extract_double (rp, d)
                     42:      mp_ptr rp;
                     43:      double d;
                     44: #endif
                     45: {
                     46:   long exp;
                     47:   unsigned sc;
                     48:   mp_limb_t manh, manl;
                     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 BITS_PER_MP_LIMB >= 32.
                     55:      4. This lits is incomplete and misspelled.
                     56:    */
                     57:
                     58:   if (d == 0.0)
                     59:     {
                     60:       rp[0] = 0;
                     61:       rp[1] = 0;
                     62: #if BITS_PER_MP_LIMB == 32
                     63:       rp[2] = 0;
                     64: #endif
                     65:       return 0;
                     66:     }
                     67:
                     68: #if _GMP_IEEE_FLOATS
                     69:   {
                     70:     union ieee_double_extract x;
                     71:     x.d = d;
                     72:
                     73:     exp = x.s.exp;
                     74:     sc = (unsigned) (exp + 2) % BITS_PER_MP_LIMB;
                     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: #else
                     79:     manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
                     80:     manl = x.s.manl << 11;
                     81: #endif
                     82:   }
                     83: #else
                     84:   {
                     85:     /* Unknown (or known to be non-IEEE) double format.  */
                     86:     exp = 0;
                     87:     if (d >= 1.0)
                     88:       {
                     89:        if (d * 0.5 == d)
                     90:          abort ();
                     91:
                     92:        while (d >= 32768.0)
                     93:          {
                     94:            d *= (1.0 / 65536.0);
                     95:            exp += 16;
                     96:          }
                     97:        while (d >= 1.0)
                     98:          {
                     99:            d *= 0.5;
                    100:            exp += 1;
                    101:          }
                    102:       }
                    103:     else if (d < 0.5)
                    104:       {
                    105:        while (d < (1.0 / 65536.0))
                    106:          {
                    107:            d *=  65536.0;
                    108:            exp -= 16;
                    109:          }
                    110:        while (d < 0.5)
                    111:          {
                    112:            d *= 2.0;
                    113:            exp -= 1;
                    114:          }
                    115:       }
                    116:
                    117:     sc = (unsigned) exp % BITS_PER_MP_LIMB;
                    118:
                    119:     d *= MP_BASE_AS_DOUBLE;
                    120: #if BITS_PER_MP_LIMB == 64
                    121:     manl = d;
                    122: #else
                    123:     manh = d;
                    124:     manl = (d - manh) * MP_BASE_AS_DOUBLE;
                    125: #endif
                    126:
                    127:     exp += 1022;
                    128:   }
                    129: #endif
                    130:
                    131:   exp = (unsigned) (exp + 1) / BITS_PER_MP_LIMB - 1024 / BITS_PER_MP_LIMB + 1;
                    132:
                    133: #if BITS_PER_MP_LIMB == 64
                    134:   if (sc != 0)
                    135:     {
                    136:       rp[1] = manl >> (BITS_PER_MP_LIMB - sc);
                    137:       rp[0] = manl << sc;
                    138:     }
                    139:   else
                    140:     {
                    141:       rp[1] = manl;
                    142:       rp[0] = 0;
                    143:     }
                    144: #else
                    145:   if (sc != 0)
                    146:     {
                    147:       rp[2] = manh >> (BITS_PER_MP_LIMB - sc);
                    148:       rp[1] = (manl >> (BITS_PER_MP_LIMB - sc)) | (manh << sc);
                    149:       rp[0] = manl << sc;
                    150:     }
                    151:   else
                    152:     {
                    153:       rp[2] = manh;
                    154:       rp[1] = manl;
                    155:       rp[0] = 0;
                    156:     }
                    157: #endif
                    158:
                    159:   return exp;
                    160: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>