From b2f4e6bf9660bc322f5f1145d77d26ba41a153fd Mon Sep 17 00:00:00 2001 From: pelissip Date: Mon, 1 Dec 2003 10:25:01 +0000 Subject: Port all the IEEE dependent remaining functions so that it uses a generic way to compute the result if it detects that the double is not in IEEE format ( _GMP_IEEE_FLOATS == 0). MPFR now should work well on non-IEEE machines. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2567 280ebfd0-de03-0410-8827-d642c229c3f4 --- exceptions.c | 4 ++-- get_d.c | 8 -------- mpfr-impl.h | 22 +++++++++++++++++++--- neg.c | 2 +- powerof2.c | 2 +- set_d.c | 10 ---------- set_ld.c | 58 +++++++++++++++++++++++++++++----------------------------- uceil_exp2.c | 19 +++++++++++++++++++ uceil_log2.c | 27 ++++++++++++++++++++++++++- ufloor_log2.c | 27 ++++++++++++++++++++++++++- 10 files changed, 123 insertions(+), 56 deletions(-) diff --git a/exceptions.c b/exceptions.c index d35e6c10d..93b37d33a 100644 --- a/exceptions.c +++ b/exceptions.c @@ -179,7 +179,7 @@ mpfr_inexflag_p (void) return __gmpfr_flags & MPFR_FLAGS_INEXACT; } -#undef mpfr_set_underflow +/* #undef mpfr_set_underflow */ /* Note: In the rounding to the nearest mode, mpfr_set_underflow always rounds away from 0. In this rounding mode, you must call @@ -209,7 +209,7 @@ mpfr_set_underflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign) return sign > 0 ? inex : -inex; } -#undef mpfr_set_overflow +/* #undef mpfr_set_overflow */ int mpfr_set_overflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign) diff --git a/get_d.c b/get_d.c index adca12929..0dcdf0358 100644 --- a/get_d.c +++ b/get_d.c @@ -28,14 +28,6 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -#ifdef XDEBUG -#undef _GMP_IEEE_FLOATS -#endif - -#ifndef _GMP_IEEE_FLOATS -#define _GMP_IEEE_FLOATS 0 -#endif - /* "double" NaN and infinities are written as explicit bytes to be sure of getting what we want, and to be sure of not depending on libm. diff --git a/mpfr-impl.h b/mpfr-impl.h index 0a8878499..03a5293ba 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -264,6 +264,14 @@ long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) ATTRIBUTE_CO # endif #endif +/* Debug non IEEE floats */ +#ifdef XDEBUG +# undef _GMP_IEEE_FLOATS +#endif +#ifndef _GMP_IEEE_FLOATS +# define _GMP_IEEE_FLOATS 0 +#endif + /* We want to test this : * (rnd == GMP_RNDU && test) || (rnd == RNDD && !test) * This macro does this test faster*/ @@ -421,28 +429,36 @@ void mpfr_inits2 _MPFR_PROTO ((mp_prec_t, mpfr_ptr, ...)); void mpfr_inits _MPFR_PROTO ((mpfr_ptr, ...)); void mpfr_clears _MPFR_PROTO ((mpfr_ptr, ...)); + int mpfr_set_underflow _MPFR_PROTO ((mpfr_ptr, mp_rnd_t, int)); int mpfr_set_overflow _MPFR_PROTO ((mpfr_ptr, mp_rnd_t, int)); void mpfr_save_emin_emax _MPFR_PROTO ((void)); void mpfr_restore_emin_emax _MPFR_PROTO ((void)); + int mpfr_add1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_sub1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)); int mpfr_can_round_raw _MPFR_PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t, mp_rnd_t, mp_rnd_t, mp_prec_t)); + double mpfr_get_d3 _MPFR_PROTO ((mpfr_srcptr, mp_exp_t, mp_rnd_t)); int mpfr_cmp2 _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr, mp_prec_t *)); -long __gmpfr_ceil_log2 _MPFR_PROTO ((double)); -long __gmpfr_floor_log2 _MPFR_PROTO ((double)); -double __gmpfr_ceil_exp2 _MPFR_PROTO ((double)); + +long __gmpfr_ceil_log2 _MPFR_PROTO ((double)); +long __gmpfr_floor_log2 _MPFR_PROTO ((double)); +double __gmpfr_ceil_exp2 _MPFR_PROTO ((double)); unsigned long __gmpfr_isqrt _MPFR_PROTO ((unsigned long)); unsigned long __gmpfr_cuberoot _MPFR_PROTO ((unsigned long)); + int mpfr_exp_2 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); int mpfr_exp3 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); int mpfr_powerof2_raw _MPFR_PROTO ((mpfr_srcptr)); + void mpfr_setmax _MPFR_PROTO ((mpfr_ptr, mp_exp_t)); void mpfr_setmin _MPFR_PROTO ((mpfr_ptr, mp_exp_t)); + long mpn_exp _MPFR_PROTO ((mp_limb_t *, mp_exp_t *, int, mp_exp_t, size_t)); + void mpfr_print_binary _MPFR_PROTO ((mpfr_srcptr)); void mpfr_set_str_binary _MPFR_PROTO ((mpfr_ptr, __gmp_const char *)); diff --git a/neg.c b/neg.c index 185a836c1..a4f6c40cf 100644 --- a/neg.c +++ b/neg.c @@ -27,7 +27,7 @@ MA 02111-1307, USA. */ int mpfr_neg (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode) { - if (a != b) + if (MPFR_UNLIKELY(a != b)) return mpfr_set4 (a, b, rnd_mode, -MPFR_SIGN(b)); else { diff --git a/powerof2.c b/powerof2.c index f4bf0d3c6..f549fdfa3 100644 --- a/powerof2.c +++ b/powerof2.c @@ -32,7 +32,7 @@ mpfr_powerof2_raw (mpfr_srcptr x) mp_limb_t *xp; mp_size_t xn; - MPFR_ASSERTN(MPFR_IS_FP(x)); + MPFR_ASSERTN(MPFR_IS_PURE_FP(x)); xp = MPFR_MANT(x); xn = (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; if (NOT_POW2(xp[xn])) diff --git a/set_d.c b/set_d.c index 7e5672b14..79cd23af9 100644 --- a/set_d.c +++ b/set_d.c @@ -34,16 +34,6 @@ MA 02111-1307, USA. */ #error "Unsupported value of BITS_PER_MP_LIMB" #endif -/* Included from gmp-2.0.2, patched to support denorms */ - -#ifdef XDEBUG -#undef _GMP_IEEE_FLOATS -#endif - -#ifndef _GMP_IEEE_FLOATS -#define _GMP_IEEE_FLOATS 0 -#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 */ diff --git a/set_ld.c b/set_ld.c index 055be494d..b99d7188d 100644 --- a/set_ld.c +++ b/set_ld.c @@ -120,35 +120,35 @@ mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode) /* since -DBL_MAX <= d <= DBL_MAX, the cast to double should not overflow here */ mpfr_set_d (u, (double) d, GMP_RNDN); - /* warning: using MPFR_IS_ZERO will cause a READ_UNINIT_MEM if u=Inf */ - if (mpfr_cmp_ui (u, 0) == 0 && (d != (long double) 0.0)) /* underflow */ - { - long double div10, div11, div12, div13; - div10 = (long double) (double) 5.5626846462680034577255e-309; /* 2^(-2^10) */ - div11 = div10 * div10; /* 2^(-2^11) */ - div12 = div11 * div11; /* 2^(-2^12) */ - div13 = div12 * div12; /* 2^(-2^13) */ - if (ABS(d) <= div13) - { - d = d / div13; /* exact */ - shift_exp -= 8192; - } - if (ABS(d) <= div12) - { - d = d / div12; /* exact */ - shift_exp -= 4096; - } - if (ABS(d) <= div11) - { - d = d / div11; /* exact */ - shift_exp -= 2048; - } - if (ABS(d) <= div10) - { - d = d / div10; /* exact */ - shift_exp -= 1024; - } - } + /* warning: using MPFR_IS_ZERO will cause a READ_UNINIT_MEM if u=Inf */ + if (mpfr_cmp_ui (u, 0) == 0 && (d != (long double) 0.0)) /* underflow */ + { + long double div10, div11, div12, div13; + div10 = (long double) (double) 5.5626846462680034577255e-309; /* 2^(-2^10) */ + div11 = div10 * div10; /* 2^(-2^11) */ + div12 = div11 * div11; /* 2^(-2^12) */ + div13 = div12 * div12; /* 2^(-2^13) */ + if (ABS(d) <= div13) + { + d = d / div13; /* exact */ + shift_exp -= 8192; + } + if (ABS(d) <= div12) + { + d = d / div12; /* exact */ + shift_exp -= 4096; + } + if (ABS(d) <= div11) + { + d = d / div11; /* exact */ + shift_exp -= 2048; + } + if (ABS(d) <= div10) + { + d = d / div10; /* exact */ + shift_exp -= 1024; + } + } } mpfr_add (t, t, u, GMP_RNDN); /* exact */ if (!mpfr_number_p (t)) diff --git a/uceil_exp2.c b/uceil_exp2.c index 8003a5164..16ad030d7 100644 --- a/uceil_exp2.c +++ b/uceil_exp2.c @@ -37,6 +37,25 @@ __gmpfr_ceil_exp2 (double d) exp++; /* now exp = ceil(d) */ x.d = 1.0; +#if _GMP_IEEE_FLOATS x.s.exp = exp <= -1022 ? 1 : 1023 + exp; +#else + if (exp >= 0) + { + while (exp != 0) + { + x.d *= 2.0; + exp--; + } + } + else + { + while (exp != 0) + { + x.d *= (1.0 / 2.0); + exp++; + } + } +#endif return x.d; } diff --git a/uceil_log2.c b/uceil_log2.c index f7ba603cd..ab3bcd8ab 100644 --- a/uceil_log2.c +++ b/uceil_log2.c @@ -24,11 +24,14 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -/* returns ceil(log(d)/log(2)) */ +/* returns ceil(log(d)/log(2)) if d > 0, + -1023 if d = +0, + and floor(log(-d)/log(2))+1 if d < 0*/ long __gmpfr_ceil_log2 (double d) { long exp; +#if _GMP_IEEE_FLOATS union ieee_double_extract x; x.d = d; @@ -37,4 +40,26 @@ __gmpfr_ceil_log2 (double d) if (x.d != 1.0) /* d: not a power of two? */ exp++; return exp; +#else + double m; + + if (d < 0.0) + return __gmpfr_floor_log2(-d)+1; + else if (d == 0.0) + return -1023; + else if (d >= 1.0) + { + exp = 0; + for( m= 1.0 ; m < d ; m *=2.0 ) + exp++; + } + else + { + exp = 1; + for( m= 1.0 ; m >= d ; m *= (1.0/2.0) ) + exp--; + } +#endif + return exp; } + diff --git a/ufloor_log2.c b/ufloor_log2.c index 7d8e05494..4b9c4a8a5 100644 --- a/ufloor_log2.c +++ b/ufloor_log2.c @@ -24,12 +24,37 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -/* returns floor(log(d)/log(2)) */ +/* returns floor(log(abs(d))/log(2)) */ long __gmpfr_floor_log2 (double d) { +#if _GMP_IEEE_FLOATS union ieee_double_extract x; x.d = d; return (long) x.s.exp - 1023; +#else + long exp; + double m; + + /* Get Abs */ + if (d < 0.0) + d = -d; + + if (d == 0.0) + return -1023; + else if (d >= 1.0) + { + exp = -1; + for( m= 1.0 ; m <= d ; m *=2.0 ) + exp++; + } + else + { + exp = 0; + for( m= 1.0 ; m > d ; m *= (1.0/2.0) ) + exp--; + } + return exp; +#endif } -- cgit v1.2.1