[BACK]Return to set_d.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/set_d.c, Revision 1.1.1.1

1.1       maekawa     1: /* mpfr_set_d, mpfr_get_d -- convert a multiple precision floating-point number
                      2:                              from/to a machine double precision float
                      3:
                      4: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
                      5:
                      6: This file is part of the MPFR Library.
                      7:
                      8: The MPFR Library is free software; you can redistribute it and/or modify
                      9: it under the terms of the GNU Library General Public License as published by
                     10: the Free Software Foundation; either version 2 of the License, or (at your
                     11: option) any later version.
                     12:
                     13: The MPFR Library is distributed in the hope that it will be useful, but
                     14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     15: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
                     16: License for more details.
                     17:
                     18: You should have received a copy of the GNU Library General Public License
                     19: along with the MPFR Library; see the file COPYING.LIB.  If not, write to
                     20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     21: MA 02111-1307, USA. */
                     22:
                     23: #if __GNUC__ /* gcc "patched" headers seem to omit isnan... */
                     24: extern int isnan(double);
                     25: #endif
                     26: #include <math.h> /* for isnan and NaN */
                     27:
                     28: #include "gmp.h"
                     29: #include "gmp-impl.h"
                     30: #include "longlong.h"
                     31: #include "mpfr.h"
                     32:
                     33: #define NaN sqrt(-1) /* ensures a machine-independent NaN */
                     34:
                     35: /* Included from gmp-2.0.2, patched to support denorms */
                     36:
                     37: #ifdef XDEBUG
                     38: #undef _GMP_IEEE_FLOATS
                     39: #endif
                     40:
                     41: #ifndef _GMP_IEEE_FLOATS
                     42: #define _GMP_IEEE_FLOATS 0
                     43: #endif
                     44:
                     45: int
                     46: #if __STDC__
                     47: __mpfr_extract_double (mp_ptr rp, double d, int e)
                     48: #else
                     49: __mpfr_extract_double (rp, d, e)
                     50:      mp_ptr rp;
                     51:      double d;
                     52:      int e;
                     53: #endif
                     54:      /* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */
                     55: {
                     56:   long exp;
                     57:   mp_limb_t manl;
                     58: #if BITS_PER_MP_LIMB == 32
                     59:   mp_limb_t manh;
                     60: #endif
                     61:
                     62:   /* BUGS
                     63:
                     64:      1. Should handle Inf and NaN in IEEE specific code.
                     65:      2. Handle Inf and NaN also in default code, to avoid hangs.
                     66:      3. Generalize to handle all BITS_PER_MP_LIMB >= 32.
                     67:      4. This lits is incomplete and misspelled.
                     68:    */
                     69:
                     70:   if (d == 0.0)
                     71:     {
                     72:       rp[0] = 0;
                     73: #if BITS_PER_MP_LIMB == 32
                     74:       if (e) rp[1] = 0;
                     75: #endif
                     76:       return 0;
                     77:     }
                     78:
                     79: #if _GMP_IEEE_FLOATS
                     80:   {
                     81:     union ieee_double_extract x;
                     82:     x.d = d;
                     83:
                     84:     exp = x.s.exp;
                     85:     if (exp)
                     86:       {
                     87: #if BITS_PER_MP_LIMB == 64
                     88:        manl = (((mp_limb_t) 1 << 63)
                     89:                | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
                     90: #else
                     91:        manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
                     92:        manl = x.s.manl << 11;
                     93: #endif
                     94:       }
                     95:     else
                     96:       {
                     97: #if BITS_PER_MP_LIMB == 64
                     98:        manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11);
                     99: #else
                    100:     manh = (x.s.manh << 11) | (x.s.manl >> 21);
                    101:        manl = x.s.manl << 11;
                    102: #endif
                    103:       }
                    104:   }
                    105: #else
                    106:   {
                    107:     /* Unknown (or known to be non-IEEE) double format.  */
                    108:     exp = 0;
                    109:     if (d >= 1.0)
                    110:       {
                    111:         if (d * 0.5 == d)
                    112:           abort ();
                    113:
                    114:         while (d >= 32768.0)
                    115:           {
                    116:             d *= (1.0 / 65536.0);
                    117:             exp += 16;
                    118:           }
                    119:         while (d >= 1.0)
                    120:           {
                    121:             d *= 0.5;
                    122:             exp += 1;
                    123:           }
                    124:       }
                    125:     else if (d < 0.5)
                    126:       {
                    127:         while (d < (1.0 / 65536.0))
                    128:           {
                    129:             d *=  65536.0;
                    130:             exp -= 16;
                    131:           }
                    132:         while (d < 0.5)
                    133:           {
                    134:             d *= 2.0;
                    135:             exp -= 1;
                    136:           }
                    137:       }
                    138:
                    139:     d *= MP_BASE_AS_DOUBLE;
                    140: #if BITS_PER_MP_LIMB == 64
                    141:     manl = d;
                    142: #else
                    143:     manh = d;
                    144:     manl = (d - manh) * MP_BASE_AS_DOUBLE;
                    145: #endif
                    146:
                    147:     exp += 1022;
                    148:   }
                    149: #endif
                    150:
                    151:   if (exp) exp = (unsigned) exp - 1022; else exp = -1021;
                    152:
                    153: #if BITS_PER_MP_LIMB == 64
                    154:       rp[0] = manl;
                    155: #else
                    156:       if (e) {
                    157:        rp[1] = manh;
                    158:        rp[0] = manl;
                    159:       }
                    160:       else {
                    161:        rp[0] = manh;
                    162:       }
                    163: #endif
                    164:
                    165:   return exp;
                    166: }
                    167:
                    168: /* End of part included from gmp-2.0.2 */
                    169: /* Part included from gmp temporary releases */
                    170: double
                    171: #if __STDC__
                    172: __mpfr_scale2 (double d, int exp)
                    173: #else
                    174: __mpfr_scale2 (d, exp)
                    175:      double d;
                    176:      int exp;
                    177: #endif
                    178: {
                    179: #if _GMP_IEEE_FLOATS
                    180:   {
                    181:     union ieee_double_extract x;
                    182:     x.d = d;
                    183:     exp += x.s.exp;
                    184:     x.s.exp = exp;
                    185:     if (exp >= 2047)
                    186:       {
                    187:         /* Return +-infinity */
                    188:         x.s.exp = 2047;
                    189:         x.s.manl = x.s.manh = 0;
                    190:       }
                    191:     else if (exp < 1)
                    192:       {
                    193:         x.s.exp = 1;            /* smallest exponent (biased) */
                    194:         /* Divide result by 2 until we have scaled it to the right IEEE
                    195:            denormalized number, but stop if it becomes zero.  */
                    196:         while (exp < 1 && x.d != 0)
                    197:           {
                    198:             x.d *= 0.5;
                    199:             exp++;
                    200:           }
                    201:       }
                    202:     return x.d;
                    203:   }
                    204: #else
                    205:   {
                    206:     double factor, r;
                    207:
                    208:     factor = 2.0;
                    209:     if (exp < 0)
                    210:       {
                    211:         factor = 0.5;
                    212:         exp = -exp;
                    213:       }
                    214:     r = d;
                    215:     if (exp != 0)
                    216:       {
                    217:         if ((exp & 1) != 0)
                    218:           r *= factor;
                    219:         exp >>= 1;
                    220:         while (exp != 0)
                    221:           {
                    222:             factor *= factor;
                    223:             if ((exp & 1) != 0)
                    224:               r *= factor;
                    225:             exp >>= 1;
                    226:           }
                    227:       }
                    228:     return r;
                    229:   }
                    230: #endif
                    231: }
                    232:
                    233:
                    234: /* End of part included from gmp */
                    235:
                    236: void
                    237: #if __STDC__
                    238: mpfr_set_d(mpfr_t r, double d, unsigned char rnd_mode)
                    239: #else
                    240: mpfr_set_d(r, d, rnd_mode)
                    241:      mpfr_t r;
                    242:      double d;
                    243:      unsigned char rnd_mode;
                    244: #endif
                    245: {
                    246:   int signd, sizer; unsigned int cnt;
                    247:
                    248:   if (d == 0) { SET_ZERO(r); return; }
                    249:   else if (isnan(d)) { SET_NAN(r); return; }
                    250:
                    251:   signd = (d < 0) ? -1 : 1;
                    252:   d = ABS (d);
                    253:   sizer = (PREC(r)-1)/BITS_PER_MP_LIMB + 1;
                    254:
                    255:   /* warning: __mpfr_extract_double requires at least two limbs */
                    256:   if (sizer < MPFR_LIMBS_PER_DOUBLE)
                    257:     EXP(r) = __mpfr_extract_double (MANT(r), d, 0);
                    258:   else
                    259:     EXP(r) = __mpfr_extract_double (MANT(r) + sizer - MPFR_LIMBS_PER_DOUBLE, d, 1);
                    260:
                    261:   if (sizer > MPFR_LIMBS_PER_DOUBLE)
                    262:     MPN_ZERO(MANT(r), sizer - MPFR_LIMBS_PER_DOUBLE);
                    263:
                    264:   count_leading_zeros(cnt, MANT(r)[sizer-1]);
                    265:   if (cnt) mpn_lshift(MANT(r), MANT(r), sizer, cnt);
                    266:
                    267:   EXP(r) -= cnt;
                    268:   if (SIZE(r)*signd<0) CHANGE_SIGN(r);
                    269:
                    270:   mpfr_round(r, rnd_mode, PREC(r));
                    271:   return;
                    272: }
                    273:
                    274: double
                    275: #if __STDC__
                    276: mpfr_get_d2(mpfr_srcptr src, long e)
                    277: #else
                    278: mpfr_get_d2(src, e)
                    279:      mpfr_srcptr(src);
                    280:      long e;
                    281: #endif
                    282: {
                    283:   double res;
                    284:   mp_size_t size, i, n_limbs_to_use;
                    285:   mp_ptr qp;
                    286:   int negative;
                    287:
                    288:   if (FLAG_NAN(src)) {
                    289: #ifdef DEBUG
                    290:     printf("recognized NaN\n");
                    291: #endif
                    292:     return NaN; }
                    293:   if (NOTZERO(src)==0) return 0.0;
                    294:   size = 1+(PREC(src)-1)/BITS_PER_MP_LIMB;
                    295:   qp = MANT(src);
                    296:   negative = (SIGN(src)==-1);
                    297:
                    298:   /* Warning: don't compute the abs(res) and set the sign afterwards,
                    299:      otherwise the current machine rounding mode will not be taken
                    300:      correctly into account. */
                    301:   /* res = (negative) ? -(double)qp[size - 1] : qp[size - 1]; */
                    302:   res = 0.0;
                    303:   /* Warning: an arbitrary number of limbs may be required for an exact
                    304:      rounding. The following code is correct but not optimal since one
                    305:      may be able to decide without considering all limbs. */
                    306:   /* n_limbs_to_use = MIN (MPFR_LIMBS_PER_DOUBLE, size); */
                    307:   n_limbs_to_use = size;
                    308:   /* Accumulate the limbs from less significant to most significant
                    309:      otherwise due to rounding we may accumulate several ulps,
                    310:      especially in rounding towards -/+infinity. */
                    311:   for (i = n_limbs_to_use; i>=1; i--)
                    312:     res = res / MP_BASE_AS_DOUBLE +
                    313:       ((negative) ? -(double)qp[size - i] : qp[size - i]);
                    314:   res = __mpfr_scale2 (res, e - BITS_PER_MP_LIMB);
                    315:
                    316:   return res;
                    317: }
                    318:
                    319: double
                    320: #if __STDC__
                    321: mpfr_get_d(mpfr_srcptr src)
                    322: #else
                    323: mpfr_get_d(src)
                    324:      mpfr_srcptr src;
                    325: #endif
                    326: {
                    327:   return mpfr_get_d2(src, EXP(src));
                    328: }
                    329:

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