=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpfr/Attic/set_d.c,v retrieving revision 1.1 retrieving revision 1.1.1.2 diff -u -p -r1.1 -r1.1.1.2 --- OpenXM_contrib/gmp/mpfr/Attic/set_d.c 2000/09/09 14:12:19 1.1 +++ OpenXM_contrib/gmp/mpfr/Attic/set_d.c 2003/08/25 16:06:07 1.1.1.2 @@ -1,37 +1,41 @@ -/* mpfr_set_d, mpfr_get_d -- convert a multiple precision floating-point number - from/to a machine double precision float +/* mpfr_set_d -- convert a machine double precision float to + a multiple precision floating-point number -Copyright (C) 1999 PolKA project, Inria Lorraine and Loria +Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. This file is part of the MPFR Library. The MPFR Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. -You should have received a copy of the GNU Library General Public License +You should have received a copy of the GNU Lesser General Public License along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#if __GNUC__ /* gcc "patched" headers seem to omit isnan... */ -extern int isnan(double); -#endif -#include /* for isnan and NaN */ - #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" #include "mpfr.h" +#include "mpfr-impl.h" -#define NaN sqrt(-1) /* ensures a machine-independent NaN */ +#if (BITS_PER_MP_LIMB==32) +#define MPFR_LIMBS_PER_DOUBLE 2 +#elif (BITS_PER_MP_LIMB >= 64) +#define MPFR_LIMBS_PER_DOUBLE 1 +#elif (BITS_PER_MP_LIMB == 16) +#define MPFR_LIMBS_PER_DOUBLE 4 +#endif +static int __mpfr_extract_double _PROTO ((mp_ptr, double)); + /* Included from gmp-2.0.2, patched to support denorms */ #ifdef XDEBUG @@ -42,15 +46,8 @@ extern int isnan(double); #define _GMP_IEEE_FLOATS 0 #endif -int -#if __STDC__ -__mpfr_extract_double (mp_ptr rp, double d, int e) -#else -__mpfr_extract_double (rp, d, e) - mp_ptr rp; - double d; - int e; -#endif +static int +__mpfr_extract_double (mp_ptr rp, double d) /* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */ { long exp; @@ -70,9 +67,6 @@ __mpfr_extract_double (rp, d, e) if (d == 0.0) { rp[0] = 0; -#if BITS_PER_MP_LIMB == 32 - if (e) rp[1] = 0; -#endif return 0; } @@ -82,23 +76,24 @@ __mpfr_extract_double (rp, d, e) x.d = d; exp = x.s.exp; - if (exp) + if (exp) { #if BITS_PER_MP_LIMB == 64 - manl = (((mp_limb_t) 1 << 63) + manl = ((MP_LIMB_T_ONE << 63) | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11)); #else - manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21); - manl = x.s.manl << 11; + manh = (MP_LIMB_T_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21); + manl = x.s.manl << 11; #endif } - else + else /* denormalized number */ { #if BITS_PER_MP_LIMB == 64 manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11); #else - manh = (x.s.manh << 11) | (x.s.manl >> 21); - manl = x.s.manl << 11; + manh = (x.s.manh << 11) /* high 21 bits */ + | (x.s.manl >> 21); /* middle 11 bits */ + manl = x.s.manl << 11; /* low 21 bits */ #endif } } @@ -148,182 +143,87 @@ __mpfr_extract_double (rp, d, e) } #endif - if (exp) exp = (unsigned) exp - 1022; else exp = -1021; + if (exp) exp = (unsigned) exp - 1022; else exp = -1021; #if BITS_PER_MP_LIMB == 64 - rp[0] = manl; + rp[0] = manl; #else - if (e) { - rp[1] = manh; - rp[0] = manl; - } - else { - rp[0] = manh; - } + rp[1] = manh; + rp[0] = manl; #endif return exp; } /* End of part included from gmp-2.0.2 */ -/* Part included from gmp temporary releases */ -double -#if __STDC__ -__mpfr_scale2 (double d, int exp) -#else -__mpfr_scale2 (d, exp) - double d; - int exp; -#endif + +int +mpfr_set_d (mpfr_ptr r, double d, mp_rnd_t rnd_mode) { -#if _GMP_IEEE_FLOATS - { - union ieee_double_extract x; - x.d = d; - exp += x.s.exp; - x.s.exp = exp; - if (exp >= 2047) - { - /* Return +-infinity */ - x.s.exp = 2047; - x.s.manl = x.s.manh = 0; - } - else if (exp < 1) - { - x.s.exp = 1; /* smallest exponent (biased) */ - /* Divide result by 2 until we have scaled it to the right IEEE - denormalized number, but stop if it becomes zero. */ - while (exp < 1 && x.d != 0) - { - x.d *= 0.5; - exp++; - } - } - return x.d; - } -#else - { - double factor, r; + int signd, sizetmp, inexact; + unsigned int cnt, k; + mpfr_ptr tmp; + TMP_DECL(marker); - factor = 2.0; - if (exp < 0) - { - factor = 0.5; - exp = -exp; - } - r = d; - if (exp != 0) - { - if ((exp & 1) != 0) - r *= factor; - exp >>= 1; - while (exp != 0) - { - factor *= factor; - if ((exp & 1) != 0) - r *= factor; - exp >>= 1; - } - } - return r; - } -#endif -} + TMP_MARK(marker); + MPFR_CLEAR_FLAGS(r); + if (d == 0) + { + union ieee_double_extract x; -/* End of part included from gmp */ + MPFR_SET_ZERO(r); + /* set correct sign */ + x.d = d; + if (((x.s.sig == 1) && (MPFR_SIGN(r) > 0)) + || ((x.s.sig == 0) && (MPFR_SIGN(r) < 0))) + MPFR_CHANGE_SIGN(r); + return 0; /* 0 is exact */ + } + else if (DOUBLE_ISNAN(d)) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + else if (DOUBLE_ISINF(d)) + { + MPFR_SET_INF(r); + if ((d > 0 && (MPFR_SIGN(r) == -1)) || (d < 0 && (MPFR_SIGN(r) == 1))) + MPFR_CHANGE_SIGN(r); + return 0; /* infinity is exact */ + } -void -#if __STDC__ -mpfr_set_d(mpfr_t r, double d, unsigned char rnd_mode) -#else -mpfr_set_d(r, d, rnd_mode) - mpfr_t r; - double d; - unsigned char rnd_mode; -#endif -{ - int signd, sizer; unsigned int cnt; + /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE, + since PREC(r) may be different from PREC(tmp), and then both variables + would have same precision in the mpfr_set4 call below. */ + tmp = (mpfr_ptr) TMP_ALLOC(sizeof(mpfr_t)); + MPFR_MANT(tmp) = TMP_ALLOC(MPFR_LIMBS_PER_DOUBLE * BYTES_PER_MP_LIMB); + MPFR_PREC(tmp) = 53; + MPFR_SIZE(tmp) = MPFR_LIMBS_PER_DOUBLE; + sizetmp = MPFR_LIMBS_PER_DOUBLE; - if (d == 0) { SET_ZERO(r); return; } - else if (isnan(d)) { SET_NAN(r); return; } - signd = (d < 0) ? -1 : 1; d = ABS (d); - sizer = (PREC(r)-1)/BITS_PER_MP_LIMB + 1; - /* warning: __mpfr_extract_double requires at least two limbs */ - if (sizer < MPFR_LIMBS_PER_DOUBLE) - EXP(r) = __mpfr_extract_double (MANT(r), d, 0); - else - EXP(r) = __mpfr_extract_double (MANT(r) + sizer - MPFR_LIMBS_PER_DOUBLE, d, 1); - - if (sizer > MPFR_LIMBS_PER_DOUBLE) - MPN_ZERO(MANT(r), sizer - MPFR_LIMBS_PER_DOUBLE); + MPFR_EXP(tmp) = __mpfr_extract_double (MPFR_MANT(tmp), d); - count_leading_zeros(cnt, MANT(r)[sizer-1]); - if (cnt) mpn_lshift(MANT(r), MANT(r), sizer, cnt); - - EXP(r) -= cnt; - if (SIZE(r)*signd<0) CHANGE_SIGN(r); + /* determine number k of zero high limbs */ + for (k = 0; k < sizetmp && MPFR_MANT(tmp)[sizetmp - 1 - k] == 0; k++); - mpfr_round(r, rnd_mode, PREC(r)); - return; -} + count_leading_zeros (cnt, MPFR_MANT(tmp)[sizetmp - 1 - k]); -double -#if __STDC__ -mpfr_get_d2(mpfr_srcptr src, long e) -#else -mpfr_get_d2(src, e) - mpfr_srcptr(src); - long e; -#endif -{ - double res; - mp_size_t size, i, n_limbs_to_use; - mp_ptr qp; - int negative; + if (cnt) + mpn_lshift (MPFR_MANT(tmp) + k, MPFR_MANT(tmp), sizetmp - k, cnt); + else if (k) + MPN_COPY (MPFR_MANT(tmp) + k, MPFR_MANT(tmp), sizetmp - k); + if (k) + MPN_ZERO (MPFR_MANT(tmp), k); - if (FLAG_NAN(src)) { -#ifdef DEBUG - printf("recognized NaN\n"); -#endif - return NaN; } - if (NOTZERO(src)==0) return 0.0; - size = 1+(PREC(src)-1)/BITS_PER_MP_LIMB; - qp = MANT(src); - negative = (SIGN(src)==-1); + MPFR_EXP(tmp) -= cnt + k * BITS_PER_MP_LIMB; - /* Warning: don't compute the abs(res) and set the sign afterwards, - otherwise the current machine rounding mode will not be taken - correctly into account. */ - /* res = (negative) ? -(double)qp[size - 1] : qp[size - 1]; */ - res = 0.0; - /* Warning: an arbitrary number of limbs may be required for an exact - rounding. The following code is correct but not optimal since one - may be able to decide without considering all limbs. */ - /* n_limbs_to_use = MIN (MPFR_LIMBS_PER_DOUBLE, size); */ - n_limbs_to_use = size; - /* Accumulate the limbs from less significant to most significant - otherwise due to rounding we may accumulate several ulps, - especially in rounding towards -/+infinity. */ - for (i = n_limbs_to_use; i>=1; i--) - res = res / MP_BASE_AS_DOUBLE + - ((negative) ? -(double)qp[size - i] : qp[size - i]); - res = __mpfr_scale2 (res, e - BITS_PER_MP_LIMB); + /* tmp is exact since PREC(tmp)=53 */ + inexact = mpfr_set4 (r, tmp, rnd_mode, signd); - return res; + TMP_FREE(marker); + return inexact; } - -double -#if __STDC__ -mpfr_get_d(mpfr_srcptr src) -#else -mpfr_get_d(src) - mpfr_srcptr src; -#endif -{ - return mpfr_get_d2(src, EXP(src)); -} -