-- cgit v1.2.1 From 2ad3badd302e4ce502e1accbd8efe0763f9ffbc6 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 12:02:22 +0000 Subject: First, reverse-merge r9975 and r9957 (whose only purpose was for fmma, but changed the exponent range and still had various issues). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10315 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/fmma.c | 18 ++++++++---------- src/mpfr-impl.h | 13 ++++++------- tests/tsprintf.c | 8 -------- 3 files changed, 14 insertions(+), 25 deletions(-) diff --git a/src/fmma.c b/src/fmma.c index 9f0af2faf..0fda41bbe 100644 --- a/src/fmma.c +++ b/src/fmma.c @@ -79,16 +79,12 @@ mpfr_fmma_fast (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, if ((up[an + bn - 1] & MPFR_LIMB_HIGHBIT) == 0) { mpn_lshift (up, up, an + bn, 1); - /* EXP(a) and EXP(b) are in [1-2^(n-2), 2^(n-2)-1] where - mpfr_exp_t has is n-bit wide, thus EXP(a)+EXP(b) is in - [2-2^(n-1), 2^(n-1)-2]. We use the fact here that 2*MPFR_EMIN_MIN-1 - is a valid exponent (see mpfr-impl.h). However we don't use - MPFR_SET_EXP() which only allows the restricted exponent range - [1-2^(n-2), 2^(n-2)-1]. */ - MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b) - 1; + MPFR_SET_EXP (u, MPFR_EXP(a) + MPFR_EXP(b) - 1); + if (MPFR_UNLIKELY(MPFR_EXP(u) == __MPFR_EXP_INF)) + goto failure; } else - MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b); + MPFR_SET_EXP (u, MPFR_EXP(a) + MPFR_EXP(b)); /* v <- c*d */ if (cn >= dn) @@ -98,10 +94,12 @@ mpfr_fmma_fast (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, if ((vp[cn + dn - 1] & MPFR_LIMB_HIGHBIT) == 0) { mpn_lshift (vp, vp, cn + dn, 1); - MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d) - 1; + MPFR_SET_EXP (v, MPFR_EXP(c) + MPFR_EXP(d) - 1); + if (MPFR_UNLIKELY(MPFR_EXP(v) == __MPFR_EXP_INF)) + goto failure; } else - MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d); + MPFR_SET_EXP (v, MPFR_EXP(c) + MPFR_EXP(d)); MPFR_PREC(u) = (an + bn) * GMP_NUMB_BITS; MPFR_PREC(v) = (cn + dn) * GMP_NUMB_BITS; diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index 3f7eb1b71..5d989f489 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -830,18 +830,17 @@ typedef intmax_t mpfr_eexp_t; /* Definition of the exponent limits for MPFR numbers. * These limits are chosen so that if e is such an exponent, then 2e-1 and - * 2e+1 are valid exponents. This is useful for intermediate computations, - * in particular the multiplication. We must have MPFR_EMIN_MIN >= 3-2^(n-2) - * = 3-MPFR_EXP_INVALID so that 2*MPFR_EMIN_MIN-1 > __MPFR_EXP_INF = 3-2^(n-1). + * 2e+1 are representable. This is useful for intermediate computations, + * in particular the multiplication. */ #undef MPFR_EMIN_MIN #undef MPFR_EMIN_MAX #undef MPFR_EMAX_MIN #undef MPFR_EMAX_MAX -#define MPFR_EMIN_MIN (3-MPFR_EXP_INVALID) -#define MPFR_EMIN_MAX (MPFR_EXP_INVALID-3) -#define MPFR_EMAX_MIN (3-MPFR_EXP_INVALID) -#define MPFR_EMAX_MAX (MPFR_EXP_INVALID-3) +#define MPFR_EMIN_MIN (1-MPFR_EXP_INVALID) +#define MPFR_EMIN_MAX (MPFR_EXP_INVALID-1) +#define MPFR_EMAX_MIN (1-MPFR_EXP_INVALID) +#define MPFR_EMAX_MAX (MPFR_EXP_INVALID-1) /* Use MPFR_GET_EXP and MPFR_SET_EXP instead of MPFR_EXP directly, unless when the exponent may be out-of-range, for instance when diff --git a/tests/tsprintf.c b/tests/tsprintf.c index fba32511e..b33355c88 100644 --- a/tests/tsprintf.c +++ b/tests/tsprintf.c @@ -1128,7 +1128,6 @@ bug20111102 (void) * for %Ra and %Rb is not done on the MPFR number itself (as it * would overflow). Note: it has been reported on comp.std.c that * some C libraries behave differently on %a, but this is a bug. - * FIXME: this function assumes e = 3 mod 4. */ static void check_emax_aux (mpfr_exp_t e) @@ -1138,9 +1137,6 @@ check_emax_aux (mpfr_exp_t e) int i; mpfr_exp_t emax; - if ((e % 4) != 3) - return; - MPFR_ASSERTN (e <= LONG_MAX); emax = mpfr_get_emax (); set_emax (e); @@ -1199,7 +1195,6 @@ check_emax (void) check_emax_aux (MPFR_EMAX_MAX); } -/* FIXME: this function assumes e = 1 mod 4 */ static void check_emin_aux (mpfr_exp_t e) { @@ -1209,9 +1204,6 @@ check_emin_aux (mpfr_exp_t e) mpfr_exp_t emin; mpz_t ee; - if ((e % 4) != 1) - return; - MPFR_ASSERTN (e >= LONG_MIN); emin = mpfr_get_emin (); set_emin (e); -- cgit v1.2.1 From 6d57636b1661601b00f4c707443b0be00b9e02ef Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 12:23:14 +0000 Subject: Also reverse-merge r9958, r9961 and r9962 (due to the previous ones). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10316 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/mpfr.h | 2 +- tests/tacosh.c | 11 ----------- tests/texp.c | 31 +++++++++++++------------------ tests/tpow.c | 7 ------- tests/tpow_z.c | 30 ++++++++++++------------------ tests/tsub1sp.c | 7 ------- 6 files changed, 26 insertions(+), 62 deletions(-) diff --git a/src/mpfr.h b/src/mpfr.h index 609807b2b..52c846283 100644 --- a/src/mpfr.h +++ b/src/mpfr.h @@ -187,7 +187,7 @@ typedef uintmax_t mpfr_uexp_t; #endif /* Definition of the standard exponent limits */ -#define MPFR_EMAX_DEFAULT ((mpfr_exp_t) (((mpfr_ulong) 1 << 30) - 3)) +#define MPFR_EMAX_DEFAULT ((mpfr_exp_t) (((mpfr_ulong) 1 << 30) - 1)) #define MPFR_EMIN_DEFAULT (-(MPFR_EMAX_DEFAULT)) /* DON'T USE THIS! (For MPFR-public macros only, see below.) diff --git a/tests/tacosh.c b/tests/tacosh.c index d10dbe4f2..fbabc3098 100644 --- a/tests/tacosh.c +++ b/tests/tacosh.c @@ -176,17 +176,7 @@ huge (void) { mpfr_t x, y, z; int inex; - mpfr_exp_t emax; - if (mpfr_get_emax_max () < 1073741823) - return; - - emax = mpfr_get_emax (); - mpfr_set_emax (1073741823); - - /* FIXME: The main purpose of this test was on 32-bit ABI, - but it is no longer run there. Solution: implement the - TODO below. */ /* TODO: extend the exponent range and use mpfr_get_emax (). */ mpfr_inits2 (32, x, y, z, (mpfr_ptr) 0); mpfr_set_ui_2exp (x, 1, 1073741822, MPFR_RNDN); @@ -203,7 +193,6 @@ huge (void) MPFR_ASSERTN (inex < 0); mpfr_clears (x, y, z, (mpfr_ptr) 0); - mpfr_set_emax (emax); } int diff --git a/tests/texp.c b/tests/texp.c index 0e85589b4..e0c9e26f2 100644 --- a/tests/texp.c +++ b/tests/texp.c @@ -273,26 +273,21 @@ check_special (void) } /* Check overflow. Corner case of mpfr_exp_2 */ - /* FIXME: The main purpose of the test below was on 32-bit ABI, - but it is no longer run there. */ mpfr_set_prec (x, 64); - if (mpfr_set_emax (1073741823) == 0) - { /* 1073741823 is in the allowed exponent range */ - mpfr_set_emin (MPFR_EMIN_DEFAULT); - mpfr_set_str (x, - "0.1011000101110010000101111111010100001100000001110001100111001101E30", - 2, MPFR_RNDN); - mpfr_exp (x, x, MPFR_RNDD); - if (mpfr_cmp_str (x, - ".1111111111111111111111111111111111111111111111111111111111111111E1073741823", - 2, MPFR_RNDN) != 0) - { - printf ("Wrong overflow detection in mpfr_exp\n"); - mpfr_dump (x); - exit (1); - } + mpfr_set_emax (MPFR_EMAX_DEFAULT); + mpfr_set_emin (MPFR_EMIN_DEFAULT); + mpfr_set_str (x, + "0.1011000101110010000101111111010100001100000001110001100111001101E30", + 2, MPFR_RNDN); + mpfr_exp (x, x, MPFR_RNDD); + if (mpfr_cmp_str (x, +".1111111111111111111111111111111111111111111111111111111111111111E1073741823", + 2, MPFR_RNDN) != 0) + { + printf ("Wrong overflow detection in mpfr_exp\n"); + mpfr_dump (x); + exit (1); } - /* Check underflow. Corner case of mpfr_exp_2 */ mpfr_set_str (x, "-0.1011000101110010000101111111011111010001110011110111100110101100E30", diff --git a/tests/tpow.c b/tests/tpow.c index 3bcf7a9a7..13819ff07 100644 --- a/tests/tpow.c +++ b/tests/tpow.c @@ -1406,13 +1406,7 @@ bug20071218 (void) { mpfr_t x, y, z, t; int tern; - mpfr_exp_t emin; - - if (mpfr_get_emin_min () > -1073741823) - return; - emin = mpfr_get_emin (); - mpfr_set_emin (-1073741823); mpfr_inits2 (64, x, y, z, t, (mpfr_ptr) 0); mpfr_set_str (x, "0x.80000000000002P-1023", 0, MPFR_RNDN); mpfr_set_str (y, "100000.000000002", 16, MPFR_RNDN); @@ -1448,7 +1442,6 @@ bug20071218 (void) exit (1); } mpfr_clears (x, y, z, t, (mpfr_ptr) 0); - mpfr_set_emin (emin); } /* With revision 5429, this gives: diff --git a/tests/tpow_z.c b/tests/tpow_z.c index 1514bd041..004435473 100644 --- a/tests/tpow_z.c +++ b/tests/tpow_z.c @@ -302,28 +302,22 @@ static void bug20080223 (void) { mpfr_t a, exp, answer; - mpfr_exp_t emin = mpfr_get_emin (); - mpfr_set_emin (mpfr_get_emin_min ()); - if (mpfr_get_emin () <= -1073741823) - { - mpfr_init2 (a, 53); - mpfr_init2 (exp, 53); - mpfr_init2 (answer, 53); + mpfr_init2 (a, 53); + mpfr_init2 (exp, 53); + mpfr_init2 (answer, 53); - mpfr_set_si (exp, -1073741824, MPFR_RNDN); + mpfr_set_si (exp, -1073741824, MPFR_RNDN); - mpfr_set_str (a, "1.999999999", 10, MPFR_RNDN); - /* a = 562949953139837/2^48 */ - mpfr_pow (answer, a, exp, MPFR_RNDN); - mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823"); - MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0); + mpfr_set_str (a, "1.999999999", 10, MPFR_RNDN); + /* a = 562949953139837/2^48 */ + mpfr_pow (answer, a, exp, MPFR_RNDN); + mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823"); + MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0); - mpfr_clear (a); - mpfr_clear (exp); - mpfr_clear (answer); - } - mpfr_set_emin (emin); + mpfr_clear (a); + mpfr_clear (exp); + mpfr_clear (answer); } static void diff --git a/tests/tsub1sp.c b/tests/tsub1sp.c index 460a5dd6f..b191e350f 100644 --- a/tests/tsub1sp.c +++ b/tests/tsub1sp.c @@ -351,19 +351,12 @@ check_special (void) MPFR_SET_NAN(x); MPFR_SET_NAN(x2); -#if 0 mpfr_set_str_binary (y, "0.1000000000000000000000000000000000000000000000000000000000000000" "E-1073741823"); mpfr_set_str_binary (z, "0.1100000000000000000000000000000000000000000000000000000000000000" "E-1073741823"); -#else - mpfr_set_ui_2exp (y, 1, MPFR_EMIN_DEFAULT-1, MPFR_RNDN); - /* y = 0.5*2^EMIN_DEFAULT */ - mpfr_set_ui_2exp (z, 3, MPFR_EMIN_DEFAULT-2, MPFR_RNDN); - /* z = 0.75*2^EMIN_DEFAULT */ -#endif inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); if (mpfr_cmp(x, x2)) -- cgit v1.2.1 From 4096945ad4bdf074618ceccd7b1fcb19ae40ec8c Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 12:24:44 +0000 Subject: Also reverse-merge r9960 (due to the previous ones). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10317 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/fmma.c | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/fmma.c b/src/fmma.c index 0fda41bbe..49fd9742e 100644 --- a/src/fmma.c +++ b/src/fmma.c @@ -57,6 +57,7 @@ mpfr_fmma_fast (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_t u, v; mp_size_t an, bn, cn, dn; mpfr_limb_ptr up, vp; + unsigned int saved_flags = __gmpfr_flags; MPFR_TMP_DECL(marker); MPFR_SAVE_EXPO_DECL (expo); @@ -118,6 +119,12 @@ mpfr_fmma_fast (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, MPFR_SAVE_EXPO_FREE (expo); return mpfr_check_range (z, inex, rnd); + + failure: + __gmpfr_flags = saved_flags; + MPFR_TMP_FREE(marker); + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_fmma_slow (z, a, b, c, d, rnd); } /* z <- a*b + c*d */ -- cgit v1.2.1 From 06e0bc4ea4af87cd171ef562e61b8efbf4f7b174 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 12:33:11 +0000 Subject: Started to implement unbounded floats (UBF) and added support in some existing functions. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10318 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/Makefile.am | 2 +- src/add.c | 8 +-- src/fmma.c | 144 +++++++++++---------------------------------- src/mpfr-impl.h | 59 +++++++++++++++++++ src/print_raw.c | 15 ++--- src/sub.c | 6 +- src/ubf.c | 177 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/vasprintf.c | 15 +++++ 8 files changed, 300 insertions(+), 126 deletions(-) create mode 100644 src/ubf.c diff --git a/src/Makefile.am b/src/Makefile.am index fae6e76b5..a9114ef98 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -59,7 +59,7 @@ buildopt.c digamma.c bernoulli.c isregular.c set_flt.c get_flt.c \ scale2.c set_z_exp.c ai.c gammaonethird.c ieee_floats.h \ grandom.c fpif.c set_float128.c get_float128.c rndna.c nrandom.c \ random_deviate.h random_deviate.c erandom.c mpfr-mini-gmp.c \ -mpfr-mini-gmp.h fmma.c log_ui.c gamma_inc.c +mpfr-mini-gmp.h fmma.c log_ui.c gamma_inc.c ubf.c libmpfr_la_LIBADD = @LIBOBJS@ diff --git a/src/add.c b/src/add.c index a953952d2..43a4507d1 100644 --- a/src/add.c +++ b/src/add.c @@ -85,8 +85,8 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } } - MPFR_ASSERTD (MPFR_IS_PURE_FP (b)); - MPFR_ASSERTD (MPFR_IS_PURE_FP (c)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); if (MPFR_UNLIKELY(MPFR_SIGN(b) != MPFR_SIGN(c))) { /* signs differ, it is a subtraction */ @@ -100,12 +100,12 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) { /* signs are equal, it's an addition */ if (MPFR_LIKELY(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c))) - if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) + if (MPFR_EXP_LESS_P (b, c)) return mpfr_add1sp(a, c, b, rnd_mode); else return mpfr_add1sp(a, b, c, rnd_mode); else - if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) + if (MPFR_EXP_LESS_P (b, c)) return mpfr_add1(a, c, b, rnd_mode); else return mpfr_add1(a, b, c, rnd_mode); diff --git a/src/fmma.c b/src/fmma.c index 49fd9742e..fbd94ff6d 100644 --- a/src/fmma.c +++ b/src/fmma.c @@ -23,108 +23,43 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" static int -mpfr_fmma_slow (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, - mpfr_srcptr d, mpfr_rnd_t rnd) +mpfr_fmma_aux (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, + mpfr_srcptr d, mpfr_rnd_t rnd, int neg) { - mpfr_t ab, cd; + mpfr_ubf_t u, v; + mp_size_t un, vn; + mpfr_limb_ptr up, vp; int inex; + MPFR_TMP_DECL(marker); - mpfr_init2 (ab, MPFR_PREC(a) + MPFR_PREC(b)); - mpfr_init2 (cd, MPFR_PREC(c) + MPFR_PREC(d)); - /* FIXME: The following multiplications may overflow or underflow - (even more often with the fact that the exponent range is not - extended), in which case the result is not exact. This should - be solved with future unbounded floats. */ - mpfr_mul (ab, a, b, MPFR_RNDZ); /* exact */ - mpfr_mul (cd, c, d, MPFR_RNDZ); /* exact */ - inex = mpfr_add (z, ab, cd, rnd); - mpfr_clear (ab); - mpfr_clear (cd); - return inex; -} + MPFR_LOG_FUNC + (("a[%Pu]=%.*Rg b[%Pu]=%.*Rg c[%Pu]=%.*Rg d[%Pu]=%.*Rg rnd=%d neg=%d", + mpfr_get_prec (a), mpfr_log_prec, a, + mpfr_get_prec (b), mpfr_log_prec, b, + mpfr_get_prec (c), mpfr_log_prec, c, + mpfr_get_prec (d), mpfr_log_prec, d, rnd, neg), + ("z[%Pu]=%.*Rg inex=%d", + mpfr_get_prec (z), mpfr_log_prec, z, inex)); -/* z <- a*b + c*d */ -static int -mpfr_fmma_fast (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, - mpfr_srcptr d, mpfr_rnd_t rnd) -{ - /* Assumes that a, b, c, d are finite and non-zero; so any multiplication - of two of them yielding an infinity is an overflow, and a - multiplication yielding 0 is an underflow. - Assumes further that z is distinct from a, b, c, d. */ - - int inex; - mpfr_t u, v; - mp_size_t an, bn, cn, dn; - mpfr_limb_ptr up, vp; - unsigned int saved_flags = __gmpfr_flags; - MPFR_TMP_DECL(marker); - MPFR_SAVE_EXPO_DECL (expo); - - MPFR_TMP_MARK(marker); - MPFR_SAVE_EXPO_MARK (expo); - - /* u=a*b, v=c*d exactly */ - an = MPFR_LIMB_SIZE(a); - bn = MPFR_LIMB_SIZE(b); - cn = MPFR_LIMB_SIZE(c); - dn = MPFR_LIMB_SIZE(d); - MPFR_TMP_INIT (up, u, (an + bn) * GMP_NUMB_BITS, an + bn); - MPFR_TMP_INIT (vp, v, (cn + dn) * GMP_NUMB_BITS, cn + dn); - - /* u <- a*b */ - if (an >= bn) - mpn_mul (up, MPFR_MANT(a), an, MPFR_MANT(b), bn); - else - mpn_mul (up, MPFR_MANT(b), bn, MPFR_MANT(a), an); - if ((up[an + bn - 1] & MPFR_LIMB_HIGHBIT) == 0) - { - mpn_lshift (up, up, an + bn, 1); - MPFR_SET_EXP (u, MPFR_EXP(a) + MPFR_EXP(b) - 1); - if (MPFR_UNLIKELY(MPFR_EXP(u) == __MPFR_EXP_INF)) - goto failure; - } - else - MPFR_SET_EXP (u, MPFR_EXP(a) + MPFR_EXP(b)); - - /* v <- c*d */ - if (cn >= dn) - mpn_mul (vp, MPFR_MANT(c), cn, MPFR_MANT(d), dn); - else - mpn_mul (vp, MPFR_MANT(d), dn, MPFR_MANT(c), cn); - if ((vp[cn + dn - 1] & MPFR_LIMB_HIGHBIT) == 0) - { - mpn_lshift (vp, vp, cn + dn, 1); - MPFR_SET_EXP (v, MPFR_EXP(c) + MPFR_EXP(d) - 1); - if (MPFR_UNLIKELY(MPFR_EXP(v) == __MPFR_EXP_INF)) - goto failure; - } - else - MPFR_SET_EXP (v, MPFR_EXP(c) + MPFR_EXP(d)); - - MPFR_PREC(u) = (an + bn) * GMP_NUMB_BITS; - MPFR_PREC(v) = (cn + dn) * GMP_NUMB_BITS; - MPFR_SIGN(u) = MPFR_MULT_SIGN(MPFR_SIGN(a), MPFR_SIGN(b)); - MPFR_SIGN(v) = MPFR_MULT_SIGN(MPFR_SIGN(c), MPFR_SIGN(d)); - - /* tentatively compute z as u+v; here we need z to be distinct - from a, b, c, d to avoid losing the input values in case we - need to call mpfr_fmma_slow */ - /* FIXME: The above comment is no longer valid. Anyway, with - unbounded floats (based on an exact multiplication like above), - it will no longer be necessary to distinguish fast and slow. */ - inex = mpfr_add (z, u, v, rnd); - - MPFR_TMP_FREE(marker); - MPFR_SAVE_EXPO_FREE (expo); - - return mpfr_check_range (z, inex, rnd); - - failure: - __gmpfr_flags = saved_flags; - MPFR_TMP_FREE(marker); - MPFR_SAVE_EXPO_FREE (expo); - return mpfr_fmma_slow (z, a, b, c, d, rnd); + MPFR_TMP_MARK (marker); + + un = MPFR_LIMB_SIZE (a) + MPFR_LIMB_SIZE (b); + vn = MPFR_LIMB_SIZE (c) + MPFR_LIMB_SIZE (d); + MPFR_TMP_INIT (up, u, (mpfr_prec_t) un * GMP_NUMB_BITS, un); + MPFR_TMP_INIT (vp, v, (mpfr_prec_t) vn * GMP_NUMB_BITS, vn); + + mpfr_ubf_mul_exact (u, a, b); + mpfr_ubf_mul_exact (v, c, d); + if (neg) + MPFR_CHANGE_SIGN (v); + inex = mpfr_add (z, (mpfr_srcptr) u, (mpfr_srcptr) v, rnd); + + MPFR_UBF_CLEAR_EXP (u); + MPFR_UBF_CLEAR_EXP (v); + + MPFR_TMP_FREE (marker); + + return inex; } /* z <- a*b + c*d */ @@ -132,13 +67,7 @@ int mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_srcptr d, mpfr_rnd_t rnd) { - mpfr_limb_ptr zp = MPFR_MANT(z); - - return (mpfr_regular_p (a) && mpfr_regular_p (b) && mpfr_regular_p (c) && - mpfr_regular_p (d) && zp != MPFR_MANT(a) && zp != MPFR_MANT(b) && - zp != MPFR_MANT(c) && zp != MPFR_MANT(d)) - ? mpfr_fmma_fast (z, a, b, c, d, rnd) - : mpfr_fmma_slow (z, a, b, c, d, rnd); + return mpfr_fmma_aux (z, a, b, c, d, rnd, 0); } /* z <- a*b - c*d */ @@ -146,8 +75,5 @@ int mpfr_fmms (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_srcptr d, mpfr_rnd_t rnd) { - mpfr_t minus_c; - - MPFR_ALIAS (minus_c, c, -MPFR_SIGN(c), MPFR_EXP(c)); - return mpfr_fmma (z, a, b, minus_c, d, rnd); + return mpfr_fmma_aux (z, a, b, c, d, rnd, 1); } diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index 5d989f489..67997c010 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -869,6 +869,10 @@ typedef intmax_t mpfr_eexp_t; # define MPFR_SET_INVALID_EXP(x) ((void) 0) #endif +#define MPFR_EXP_LESS_P(x,y) \ + (MPFR_UNLIKELY (MPFR_IS_UBF (x) || MPFR_IS_UBF (y)) ? \ + mpfr_ubf_exp_less_p (x, y) : MPFR_GET_EXP (x) < MPFR_GET_EXP (y)) + /****************************************************** ********* Singular values (NAN, INF, ZERO) ********* @@ -878,6 +882,7 @@ typedef intmax_t mpfr_eexp_t; # define MPFR_EXP_ZERO (MPFR_EXP_MIN+1) # define MPFR_EXP_NAN (MPFR_EXP_MIN+2) # define MPFR_EXP_INF (MPFR_EXP_MIN+3) +# define MPFR_EXP_UBF (MPFR_EXP_MIN+4) #define MPFR_IS_NAN(x) (MPFR_EXP(x) == MPFR_EXP_NAN) #define MPFR_SET_NAN(x) (MPFR_EXP(x) = MPFR_EXP_NAN) @@ -886,6 +891,8 @@ typedef intmax_t mpfr_eexp_t; #define MPFR_IS_ZERO(x) (MPFR_EXP(x) == MPFR_EXP_ZERO) #define MPFR_SET_ZERO(x) (MPFR_EXP(x) = MPFR_EXP_ZERO) #define MPFR_NOTZERO(x) (MPFR_EXP(x) != MPFR_EXP_ZERO) +#define MPFR_IS_UBF(x) (MPFR_EXP(x) == MPFR_EXP_UBF) +#define MPFR_SET_UBF(x) (MPFR_EXP(x) = MPFR_EXP_UBF) #define MPFR_IS_NORMALIZED(x) \ (MPFR_LIMB_MSB (MPFR_MANT(x)[MPFR_LAST_LIMB(x)]) != 0) @@ -897,6 +904,12 @@ typedef intmax_t mpfr_eexp_t; (MPFR_ASSERTD (MPFR_EXP (x) >= MPFR_EMIN_MIN && \ MPFR_EXP (x) <= MPFR_EMAX_MAX && \ MPFR_IS_NORMALIZED (x)), 1)) +#define MPFR_IS_PURE_UBF(x) \ + (!MPFR_IS_SINGULAR(x) && \ + (MPFR_ASSERTD ((MPFR_IS_UBF (x) || \ + (MPFR_EXP (x) >= MPFR_EMIN_MIN && \ + MPFR_EXP (x) <= MPFR_EMAX_MAX)) && \ + MPFR_IS_NORMALIZED (x)), 1)) #define MPFR_ARE_SINGULAR(x,y) \ (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)) || MPFR_UNLIKELY(MPFR_IS_SINGULAR(y))) @@ -2201,4 +2214,50 @@ __MPFR_DECLSPEC extern int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2]; #endif /* MPFR_COV_CHECK */ +/****************************************************** + ***************** Unbounded Floats ***************** + ******************************************************/ + +#if defined (__cplusplus) +extern "C" { +#endif + +/* An UBF is like a MPFR number, but with an additional mpz_t member, + which is assumed to be present (with a value in it) when the usual + exponent field has the value MPFR_EXP_UBF. The goal of this compatible + representation is to easily be able to support UBF in "normal" code + and hopefully avoid aliasing issues at the same time. And code that + accepts UBF in input should also accept mpfr_t as a consequence; this + makes mpfr_t to UBF conversion unnecessary. + When an input of a public function is a UBF, the semantic remains + internal to MPFR and can change in the future. + Note that functions used for logging need to support UBF. *FIXME* */ + +typedef struct { + mpfr_prec_t _mpfr_prec; + mpfr_sign_t _mpfr_sign; + mpfr_exp_t _mpfr_exp; + mp_limb_t *_mpfr_d; + mpz_t _mpfr_zexp; +} __mpfr_ubf_struct; + +typedef __mpfr_ubf_struct mpfr_ubf_t[1]; +typedef __mpfr_ubf_struct *mpfr_ubf_ptr; + +__MPFR_DECLSPEC void + mpfr_ubf_mul_exact _MPFR_PROTO ((mpfr_ubf_ptr, mpfr_srcptr, mpfr_srcptr)); +__MPFR_DECLSPEC int + mpfr_ubf_exp_less_p _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr)); + +#if defined (__cplusplus) +} +#endif + +#define MPFR_ZEXP(x) \ + ((void) (x)->_mpfr_exp /* to check that x has a correct type */, \ + ((mpfr_ubf_ptr) (x))->_mpfr_zexp) + +#define MPFR_UBF_CLEAR_EXP(x) \ + ((void) (MPFR_IS_UBF (u) && (mpz_clear (MPFR_ZEXP (x)), 0))) + #endif /* __MPFR_IMPL_H__ */ diff --git a/src/print_raw.c b/src/print_raw.c index af9c6c1b9..cd19ddee7 100644 --- a/src/print_raw.c +++ b/src/print_raw.c @@ -49,25 +49,22 @@ mpfr_fprint_binary (FILE *stream, mpfr_srcptr x) px = MPFR_PREC (x); fprintf (stream, "0."); - for (n = (px - 1) / GMP_NUMB_BITS; ; n--) + for (n = (px - 1) / GMP_NUMB_BITS; n >= 0; n--) { mp_limb_t wd, t; - MPFR_ASSERTN (n >= 0); wd = mx[n]; for (t = MPFR_LIMB_HIGHBIT; t != 0; t >>= 1) { putc ((wd & t) == 0 ? '0' : '1', stream); if (--px == 0) - { - mpfr_exp_t ex; - - ex = MPFR_EXP (x); - fprintf (stream, "E%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) ex); - return; - } + break; } } + if (MPFR_IS_UBF (x)) + gmp_fprintf (stream, "E%Zd", MPFR_ZEXP (x)); + else + fprintf (stream, "E%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) MPFR_EXP (x)); } } diff --git a/src/sub.c b/src/sub.c index 457f67d9d..f6a4e16d5 100644 --- a/src/sub.c +++ b/src/sub.c @@ -79,8 +79,8 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } } - MPFR_ASSERTD (MPFR_IS_PURE_FP (b)); - MPFR_ASSERTD (MPFR_IS_PURE_FP (c)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); if (MPFR_LIKELY (MPFR_SIGN (b) == MPFR_SIGN (c))) { /* signs are equal, it's a real subtraction */ @@ -92,7 +92,7 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } else { /* signs differ, it's an addition */ - if (MPFR_GET_EXP (b) < MPFR_GET_EXP (c)) + if (MPFR_EXP_LESS_P (b, c)) { /* exchange rounding modes toward +/- infinity */ int inexact; rnd_mode = MPFR_INVERT_RND (rnd_mode); diff --git a/src/ubf.c b/src/ubf.c new file mode 100644 index 000000000..3f1362f19 --- /dev/null +++ b/src/ubf.c @@ -0,0 +1,177 @@ +/* Functions to work with unbounded floats (limited low-level interface). + +Copyright 2016 Free Software Foundation, Inc. +Contributed by the AriC and Caramba projects, INRIA. + +This file is part of the GNU MPFR Library. + +The GNU MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU 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 Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see +http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., +51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ + +#include "mpfr-impl.h" + +/* Note: In MPFR math functions, even if UBF code is not called first, + we may still need to handle special values NaN and infinities here. + Indeed, for MAXR * MAXR + (-inf), even though this is a special case, + the computation will also generate an overflow due to MAXR * MAXR, + so that UBF code will be called anyway (except if special cases are + detected and handled separately, but for polynomial, this should not + be needed). */ + +/* This function does not change the flags. */ +static void +mpfr_get_zexp (mpz_ptr ez, mpfr_srcptr x) +{ + mpz_init (ez); + + if (MPFR_IS_UBF (x)) + mpz_set (ez, MPFR_ZEXP (x)); + else + { + mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE]; + mpfr_t e; + int inex; + + /* TODO: Once this has been tested, optimize based on whether + _MPFR_EXP_FORMAT <= 3. */ + MPFR_TMP_INIT1 (e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT); + MPFR_DBGRES (inex = mpfr_set_exp_t (e, MPFR_GET_EXP (x), MPFR_RNDN)); + MPFR_ASSERTD (inex == 0); + MPFR_DBGRES (inex = mpfr_get_z (ez, e, MPFR_RNDN)); + MPFR_ASSERTD (inex == 0); + } +} + +/* Exact product. The number a is assumed to have enough allocated memory, + where the trailing bits are regarded as being part of the input numbers + (no reallocation is attempted and no check is performed as MPFR_TMP_INIT + could have been used). The arguments b and c may actually be UBF numbers + (mpfr_srcptr can be seen a bit like void *, but is stronger). + This function does not change the flags, except in case of NaN. */ +void +mpfr_ubf_mul_exact (mpfr_ubf_ptr a, mpfr_srcptr b, mpfr_srcptr c) +{ + MPFR_LOG_FUNC + (("b[%Pu]=%.*Rg c[%Pu]=%.*Rg", + mpfr_get_prec (b), mpfr_log_prec, b, + mpfr_get_prec (c), mpfr_log_prec, c), + ("a[%Pu]=%.*Rg", + mpfr_get_prec (a), mpfr_log_prec, a)); + + MPFR_ASSERTD ((mpfr_ptr) a != b); + MPFR_ASSERTD ((mpfr_ptr) a != c); + MPFR_SIGN (a) = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); + + if (MPFR_ARE_SINGULAR (b, c)) + { + if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c)) + MPFR_SET_NAN (a); + else if (MPFR_IS_INF (b)) + { + if (MPFR_NOTZERO (c)) + MPFR_SET_INF (a); + else + MPFR_SET_NAN (a); + } + else if (MPFR_IS_INF (c)) + { + if (!MPFR_IS_ZERO (b)) + MPFR_SET_INF (a); + else + MPFR_SET_NAN (a); + } + else + { + MPFR_ASSERTD (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c)); + MPFR_SET_ZERO (a); + } + } + else + { + mpfr_exp_t e; + mp_size_t bn, cn; + mpfr_limb_ptr ap; + mp_limb_t u, v; + int m; + + /* Note about the code below: For the choice of the precision of + * the result a, one could choose PREC(b) + PREC(c), instead of + * taking whole limbs into account, but in most cases where one + * would gain one limb, one would need to copy the significand + * instead of a no-op (see the mul.c code). + * But in the case MPFR_LIMB_MSB (u) == 0, if the result fits in + * an-1 limbs, one could actually do + * mpn_rshift (ap, ap, k, GMP_NUMB_BITS - 1) + * instead of + * mpn_lshift (ap, ap, k, 1) + * to gain one limb (and reduce the precision), replacing a shift + * by another one. Would this be interesting? + */ + + bn = MPFR_LIMB_SIZE (b); + cn = MPFR_LIMB_SIZE (c); + + ap = MPFR_MANT (a); + + u = (bn >= cn) ? + mpn_mul (ap, MPFR_MANT (b), bn, MPFR_MANT (c), cn) : + mpn_mul (ap, MPFR_MANT (c), cn, MPFR_MANT (b), bn); + if (MPFR_UNLIKELY (MPFR_LIMB_MSB (u) == 0)) + { + m = 1; + MPFR_DBGRES (v = mpn_lshift (ap, ap, bn + cn, 1)); + MPFR_ASSERTD (v == 0); + } + else + m = 0; + + if (! MPFR_IS_UBF (b) && ! MPFR_IS_UBF (c) && + (e = MPFR_GET_EXP (b) + MPFR_GET_EXP (c) - m, + MPFR_EXP_IN_RANGE (e))) + { + MPFR_SET_EXP (a, e); + } + else + { + mpz_t be, ce; + + mpz_init (MPFR_ZEXP (a)); + + /* This may involve copies of mpz_t, but exponents should not be + very large integers anyway. */ + mpfr_get_zexp (be, b); + mpfr_get_zexp (ce, c); + mpz_add (MPFR_ZEXP (a), be, ce); + mpz_clear (be); + mpz_clear (ce); + mpz_sub_ui (MPFR_ZEXP (a), MPFR_ZEXP (a), m); + MPFR_SET_UBF (a); + } + } +} + +int +mpfr_ubf_exp_less_p (mpfr_srcptr x, mpfr_srcptr y) +{ + mpz_t xe, ye; + int c; + + mpfr_get_zexp (xe, x); + mpfr_get_zexp (ye, y); + c = mpz_cmp (xe, ye) < 0; + mpz_clear (xe); + mpz_clear (ye); + return c; +} diff --git a/src/vasprintf.c b/src/vasprintf.c index 634b5b189..3238ba530 100644 --- a/src/vasprintf.c +++ b/src/vasprintf.c @@ -1552,6 +1552,21 @@ partition_number (struct number_parts *np, mpfr_srcptr p, } } } + else if (MPFR_UNLIKELY (MPFR_IS_UBF (p))) + { + /* mpfr_get_str does not support UBF, so that UBF numbers are regarded + as special cases here. This is not much a problem since UBF numbers + are internal to MPFR and here, they only for logging. */ + if (np->pad_type == LEADING_ZEROS) + /* change to right justification padding with left spaces */ + np->pad_type = LEFT; + + np->ip_size = 3; + str = (char *) (*__gmp_allocate_func) (1 + np->ip_size); + strcpy (str, uppercase ? "UBF" : "ubf"); + np->ip_ptr = register_string (np->sl, str); + /* TODO: output more information (e.g. the exponent) if need be. */ + } else { MPFR_ASSERTD (MPFR_IS_PURE_FP (p)); -- cgit v1.2.1 From d7efcf0fed1ff92af216e8b3f7b00707bd2375dd Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 13:15:49 +0000 Subject: [src/ubf.c] Support reduced exponent range in mpfr_get_zexp. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10320 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/ubf.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ubf.c b/src/ubf.c index 3f1362f19..1aa67d759 100644 --- a/src/ubf.c +++ b/src/ubf.c @@ -43,14 +43,17 @@ mpfr_get_zexp (mpz_ptr ez, mpfr_srcptr x) mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE]; mpfr_t e; int inex; + MPFR_SAVE_EXPO_DECL (expo); /* TODO: Once this has been tested, optimize based on whether _MPFR_EXP_FORMAT <= 3. */ MPFR_TMP_INIT1 (e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT); + MPFR_SAVE_EXPO_MARK (expo); MPFR_DBGRES (inex = mpfr_set_exp_t (e, MPFR_GET_EXP (x), MPFR_RNDN)); MPFR_ASSERTD (inex == 0); MPFR_DBGRES (inex = mpfr_get_z (ez, e, MPFR_RNDN)); MPFR_ASSERTD (inex == 0); + MPFR_SAVE_EXPO_FREE (expo); } } -- cgit v1.2.1 From 912aa998b518869863f8b93e717dbb4e6edd5f1d Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 14:23:11 +0000 Subject: Added UBF support for mpfr_cmp2. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10322 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/cmp2.c | 31 +++++++++++++++++++++---------- src/mpfr-impl.h | 2 ++ src/ubf.c | 34 ++++++++++++++++++++++++++++++++++ 3 files changed, 57 insertions(+), 10 deletions(-) diff --git a/src/cmp2.c b/src/cmp2.c index ee2a8d138..ccf8175ab 100644 --- a/src/cmp2.c +++ b/src/cmp2.c @@ -38,6 +38,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) { mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0; mp_size_t bn, cn; + mpfr_exp_t sdiff_exp; mpfr_uexp_t diff_exp; mpfr_prec_t res = 0; int sign; @@ -48,13 +49,21 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) /* the cases b=0 or c=0 are also treated apart in agm and sub (which calls sub1) */ - MPFR_ASSERTD (MPFR_IS_PURE_FP(b)); - MPFR_ASSERTD (MPFR_IS_PURE_FP(c)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); - if (MPFR_GET_EXP (b) >= MPFR_GET_EXP (c)) + sdiff_exp = MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)) ? + mpfr_ubf_diff_exp (b, c) : MPFR_GET_EXP (b) - MPFR_GET_EXP (c); + + /* The returned result is restricted to [MPFR_EXP_MIN,MPFR_EXP_MAX], + which is the range of the mpfr_exp_t type. But under the condition + below, the value of cancel will not be affected. */ + MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX); + + if (sdiff_exp >= 0) { sign = 1; - diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c); + diff_exp = sdiff_exp; bp = MPFR_MANT(b); cp = MPFR_MANT(c); @@ -85,7 +94,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) if (MPFR_UNLIKELY (cn < 0)) /* c discards exactly the upper part of b */ { - unsigned int z; + int z; MPFR_ASSERTD (bn >= 0); @@ -96,7 +105,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) res += GMP_NUMB_BITS; } - count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */ + count_leading_zeros (z, bp[bn]); /* bp[bn] <> 0 */ *cancel = res + z; return sign; } @@ -118,7 +127,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) else /* MPFR_EXP(b) < MPFR_EXP(c) */ { sign = -1; - diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (c) - MPFR_GET_EXP (b); + diff_exp = - (mpfr_uexp_t) sdiff_exp; bp = MPFR_MANT(c); cp = MPFR_MANT(b); @@ -152,6 +161,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0) && (high_dif == 0) && (dif == 1))) { /* dif=1 implies diff_exp = 0 or 1 */ + MPFR_ASSERTD (diff_exp <= 1); bb = (bn >= 0) ? bp[bn] : 0; cc = lastc; if (cn >= 0) @@ -160,8 +170,9 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) { cc += cp[cn]; } - else /* diff_exp = 1 */ + else { + MPFR_ASSERTD (diff_exp == 1); cc += cp[cn] >> 1; lastc = cp[cn] << (GMP_NUMB_BITS - 1); } @@ -188,9 +199,9 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) } else /* high_dif == 0 */ { - unsigned int z; + int z; - count_leading_zeros(z, dif); /* dif > 1 here */ + count_leading_zeros (z, dif); /* dif > 1 here */ res += z; if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (GMP_NUMB_BITS - z - 1)))) { /* dif is not a power of two */ diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index 67997c010..cf4905295 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -2248,6 +2248,8 @@ __MPFR_DECLSPEC void mpfr_ubf_mul_exact _MPFR_PROTO ((mpfr_ubf_ptr, mpfr_srcptr, mpfr_srcptr)); __MPFR_DECLSPEC int mpfr_ubf_exp_less_p _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr)); +__MPFR_DECLSPEC mpfr_exp_t + mpfr_ubf_diff_exp _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr)); #if defined (__cplusplus) } diff --git a/src/ubf.c b/src/ubf.c index 1aa67d759..fe596c32e 100644 --- a/src/ubf.c +++ b/src/ubf.c @@ -178,3 +178,37 @@ mpfr_ubf_exp_less_p (mpfr_srcptr x, mpfr_srcptr y) mpz_clear (ye); return c; } + +/* Return the difference of the exponents of x and y, restricted to + the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */ +mpfr_exp_t +mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y) +{ + mpz_t xe, ye; + mp_size_t n; + mpfr_eexp_t e; + mpfr_t d; + int inex; + MPFR_SAVE_EXPO_DECL (expo); + + mpfr_get_zexp (xe, x); + mpfr_get_zexp (ye, y); + mpz_sub (xe, xe, ye); + mpz_clear (ye); + n = ABSIZ(xe); /* limb size of xe */ + if (n == 0) + return 0; + MPFR_SAVE_EXPO_MARK (expo); + mpfr_init2 (d, n * GMP_NUMB_BITS); + MPFR_DBGRES (inex = mpfr_set_z (d, xe, MPFR_RNDN)); + MPFR_ASSERTD (inex == 0); + mpz_clear (xe); + e = mpfr_get_exp_t (d, MPFR_RNDZ); + mpfr_clear (d); + MPFR_SAVE_EXPO_FREE (expo); + if (MPFR_UNLIKELY (e < MPFR_EXP_MIN)) + return MPFR_EXP_MIN; + if (MPFR_UNLIKELY (e > MPFR_EXP_MAX)) + return MPFR_EXP_MAX; + return e; +} -- cgit v1.2.1 From 6550eae6e9f3125fd2aa3deb1f9d7301ae8f5e35 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 14:26:17 +0000 Subject: [src/ubf.c] Memory leak in some case. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10323 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/ubf.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/ubf.c b/src/ubf.c index fe596c32e..8d897c61d 100644 --- a/src/ubf.c +++ b/src/ubf.c @@ -197,7 +197,10 @@ mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y) mpz_clear (ye); n = ABSIZ(xe); /* limb size of xe */ if (n == 0) - return 0; + { + mpz_clear (xe); + return 0; + } MPFR_SAVE_EXPO_MARK (expo); mpfr_init2 (d, n * GMP_NUMB_BITS); MPFR_DBGRES (inex = mpfr_set_z (d, xe, MPFR_RNDN)); -- cgit v1.2.1 From 1bd2091c3dd82bfdac476a6243edd06f3616300b Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 14:46:57 +0000 Subject: [tests/tfmma.c] Added max_tests. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10324 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/tfmma.c b/tests/tfmma.c index b825de646..4140ae05a 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -177,6 +177,51 @@ zero_tests (void) set_emax (emax); } +static void +max_tests (void) +{ + mpfr_exp_t emin, emax; + mpfr_t x, y1, y2; + int r; + int inex1, inex2; + mpfr_flags_t flags1, flags2; + + emin = mpfr_get_emin (); + emax = mpfr_get_emax (); + set_emin (MPFR_EMIN_MIN); + set_emax (MPFR_EMAX_MAX); + + mpfr_inits2 (64, x, y1, y2, (mpfr_ptr) 0); + mpfr_setmax (x, MPFR_EMAX_MAX); + flags1 = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; + RND_LOOP (r) + { + inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r); + mpfr_clear_flags (); + inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r); + flags2 = __gmpfr_flags; + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (y1, y2))) + { + printf ("Error in max_tests for %s", + mpfr_print_rnd_mode ((mpfr_rnd_t) r)); + printf ("Expected "); + mpfr_dump (y1); + printf (" with inex = %d, flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (y2); + printf (" with inex = %d, flags =", inex2); + flags_out (flags2); + exit (1); + } + } + mpfr_clears (x, y1, y2, (mpfr_ptr) 0); + + set_emin (emin); + set_emax (emax); +} + int main (int argc, char *argv[]) { @@ -184,6 +229,7 @@ main (int argc, char *argv[]) random_tests (); zero_tests (); + max_tests (); tests_end_mpfr (); return 0; -- cgit v1.2.1 From c3d150bd5da99585884e1733909d32a61ed87078 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 23 May 2016 15:08:40 +0000 Subject: More UBF support: * mpfr-impl.h, ubf.c: added mpfr_ubf_zexp2exp function. * add1.c: support the case where b (the first input) is an UBF. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10325 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/add1.c | 14 +++++++++++--- src/mpfr-impl.h | 2 ++ src/ubf.c | 37 +++++++++++++++++++++++-------------- 3 files changed, 36 insertions(+), 17 deletions(-) diff --git a/src/add1.c b/src/add1.c index 2c46eaaf3..497b3e989 100644 --- a/src/add1.c +++ b/src/add1.c @@ -35,8 +35,17 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) mpfr_uexp_t diff_exp; MPFR_TMP_DECL(marker); - MPFR_ASSERTD(MPFR_IS_PURE_FP(b)); - MPFR_ASSERTD(MPFR_IS_PURE_FP(c)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); + + if (MPFR_UNLIKELY (MPFR_IS_UBF (b))) + { + exp = mpfr_ubf_zexp2exp (MPFR_ZEXP (b)); + if (exp > __gmpfr_emax) + return mpfr_overflow (a, rnd_mode, MPFR_SIGN (b));; + } + else + exp = MPFR_GET_EXP (b); MPFR_TMP_MARK(marker); @@ -68,7 +77,6 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) MPN_COPY(cp, ap, cn); } - exp = MPFR_GET_EXP (b); MPFR_SET_SAME_SIGN(a, b); MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN(b)); /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */ diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index cf4905295..b9ed7d159 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -2248,6 +2248,8 @@ __MPFR_DECLSPEC void mpfr_ubf_mul_exact _MPFR_PROTO ((mpfr_ubf_ptr, mpfr_srcptr, mpfr_srcptr)); __MPFR_DECLSPEC int mpfr_ubf_exp_less_p _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr)); +__MPFR_DECLSPEC mpfr_exp_t + mpfr_ubf_zexp2exp _MPFR_PROTO ((mpz_ptr)); __MPFR_DECLSPEC mpfr_exp_t mpfr_ubf_diff_exp _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr)); diff --git a/src/ubf.c b/src/ubf.c index 8d897c61d..e43e0e41d 100644 --- a/src/ubf.c +++ b/src/ubf.c @@ -179,33 +179,25 @@ mpfr_ubf_exp_less_p (mpfr_srcptr x, mpfr_srcptr y) return c; } -/* Return the difference of the exponents of x and y, restricted to +/* Convert an mpz_t to an mpfr_exp_t, restricted to the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */ mpfr_exp_t -mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y) +mpfr_ubf_zexp2exp (mpz_ptr ez) { - mpz_t xe, ye; mp_size_t n; mpfr_eexp_t e; mpfr_t d; int inex; MPFR_SAVE_EXPO_DECL (expo); - mpfr_get_zexp (xe, x); - mpfr_get_zexp (ye, y); - mpz_sub (xe, xe, ye); - mpz_clear (ye); - n = ABSIZ(xe); /* limb size of xe */ + n = ABSIZ (ez); /* limb size of ez */ if (n == 0) - { - mpz_clear (xe); - return 0; - } + return 0; + MPFR_SAVE_EXPO_MARK (expo); mpfr_init2 (d, n * GMP_NUMB_BITS); - MPFR_DBGRES (inex = mpfr_set_z (d, xe, MPFR_RNDN)); + MPFR_DBGRES (inex = mpfr_set_z (d, ez, MPFR_RNDN)); MPFR_ASSERTD (inex == 0); - mpz_clear (xe); e = mpfr_get_exp_t (d, MPFR_RNDZ); mpfr_clear (d); MPFR_SAVE_EXPO_FREE (expo); @@ -215,3 +207,20 @@ mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y) return MPFR_EXP_MAX; return e; } + +/* Return the difference of the exponents of x and y, restricted to + the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */ +mpfr_exp_t +mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y) +{ + mpz_t xe, ye; + mpfr_exp_t e; + + mpfr_get_zexp (xe, x); + mpfr_get_zexp (ye, y); + mpz_sub (xe, xe, ye); + mpz_clear (ye); + e = mpfr_ubf_zexp2exp (xe); + mpz_clear (xe); + return e; +} -- cgit v1.2.1 From 491a95a1f4ea6b1c6304092caec65735bee8dc61 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 24 May 2016 09:41:49 +0000 Subject: [tests/tfmma.c] Added near_overflow_tests, which crashes. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10340 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 81 insertions(+), 2 deletions(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index 4140ae05a..d6fcafbaa 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -22,8 +22,6 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-test.h" -/* TODO: add more tests, with special values and exception checking. */ - /* check both mpfr_fmma and mpfr_fmms */ static void random_test (mpfr_t a, mpfr_t b, mpfr_t c, mpfr_t d, mpfr_rnd_t rnd) @@ -222,6 +220,85 @@ max_tests (void) set_emax (emax); } +/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't. + * With MPFR r10337 and gcc -fsanitize=undefined -fno-sanitize-recover, + * this triggers the following error on a 64-bit machine: + * /usr/local/gmp-6.0.0-debug/include/gmp.h:2142:3: + * runtime error: store to null pointer of type 'mp_limb_t' + * (mpn_add_1 is called on a null pointer). + */ +static void +near_overflow_tests (void) +{ + mpfr_exp_t old_emax; + int i; + + old_emax = mpfr_get_emax (); + + for (i = 0; i < 40; i++) + { + mpfr_exp_t emax; + mpfr_t a, k, c, d, z1, z2; + mpfr_rnd_t rnd; + mpfr_prec_t prec_a; + int add = i & 1, inex, inex1, inex2; + mpfr_flags_t flags1, flags2; + + /* In most cases, we do the test with the maximum exponent. */ + emax = i % 8 != 0 ? MPFR_EMAX_MAX : 64 + (randlimb () % 1); + set_emax (emax); + + rnd = RND_RAND (); + prec_a = 64 + randlimb () % 100; + + mpfr_init2 (a, prec_a); + mpfr_urandom (a, RANDS, MPFR_RNDN); + mpfr_set_exp (a, emax/2 + 32); + + mpfr_init2 (k, prec_a - 32); + mpfr_urandom (k, RANDS, MPFR_RNDN); + mpfr_set_exp (k, emax/2); + + mpfr_init2 (c, prec_a + 1); + inex = mpfr_add (c, a, k, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + + mpfr_init2 (d, prec_a); + inex = mpfr_sub (d, a, k, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + if (add) + mpfr_neg (d, d, MPFR_RNDN); + + mpfr_clear_flags (); + inex1 = mpfr_sqr (z1, k, rnd); + flags1 = __gmpfr_flags; + + mpfr_clear_flags (); + inex2 = add ? + mpfr_fmma (z2, a, a, c, d, rnd) : + mpfr_fmms (z2, a, a, c, d, rnd); + flags2 = __gmpfr_flags; + + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (z1, z2))) + { + printf ("Error in near_overflow_tests for %s", + mpfr_print_rnd_mode (rnd)); + printf ("Expected "); + mpfr_dump (z1); + printf (" with inex = %d, flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (z2); + printf (" with inex = %d, flags =", inex2); + flags_out (flags2); + exit (1); + } + } + + set_emax (old_emax); +} + int main (int argc, char *argv[]) { @@ -230,6 +307,8 @@ main (int argc, char *argv[]) random_tests (); zero_tests (); max_tests (); + near_overflow_tests (); + /* TODO: near_underflow_tests (); */ tests_end_mpfr (); return 0; -- cgit v1.2.1 From 7421242f3a086fcb6cdd3c89794ab95b3c62b29c Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 24 May 2016 09:48:15 +0000 Subject: [tests/tfmma.c] Completed near_overflow_tests (there's still a crash, this time really in fmma.c). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10341 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index d6fcafbaa..e2c7e7530 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -220,13 +220,7 @@ max_tests (void) set_emax (emax); } -/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't. - * With MPFR r10337 and gcc -fsanitize=undefined -fno-sanitize-recover, - * this triggers the following error on a 64-bit machine: - * /usr/local/gmp-6.0.0-debug/include/gmp.h:2142:3: - * runtime error: store to null pointer of type 'mp_limb_t' - * (mpn_add_1 is called on a null pointer). - */ +/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't. */ static void near_overflow_tests (void) { @@ -240,7 +234,7 @@ near_overflow_tests (void) mpfr_exp_t emax; mpfr_t a, k, c, d, z1, z2; mpfr_rnd_t rnd; - mpfr_prec_t prec_a; + mpfr_prec_t prec_a, prec_z; int add = i & 1, inex, inex1, inex2; mpfr_flags_t flags1, flags2; @@ -250,6 +244,7 @@ near_overflow_tests (void) rnd = RND_RAND (); prec_a = 64 + randlimb () % 100; + prec_z = MPFR_PREC_MIN + randlimb () % 160; mpfr_init2 (a, prec_a); mpfr_urandom (a, RANDS, MPFR_RNDN); @@ -269,10 +264,12 @@ near_overflow_tests (void) if (add) mpfr_neg (d, d, MPFR_RNDN); + mpfr_init2 (z1, prec_z); mpfr_clear_flags (); inex1 = mpfr_sqr (z1, k, rnd); flags1 = __gmpfr_flags; + mpfr_init2 (z2, prec_z); mpfr_clear_flags (); inex2 = add ? mpfr_fmma (z2, a, a, c, d, rnd) : @@ -294,6 +291,8 @@ near_overflow_tests (void) flags_out (flags2); exit (1); } + + mpfr_clears (a, k, c, d, z1, z2, (mpfr_ptr) 0); } set_emax (old_emax); -- cgit v1.2.1 From dd5c6275d50efc8c669443f1614ac88f63f3ae66 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 24 May 2016 10:09:58 +0000 Subject: [src/vasprintf.c] Output the sign of UBF like with infinities. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10343 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/vasprintf.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/vasprintf.c b/src/vasprintf.c index 3238ba530..3e5fb0672 100644 --- a/src/vasprintf.c +++ b/src/vasprintf.c @@ -1561,6 +1561,9 @@ partition_number (struct number_parts *np, mpfr_srcptr p, /* change to right justification padding with left spaces */ np->pad_type = LEFT; + if (MPFR_IS_NEG (p)) + np->sign = '-'; + np->ip_size = 3; str = (char *) (*__gmp_allocate_func) (1 + np->ip_size); strcpy (str, uppercase ? "UBF" : "ubf"); -- cgit v1.2.1 From 1c5bb2a14facbcfbe736753fdc572507e3cbc867 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 24 May 2016 16:43:03 +0000 Subject: [src/sub1.c] Started to add UBF support (still incomplete). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10358 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/sub1.c | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/sub1.c b/src/sub1.c index a0bee3d6a..d9ced1646 100644 --- a/src/sub1.c +++ b/src/sub1.c @@ -28,11 +28,13 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., a positive value otherwise. */ +/* TODO: check the code in case exp_b == MPFR_EXP_MAX. */ + int mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) { int sign; - mpfr_exp_t diff_exp; + mpfr_exp_t diff_exp, exp_b; mpfr_prec_t cancel, cancel1; mp_size_t cancel2, an, bn, cn, cn0; mp_limb_t *ap, *bp, *cp; @@ -87,7 +89,17 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) else MPFR_SET_SAME_SIGN (a,b); - diff_exp = MPFR_GET_EXP (b) - MPFR_GET_EXP (c); + if (MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c))) + { + exp_b = MPFR_IS_UBF (b) ? + mpfr_ubf_zexp2exp (MPFR_ZEXP (b)) : MPFR_GET_EXP (b); + diff_exp = mpfr_ubf_diff_exp (b, c); + } + else + { + exp_b = MPFR_GET_EXP (b); + diff_exp = exp_b - MPFR_GET_EXP (c); + } MPFR_ASSERTD (diff_exp >= 0); aq = MPFR_GET_PREC (a); @@ -104,7 +116,9 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) = B.BBBBBBBBBBBBBBB - C.CCCCCCCCCCCCC */ /* A = S*ABS(B) +/- ulp(a) */ - MPFR_SET_EXP (a, MPFR_GET_EXP (b)); + /* FIXME: The following is incorrect when exp_b > __gmpfr_emax, + but we need to add a test first for complete coverage. */ + MPFR_SET_EXP (a, exp_b); MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq, rnd_mode, MPFR_SIGN (a), if (MPFR_UNLIKELY (++ MPFR_EXP (a) > __gmpfr_emax)) @@ -614,7 +628,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) mpfr_exp_t exp_a; cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */ - exp_a = MPFR_GET_EXP (b) - cancel; + exp_a = exp_b - cancel; if (MPFR_UNLIKELY(exp_a < __gmpfr_emin)) { MPFR_TMP_FREE(marker); @@ -631,10 +645,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) /* in case cancel = 0, add_exp can still be 1, in case b is just below a power of two, c is very small, prec(a) < prec(b), and rnd=away or nearest */ - mpfr_exp_t exp_b; - - exp_b = MPFR_GET_EXP (b); - if (MPFR_UNLIKELY(add_exp && exp_b == __gmpfr_emax)) + if (MPFR_UNLIKELY(add_exp && exp_b >= __gmpfr_emax)) { MPFR_TMP_FREE(marker); return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); -- cgit v1.2.1 From 53ea8ba36880bbd8c6eae7333f437ecf6bdbbded Mon Sep 17 00:00:00 2001 From: vlefevre Date: Sun, 29 May 2016 16:19:50 +0000 Subject: [tests/tfmma.c] Added overflow tests that trigger an assertion failure in sub1.c (because UBF support is not complete yet). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10381 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 59 ++++++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 50 insertions(+), 9 deletions(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index e2c7e7530..6b24a25fc 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -220,9 +220,12 @@ max_tests (void) set_emax (emax); } -/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't. */ +/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't. + * a^2 + cd where a^2 overflows and cd doesn't affect the overflow + * (and cd may even underflow). + */ static void -near_overflow_tests (void) +overflow_tests (void) { mpfr_exp_t old_emax; int i; @@ -231,7 +234,7 @@ near_overflow_tests (void) for (i = 0; i < 40; i++) { - mpfr_exp_t emax; + mpfr_exp_t emax, exp_a; mpfr_t a, k, c, d, z1, z2; mpfr_rnd_t rnd; mpfr_prec_t prec_a, prec_z; @@ -241,18 +244,19 @@ near_overflow_tests (void) /* In most cases, we do the test with the maximum exponent. */ emax = i % 8 != 0 ? MPFR_EMAX_MAX : 64 + (randlimb () % 1); set_emax (emax); + exp_a = emax/2 + 32; rnd = RND_RAND (); prec_a = 64 + randlimb () % 100; prec_z = MPFR_PREC_MIN + randlimb () % 160; mpfr_init2 (a, prec_a); - mpfr_urandom (a, RANDS, MPFR_RNDN); - mpfr_set_exp (a, emax/2 + 32); + mpfr_urandom (a, RANDS, MPFR_RNDU); + mpfr_set_exp (a, exp_a); mpfr_init2 (k, prec_a - 32); - mpfr_urandom (k, RANDS, MPFR_RNDN); - mpfr_set_exp (k, emax/2); + mpfr_urandom (k, RANDS, MPFR_RNDU); + mpfr_set_exp (k, exp_a - 32); mpfr_init2 (c, prec_a + 1); inex = mpfr_add (c, a, k, MPFR_RNDN); @@ -279,7 +283,44 @@ near_overflow_tests (void) if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && mpfr_equal_p (z1, z2))) { - printf ("Error in near_overflow_tests for %s", + printf ("Error 1 in overflow_tests for %s", + mpfr_print_rnd_mode (rnd)); + printf ("Expected "); + mpfr_dump (z1); + printf (" with inex = %d, flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (z2); + printf (" with inex = %d, flags =", inex2); + flags_out (flags2); + exit (1); + } + + /* c and d such that a^2 +/- cd ~= a^2 (overflow) */ + mpfr_urandom (c, RANDS, MPFR_RNDU); + mpfr_set_exp (c, randlimb () % 1 ? exp_a - 2 : __gmpfr_emin); + if (randlimb () % 1) + mpfr_neg (c, c, MPFR_RNDN); + mpfr_urandom (d, RANDS, MPFR_RNDU); + mpfr_set_exp (d, randlimb () % 1 ? exp_a - 2 : __gmpfr_emin); + if (randlimb () % 1) + mpfr_neg (d, d, MPFR_RNDN); + + mpfr_clear_flags (); + inex1 = mpfr_sqr (z1, a, rnd); + flags1 = __gmpfr_flags; + MPFR_ASSERTN (flags1 == (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)); + + mpfr_clear_flags (); + inex2 = add ? + mpfr_fmma (z2, a, a, c, d, rnd) : + mpfr_fmms (z2, a, a, c, d, rnd); + flags2 = __gmpfr_flags; + + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (z1, z2))) + { + printf ("Error 2 in overflow_tests for %s", mpfr_print_rnd_mode (rnd)); printf ("Expected "); mpfr_dump (z1); @@ -306,7 +347,7 @@ main (int argc, char *argv[]) random_tests (); zero_tests (); max_tests (); - near_overflow_tests (); + overflow_tests (); /* TODO: near_underflow_tests (); */ tests_end_mpfr (); -- cgit v1.2.1 From b635ceae1de3f27dcce172a216968b2cca1de4fa Mon Sep 17 00:00:00 2001 From: vlefevre Date: Sun, 29 May 2016 17:36:39 +0000 Subject: [tests/tfmma.c] Forgot a \n in printf. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10382 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index 6b24a25fc..fcd291557 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -201,7 +201,7 @@ max_tests (void) if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && mpfr_equal_p (y1, y2))) { - printf ("Error in max_tests for %s", + printf ("Error in max_tests for %s\n", mpfr_print_rnd_mode ((mpfr_rnd_t) r)); printf ("Expected "); mpfr_dump (y1); @@ -283,7 +283,7 @@ overflow_tests (void) if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && mpfr_equal_p (z1, z2))) { - printf ("Error 1 in overflow_tests for %s", + printf ("Error 1 in overflow_tests for %s\n", mpfr_print_rnd_mode (rnd)); printf ("Expected "); mpfr_dump (z1); @@ -320,7 +320,7 @@ overflow_tests (void) if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && mpfr_equal_p (z1, z2))) { - printf ("Error 2 in overflow_tests for %s", + printf ("Error 2 in overflow_tests for %s\n", mpfr_print_rnd_mode (rnd)); printf ("Expected "); mpfr_dump (z1); -- cgit v1.2.1 From 8a96712b7005b9a561bf6bf3144678217bb8c67b Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 30 May 2016 13:44:49 +0000 Subject: [src/sub1.c] Added UBF support. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10389 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/sub1.c | 100 +++++++++++++++++++++++++++++++++++-------------------------- 1 file changed, 58 insertions(+), 42 deletions(-) diff --git a/src/sub1.c b/src/sub1.c index 8733fe4dc..6a5a910a6 100644 --- a/src/sub1.c +++ b/src/sub1.c @@ -111,64 +111,74 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) but it is more expensive and not very useful */ if (MPFR_UNLIKELY (MAX (aq, bq) + 2 <= diff_exp)) { + MPFR_LOG_MSG (("case c small\n", 0)); + /* Remember, we can't have an exact result! */ /* A.AAAAAAAAAAAAAAAAA = B.BBBBBBBBBBBBBBB - C.CCCCCCCCCCCCC */ /* A = S*ABS(B) +/- ulp(a) */ - /* FIXME: The following is incorrect when exp_b > __gmpfr_emax, - but we need to add a test first for complete coverage. */ - MPFR_SET_EXP (a, exp_b); + MPFR_EXP (a) = exp_b; /* may be up to MPFR_EXP_MAX */ MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq, - rnd_mode, MPFR_SIGN (a), ++ MPFR_EXP (a)); + rnd_mode, MPFR_SIGN (a), + if (MPFR_EXP (a) != MPFR_EXP_MAX) + ++ MPFR_EXP (a)); + MPFR_LOG_MSG (("inexact=%d\n", inexact)); if (inexact == 0) { - /* a = b (Exact) - But we know it isn't (Since we have to remove `c') - So if we round to Zero, we have to remove one ulp. - Otherwise the result is correctly rounded. */ - /* An overflow is not possible. */ - MPFR_ASSERTD (MPFR_EXP (a) <= __gmpfr_emax); - if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a))) + /* a = b, but the exact value of b - c is a bit below. Then, + except for directed rounding similar to toward zero and + before overflow checking: a is the correctly rounded value + and since |b| - |c| < |a|, the ternary value value is given + by the sign of a. */ + if (! MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a))) { - mpfr_nexttozero (a); - MPFR_RET (- MPFR_INT_SIGN (a)); + inexact = MPFR_INT_SIGN (a); + goto check_overflow; } - MPFR_RET (MPFR_INT_SIGN (a)); } - else + else /* inexact != 0 */ { /* A.AAAAAAAAAAAAAA = B.BBBBBBBBBBBBBBB - C.CCCCCCCCCCCCC */ - /* It isn't exact so Prec(b) > Prec(a) and the last - Prec(b)-Prec(a) bits of `b' are not zeros. - Which means that removing c from b can't generate a carry - execpt in case of even rounding. - In all other case the result and the inexact flag should be - correct (We can't have an exact result). - In case of EVEN rounding: + /* It isn't exact, so PREC(b) > PREC(a) and the last + PREC(b)-PREC(a) bits of b are not all zeros. + Subtracting c from b will not have an effect on the rounding + except in case of a midpoint in the round-to-nearest mode, + when the even rounding was done away from zero instead of + toward zero. + In case of even rounding: 1.BBBBBBBBBBBBBx10 - 1.CCCCCCCCCCCC - = 1.BBBBBBBBBBBBBx01 Rounded to Prec(b) - = 1.BBBBBBBBBBBBBx Nearest / Rounded to Prec(a) + = 1.BBBBBBBBBBBBBx01 Rounded to PREC(b) + = 1.BBBBBBBBBBBBBx Nearest / Rounded to PREC(a) Set gives: 1.BBBBBBBBBBBBB0 if inexact == EVEN_INEX (x == 0) 1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1) - which means we get a wrong rounded result if x==1, - i.e. inexact= MPFR_EVEN_INEX */ - if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX * MPFR_INT_SIGN (a))) - { - if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax)) - mpfr_setmax (a, __gmpfr_emax); - else - mpfr_nexttozero (a); - inexact = -MPFR_INT_SIGN (a); - } - else if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax)) - inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)); - MPFR_RET (inexact); + which means we get a wrong rounded result if x == 1, + i.e. inexact == MPFR_EVEN_INEX (for positive numbers). */ + if (MPFR_LIKELY (inexact != MPFR_EVEN_INEX * MPFR_INT_SIGN (a))) + goto check_overflow; } + /* We need to take the value preceding |a|. We can't use + mpfr_nexttozero due to a possible out-of-range exponent. + But this will allow us to have more specific code. */ + MPFR_LOG_MSG (("correcting the value of a", 0)); + sh = (mpfr_prec_t) an * GMP_NUMB_BITS - aq; + mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh); + if (MPFR_UNLIKELY (MPFR_LIMB_MSB (ap[an-1]) == 0)) + { + MPFR_EXP (a) --; + /* The following is valid whether an = 1 or an > 1. */ + ap[an-1] |= MPFR_LIMB_HIGHBIT; + } + inexact = - MPFR_INT_SIGN (a); + check_overflow: + if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax)) + return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)); + else + MPFR_RET (inexact); } /* reserve a space to store b aligned with the result, i.e. shifted by @@ -633,15 +643,20 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */ exp_a = exp_b - cancel; - if (MPFR_UNLIKELY(exp_a < __gmpfr_emin)) + if (MPFR_UNLIKELY (exp_a < __gmpfr_emin)) { - MPFR_TMP_FREE(marker); + MPFR_TMP_FREE (marker); if (rnd_mode == MPFR_RNDN && (exp_a < __gmpfr_emin - 1 || (inexact >= 0 && mpfr_powerof2_raw (a)))) rnd_mode = MPFR_RNDZ; return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); } + if (MPFR_UNLIKELY (exp_a > __gmpfr_emax)) + { + MPFR_TMP_FREE (marker); + return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)); + } MPFR_SET_EXP (a, exp_a); } else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */ @@ -649,10 +664,11 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) /* in case cancel = 0, add_exp can still be 1, in case b is just below a power of two, c is very small, prec(a) < prec(b), and rnd=away or nearest */ - if (MPFR_UNLIKELY(add_exp && exp_b >= __gmpfr_emax)) + MPFR_ASSERTD (add_exp == 0 || add_exp == 1); + if (MPFR_UNLIKELY (add_exp && exp_b >= __gmpfr_emax)) { - MPFR_TMP_FREE(marker); - return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); + MPFR_TMP_FREE (marker); + return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)); } MPFR_SET_EXP (a, exp_b + add_exp); } -- cgit v1.2.1 From ad988f2360ad2b59b54769b283b8c2f7bf04a13f Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 30 May 2016 14:26:34 +0000 Subject: [src/sub1.c] Forgot a \n in a log message. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10391 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/sub1.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sub1.c b/src/sub1.c index 6a5a910a6..d81b2617d 100644 --- a/src/sub1.c +++ b/src/sub1.c @@ -164,7 +164,7 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) /* We need to take the value preceding |a|. We can't use mpfr_nexttozero due to a possible out-of-range exponent. But this will allow us to have more specific code. */ - MPFR_LOG_MSG (("correcting the value of a", 0)); + MPFR_LOG_MSG (("correcting the value of a\n", 0)); sh = (mpfr_prec_t) an * GMP_NUMB_BITS - aq; mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh); if (MPFR_UNLIKELY (MPFR_LIMB_MSB (ap[an-1]) == 0)) -- cgit v1.2.1 From 64b6a91bfdedfb9d031917298588afe9012266dc Mon Sep 17 00:00:00 2001 From: vlefevre Date: Thu, 2 Jun 2016 15:36:44 +0000 Subject: [tests/tfmma.c] Added test cases where the precision of the result is twice the precision of each input, which can currently involve add1sp.c and sub1sp.c code. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10403 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 109 +++++++++++++++++++++++++++++++++------------------------- 1 file changed, 63 insertions(+), 46 deletions(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index fcd291557..7d4588a3a 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -121,7 +121,7 @@ zero_tests (void) set_emin (MPFR_EMIN_MIN); set_emax (MPFR_EMAX_MAX); - mpfr_inits2 (64, a, b, c, d, res, (mpfr_ptr) 0); + mpfr_inits2 (64, a, b, c, d, (mpfr_ptr) 0); for (i = 0; i <= 4; i++) { switch (i) @@ -141,35 +141,44 @@ zero_tests (void) } RND_LOOP (r) { - int inex; + int j, inex; mpfr_flags_t flags; mpfr_set (b, a, MPFR_RNDN); mpfr_set (c, a, MPFR_RNDN); mpfr_neg (d, a, MPFR_RNDN); - mpfr_clear_flags (); - inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r); - flags = __gmpfr_flags; - if (flags != 0 || inex != 0 || ! mpfr_zero_p (res) || - (r == MPFR_RNDD ? MPFR_IS_POS (res) : MPFR_IS_NEG (res))) + /* We also want to test cases where the precision of the + result is twice the precision of each input, in case + add1sp.c and/or sub1sp.c could be involved. */ + for (j = 1; j <= 2; j++) { - printf ("Error in zero_tests on i = %d, %s\n", - i, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); - printf ("Expected %c0, inex = 0\n", r == MPFR_RNDD ? '-' : '+'); - printf ("Got "); - if (MPFR_IS_POS (res)) - printf ("+"); - mpfr_out_str (stdout, 16, 0, res, MPFR_RNDN); - printf (", inex = %d\n", inex); - printf ("Expected flags:"); - flags_out (0); - printf ("Got flags: "); - flags_out (flags); - exit (1); - } - } - } - mpfr_clears (a, b, c, d, res, (mpfr_ptr) 0); + mpfr_init2 (res, 64 * j); + mpfr_clear_flags (); + inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r); + flags = __gmpfr_flags; + if (flags != 0 || inex != 0 || ! mpfr_zero_p (res) || + (r == MPFR_RNDD ? MPFR_IS_POS (res) : MPFR_IS_NEG (res))) + { + printf ("Error in zero_tests on i = %d, %s\n", + i, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); + printf ("Expected %c0, inex = 0\n", + r == MPFR_RNDD ? '-' : '+'); + printf ("Got "); + if (MPFR_IS_POS (res)) + printf ("+"); + mpfr_out_str (stdout, 16, 0, res, MPFR_RNDN); + printf (", inex = %d\n", inex); + printf ("Expected flags:"); + flags_out (0); + printf ("Got flags: "); + flags_out (flags); + exit (1); + } + mpfr_clear (res); + } /* j */ + } /* r */ + } /* i */ + mpfr_clears (a, b, c, d, (mpfr_ptr) 0); set_emin (emin); set_emax (emax); @@ -181,7 +190,7 @@ max_tests (void) mpfr_exp_t emin, emax; mpfr_t x, y1, y2; int r; - int inex1, inex2; + int i, inex1, inex2; mpfr_flags_t flags1, flags2; emin = mpfr_get_emin (); @@ -189,32 +198,40 @@ max_tests (void) set_emin (MPFR_EMIN_MIN); set_emax (MPFR_EMAX_MAX); - mpfr_inits2 (64, x, y1, y2, (mpfr_ptr) 0); + mpfr_init2 (x, 64); mpfr_setmax (x, MPFR_EMAX_MAX); flags1 = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; RND_LOOP (r) { - inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r); - mpfr_clear_flags (); - inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r); - flags2 = __gmpfr_flags; - if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && - mpfr_equal_p (y1, y2))) + /* We also want to test cases where the precision of the + result is twice the precision of each input, in case + add1sp.c and/or sub1sp.c could be involved. */ + for (i = 1; i <= 2; i++) { - printf ("Error in max_tests for %s\n", - mpfr_print_rnd_mode ((mpfr_rnd_t) r)); - printf ("Expected "); - mpfr_dump (y1); - printf (" with inex = %d, flags =", inex1); - flags_out (flags1); - printf ("Got "); - mpfr_dump (y2); - printf (" with inex = %d, flags =", inex2); - flags_out (flags2); - exit (1); - } - } - mpfr_clears (x, y1, y2, (mpfr_ptr) 0); + mpfr_inits2 (64 * i, y1, y2, (mpfr_ptr) 0); + inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r); + mpfr_clear_flags (); + inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r); + flags2 = __gmpfr_flags; + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (y1, y2))) + { + printf ("Error in max_tests for %s\n", + mpfr_print_rnd_mode ((mpfr_rnd_t) r)); + printf ("Expected "); + mpfr_dump (y1); + printf (" with inex = %d, flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (y2); + printf (" with inex = %d, flags =", inex2); + flags_out (flags2); + exit (1); + } + mpfr_clears (y1, y2, (mpfr_ptr) 0); + } /* i */ + } /* r */ + mpfr_clear (x); set_emin (emin); set_emax (emax); -- cgit v1.2.1 From ce6db01917269702e951ad9310f68dce26f4f3db Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 3 Jun 2016 09:05:52 +0000 Subject: [src/mpfr-impl.h] UBF support: added MPFR_IS_SINGULAR_OR_UBF and MPFR_ARE_SINGULAR_OR_UBF macros. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10415 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/mpfr-impl.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index e0167eff6..d215f5242 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -909,6 +909,7 @@ typedef intmax_t mpfr_eexp_t; #define MPFR_IS_FP(x) (!MPFR_IS_NAN(x) && !MPFR_IS_INF(x)) #define MPFR_IS_SINGULAR(x) (MPFR_EXP(x) <= MPFR_EXP_INF) +#define MPFR_IS_SINGULAR_OR_UBF(x) (MPFR_EXP(x) <= MPFR_EXP_UBF) #define MPFR_IS_PURE_FP(x) \ (!MPFR_IS_SINGULAR(x) && \ (MPFR_ASSERTD (MPFR_EXP (x) >= MPFR_EMIN_MIN && \ @@ -924,6 +925,10 @@ typedef intmax_t mpfr_eexp_t; #define MPFR_ARE_SINGULAR(x,y) \ (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)) || MPFR_UNLIKELY(MPFR_IS_SINGULAR(y))) +#define MPFR_ARE_SINGULAR_OR_UBF(x,y) \ + (MPFR_UNLIKELY(MPFR_IS_SINGULAR_OR_UBF(x)) || \ + MPFR_UNLIKELY(MPFR_IS_SINGULAR_OR_UBF(y))) + #define MPFR_IS_POWER_OF_2(x) \ (mpfr_cmp_ui_2exp ((x), 1, MPFR_GET_EXP (x) - 1) == 0) -- cgit v1.2.1 From 48d7357b564fb208526641a4b415ba20b17e86cb Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 3 Jun 2016 09:21:00 +0000 Subject: [src/{add,sub}.c] Consider UBF numbers as special cases so that mpfr_sub1sp and mpfr_add1sp, which do not support UBF, are never called on UBF numbers. This should also (very slightly) speed up the normal cases. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10416 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/add.c | 28 ++++++++++++++++++++-------- src/sub.c | 30 ++++++++++++++++++++++++------ 2 files changed, 44 insertions(+), 14 deletions(-) diff --git a/src/add.c b/src/add.c index 43a4507d1..969a64d0a 100644 --- a/src/add.c +++ b/src/add.c @@ -31,7 +31,7 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) mpfr_get_prec (c), mpfr_log_prec, c, rnd_mode), ("a[%Pu]=%.*Rg", mpfr_get_prec (a), mpfr_log_prec, a)); - if (MPFR_ARE_SINGULAR(b,c)) + if (MPFR_ARE_SINGULAR_OR_UBF (b, c)) { if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { @@ -59,7 +59,7 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) MPFR_SET_SAME_SIGN(a, c); MPFR_RET(0); /* exact */ } - /* now either b or c is zero */ + /* now both b and c are finite numbers */ else if (MPFR_IS_ZERO(b)) { if (MPFR_IS_ZERO(c)) @@ -78,15 +78,27 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } return mpfr_set (a, c, rnd_mode); } - else + else if (MPFR_IS_ZERO(c)) { - MPFR_ASSERTD(MPFR_IS_ZERO(c)); return mpfr_set (a, b, rnd_mode); } + else + { + MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); + /* mpfr_sub1sp and mpfr_add1sp are not intended to support UBF, + for which optimization is less important. */ + if (MPFR_SIGN(b) != MPFR_SIGN(c)) + return mpfr_sub1 (a, b, c, rnd_mode); + else if (MPFR_EXP_LESS_P (b, c)) + return mpfr_add1 (a, c, b, rnd_mode); + else + return mpfr_add1 (a, b, c, rnd_mode); + } } - MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); - MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); + MPFR_ASSERTD (MPFR_IS_PURE_FP (b)); + MPFR_ASSERTD (MPFR_IS_PURE_FP (c)); if (MPFR_UNLIKELY(MPFR_SIGN(b) != MPFR_SIGN(c))) { /* signs differ, it is a subtraction */ @@ -100,12 +112,12 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) { /* signs are equal, it's an addition */ if (MPFR_LIKELY(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c))) - if (MPFR_EXP_LESS_P (b, c)) + if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) return mpfr_add1sp(a, c, b, rnd_mode); else return mpfr_add1sp(a, b, c, rnd_mode); else - if (MPFR_EXP_LESS_P (b, c)) + if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) return mpfr_add1(a, c, b, rnd_mode); else return mpfr_add1(a, b, c, rnd_mode); diff --git a/src/sub.c b/src/sub.c index f6a4e16d5..c95cfdb0a 100644 --- a/src/sub.c +++ b/src/sub.c @@ -31,7 +31,7 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) mpfr_get_prec (c), mpfr_log_prec, c, rnd_mode), ("a[%Pu]=%.*Rg", mpfr_get_prec (a), mpfr_log_prec, a)); - if (MPFR_ARE_SINGULAR (b,c)) + if (MPFR_ARE_SINGULAR_OR_UBF (b,c)) { if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c)) { @@ -72,15 +72,33 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) else return mpfr_neg (a, c, rnd_mode); } - else + else if (MPFR_IS_ZERO (c)) { - MPFR_ASSERTD (MPFR_IS_ZERO (c)); return mpfr_set (a, b, rnd_mode); } + else + { + MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); + /* mpfr_sub1sp and mpfr_add1sp are not intended to support UBF, + for which optimization is less important. */ + if (MPFR_SIGN(b) == MPFR_SIGN(c)) + return mpfr_sub1 (a, b, c, rnd_mode); + else if (MPFR_EXP_LESS_P (b, c)) + { + int inexact; + rnd_mode = MPFR_INVERT_RND (rnd_mode); + inexact = mpfr_add1 (a, c, b, rnd_mode); + MPFR_CHANGE_SIGN (a); + return -inexact; + } + else + return mpfr_add1 (a, b, c, rnd_mode); + } } - MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); - MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); + MPFR_ASSERTD (MPFR_IS_PURE_FP (b)); + MPFR_ASSERTD (MPFR_IS_PURE_FP (c)); if (MPFR_LIKELY (MPFR_SIGN (b) == MPFR_SIGN (c))) { /* signs are equal, it's a real subtraction */ @@ -92,7 +110,7 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } else { /* signs differ, it's an addition */ - if (MPFR_EXP_LESS_P (b, c)) + if (MPFR_GET_EXP (b) < MPFR_GET_EXP (c)) { /* exchange rounding modes toward +/- infinity */ int inexact; rnd_mode = MPFR_INVERT_RND (rnd_mode); -- cgit v1.2.1 From 5ce3d5abeedd9513022dc9ddb62af75adca99e6d Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 3 Jun 2016 12:44:20 +0000 Subject: [tests/tfmma.c] Replaced precision 64 by GMP_NUMB_BITS (to make sure that the add1sp and sub1sp conditions are satisfied with the current src code, these tests may rely on the fact that there are no trailing bits, i.e. that the precision is a multiple of GMP_NUMB_BITS). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10418 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index 7d4588a3a..a5b473ffe 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -121,7 +121,7 @@ zero_tests (void) set_emin (MPFR_EMIN_MIN); set_emax (MPFR_EMAX_MAX); - mpfr_inits2 (64, a, b, c, d, (mpfr_ptr) 0); + mpfr_inits2 (GMP_NUMB_BITS, a, b, c, d, (mpfr_ptr) 0); for (i = 0; i <= 4; i++) { switch (i) @@ -152,7 +152,7 @@ zero_tests (void) add1sp.c and/or sub1sp.c could be involved. */ for (j = 1; j <= 2; j++) { - mpfr_init2 (res, 64 * j); + mpfr_init2 (res, GMP_NUMB_BITS * j); mpfr_clear_flags (); inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r); flags = __gmpfr_flags; @@ -198,7 +198,7 @@ max_tests (void) set_emin (MPFR_EMIN_MIN); set_emax (MPFR_EMAX_MAX); - mpfr_init2 (x, 64); + mpfr_init2 (x, GMP_NUMB_BITS); mpfr_setmax (x, MPFR_EMAX_MAX); flags1 = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; RND_LOOP (r) @@ -208,7 +208,7 @@ max_tests (void) add1sp.c and/or sub1sp.c could be involved. */ for (i = 1; i <= 2; i++) { - mpfr_inits2 (64 * i, y1, y2, (mpfr_ptr) 0); + mpfr_inits2 (GMP_NUMB_BITS * i, y1, y2, (mpfr_ptr) 0); inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r); mpfr_clear_flags (); inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r); -- cgit v1.2.1 From 675f1ca65aa8099fc78f4cee13ed17436652cf3c Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 3 Jun 2016 12:47:12 +0000 Subject: [tests/tfmma.c] Test (1/2)x + (1/2)x = x near underflow. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10419 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index a5b473ffe..b2de7722d 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -356,6 +356,65 @@ overflow_tests (void) set_emax (old_emax); } +/* (1/2)x + (1/2)x = x tested near underflow */ +static void +half_plus_half (void) +{ + mpfr_exp_t emin; + mpfr_t h, x1, x2, y; + int neg, r, i, inex; + mpfr_flags_t flags; + + emin = mpfr_get_emin (); + set_emin (MPFR_EMIN_MIN); + mpfr_inits2 (4, h, x1, (mpfr_ptr) 0); + mpfr_init2 (x2, GMP_NUMB_BITS); + mpfr_set_ui_2exp (h, 1, -1, MPFR_RNDN); + + for (mpfr_setmin (x1, __gmpfr_emin); + MPFR_GET_EXP (x1) < __gmpfr_emin + 2; + mpfr_nextabove (x1)) + { + inex = mpfr_set (x2, x1, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + for (neg = 0; neg < 2; neg++) + { + RND_LOOP (r) + { + /* We also want to test cases where the precision of the + result is twice the precision of each input, in case + add1sp.c and/or sub1sp.c could be involved. */ + for (i = 1; i <= 2; i++) + { + mpfr_init2 (y, GMP_NUMB_BITS * i); + mpfr_clear_flags (); + inex = mpfr_fmma (y, h, x2, h, x2, (mpfr_rnd_t) r); + flags = __gmpfr_flags; + if (! (flags == 0 && inex == 0 && mpfr_equal_p (y, x2))) + { + printf ("Error in half_plus_half for %s\n", + mpfr_print_rnd_mode (r)); + printf ("Expected "); + mpfr_dump (x2); + printf (" with inex = 0, flags ="); + flags_out (0); + printf ("Got "); + mpfr_dump (y); + printf (" with inex = %d, flags =", inex); + flags_out (flags); + exit (1); + } + mpfr_clear (y); + } + } + mpfr_neg (x2, x2, MPFR_RNDN); + } + } + + mpfr_clears (h, x1, x2, y, (mpfr_ptr) 0); + set_emin (emin); +} + int main (int argc, char *argv[]) { @@ -365,7 +424,7 @@ main (int argc, char *argv[]) zero_tests (); max_tests (); overflow_tests (); - /* TODO: near_underflow_tests (); */ + half_plus_half (); tests_end_mpfr (); return 0; -- cgit v1.2.1 From c9c5c377ce687f0c1dca19aea9f5464719bd8b06 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 3 Jun 2016 13:03:52 +0000 Subject: [tests/tfmma.c] Correction. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10420 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index b2de7722d..8939fa0f0 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -411,7 +411,7 @@ half_plus_half (void) } } - mpfr_clears (h, x1, x2, y, (mpfr_ptr) 0); + mpfr_clears (h, x1, x2, (mpfr_ptr) 0); set_emin (emin); } -- cgit v1.2.1 From 6344605a1d9860deeedea9773081173351fa82ce Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 3 Jun 2016 13:05:57 +0000 Subject: [src/add1.c] Completed UBF support. Note: due to the restriction on the exponent values, diff_exp does not need to be the unsigned integer type mpfr_uexp_t (mpfr_exp_t is sufficient). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10421 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/add1.c | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/add1.c b/src/add1.c index 437a1d711..37284a9b1 100644 --- a/src/add1.c +++ b/src/add1.c @@ -30,9 +30,8 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) mp_limb_t *ap, *bp, *cp; mpfr_prec_t aq, bq, cq, aq2; mp_size_t an, bn, cn; - mpfr_exp_t difw, exp; + mpfr_exp_t difw, exp, diff_exp; int sh, rb, fb, inex; - mpfr_uexp_t diff_exp; MPFR_TMP_DECL(marker); MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); @@ -80,9 +79,15 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) MPFR_SET_SAME_SIGN(a, b); MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN(b)); /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */ - /* Note: exponents can be negative, but the unsigned subtraction is - a modular subtraction, so that one gets the correct result. */ - diff_exp = (mpfr_uexp_t) exp - MPFR_GET_EXP(c); + if (MPFR_UNLIKELY (MPFR_IS_UBF (c))) + { + MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX); + diff_exp = mpfr_ubf_diff_exp (b, c); + } + else + diff_exp = exp - MPFR_GET_EXP (c); + + MPFR_ASSERTD (diff_exp >= 0); /* * 1. Compute the significant part A', the non-significant bits of A -- cgit v1.2.1 From 7309f2c16032efafedaaa22dbf2b6a0a396d0876 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 3 Jun 2016 13:28:45 +0000 Subject: [tests/tfmma.c] Forgot a cast to mpfr_rnd_t for C++ compatibility. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10422 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tfmma.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tfmma.c b/tests/tfmma.c index 8939fa0f0..c45995a57 100644 --- a/tests/tfmma.c +++ b/tests/tfmma.c @@ -393,7 +393,7 @@ half_plus_half (void) if (! (flags == 0 && inex == 0 && mpfr_equal_p (y, x2))) { printf ("Error in half_plus_half for %s\n", - mpfr_print_rnd_mode (r)); + mpfr_print_rnd_mode ((mpfr_rnd_t) r)); printf ("Expected "); mpfr_dump (x2); printf (" with inex = 0, flags ="); -- cgit v1.2.1 From 447f45d74a98367210e108ce7bde60c2d945c840 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 6 Jun 2016 16:19:21 +0000 Subject: [acinclude.m4] When checking if __float128 is available, we now also check whether C99 constants (in particular the __float128 ones, such as 0x1.fp+16383q) are supported since this is now required with the __float128 support. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10440 280ebfd0-de03-0410-8827-d642c229c3f4 --- acinclude.m4 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/acinclude.m4 b/acinclude.m4 index be917becb..1fefe840e 100644 --- a/acinclude.m4 +++ b/acinclude.m4 @@ -592,15 +592,17 @@ Please use another compiler or build MPFR without --enable-decimal-float.]) fi]) fi -dnl Check if __float128 is available +dnl Check if __float128 is available. We also require the compiler +dnl to support C99 constants (this prevents the __float128 support +dnl with GCC's -std=c90, but who cares?). if test "$enable_float128" != no; then - AC_MSG_CHECKING(if compiler knows __float128) - AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[__float128 x;]])], + AC_MSG_CHECKING(if compiler knows __float128 with C99 constants) + AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[__float128 x = 0x1.fp+16383q;]])], [AC_MSG_RESULT(yes) AC_DEFINE([MPFR_WANT_FLOAT128],1,[Build float128 functions])], [AC_MSG_RESULT(no) if test "$enable_float128" = yes; then - AC_MSG_ERROR([compiler doesn't know __float128 + AC_MSG_ERROR([compiler doesn't know __float128 with C99 constants Please use another compiler or build MPFR without --enable-float128.]) fi]) fi -- cgit v1.2.1 From b91d95a975f59c200fd8f455064a19dd2c588173 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 6 Jun 2016 23:09:17 +0000 Subject: [tests/tlgamma.c] Added a test causing a "too much memory" error with tcc 0.9.27~git20151227.933c223-1 (there's already one in special(), but this one is a simpler, standalone test). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10442 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tlgamma.c | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/tests/tlgamma.c b/tests/tlgamma.c index 9d3689e95..a97daee20 100644 --- a/tests/tlgamma.c +++ b/tests/tlgamma.c @@ -385,11 +385,30 @@ mpfr_lgamma1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) return mpfr_lgamma (y, &sign, x, r); } +/* Since r10377, the following test causes a "too much memory" error + when MPFR is built with Debian's tcc 0.9.27~git20151227.933c223-1 + on x86_64. */ +static void +tcc_bug20160606 (void) +{ + mpfr_t x, y; + int sign; + + mpfr_init2 (x, 53); + mpfr_init2 (y, 53); + mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN); + mpfr_lgamma (y, &sign, x, MPFR_RNDN); + mpfr_clear (x); + mpfr_clear (y); +} + int main (void) { tests_start_mpfr (); + tcc_bug20160606 (); + special (); test_generic (MPFR_PREC_MIN, 100, 2); -- cgit v1.2.1 From 11851de31420e372e9885ffa663d5557599c3a48 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 6 Jun 2016 23:43:44 +0000 Subject: [src/uceil_log2.c] Correction in __gmpfr_ceil_log2, avoiding an incorrect result with tcc: x.s.exp is declared as an unsigned bit-field, so that tcc considers that x.s.exp - 1023 is unsigned. However, since all the values of x.s.exp are representable in an int, according to the integer promotion rules, x.s.exp should be converted to an int, so that the subtraction is signed. So, this appears to be a bug in tcc. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10443 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/uceil_log2.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/uceil_log2.c b/src/uceil_log2.c index bfb799cae..27aaca353 100644 --- a/src/uceil_log2.c +++ b/src/uceil_log2.c @@ -33,7 +33,7 @@ __gmpfr_ceil_log2 (double d) union mpfr_ieee_double_extract x; x.d = d; - exp = x.s.exp - 1023; + exp = (long) x.s.exp - 1023; x.s.exp = 1023; /* value for 1 <= d < 2 */ if (x.d != 1.0) /* d: not a power of two? */ exp++; -- cgit v1.2.1 From 8ea8f53d9622c5929cca06513318efcb8acd68b3 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 6 Jun 2016 23:45:31 +0000 Subject: [tests/tlgamma.c] Updated comment of the test added in r10442. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10444 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tlgamma.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/tlgamma.c b/tests/tlgamma.c index a97daee20..bebe85848 100644 --- a/tests/tlgamma.c +++ b/tests/tlgamma.c @@ -387,7 +387,8 @@ mpfr_lgamma1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) /* Since r10377, the following test causes a "too much memory" error when MPFR is built with Debian's tcc 0.9.27~git20151227.933c223-1 - on x86_64. */ + on x86_64. This was due to a bug in __gmpfr_ceil_log2, now fixed + in r10443. */ static void tcc_bug20160606 (void) { -- cgit v1.2.1 From 8d28bfa8394b1b96703ba2e29453ea7092ef7922 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 7 Jun 2016 00:18:47 +0000 Subject: [tests/tlgamma.c] Corrected updated comment from r10444. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10445 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tlgamma.c | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/tlgamma.c b/tests/tlgamma.c index bebe85848..eca1e4b27 100644 --- a/tests/tlgamma.c +++ b/tests/tlgamma.c @@ -387,8 +387,10 @@ mpfr_lgamma1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) /* Since r10377, the following test causes a "too much memory" error when MPFR is built with Debian's tcc 0.9.27~git20151227.933c223-1 - on x86_64. This was due to a bug in __gmpfr_ceil_log2, now fixed - in r10443. */ + on x86_64. The problem came from __gmpfr_ceil_log2, now fixed in + r10443 (according to the integer promotion rules, this appeared to + be a bug in tcc, not in MPFR; however replying on such an obscure + rule was not a good idea). */ static void tcc_bug20160606 (void) { -- cgit v1.2.1 From a0e724e6055d10541294e97632ba01db71e81a55 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 7 Jun 2016 00:23:13 +0000 Subject: [src/uceil_log2.c] Added a comment. Cosmetic changes. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10446 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/uceil_log2.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/uceil_log2.c b/src/uceil_log2.c index 27aaca353..7ee5166ae 100644 --- a/src/uceil_log2.c +++ b/src/uceil_log2.c @@ -33,6 +33,8 @@ __gmpfr_ceil_log2 (double d) union mpfr_ieee_double_extract x; x.d = d; + /* The cast below is useless in theory, but let us not depend on the + integer promotion rules (for instance, tcc is currently wrong). */ exp = (long) x.s.exp - 1023; x.s.exp = 1023; /* value for 1 <= d < 2 */ if (x.d != 1.0) /* d: not a power of two? */ @@ -42,19 +44,19 @@ __gmpfr_ceil_log2 (double d) double m; if (d < 0.0) - return __gmpfr_floor_log2(-d)+1; + 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 ) + for (m = 1.0; m < d; m *= 2.0) exp++; } else { exp = 1; - for( m= 1.0 ; m >= d ; m *= (1.0/2.0) ) + for (m = 1.0; m >= d; m *= 0.5) exp--; } #endif /* _MPFR_IEEE_FLOATS */ -- cgit v1.2.1 From 45abf6eaaf2b4986c138cc1a961466ee5d239ba2 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 7 Jun 2016 00:44:17 +0000 Subject: [tests/tlgamma.c] Typo in comment. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10447 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tlgamma.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tlgamma.c b/tests/tlgamma.c index eca1e4b27..3d686a6fc 100644 --- a/tests/tlgamma.c +++ b/tests/tlgamma.c @@ -389,7 +389,7 @@ mpfr_lgamma1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) when MPFR is built with Debian's tcc 0.9.27~git20151227.933c223-1 on x86_64. The problem came from __gmpfr_ceil_log2, now fixed in r10443 (according to the integer promotion rules, this appeared to - be a bug in tcc, not in MPFR; however replying on such an obscure + be a bug in tcc, not in MPFR; however relying on such an obscure rule was not a good idea). */ static void tcc_bug20160606 (void) -- cgit v1.2.1 From e31290a255dd7c7e4ef5dac84d6a38e215df66f3 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 7 Jun 2016 07:26:52 +0000 Subject: [tests/tprintf.c] Removed tests of native %'g and %'f (from r8292) as the ' flag is an extension from Single UNIX Specification and in particular, they fail with MinGW under Wine. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10449 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tprintf.c | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/tests/tprintf.c b/tests/tprintf.c index aa5a540c6..fbc32734f 100644 --- a/tests/tprintf.c +++ b/tests/tprintf.c @@ -428,7 +428,6 @@ test_locale (void) int i; char *s = NULL; mpfr_t x; - double y; int count; for(i = 0; i < numberof(tab_locale) && s == NULL; i++) @@ -439,12 +438,11 @@ test_locale (void) mpfr_init2 (x, 113); mpfr_set_ui (x, 10000, MPFR_RNDN); - y = 100000; - count = mpfr_printf ("(1) 10000=%'Rg 100000=%'g \n", x, y); - check_length (10000, count, 33, d); - count = mpfr_printf ("(2) 10000=%'Rf 100000=%'f \n", x, y); - check_length (10001, count, 47, d); + count = mpfr_printf ("(1) 10000=%'Rg \n", x); + check_length (10000, count, 18, d); + count = mpfr_printf ("(2) 10000=%'Rf \n", x); + check_length (10001, count, 25, d); mpfr_clear (x); } -- cgit v1.2.1 From c22098e9b90dd2c7ef6ac7040e887b6b804cc603 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 7 Jun 2016 08:31:31 +0000 Subject: [src/mpfr-impl.h] Updated a comment (removing a FIXME). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/ubf@10451 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/mpfr-impl.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index d215f5242..cc035d346 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -2268,7 +2268,9 @@ extern "C" { makes mpfr_t to UBF conversion unnecessary. When an input of a public function is a UBF, the semantic remains internal to MPFR and can change in the future. - Note that functions used for logging need to support UBF. *FIXME* */ + Note that functions used for logging need to support UBF (currently + done by printing that a number is a UBF, as it may be difficult to + do more without significant changes). */ typedef struct { mpfr_prec_t _mpfr_prec; -- cgit v1.2.1 From b81607df91962ec18d3adb313e1e9394ba9c2c91 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Tue, 7 Jun 2016 10:22:03 +0000 Subject: [tools/mbench/Makefile] Added multiarch support for GMP. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10452 280ebfd0-de03-0410-8827-d642c229c3f4 --- tools/mbench/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/mbench/Makefile b/tools/mbench/Makefile index 94dcc9314..06092c18e 100644 --- a/tools/mbench/Makefile +++ b/tools/mbench/Makefile @@ -60,7 +60,7 @@ INCLUDES=-I. -I$(MPFR_INCLUDE) -I$(MPFR) -I$(GMP)/include/ -I$(GMP) endif # search first for real install, then for build directory -LIBS=`(test -f $(MPFR_LIB)/libmpfr.a && echo $(MPFR_LIB)/libmpfr.a)` `(test -f $(MPFR)/.libs/libmpfr.a && echo $(MPFR)/.libs/libmpfr.a)` `(test -f $(MPFR)/src/.libs/libmpfr.a && echo $(MPFR)/src/.libs/libmpfr.a)` `(test -f $(GMP)/lib/libgmp.a && echo $(GMP)/lib/libgmp.a)` `(test -f $(GMP)/.libs/libgmp.a && echo $(GMP)/.libs/libgmp.a)` +LIBS=`(test -f $(MPFR_LIB)/libmpfr.a && echo $(MPFR_LIB)/libmpfr.a)` `(test -f $(MPFR)/.libs/libmpfr.a && echo $(MPFR)/.libs/libmpfr.a)` `(test -f $(MPFR)/src/.libs/libmpfr.a && echo $(MPFR)/src/.libs/libmpfr.a)` `(test -f $(GMP)/lib/libgmp.a && echo $(GMP)/lib/libgmp.a)` `(test -f $(GMP)/lib/$$(gcc -print-multiarch 2>/dev/null)/libgmp.a && echo $(GMP)/lib/$$(gcc -print-multiarch 2>/dev/null)/libgmp.a)` `(test -f $(GMP)/.libs/libgmp.a && echo $(GMP)/.libs/libgmp.a)` XLIBS=`test -f $(PARI)/lib/libpari.a && echo $(PARI)/lib/libpari.a` \ `test -f $(NTL)/lib/libntl.a && echo $(NTL)/lib/libntl.a` \ -- cgit v1.2.1 From 45c19367581741b348bee5cd8deb9ce7b60e0dc7 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Wed, 8 Jun 2016 12:06:57 +0000 Subject: Value coverage for tsum: After r9984 to differentiate sq > MPFR_PREC_MIN and sq == MPFR_PREC_MIN, 42 tests were not done for sq == MPFR_PREC_MIN on a 64-bit machine. In the check4 test, changed a 2 to MPFR_PREC_MIN. This reduces to 12 tests that are not done. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10456 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tsum.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tsum.c b/tests/tsum.c index 78946904c..f35ba314a 100644 --- a/tests/tsum.c +++ b/tests/tsum.c @@ -630,7 +630,7 @@ check4 (void) /* No GNU style for the many nested loops... */ for (k = 1; k <= 64; k++) { mpfr_set_si_2exp (t[0], -1, k, MPFR_RNDN); - for (n = k + 2; n <= k + 65; n++) { + for (n = k + MPFR_PREC_MIN; n <= k + 65; n++) { prec = n - k; mpfr_set_prec (sum1, prec); mpfr_set_prec (sum2, prec); -- cgit v1.2.1 From bbd084c4b3e5c4802bf872afe11371a682713fdb Mon Sep 17 00:00:00 2001 From: vlefevre Date: Wed, 8 Jun 2016 12:28:43 +0000 Subject: [tests/tsum.c] Improved check3: * Also do the tests with output precision MPFR_PREC_MIN; this completes the value coverage for tsum. * Also compare the flags. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10457 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tsum.c | 69 +++++++++++++++++++++++++++++++++++------------------------- 1 file changed, 40 insertions(+), 29 deletions(-) diff --git a/tests/tsum.c b/tests/tsum.c index f35ba314a..8c65d38ec 100644 --- a/tests/tsum.c +++ b/tests/tsum.c @@ -518,18 +518,17 @@ check2 (void) * t[17] = 2^(17*9+1) * j for -4 <= j <= 4. * t[18] = 2^(-1) * k for -1 <= k <= 1. * t[19] = 2^(-17*8) * m for -3 <= m <= 3. - * prec = 17*9+4 + * prec = MPFR_PREC_MIN and 17*9+4 */ static void check3 (void) { mpfr_t sum1, sum2, s1, s2, s3, s4, t[20]; mpfr_ptr p[20]; - int i, s, j, k, m, r, inex1, inex2; - int prec = 17*9+4; + mpfr_flags_t flags1, flags2; + int i, s, j, k, m, q, r, inex1, inex2; + int prec[2] = { MPFR_PREC_MIN, 17*9+4 }; - mpfr_init2 (sum1, prec); - mpfr_init2 (sum2, prec); mpfr_init2 (s1, 17*17); mpfr_init2 (s2, 17*17+4); mpfr_init2 (s3, 17*17+4); @@ -564,39 +563,51 @@ check3 (void) mpfr_set_si_2exp (t[19], m, -17*8, MPFR_RNDN); inex1 = mpfr_add (s4, s3, t[19], MPFR_RNDN); MPFR_ASSERTN (inex1 == 0); - RND_LOOP (r) + for (q = 0; q < 2; q++) { - inex1 = mpfr_set (sum1, s4, (mpfr_rnd_t) r); - inex2 = mpfr_sum (sum2, p, 20, (mpfr_rnd_t) r); - MPFR_ASSERTN (mpfr_check (sum1)); - MPFR_ASSERTN (mpfr_check (sum2)); - if (!(mpfr_equal_p (sum1, sum2) && - SAME_SIGN (inex1, inex2))) + mpfr_inits2 (prec[q], sum1, sum2, (mpfr_ptr) 0); + RND_LOOP (r) { - printf ("Error in check3 on %s, " - "s = %d, j = %d, k = %d, m = %d\n", - mpfr_print_rnd_mode ((mpfr_rnd_t) r), - s, j, k, m); - printf ("Expected "); - mpfr_dump (sum1); - printf ("with inex = %d\n", inex1); - printf ("Got "); - mpfr_dump (sum2); - printf ("with inex = %d\n", inex2); - exit (1); + mpfr_clear_flags (); + inex1 = mpfr_set (sum1, s4, (mpfr_rnd_t) r); + flags1 = __gmpfr_flags; + mpfr_clear_flags (); + inex2 = mpfr_sum (sum2, p, 20, (mpfr_rnd_t) r); + flags2 = __gmpfr_flags; + MPFR_ASSERTN (mpfr_check (sum1)); + MPFR_ASSERTN (mpfr_check (sum2)); + if (!(mpfr_equal_p (sum1, sum2) && + SAME_SIGN (inex1, inex2) && + flags1 == flags2)) + { + printf ("Error in check3 on %s, " + "s = %d, j = %d, k = %d, m = %d\n", + mpfr_print_rnd_mode ((mpfr_rnd_t) r), + s, j, k, m); + printf ("Expected "); + mpfr_dump (sum1); + printf ("with inex = %d and flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (sum2); + printf ("with inex = %d and flags =", inex2); + flags_out (flags2); + exit (1); + } } - } - } - } - } + mpfr_clears (sum1, sum2, (mpfr_ptr) 0); + } /* q */ + } /* m */ + } /* k */ + } /* j */ for (i = 0; i < 17; i++) mpfr_neg (t[i], t[i], MPFR_RNDN); mpfr_neg (s1, s1, MPFR_RNDN); - } + } /* s */ for (i = 0; i < 20; i++) mpfr_clear (t[i]); - mpfr_clears (sum1, sum2, s1, s2, s3, s4, (mpfr_ptr) 0); + mpfr_clears (s1, s2, s3, s4, (mpfr_ptr) 0); } /* Test of s * (q * 2^(n-1) - 2^k) + h + i * 2^(-2) + j * 2^(-2) -- cgit v1.2.1 From 02571b05f694673c75ce656572685f4ac422c0d5 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Thu, 9 Jun 2016 11:51:52 +0000 Subject: [tests/memory.c] The MPFR_TESTS_MEMORY_LIMIT environment variable can now contain an integer specifying the memory limit for the tests, or 0 for unlimited, the default still being 2^22 = 4 MB. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10459 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/memory.c | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/tests/memory.c b/tests/memory.c index ac478117d..0a52b7799 100644 --- a/tests/memory.c +++ b/tests/memory.c @@ -42,6 +42,7 @@ struct header { static struct header *tests_memory_list; static size_t tests_total_size = 0; +static size_t tests_memory_limit = 1UL << 22; /* default memory limit */ MPFR_LOCK_DECL(mpfr_lock_memory) static void * @@ -103,14 +104,11 @@ tests_memory_valid (void *ptr) } */ -/* FIXME: we might have an environment variable MPFR_TESTS_NO_MEMORY_LIMIT - to disable the memory limit check, for example if we want to compute - zillions of digits of Pi with ./tconst_pi 1000000000000 */ static void tests_addsize (size_t size) { tests_total_size += size; - if (tests_total_size > 1UL << 22) + if (tests_total_size > tests_memory_limit) { /* The total size taken by MPFR on the heap is more than 4 MB: either a bug or a huge inefficiency. */ @@ -244,8 +242,23 @@ tests_free (void *ptr, size_t size) void tests_memory_start (void) { + char *p; + tests_memory_list = NULL; mp_set_memory_functions (tests_allocate, tests_reallocate, tests_free); + + /* The memory limit can be changed with the MPFR_TESTS_MEMORY_LIMIT + environment variable. This is normally not necessary (a failure + would mean a bug), thus not recommended, for "make check". But + some test programs can take arguments for particular tests, which + may need more memory. */ + p = getenv ("MPFR_TESTS_MEMORY_LIMIT"); + if (p != NULL) + { + tests_memory_limit = strtoul (p, NULL, 0); + if (tests_memory_limit == 0) + tests_memory_limit = (size_t) -1; /* no memory limit */ + } } void -- cgit v1.2.1 From 73fc3d48eb2ec842b670bac82199912fd382a6b5 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Thu, 9 Jun 2016 12:15:14 +0000 Subject: Other changes concerning the memory limit for the tests. * Export the tests_memory_limit variable so that it can be accessed in test programs (e.g. read by tversion as mentioned below, or modified by a test program when executed with particular arguments). * In tversion, print a warning when the memory limit has been modified (with the MPFR_TESTS_MEMORY_LIMIT environment variable). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10460 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/memory.c | 14 ++++++++------ tests/mpfr-test.h | 3 +++ tests/tversion.c | 5 +++++ 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/tests/memory.c b/tests/memory.c index 0a52b7799..640b2b5d9 100644 --- a/tests/memory.c +++ b/tests/memory.c @@ -40,9 +40,16 @@ struct header { struct header *next; }; +/* The memory limit can be changed with the MPFR_TESTS_MEMORY_LIMIT + environment variable. This is normally not necessary (a failure + would mean a bug), thus not recommended, for "make check". But + some test programs can take arguments for particular tests, which + may need more memory. This variable is exported, so that such + programs may also change the memory limit. */ +size_t tests_memory_limit = DEFAULT_MEMORY_LIMIT; + static struct header *tests_memory_list; static size_t tests_total_size = 0; -static size_t tests_memory_limit = 1UL << 22; /* default memory limit */ MPFR_LOCK_DECL(mpfr_lock_memory) static void * @@ -247,11 +254,6 @@ tests_memory_start (void) tests_memory_list = NULL; mp_set_memory_functions (tests_allocate, tests_reallocate, tests_free); - /* The memory limit can be changed with the MPFR_TESTS_MEMORY_LIMIT - environment variable. This is normally not necessary (a failure - would mean a bug), thus not recommended, for "make check". But - some test programs can take arguments for particular tests, which - may need more memory. */ p = getenv ("MPFR_TESTS_MEMORY_LIMIT"); if (p != NULL) { diff --git a/tests/mpfr-test.h b/tests/mpfr-test.h index 221375386..616ef4119 100644 --- a/tests/mpfr-test.h +++ b/tests/mpfr-test.h @@ -79,6 +79,9 @@ extern "C" { int test_version (void); +/* Memory handling */ +#define DEFAULT_MEMORY_LIMIT (1UL << 22) +extern size_t tests_memory_limit; void tests_memory_start (void); void tests_memory_end (void); diff --git a/tests/tversion.c b/tests/tversion.c index b49eba0cd..6bd87d47b 100644 --- a/tests/tversion.c +++ b/tests/tversion.c @@ -219,6 +219,11 @@ main (void) tests_start_mpfr (); if (locale != NULL) printf ("[tversion] Locale: %s\n", locale); + /* The memory limit should not be changed for "make check". + The warning below signals a possible user mistake. */ + if (tests_memory_limit != DEFAULT_MEMORY_LIMIT) + printf ("[tversion] Warning! Memory limit changed to %zu\n", + tests_memory_limit); tests_end_mpfr (); return err; -- cgit v1.2.1 From 180063b3f84604c573e07d8a009e2e021ed49f71 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Thu, 9 Jun 2016 12:41:16 +0000 Subject: [doc/README.dev] Document environment variables that affect the tests. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10461 280ebfd0-de03-0410-8827-d642c229c3f4 --- doc/README.dev | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/doc/README.dev b/doc/README.dev index 4d9780df0..5f7279392 100644 --- a/doc/README.dev +++ b/doc/README.dev @@ -466,6 +466,35 @@ List of macros used for checking MPFR: =========================================================================== +Environment variables that affect the tests: + ++ GMP_CHECK_RANDOMIZE: Seed for the random functions, except for 0 or 1, + in which case a random (time based) seed is used. + By default, a fixed seed is used. Only developers + and testers should change the seed. + ++ MPFR_CHECK_LARGEMEM: Define to enable expensive tests. + ++ MPFR_CHECK_LIBC_PRINTF: Define to enable comparisons with the printf + function of the C library. These comparisons are + disabled by default as failures could be due to + the C library itself on some machines, and they + do not affect MPFR. + ++ MPFR_DEBUG_BADCASES: For debugging (see tests.c, function bad_cases). + ++ MPFR_SUSPICIOUS_OVERFLOW: Define to check suspicious overflow in the + generic tests (tgeneric.c). For developers and + testers. + ++ MPFR_TESTS_MEMORY_LIMIT: The memory limit for the tests (default is + 2^22 = 4 MB). Set to 0 for unlimited. + ++ MPFR_TESTS_TIMEOUT: When timeout in the tests is enabled, this + overrides the value of the macro. + +=========================================================================== + Before testing any macro in a .c file, one needs: #ifdef HAVE_CONFIG_H -- cgit v1.2.1 From 5080e635a3b96781ee010c1e94ccea27b8e33ef7 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Thu, 9 Jun 2016 14:09:11 +0000 Subject: [tests/tversion.c] Do not use "%zu" with printf (added in r10460). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10462 280ebfd0-de03-0410-8827-d642c229c3f4 --- tests/tversion.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/tversion.c b/tests/tversion.c index 6bd87d47b..45eed5042 100644 --- a/tests/tversion.c +++ b/tests/tversion.c @@ -220,10 +220,12 @@ main (void) if (locale != NULL) printf ("[tversion] Locale: %s\n", locale); /* The memory limit should not be changed for "make check". - The warning below signals a possible user mistake. */ + The warning below signals a possible user mistake. + Do not use "%zu" because it is not available in C90; + the type mpfr_ueexp_t should be sufficiently large. */ if (tests_memory_limit != DEFAULT_MEMORY_LIMIT) - printf ("[tversion] Warning! Memory limit changed to %zu\n", - tests_memory_limit); + printf ("[tversion] Warning! Memory limit changed to %" MPFR_EXP_FSPEC + "u\n", (mpfr_ueexp_t) tests_memory_limit); tests_end_mpfr (); return err; -- cgit v1.2.1 From 015d5c91210a13979984b778635566534ca3390c Mon Sep 17 00:00:00 2001 From: vlefevre Date: Thu, 9 Jun 2016 14:30:32 +0000 Subject: [src/mpfr-impl.h] Define mpfr_ueexp_t (needed for r10462). git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10463 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/mpfr-impl.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index cc035d346..a557a951a 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -820,11 +820,13 @@ union ieee_double_decimal64 { double d; _Decimal64 d64; }; #if _MPFR_EXP_FORMAT <= 3 typedef long int mpfr_eexp_t; +typedef unsigned long int mpfr_ueexp_t; # define mpfr_get_exp_t(x,r) mpfr_get_si((x),(r)) # define mpfr_set_exp_t(x,e,r) mpfr_set_si((x),(e),(r)) # define MPFR_EXP_FSPEC "l" #else typedef intmax_t mpfr_eexp_t; +typedef uintmax_t mpfr_ueexp_t; # define mpfr_get_exp_t(x,r) mpfr_get_sj((x),(r)) # define mpfr_set_exp_t(x,e,r) mpfr_set_sj((x),(e),(r)) # define MPFR_EXP_FSPEC "j" -- cgit v1.2.1 From ee220321ab2111cb425720920ff78b650b325646 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 10 Jun 2016 13:47:29 +0000 Subject: [src/sum.c] For the future support of reused arguments: TMD resolution is now done in a specific TMD accumulator, allocated at the same time as the main accumulator. This TMD accumulator currently just takes at most 2 limbs in practice, so that's not a problem. [doc/sum.txt] Added TODO for things that will change. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10464 280ebfd0-de03-0410-8827-d642c229c3f4 --- doc/sum.txt | 5 ++++ src/sum.c | 77 +++++++++++++++++++++++++++++++++++++++++-------------------- 2 files changed, 57 insertions(+), 25 deletions(-) diff --git a/doc/sum.txt b/doc/sum.txt index 24f1cb903..147602266 100644 --- a/doc/sum.txt +++ b/doc/sum.txt @@ -733,6 +733,8 @@ Accumulator: [--------]-----------------------------------] where cq = logn + 1, sq is the target precision, and dq ≥ logn + 2. +TODO: Update for code supporting reuse. + Memory is allocated both for the accumulator and for the temporary area needed by sum_raw. For performance reasons, the allocation is done in the stack if the size is small enough (see MPFR_TMP_LIMBS_ALLOC macro). @@ -937,6 +939,9 @@ We distinguish two cases: (−1, 0 or +1) of a secondary term with a second sum_raw invocation with a small-precision accumulator. + TODO: Update for code supporting reuse. Still explain what was done + before. + Since the uncorrected significand (which consists of the bits of exponents in ⟦max(u,minexp),e⟦ as mentioned above) has already been copied from the accumulator to the destination, we can reuse the diff --git a/src/sum.c b/src/sum.c index 280cfa63f..9c309288c 100644 --- a/src/sum.c +++ b/src/sum.c @@ -511,6 +511,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, mp_limb_t *wp; /* pointer to the accumulator */ mp_size_t ts; /* size of the temporary area, in limbs */ mp_size_t ws; /* size of the accumulator, in limbs */ + mp_size_t zs; /* size of the TMD accumulator, in limbs */ mpfr_prec_t wq; /* size of the accumulator, in bits */ int logn; /* ceil(log2(rn)) */ int cq; @@ -552,6 +553,9 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, wq = (mpfr_prec_t) ws * GMP_NUMB_BITS; MPFR_ASSERTD (wq - cq - sq >= 4); + /* TODO: timings, comparing with a larger zs. */ + zs = MPFR_PREC2LIMBS (wq - sq); + MPFR_LOG_MSG (("cq=%d sq=%Pd logn=%d wq=%Pd\n", cq, sq, logn, wq)); /* An input block will have up to wq - cq bits, and its shifted value @@ -560,7 +564,28 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPFR_TMP_MARK (marker); - tp = MPFR_TMP_LIMBS_ALLOC (ts + ws); + /* Note: If the TMD does not occur, which should be the case for most + sums, allocating zs limbs is not necessary. However, we choose to + do this now (thus in all cases) because zs is very small, so that + the difference on the memory footprint will not be noticeable. + More precisely, zs is at most 2 in practice with the current code; + we may want to increase it in order to avoid performance issues in + some unlikely corner cases, but even in this case, it will remain + small. + One will have: + [------ ts ------][------ ws ------][- zs -] + The following would probably be better: + [------ ts ------] [------ ws ------] + [- zs -] + i.e. where the TMD accumulator (partially or completely) takes + some unneeded part of the temporary area in order to improve + data locality. But + * in low precision, data locality is regarded as ensured even + with the actual choice; + * in high precision, data locality for TMD resolution may not + be that important. + */ + tp = MPFR_TMP_LIMBS_ALLOC (ts + ws + zs); wp = tp + ts; MPN_ZERO (wp, ws); /* zero the accumulator */ @@ -882,25 +907,27 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, else { mpfr_exp_t err; /* exponent of the error bound */ - mp_size_t zs; - int sst; /* sign of the secondary term */ + mp_size_t zz; /* number of limbs to zero in the TMD accumulator */ + mp_limb_t *zp; /* pointer to the TMD accumulator */ + mpfr_prec_t zq; /* size of the TMD accumulator, in bits */ + int sst; /* sign of the secondary term */ MPFR_ASSERTD (maxexp > MPFR_EXP_MIN); MPFR_ASSERTD (tmd == 1 || tmd == 2); - /* New accumulator size */ - ws = MPFR_PREC2LIMBS (wq - sq); - wq = (mpfr_prec_t) ws * GMP_NUMB_BITS; + /* TMD accumulator */ + zp = wp + ws; + zq = (mpfr_prec_t) zs * GMP_NUMB_BITS; err = maxexp + logn; MPFR_LOG_MSG (("TMD with" " maxexp=%" MPFR_EXP_FSPEC "d" " err=%" MPFR_EXP_FSPEC "d" - " ws=%Pd" - " wq=%Pd\n", + " zs=%Pd", + " zq=%Pd\n", (mpfr_eexp_t) maxexp, (mpfr_eexp_t) err, - (mpfr_prec_t) ws, wq)); + (mpfr_prec_t) zs, zq)); /* The d-1 bits from u-2 to u-d (= err) are identical. */ @@ -923,22 +950,22 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, { wi++; /* number of words with represented bits */ td = GMP_NUMB_BITS - td; - zs = ws - wi; - MPFR_ASSERTD (zs >= 0 && zs < ws); - mpn_lshift (wp + zs, wp, wi, td); + zz = zs - wi; + MPFR_ASSERTD (zz >= 0 && zz < zs); + mpn_lshift (zp + zz, wp, wi, td); } else { MPFR_ASSERTD (wi > 0); - zs = ws - wi; - MPFR_ASSERTD (zs >= 0 && zs < ws); - if (zs > 0) - MPN_COPY_DECR (wp + zs, wp, wi); + zz = zs - wi; + MPFR_ASSERTD (zz >= 0 && zz < zs); + if (zz > 0) + MPN_COPY_DECR (zp + zz, wp, wi); } /* Compute minexp = minexp - (zs * GMP_NUMB_BITS + td) safely. */ - UPDATE_MINEXP (minexp, zs * GMP_NUMB_BITS + td); - MPFR_ASSERTD (minexp == err + 2 - wq); + UPDATE_MINEXP (minexp, zz * GMP_NUMB_BITS + td); + MPFR_ASSERTD (minexp == err + 2 - zq); } else /* err < minexp */ { @@ -948,25 +975,25 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, from maxexp, with cq bits reserved to avoid an overflow (as in the early steps). */ MPFR_LOG_MSG (("[TMD] err < minexp\n", 0)); - zs = ws; + zz = zs; - /* Compute minexp = maxexp - (wq - cq) safely. */ - UPDATE_MINEXP (maxexp, wq - cq); - MPFR_ASSERTD (minexp == err + 1 - wq); + /* Compute minexp = maxexp - (zq - cq) safely. */ + UPDATE_MINEXP (maxexp, zq - cq); + MPFR_ASSERTD (minexp == err + 1 - zq); } - MPN_ZERO (wp, zs); + MPN_ZERO (zp, zz); /* We need to determine the sign sst of the secondary term. In sum_raw, since the truncated sum corresponding to this secondary term will be in [2^(e-1),2^e] and the error strictly less than 2^err, we can stop the iterations when e - err >= 1 (this bound is the 11th argument of sum_raw). */ - cancel = sum_raw (wp, ws, wq, x, n, minexp, maxexp, tp, ts, + cancel = sum_raw (zp, zs, zq, x, n, minexp, maxexp, tp, ts, logn, 1, NULL, &minexp, &maxexp); if (cancel != 0) - sst = MPFR_LIMB_MSB (wp[ws-1]) == 0 ? 1 : -1; + sst = MPFR_LIMB_MSB (zp[zs-1]) == 0 ? 1 : -1; else if (tmd == 1) sst = 0; else -- cgit v1.2.1 From 8db925017daa3bc13bc90af1dc10e736a6ba7e07 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 10 Jun 2016 13:58:34 +0000 Subject: [src/sum.c] Corrected a MPFR_LOG_MSG modified in the latest change. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10465 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/sum.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sum.c b/src/sum.c index 9c309288c..e95ffd3a4 100644 --- a/src/sum.c +++ b/src/sum.c @@ -924,7 +924,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPFR_LOG_MSG (("TMD with" " maxexp=%" MPFR_EXP_FSPEC "d" " err=%" MPFR_EXP_FSPEC "d" - " zs=%Pd", + " zs=%Pd" " zq=%Pd\n", (mpfr_eexp_t) maxexp, (mpfr_eexp_t) err, (mpfr_prec_t) zs, zq)); -- cgit v1.2.1 From 78c91733f56aab5f26b8843fbe867a7970118cbd Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 13 Jun 2016 06:40:25 +0000 Subject: [src/sum.c] Improved a comment. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10466 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/sum.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sum.c b/src/sum.c index e95ffd3a4..734e6c9ec 100644 --- a/src/sum.c +++ b/src/sum.c @@ -591,7 +591,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPN_ZERO (wp, ws); /* zero the accumulator */ { - mpfr_exp_t minexp; /* exponent of the LSB of the block */ + mpfr_exp_t minexp; /* exponent of the LSB of the block for sum_raw */ mpfr_prec_t cancel; /* number of cancelled bits */ mpfr_exp_t e; /* temporary exponent of the result */ mpfr_exp_t u; /* temporary exponent of the ulp (quantum) */ -- cgit v1.2.1 From ef2f6ef69c12a021fef23937dd6398a396c8ee33 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 13 Jun 2016 09:06:29 +0000 Subject: [src/sum.c] Copy the significand to the destination after resolving the TMD in order to support reused arguments. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10467 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/sum.c | 175 ++++++++++++++++++++++++++++++++++++-------------------------- 1 file changed, 101 insertions(+), 74 deletions(-) diff --git a/src/sum.c b/src/sum.c index 734e6c9ec..87ac8a13d 100644 --- a/src/sum.c +++ b/src/sum.c @@ -23,16 +23,15 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" -/* FIXME: In the case where one or several input pointers point to the - output variable, we need to store the significand in a new temporary - area (as usual), because these inputs may still need to be read for - the possible TMD resolution. Alternatively, since this is not - necessarily a rare case (doing s += sum(x[i],0<=i= MPFR_EXP_MIN (with an assertion failure if this is not the case); @@ -78,12 +77,12 @@ int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2] = { 0 }; easily testable due to these huge precisions. Moreover, switching to a 64-bit ABI would be a better solution for such computations. So, let's leave this unimplemented. */ -#define UPDATE_MINEXP(E,SH) \ +#define SAFE_SUB(V,E,SH) \ do \ { \ mpfr_prec_t sh = (SH); \ MPFR_ASSERTN ((E) >= MPFR_EXP_MIN + sh); \ - minexp = (E) - sh; \ + V = (E) - sh; \ } \ while (0) @@ -108,8 +107,8 @@ int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2] = { 0 }; * logn: ceil(log2(rn)), where rn is the number of regular inputs. * prec: lower bound for e - err (as described above). * ep: pointer to mpfr_exp_t (see below), or a null pointer. - * minexpp: pointer to mpfr_exp_t (see below). - * maxexpp: pointer to mpfr_exp_t (see below). + * minexpp: pointer to mpfr_exp_t (see below), or a null pointer. + * maxexpp: pointer to mpfr_exp_t (see below), or a null pointer. * * Preconditions: * prec >= 1 @@ -119,8 +118,9 @@ int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2] = { 0 }; * the exact sum for this sum_raw invocation is 0), otherwise the number * of cancelled bits (>= 1), defined as the number of identical bits on * the most significant part of the accumulator. In the latter case, it - * also returns the following data in variables passed by reference: - * - in ep: the exponent e of the computed result (unless ep is NULL); + * also returns the following data in variables passed by reference, if + * the pointers are not NULL: + * - in ep: the exponent e of the computed result; * - in minexpp: the last value of minexp; * - in maxexpp: the new value of maxexp (for the next iteration after * the first invocation of sum_raw in the main code). @@ -430,10 +430,20 @@ sum_raw (mp_limb_t *wp, mp_size_t ws, mpfr_prec_t wq, mpfr_ptr *const x, MPFR_LOG_MSG (("(err=%" MPFR_EXP_FSPEC "d) <= (e=%" MPFR_EXP_FSPEC "d) - (prec=%Pd)\n", (mpfr_eexp_t) err, (mpfr_eexp_t) e, prec)); + /* To avoid tests or copies, we consider the only two cases + that will occur in sum_aux. */ + MPFR_ASSERTD ((ep != NULL && + minexpp != NULL && + maxexpp != NULL) || + (ep == NULL && + minexpp == NULL && + maxexpp == NULL)); if (ep != NULL) - *ep = e; - *minexpp = minexp; - *maxexpp = maxexp2; + { + *ep = e; + *minexpp = minexp; + *maxexpp = maxexp2; + } MPFR_LOG_MSG (("return with minexp=%" MPFR_EXP_FSPEC "d maxexp2=%" MPFR_EXP_FSPEC "d%s\n", (mpfr_eexp_t) minexp, (mpfr_eexp_t) maxexp2, @@ -476,7 +486,7 @@ sum_raw (mp_limb_t *wp, mp_size_t ws, mpfr_prec_t wq, mpfr_ptr *const x, MPN_COPY_DECR (wp + shifts, wp, ws - shifts); MPN_ZERO (wp, shifts); /* Compute minexp = minexp - shiftq safely. */ - UPDATE_MINEXP (minexp, shiftq); + SAFE_SUB (minexp, minexp, shiftq); MPFR_ASSERTD (minexp < maxexp2); } } @@ -489,7 +499,7 @@ sum_raw (mp_limb_t *wp, mp_size_t ws, mpfr_prec_t wq, mpfr_ptr *const x, { MPFR_LOG_MSG (("accumulator = 0, reiterate\n", 0)); /* Compute minexp = maxexp2 - (wq - (logn + 1)) safely. */ - UPDATE_MINEXP (maxexp2, wq - (logn + 1)); + SAFE_SUB (minexp, maxexp2, wq - (logn + 1)); /* Note: the logn + 1 corresponds to cq in the main code. */ } } @@ -595,6 +605,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, mpfr_prec_t cancel; /* number of cancelled bits */ mpfr_exp_t e; /* temporary exponent of the result */ mpfr_exp_t u; /* temporary exponent of the ulp (quantum) */ + mp_limb_t lbit; /* last bit (useful if even rounding) */ mp_limb_t rbit; /* rounding bit (corrected in halfway case) */ int corr; /* correction term (from -1 to 2) */ int sd, sh; /* shift counts */ @@ -607,7 +618,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPFR_LOG_MSG (("Compute an approximation with sum_raw...\n", 0)); /* Compute minexp = maxexp - (wq - cq) safely. */ - UPDATE_MINEXP (maxexp, wq - cq); + SAFE_SUB (minexp, maxexp, wq - cq); MPFR_ASSERTD (wq >= logn + sq + 5); cancel = sum_raw (wp, ws, wq, x, n, minexp, maxexp, tp, ts, logn, sq + 3, &e, &minexp, &maxexp); @@ -660,7 +671,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, if (MPFR_LIKELY (u > minexp)) { mpfr_prec_t tq; - mp_size_t ei, fi, wi; + mp_size_t wi; int td; tq = u - minexp; @@ -669,24 +680,9 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, wi = tq / GMP_NUMB_BITS; - if (MPFR_LIKELY (sh != 0)) - { - ei = (e - minexp) / GMP_NUMB_BITS; - fi = ei - (sn - 1); - MPFR_ASSERTD (fi == wi || fi == wi + 1); - mpn_lshift (sump, wp + fi, sn, sh); - if (fi != wi) - sump[0] |= wp[wi] >> (GMP_NUMB_BITS - sh); - } - else - { - MPFR_ASSERTD ((mpfr_prec_t) (ws - (wi + sn)) * GMP_NUMB_BITS - == cancel); - MPN_COPY (sump, wp + wi, sn); - } - /* Determine the rounding bit, which is represented. */ td = tq % GMP_NUMB_BITS; + lbit = (wp[wi] >> td) & 1; rbit = td >= 1 ? ((wp[wi] >> (td - 1)) & MPFR_LIMB_ONE) : (MPFR_ASSERTD (wi >= 1), wp[wi-1] >> (GMP_NUMB_BITS - 1)); MPFR_ASSERTD (rbit == 0 || rbit == 1); @@ -699,8 +695,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, inex = rounding bit || sticky bit. In round to nearest, also determine the rounding direction: obtained from the rounding bit possibly except in halfway cases. */ - if (MPFR_LIKELY (rbit == 0 || - (rnd == MPFR_RNDN && ((wp[wi] >> td) & 1) == 0))) + if (MPFR_LIKELY (rbit == 0 || (rnd == MPFR_RNDN && lbit == 0))) { /* We need to determine the sticky bit, either to set inex (if the rounding bit is 0) or to possibly "correct" rbit @@ -716,10 +711,8 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, if (!inex) { - mp_size_t wj = wi; - - while (!inex && wj > 0) - inex = wp[--wj] != 0; + while (!inex && wi > 0) + inex = wp[--wi] != 0; if (!inex && rbit != 0) { /* sticky bit = 0, rounding bit = 1, @@ -835,32 +828,24 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, } else /* u <= minexp */ { - mp_size_t en; - - en = (e - minexp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS; - if (MPFR_LIKELY (sh != 0)) - mpn_lshift (sump + sn - en, wp, en, sh); - else if (MPFR_UNLIKELY (en > 0)) - MPN_COPY (sump + sn - en, wp, en); - if (sn > en) - MPN_ZERO (sump, sn - en); - - /* The exact value of the accumulator has been copied. + /* The exact value of the accumulator will be copied. * The TMD occurs if and only if there are bits still * not taken into account, and if it occurs, this is * necessarily on a machine number (-> tmd = 1). */ + lbit = u == minexp ? wp[0] & 1 : 0; rbit = 0; inex = tmd = maxexp != MPFR_EXP_MIN; } - /* Leading bit: 1 if positive, 0 if negative. */ - pos = sump[sn-1] >> (GMP_NUMB_BITS - 1); + /* pos = 1 if positive, 0 if negative. */ + /* TODO: change to neg? */ + pos = ! (wp[ws-1] >> (GMP_NUMB_BITS - 1)); MPFR_ASSERTD (rbit == 0 || rbit == 1); - MPFR_LOG_MSG (("tmd=%d rbit=%d inex=%d pos=%d\n", - tmd, (int) rbit, inex, pos)); + MPFR_LOG_MSG (("tmd=%d lbit=%d rbit=%d inex=%d pos=%d\n", + tmd, (int) lbit, (int) rbit, inex, pos)); /* Here, if the final sum is known to be exact, inex = 0, otherwise * inex = 1. We have a truncated significand, a trailing term t such @@ -906,12 +891,18 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, } else { + mpfr_exp_t minexp2; + mpfr_prec_t cancel2; mpfr_exp_t err; /* exponent of the error bound */ mp_size_t zz; /* number of limbs to zero in the TMD accumulator */ mp_limb_t *zp; /* pointer to the TMD accumulator */ mpfr_prec_t zq; /* size of the TMD accumulator, in bits */ int sst; /* sign of the secondary term */ + /* TMD case. Here we use a new variable minexp2, with the same + meaning as minexp, as we want to keep the minexp value for + the copy to the destination. */ + MPFR_ASSERTD (maxexp > MPFR_EXP_MIN); MPFR_ASSERTD (tmd == 1 || tmd == 2); @@ -963,9 +954,9 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPN_COPY_DECR (zp + zz, wp, wi); } - /* Compute minexp = minexp - (zs * GMP_NUMB_BITS + td) safely. */ - UPDATE_MINEXP (minexp, zz * GMP_NUMB_BITS + td); - MPFR_ASSERTD (minexp == err + 2 - zq); + /* Compute minexp2 = minexp - (zs * GMP_NUMB_BITS + td) safely. */ + SAFE_SUB (minexp2, minexp, zz * GMP_NUMB_BITS + td); + MPFR_ASSERTD (minexp2 == err + 2 - zq); } else /* err < minexp */ { @@ -977,9 +968,9 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPFR_LOG_MSG (("[TMD] err < minexp\n", 0)); zz = zs; - /* Compute minexp = maxexp - (zq - cq) safely. */ - UPDATE_MINEXP (maxexp, zq - cq); - MPFR_ASSERTD (minexp == err + 1 - zq); + /* Compute minexp2 = maxexp - (zq - cq) safely. */ + SAFE_SUB (minexp2, maxexp, zq - cq); + MPFR_ASSERTD (minexp2 == err + 1 - zq); } MPN_ZERO (zp, zz); @@ -989,10 +980,10 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, secondary term will be in [2^(e-1),2^e] and the error strictly less than 2^err, we can stop the iterations when e - err >= 1 (this bound is the 11th argument of sum_raw). */ - cancel = sum_raw (zp, zs, zq, x, n, minexp, maxexp, tp, ts, - logn, 1, NULL, &minexp, &maxexp); + cancel2 = sum_raw (zp, zs, zq, x, n, minexp2, maxexp, tp, ts, + logn, 1, NULL, NULL, NULL); - if (cancel != 0) + if (cancel2 != 0) sst = MPFR_LIMB_MSB (zp[zs-1]) == 0 ? 1 : -1; else if (tmd == 1) sst = 0; @@ -1003,7 +994,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, which depends on the last bit of the pre-rounded result. */ MPFR_ASSERTD (rnd == MPFR_RNDN && tmd == 2); - sst = (sump[0] & (MPFR_LIMB_ONE << sd)) ? 1 : -1; + sst = lbit != 0 ? 1 : -1; } MPFR_LOG_MSG (("[TMD] tmd=%d rbit=%d sst=%d\n", @@ -1011,7 +1002,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, /* Do not consider the corrected sst for MPFR_COV_SET */ MPFR_COV_SET (sum_tmd[(int) rnd][tmd-1][rbit] - [cancel == 0 ? 1 : sst+1][pos][sq > MPFR_PREC_MIN]); + [cancel2 == 0 ? 1 : sst+1][pos][sq > MPFR_PREC_MIN]); inex = MPFR_IS_LIKE_RNDD (rnd, pos ? 1 : -1) ? (sst ? -1 : 0) : @@ -1019,7 +1010,7 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, (MPFR_ASSERTD (rnd == MPFR_RNDN), tmd == 1 ? - sst : sst); - if (tmd == 2 && sst == (rbit ? -1 : 1)) + if (tmd == 2 && sst == (rbit != 0 ? -1 : 1)) corr = 1 - (int) rbit; else if (MPFR_IS_LIKE_RNDD (rnd, pos ? 1 : -1) && sst == -1) corr = (int) rbit - 1; @@ -1039,6 +1030,42 @@ sum_aux (mpfr_ptr sum, mpfr_ptr *const x, unsigned long n, mpfr_rnd_t rnd, MPFR_SIGN (sum) = pos ? MPFR_SIGN_POS : MPFR_SIGN_NEG; + if (MPFR_LIKELY (u > minexp)) + { + mp_size_t wi; + + /* Recompute the initial value of wi. */ + wi = (u - minexp) / GMP_NUMB_BITS; + if (MPFR_LIKELY (sh != 0)) + { + mp_size_t fi; + + fi = (e - minexp) / GMP_NUMB_BITS - (sn - 1); + MPFR_ASSERTD (fi == wi || fi == wi + 1); + mpn_lshift (sump, wp + fi, sn, sh); + if (fi != wi) + sump[0] |= wp[wi] >> (GMP_NUMB_BITS - sh); + } + else + { + MPFR_ASSERTD ((mpfr_prec_t) (ws - (wi + sn)) * GMP_NUMB_BITS + == cancel); + MPN_COPY (sump, wp + wi, sn); + } + } + else /* u <= minexp */ + { + mp_size_t en; + + en = (e - minexp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS; + if (MPFR_LIKELY (sh != 0)) + mpn_lshift (sump + sn - en, wp, en, sh); + else if (MPFR_UNLIKELY (en > 0)) + MPN_COPY (sump + sn - en, wp, en); + if (sn > en) + MPN_ZERO (sump, sn - en); + } + if (MPFR_UNLIKELY (sq == 1)) /* precision 1 */ { sump[0] = MPFR_LIMB_HIGHBIT; -- cgit v1.2.1 From c434ba967adf9a3e692d25e0caa2f4d4ca36fe66 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Wed, 15 Jun 2016 09:46:29 +0000 Subject: =?UTF-8?q?[doc/README.dev]=20Update=20(GCC=20trunk=20=E2=86=92=20?= =?UTF-8?q?GCC=205).?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10468 280ebfd0-de03-0410-8827-d642c229c3f4 --- doc/README.dev | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/README.dev b/doc/README.dev index 5f7279392..1ac6cf53c 100644 --- a/doc/README.dev +++ b/doc/README.dev @@ -266,7 +266,7 @@ To make a release (for the MPFR team): *after* the -fsanitize=undefined one. GCC 4.9 also supports "-fsanitize=undefined", but it just gives - diagnostic messages at runtime, not a failure; GCC trunk supports + diagnostic messages at runtime, not a failure; GCC 5 supports -fno-sanitize-recover like clang. Test with GCC's AddressSanitizer (-fsanitize=address). -- cgit v1.2.1 From 89cf3d8e5c6a390877df3744a7dc117fbf057de6 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 17 Jun 2016 15:17:34 +0000 Subject: [TODO] Suggest the use the keyword "static" in array indices of parameter declarations with C99 compilers (6.7.5.3p7) when the pointer is expected not to be null. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10476 280ebfd0-de03-0410-8827-d642c229c3f4 --- TODO | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/TODO b/TODO index 39fe6d47e..eda40784d 100644 --- a/TODO +++ b/TODO @@ -567,6 +567,18 @@ Table of contents: The exponent delta should be sufficient for now since it can be bounded by MPFR_PREC_MAX+1 if need be. +- use the keyword "static" in array indices of parameter declarations with + C99 compilers (6.7.5.3p7) when the pointer is expected not to be null? + For instance, if mpfr.h is changed to have: + __MPFR_DECLSPEC void mpfr_dump (const __mpfr_struct [static 1]); + and one calls + mpfr_dump (NULL); + one gets a warning with Clang. This is just an example; this needs to be + done in a clean way. + See: + http://stackoverflow.com/a/3430353/3782797 + https://hamberg.no/erlend/posts/2013-02-18-static-array-indices.html + ############################################################################## 7. Portability -- cgit v1.2.1 From 8c44705c9f06a8701666c53e476d536b6bdfea04 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Fri, 17 Jun 2016 15:38:51 +0000 Subject: [src/mpfr.h] Coding style: added spaces. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10477 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/mpfr.h | 85 +++++++++++++++++++++++++++++++------------------------------- 1 file changed, 43 insertions(+), 42 deletions(-) diff --git a/src/mpfr.h b/src/mpfr.h index 9c7a75f3f..ad14b5672 100644 --- a/src/mpfr.h +++ b/src/mpfr.h @@ -445,7 +445,7 @@ __MPFR_DECLSPEC int __MPFR_DECLSPEC int mpfr_set_si_2exp (mpfr_ptr, long, mpfr_exp_t, mpfr_rnd_t); __MPFR_DECLSPEC int - mpfr_set_ui_2exp (mpfr_ptr,unsigned long,mpfr_exp_t,mpfr_rnd_t); + mpfr_set_ui_2exp (mpfr_ptr, unsigned long, mpfr_exp_t, mpfr_rnd_t); #ifndef MPFR_USE_MINI_GMP /* mini-gmp does not provide mpq_t, we disable the following functions */ __MPFR_DECLSPEC int @@ -596,7 +596,7 @@ __MPFR_DECLSPEC int mpfr_div_d (mpfr_ptr, mpfr_srcptr, __MPFR_DECLSPEC int mpfr_d_div (mpfr_ptr, double, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_sqr (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_sqr (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_const_pi (mpfr_ptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_const_log2 (mpfr_ptr, mpfr_rnd_t); @@ -606,8 +606,8 @@ __MPFR_DECLSPEC int mpfr_const_catalan (mpfr_ptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_agm (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_log (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_log2 (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_log (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_log2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_log10 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_log1p (mpfr_ptr, mpfr_srcptr, @@ -615,14 +615,14 @@ __MPFR_DECLSPEC int mpfr_log1p (mpfr_ptr, mpfr_srcptr, __MPFR_DECLSPEC int mpfr_log_ui (mpfr_ptr, unsigned long, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_exp (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_exp2 (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_exp (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_exp2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_exp10 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_expm1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_eint (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_li2 (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_eint (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_li2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_cmp (mpfr_srcptr, mpfr_srcptr); __MPFR_DECLSPEC int mpfr_cmp3 (mpfr_srcptr, mpfr_srcptr, int); @@ -654,7 +654,7 @@ __MPFR_DECLSPEC int mpfr_mul_2si (mpfr_ptr, mpfr_srcptr, __MPFR_DECLSPEC int mpfr_div_2si (mpfr_ptr, mpfr_srcptr, long, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_rint (mpfr_ptr,mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_rint (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_roundeven (mpfr_ptr, mpfr_srcptr); __MPFR_DECLSPEC int mpfr_round (mpfr_ptr, mpfr_srcptr); __MPFR_DECLSPEC int mpfr_trunc (mpfr_ptr, mpfr_srcptr); @@ -670,7 +670,7 @@ __MPFR_DECLSPEC int mpfr_rint_ceil (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_rint_floor (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_frac (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_frac (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_modf (mpfr_ptr, mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_remquo (mpfr_ptr, long*, mpfr_srcptr, @@ -688,7 +688,7 @@ __MPFR_DECLSPEC int mpfr_fits_uint_p (mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_fits_sint_p (mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_fits_ushort_p (mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_fits_sshort_p (mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_fits_uintmax_p (mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_fits_uintmax_p (mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_fits_intmax_p (mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC void mpfr_extract (mpz_ptr, mpfr_srcptr, @@ -708,51 +708,52 @@ __MPFR_DECLSPEC int mpfr_greaterequal_p (mpfr_srcptr, mpfr_srcptr); __MPFR_DECLSPEC int mpfr_less_p (mpfr_srcptr, mpfr_srcptr); __MPFR_DECLSPEC int mpfr_lessequal_p (mpfr_srcptr, mpfr_srcptr); -__MPFR_DECLSPEC int mpfr_lessgreater_p (mpfr_srcptr,mpfr_srcptr); +__MPFR_DECLSPEC int mpfr_lessgreater_p (mpfr_srcptr, mpfr_srcptr); __MPFR_DECLSPEC int mpfr_equal_p (mpfr_srcptr, mpfr_srcptr); __MPFR_DECLSPEC int mpfr_unordered_p (mpfr_srcptr, mpfr_srcptr); -__MPFR_DECLSPEC int mpfr_atanh (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_acosh (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_asinh (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_cosh (mpfr_ptr,mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_sinh (mpfr_ptr,mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_tanh (mpfr_ptr,mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_atanh (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_acosh (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_asinh (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_cosh (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_sinh (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_tanh (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_sinh_cosh (mpfr_ptr, mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_sech (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_csch (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_coth (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_sech (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_csch (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_coth (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_acos (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_asin (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_atan (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_sin (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_acos (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_asin (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_atan (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_sin (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_sin_cos (mpfr_ptr, mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_cos (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_tan (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_atan2 (mpfr_ptr,mpfr_srcptr,mpfr_srcptr, +__MPFR_DECLSPEC int mpfr_cos (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_tan (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_atan2 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_sec (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_csc (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_cot (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_sec (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_csc (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_cot (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_hypot (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_erf (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_erfc (mpfr_ptr, mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_cbrt (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_root (mpfr_ptr,mpfr_srcptr,unsigned long,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_gamma (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_gamma_inc (mpfr_ptr,mpfr_srcptr, +__MPFR_DECLSPEC int mpfr_erf (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_erfc (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_cbrt (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_root (mpfr_ptr, mpfr_srcptr, unsigned long, + mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_gamma (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_gamma_inc (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_lngamma (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_lgamma (mpfr_ptr,int*,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_digamma (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_zeta (mpfr_ptr,mpfr_srcptr,mpfr_rnd_t); -__MPFR_DECLSPEC int mpfr_zeta_ui (mpfr_ptr,unsigned long,mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_lngamma (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_lgamma (mpfr_ptr, int *, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_digamma (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_zeta (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_zeta_ui (mpfr_ptr, unsigned long, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_fac_ui (mpfr_ptr, unsigned long int, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_j0 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); -- cgit v1.2.1