diff -Naurd mpfr-2.3.0-a/Makefile.am mpfr-2.3.0-b/Makefile.am --- mpfr-2.3.0-a/Makefile.am 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/Makefile.am 2007-09-02 23:19:32.000000000 +0000 @@ -32,8 +32,13 @@ MAKEINFOFLAGS = --enable-encoding -get_patches.c: PATCHES get_patches.sh - ./get_patches.sh > $@ || rm -f $@ +# Important note: If for some reason, srcdir is read-only at build time +# (and you use objdir != srcdir), then you need to rebuild get_patches.c +# (with "make get_patches.c") just after patching the MPFR source. This +# should not be a problem in practice, in particular because "make dist" +# automatically rebuilds get_patches.c before generating the archives. +$(srcdir)/get_patches.c: PATCHES get_patches.sh + (cd $(srcdir) && ./get_patches.sh) > $@ || rm -f $@ # Do not add get_patches.c to CLEANFILES so that this file doesn't # need to be (re)built as long as no patches are applied. Anyway the diff -Naurd mpfr-2.3.0-a/Makefile.in mpfr-2.3.0-b/Makefile.in --- mpfr-2.3.0-a/Makefile.in 2007-08-29 10:27:18.000000000 +0000 +++ mpfr-2.3.0-b/Makefile.in 2007-09-02 23:59:30.000000000 +0000 @@ -1665,8 +1665,13 @@ uninstall-info-am uninstall-libLTLIBRARIES -get_patches.c: PATCHES get_patches.sh - ./get_patches.sh > $@ || rm -f $@ +# Important note: If for some reason, srcdir is read-only at build time +# (and you use objdir != srcdir), then you need to rebuild get_patches.c +# (with "make get_patches.c") just after patching the MPFR source. This +# should not be a problem in practice, in particular because "make dist" +# automatically rebuilds get_patches.c before generating the archives. +$(srcdir)/get_patches.c: PATCHES get_patches.sh + (cd $(srcdir) && ./get_patches.sh) > $@ || rm -f $@ tune: $(MAKE) $(AM_MAKEFLAGS) tuneup$(EXEEXT) diff -Naurd mpfr-2.3.0-a/PATCHES mpfr-2.3.0-b/PATCHES --- mpfr-2.3.0-a/PATCHES 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/PATCHES 2007-09-02 23:59:50.000000000 +0000 @@ -0,0 +1 @@ +get_patches diff -Naurd mpfr-2.3.0-a/VERSION mpfr-2.3.0-b/VERSION --- mpfr-2.3.0-a/VERSION 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/VERSION 2007-09-02 23:55:55.000000000 +0000 @@ -1 +1 @@ -2.3.0 +2.3.0-p1 diff -Naurd mpfr-2.3.0-a/mpfr.h mpfr-2.3.0-b/mpfr.h --- mpfr-2.3.0-a/mpfr.h 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/mpfr.h 2007-09-02 23:55:55.000000000 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 2 #define MPFR_VERSION_MINOR 3 #define MPFR_VERSION_PATCHLEVEL 0 -#define MPFR_VERSION_STRING "2.3.0" +#define MPFR_VERSION_STRING "2.3.0-p1" /* Macros dealing with MPFR VERSION */ #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) diff -Naurd mpfr-2.3.0-a/tests/tversion.c mpfr-2.3.0-b/tests/tversion.c --- mpfr-2.3.0-a/tests/tversion.c 2007-08-29 10:18:10.000000000 +0000 +++ mpfr-2.3.0-b/tests/tversion.c 2007-09-02 23:55:55.000000000 +0000 @@ -46,7 +46,7 @@ version = mpfr_get_version (); /* This test is disabled when a suffix (e.g. -dev) has been defined. */ -#if 1 +#if 0 sprintf (buffer, "%d.%d.%d", MPFR_VERSION_MAJOR, MPFR_VERSION_MINOR, MPFR_VERSION_PATCHLEVEL); if (strcmp (buffer, version) != 0) diff -Naurd mpfr-2.3.0-a/version.c mpfr-2.3.0-b/version.c --- mpfr-2.3.0-a/version.c 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/version.c 2007-09-02 23:55:55.000000000 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "2.3.0"; + return "2.3.0-p1"; } diff -Naurd mpfr-2.3.0-a/PATCHES mpfr-2.3.0-b/PATCHES --- mpfr-2.3.0-a/PATCHES 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/PATCHES 2007-09-03 00:03:46.000000000 +0000 @@ -0,0 +1 @@ +mpfr_acosh diff -Naurd mpfr-2.3.0-a/VERSION mpfr-2.3.0-b/VERSION --- mpfr-2.3.0-a/VERSION 2007-09-02 23:55:55.000000000 +0000 +++ mpfr-2.3.0-b/VERSION 2007-09-03 00:02:12.000000000 +0000 @@ -1 +1 @@ -2.3.0-p1 +2.3.0-p2 diff -Naurd mpfr-2.3.0-a/acosh.c mpfr-2.3.0-b/acosh.c --- mpfr-2.3.0-a/acosh.c 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/acosh.c 2007-08-31 17:20:08.000000000 +0000 @@ -73,7 +73,7 @@ /* Declaration of the size variables */ mp_prec_t Ny = MPFR_PREC(y); /* Precision of output variable */ mp_prec_t Nt; /* Precision of the intermediary variable */ - mp_exp_t err, exp_te, exp_ti; /* Precision of error */ + mp_exp_t err, exp_te, d; /* Precision of error */ MPFR_ZIV_DECL (loop); /* compute the precision of intermediary variable */ @@ -91,13 +91,35 @@ mpfr_mul (t, x, x, GMP_RNDD); /* x^2 */ exp_te = MPFR_GET_EXP (t); mpfr_sub_ui (t, t, 1, GMP_RNDD); /* x^2-1 */ - exp_ti = MPFR_GET_EXP (t); - mpfr_sqrt (t, t, GMP_RNDN); /* sqrt(x^2-1) */ - mpfr_add (t, t, x, GMP_RNDN); /* sqrt(x^2-1)+x */ - mpfr_log (t, t, GMP_RNDN); /* ln(sqrt(x^2-1)+x)*/ + if (MPFR_UNLIKELY (MPFR_IS_ZERO (t))) + { + mpfr_t z; + + /* This means that x is very close to 1: x = 1 + z with + z < 2^(-Nt). Instead of increasing the precision, let's + compute x^2-1 by (x+1)(x-1) with an accuracy of about + Nt bits. */ + mpfr_init2 (z, Nt); + mpfr_add_ui (t, x, 1, GMP_RNDD); + mpfr_sub_ui (z, x, 1, GMP_RNDD); + mpfr_mul (t, t, z, GMP_RNDD); + d = 2; + mpfr_sqrt (t, t, GMP_RNDN); /* sqrt(x^2-1) */ + mpfr_add (t, t, z, GMP_RNDN); /* sqrt(x^2-1)+z */ + mpfr_clear (z); + mpfr_log1p (t, t, GMP_RNDN); /* log1p(sqrt(x^2-1)+z) */ + } + else + { + d = exp_te - MPFR_GET_EXP (t); + d = MAX (1, d); + mpfr_sqrt (t, t, GMP_RNDN); /* sqrt(x^2-1) */ + mpfr_add (t, t, x, GMP_RNDN); /* sqrt(x^2-1)+x */ + mpfr_log (t, t, GMP_RNDN); /* ln(sqrt(x^2-1)+x) */ + } /* error estimate -- see algorithms.tex */ - err = 3 + MAX (1, exp_te - exp_ti) - MPFR_GET_EXP(t); + err = 3 + d - MPFR_GET_EXP (t); /* error is bounded by 1/2 + 2^err <= 2^(1+max(-1,err)) */ err = 1 + MAX (-1, err); if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Ny, rnd_mode))) @@ -117,9 +139,3 @@ MPFR_SAVE_EXPO_FREE (expo); return mpfr_check_range (y, inexact, rnd_mode); } - - - - - - diff -Naurd mpfr-2.3.0-a/mpfr.h mpfr-2.3.0-b/mpfr.h --- mpfr-2.3.0-a/mpfr.h 2007-09-02 23:55:55.000000000 +0000 +++ mpfr-2.3.0-b/mpfr.h 2007-09-03 00:02:12.000000000 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 2 #define MPFR_VERSION_MINOR 3 #define MPFR_VERSION_PATCHLEVEL 0 -#define MPFR_VERSION_STRING "2.3.0-p1" +#define MPFR_VERSION_STRING "2.3.0-p2" /* Macros dealing with MPFR VERSION */ #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) diff -Naurd mpfr-2.3.0-a/tests/tacosh.c mpfr-2.3.0-b/tests/tacosh.c --- mpfr-2.3.0-a/tests/tacosh.c 2007-08-29 10:18:10.000000000 +0000 +++ mpfr-2.3.0-b/tests/tacosh.c 2007-08-31 17:20:00.000000000 +0000 @@ -123,12 +123,40 @@ mpfr_clear (y); } +/* With MPFR 2.3.0, this yields an assertion failure in mpfr_acosh. */ +static void +bug20070831 (void) +{ + mpfr_t x, y, z; + int inex; + + mpfr_init2 (x, 256); + mpfr_init2 (y, 32); + mpfr_init2 (z, 32); + mpfr_set_ui (x, 1, GMP_RNDN); + mpfr_nextabove (x); + inex = mpfr_acosh (y, x, GMP_RNDZ); + mpfr_set_ui_2exp (z, 1, -127, GMP_RNDN); + mpfr_nextbelow (z); + MPFR_ASSERTN (inex < 0); + if (!mpfr_equal_p (y, z)) + { + printf ("Error in bug20070831:\nexpected "); + mpfr_dump (z); + printf ("got "); + mpfr_dump (y); + exit (1); + } + mpfr_clears (x, y, z, (void *) 0); +} + int main (int argc, char *argv[]) { tests_start_mpfr (); special (); + bug20070831 (); test_generic (2, 100, 25); diff -Naurd mpfr-2.3.0-a/version.c mpfr-2.3.0-b/version.c --- mpfr-2.3.0-a/version.c 2007-09-02 23:55:55.000000000 +0000 +++ mpfr-2.3.0-b/version.c 2007-09-03 00:02:12.000000000 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "2.3.0-p1"; + return "2.3.0-p2"; } diff -Naurd mpfr-2.3.0-a/PATCHES mpfr-2.3.0-b/PATCHES --- mpfr-2.3.0-a/PATCHES 2007-10-05 12:20:54.000000000 +0000 +++ mpfr-2.3.0-b/PATCHES 2007-10-05 12:21:12.000000000 +0000 @@ -0,0 +1 @@ +mpfr_atan2 diff -Naurd mpfr-2.3.0-a/VERSION mpfr-2.3.0-b/VERSION --- mpfr-2.3.0-a/VERSION 2007-09-03 00:02:12.000000000 +0000 +++ mpfr-2.3.0-b/VERSION 2007-10-05 12:21:06.000000000 +0000 @@ -1 +1 @@ -2.3.0-p2 +2.3.0-p3 diff -Naurd mpfr-2.3.0-a/atan2.c mpfr-2.3.0-b/atan2.c --- mpfr-2.3.0-a/atan2.c 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/atan2.c 2007-10-05 12:21:01.000000000 +0000 @@ -170,6 +170,7 @@ /* use atan2(y,x) = atan(y/x) */ for (;;) { + mpfr_clear_flags (); if (mpfr_div (tmp, y, x, GMP_RNDN) == 0) { /* Result is exact. */ diff -Naurd mpfr-2.3.0-a/mpfr.h mpfr-2.3.0-b/mpfr.h --- mpfr-2.3.0-a/mpfr.h 2007-09-03 00:02:12.000000000 +0000 +++ mpfr-2.3.0-b/mpfr.h 2007-10-05 12:21:06.000000000 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 2 #define MPFR_VERSION_MINOR 3 #define MPFR_VERSION_PATCHLEVEL 0 -#define MPFR_VERSION_STRING "2.3.0-p2" +#define MPFR_VERSION_STRING "2.3.0-p3" /* Macros dealing with MPFR VERSION */ #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) diff -Naurd mpfr-2.3.0-a/tests/tatan.c mpfr-2.3.0-b/tests/tatan.c --- mpfr-2.3.0-a/tests/tatan.c 2007-08-29 10:18:10.000000000 +0000 +++ mpfr-2.3.0-b/tests/tatan.c 2007-10-05 12:21:01.000000000 +0000 @@ -405,6 +405,40 @@ mpfr_clears (a, x, y, (void *) 0); } +/* Bug found by Robert Bajema (regression in MPFR 2.3.0). + The cause is the underflow flag set before the mpfr_atan2 call. */ +static void +atan2_bug_20071003 (void) +{ + mpfr_t a, x, y, z; + + mpfr_inits (a, x, y, z, (void *) 0); + + mpfr_set_underflow (); + mpfr_set_str_binary (y, + "-0.10100110110100110111010110111111100110100010001110110E2"); + mpfr_set_str_binary (x, + "0.10100101010110010100010010111000110110011110001011110E3"); + mpfr_set_str_binary (z, + "-0.11101111001101101100111011001101000010010111101110110E-1"); + mpfr_atan2 (a, y, x, GMP_RNDN); + if (! mpfr_equal_p (a, z)) + { + printf ("mpfr_atan2 fails on:\n"); + printf (" y = "); + mpfr_dump (y); + printf (" x = "); + mpfr_dump (x); + printf ("Expected "); + mpfr_dump (z); + printf ("Got "); + mpfr_dump (a); + exit (1); + } + + mpfr_clears (a, x, y, z, (void *) 0); +} + int main (int argc, char *argv[]) { @@ -414,6 +448,7 @@ special (); special_atan2 (); smallvals_atan2 (); + atan2_bug_20071003 (); test_generic_atan (2, 200, 17); test_generic_atan2 (2, 200, 17); diff -Naurd mpfr-2.3.0-a/version.c mpfr-2.3.0-b/version.c --- mpfr-2.3.0-a/version.c 2007-09-03 00:02:12.000000000 +0000 +++ mpfr-2.3.0-b/version.c 2007-10-05 12:21:06.000000000 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "2.3.0-p2"; + return "2.3.0-p3"; } diff -Naurd mpfr-2.3.0-a/PATCHES mpfr-2.3.0-b/PATCHES --- mpfr-2.3.0-a/PATCHES 2007-10-23 03:03:12.000000000 +0000 +++ mpfr-2.3.0-b/PATCHES 2007-10-23 03:06:18.000000000 +0000 @@ -0,0 +1 @@ +mpfr_subnormalize diff -Naurd mpfr-2.3.0-a/VERSION mpfr-2.3.0-b/VERSION --- mpfr-2.3.0-a/VERSION 2007-10-05 12:21:06.000000000 +0000 +++ mpfr-2.3.0-b/VERSION 2007-10-23 03:05:51.000000000 +0000 @@ -1 +1 @@ -2.3.0-p3 +2.3.0-p4 diff -Naurd mpfr-2.3.0-a/mpfr.h mpfr-2.3.0-b/mpfr.h --- mpfr-2.3.0-a/mpfr.h 2007-10-05 12:21:06.000000000 +0000 +++ mpfr-2.3.0-b/mpfr.h 2007-10-23 03:05:51.000000000 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 2 #define MPFR_VERSION_MINOR 3 #define MPFR_VERSION_PATCHLEVEL 0 -#define MPFR_VERSION_STRING "2.3.0-p3" +#define MPFR_VERSION_STRING "2.3.0-p4" /* Macros dealing with MPFR VERSION */ #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) diff -Naurd mpfr-2.3.0-a/mpfr.texi mpfr-2.3.0-b/mpfr.texi --- mpfr-2.3.0-a/mpfr.texi 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/mpfr.texi 2007-10-23 03:05:25.000000000 +0000 @@ -2038,7 +2038,8 @@ @var{rnd} and @var{t} must be the used rounding mode for computing @var{x} and the returned ternary value when computing @var{x}. The subnormal exponent range is from @code{emin} to @code{emin+PREC(x)-1}. -This functions assumes that @code{emax-emin >= PREC(x)}. +If the result cannot be represented in the current exponent range +(due to a too small @code{emax}), the behavior is undefined. Note that unlike most functions, the result is compared to the exact one, not the input value @var{x}, i.e.@: the ternary value is propagated. This is a preliminary interface. diff -Naurd mpfr-2.3.0-a/subnormal.c mpfr-2.3.0-b/subnormal.c --- mpfr-2.3.0-a/subnormal.c 2007-08-29 10:18:11.000000000 +0000 +++ mpfr-2.3.0-b/subnormal.c 2007-10-23 03:05:25.000000000 +0000 @@ -26,10 +26,10 @@ /* For GMP_RNDN, we can have a problem of double rounding. In such a case, this table helps to conclude what to do (y positive): Rounding Bit | Sticky Bit | inexact | Action | new inexact - 0 | ? | ? | Trunc | sticky - 1 | 0 | -1 | Trunc | + 0 | ? | ? | Trunc | sticky + 1 | 0 | 1 | Trunc | 1 | 0 | 0 | Trunc if even | - 1 | 0 | 1 | AddOneUlp | + 1 | 0 | -1 | AddOneUlp | 1 | 1 | ? | AddOneUlp | For other rounding mode, there isn't such a problem. @@ -41,8 +41,6 @@ { int inexact = 0; - MPFR_ASSERTD ((mpfr_uexp_t) __gmpfr_emax - __gmpfr_emin >= MPFR_PREC (y)); - /* The subnormal exponent range are from: mpfr_emin to mpfr_emin + MPFR_PREC(y) - 1 */ if (MPFR_LIKELY (MPFR_IS_SINGULAR (y) @@ -69,23 +67,24 @@ and use the previous table to conclude. */ s = MPFR_LIMB_SIZE (y) - 1; mant = MPFR_MANT (y) + s; - rb = *mant & (MPFR_LIMB_HIGHBIT>>1); + rb = *mant & (MPFR_LIMB_HIGHBIT >> 1); if (rb == 0) goto set_min; - sb = *mant & ((MPFR_LIMB_HIGHBIT>>1)-1); + sb = *mant & ((MPFR_LIMB_HIGHBIT >> 1) - 1); while (sb == 0 && s-- != 0) sb = *--mant; if (sb != 0) goto set_min_p1; /* Rounding bit is 1 and sticky bit is 0. We need to examine old inexact flag to conclude. */ - if (old_inexact * MPFR_SIGN (y) < 0) + if ((old_inexact > 0 && MPFR_SIGN (y) > 0) || + (old_inexact < 0 && MPFR_SIGN (y) < 0)) goto set_min; - /* If inexact != 0, return 0.1*2^emin+1. + /* If inexact != 0, return 0.1*2^(emin+1). Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0 - So we have 0.1100000000000000000000000*2^emin exactly!!! - we choose to return 0.1*2^emin+1 which minimizes the relative - error. */ + So we have 0.1100000000000000000000000*2^emin exactly. + We return 0.1*2^(emin+1) according to the even-rounding + rule on subnormals. */ goto set_min_p1; } else if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y))) @@ -97,7 +96,8 @@ else { set_min_p1: - mpfr_setmin (y, __gmpfr_emin+1); + /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */ + mpfr_setmin (y, __gmpfr_emin + 1); inexact = MPFR_SIGN (y); } } @@ -120,15 +120,17 @@ MPFR_SET_EXP (dest, MPFR_GET_EXP (dest)+1)); if (MPFR_LIKELY (old_inexact != 0)) { - if (MPFR_UNLIKELY(rnd==GMP_RNDN && (inexact == MPFR_EVEN_INEX - || inexact == -MPFR_EVEN_INEX))) + if (MPFR_UNLIKELY(rnd == GMP_RNDN && (inexact == MPFR_EVEN_INEX + || inexact == -MPFR_EVEN_INEX))) { - if (old_inexact*inexact*MPFR_INT_SIGN (y) > 0) + /* if both roundings are in the same direction, we have to go + back in the other direction */ + if (SAME_SIGN (inexact, old_inexact)) { - if (inexact < 0) - mpfr_nexttoinf (dest); - else + if (SAME_SIGN (inexact, MPFR_INT_SIGN (y))) mpfr_nexttozero (dest); + else + mpfr_nexttoinf (dest); inexact = -inexact; } } @@ -136,7 +138,8 @@ inexact = old_inexact; } old_inexact = mpfr_set (y, dest, rnd); - MPFR_ASSERTD (old_inexact == 0); + MPFR_ASSERTN (old_inexact == 0); + MPFR_ASSERTN (MPFR_IS_PURE_FP (y)); mpfr_clear (dest); } return inexact; diff -Naurd mpfr-2.3.0-a/tests/tsubnormal.c mpfr-2.3.0-b/tests/tsubnormal.c --- mpfr-2.3.0-a/tests/tsubnormal.c 2007-08-29 10:18:09.000000000 +0000 +++ mpfr-2.3.0-b/tests/tsubnormal.c 2007-10-23 03:05:25.000000000 +0000 @@ -22,6 +22,7 @@ #include #include +#include #include "mpfr-test.h" @@ -41,6 +42,8 @@ {"0.11001E-10", 0, GMP_RNDZ, "0.1E-10"}, {"0.11001E-10", 0, GMP_RNDU, "0.1E-9"}, {"0.11000E-10", 0, GMP_RNDN, "0.1E-9"}, + {"0.11000E-10", -1, GMP_RNDN, "0.1E-9"}, + {"0.11000E-10", 1, GMP_RNDN, "0.1E-10"}, {"0.11111E-8", 0, GMP_RNDN, "0.10E-7"}, {"0.10111E-8", 0, GMP_RNDN, "0.11E-8"}, {"0.11110E-8", -1, GMP_RNDN, "0.10E-7"}, @@ -51,7 +54,7 @@ check1 (void) { mpfr_t x; - int i, j; + int i, j, k, s, old_inex; mpfr_set_default_prec (9); mpfr_set_emin (-10); @@ -59,20 +62,115 @@ mpfr_init (x); for (i = 0; i < (sizeof (tab) / sizeof (tab[0])); i++) - { - mpfr_set_str (x, tab[i].in, 2, GMP_RNDN); - j = mpfr_subnormalize (x, tab[i].i, tab[i].rnd); - if (mpfr_cmp_str (x, tab[i].out, 2, GMP_RNDN) != 0) + for (s = 0; s <= (tab[i].rnd == GMP_RNDN); s++) + for (k = 0; k <= 1; k++) { - printf ("Error for i=%d\nFor:%s\nExpected:%s\nGot:", - i, tab[i].in, tab[i].out); - mpfr_dump (x); - exit (1); + mpfr_set_str (x, tab[i].in, 2, GMP_RNDN); + old_inex = tab[i].i; + if (s) + { + mpfr_neg (x, x, GMP_RNDN); + old_inex = - old_inex; + } + if (k && old_inex) + old_inex = old_inex < 0 ? INT_MIN : INT_MAX; + j = mpfr_subnormalize (x, old_inex, tab[i].rnd); + /* TODO: test j. */ + if (s) + mpfr_neg (x, x, GMP_RNDN); + if (mpfr_cmp_str (x, tab[i].out, 2, GMP_RNDN) != 0) + { + char *sgn = s ? "-" : ""; + printf ("Error for i = %d (old_inex = %d), k = %d, x = %s%s\n" + "Expected: %s%s\nGot: ", i, old_inex, k, + sgn, tab[i].in, sgn, tab[i].out); + if (s) + mpfr_neg (x, x, GMP_RNDN); + mpfr_dump (x); + exit (1); + } } - } mpfr_clear (x); } +/* bug found by Kevin P. Rauch on 22 Oct 2007 */ +static void +check2 (void) +{ + mpfr_t x, y, z; + int tern; + + mpfr_init2 (x, 32); + mpfr_init2 (y, 32); + mpfr_init2 (z, 32); + + mpfr_set_ui (x, 0xC0000000U, GMP_RNDN); + mpfr_neg (x, x, GMP_RNDN); + mpfr_set_ui (y, 0xFFFFFFFEU, GMP_RNDN); + mpfr_set_exp (x, 0); + mpfr_set_exp (y, 0); + mpfr_set_emin (-29); + + tern = mpfr_mul (z, x, y, GMP_RNDN); + /* z = -0.BFFFFFFE, tern > 0 */ + + tern = mpfr_subnormalize (z, tern, GMP_RNDN); + /* z should be -0.75 */ + MPFR_ASSERTN (tern < 0 && mpfr_cmp_si_2exp (z, -3, -2) == 0); + + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); +} + +/* bug found by Kevin P. Rauch on 22 Oct 2007 */ +static void +check3 (void) +{ + mpfr_t x, y, z; + int tern; + + mpfr_init2 (x, 32); + mpfr_init2 (y, 32); + mpfr_init2 (z, 32); + + mpfr_set_ui (x, 0xBFFFFFFFU, GMP_RNDN); /* 3221225471/2^32 */ + mpfr_set_ui (y, 0x80000001U, GMP_RNDN); /* 2147483649/2^32 */ + mpfr_set_exp (x, 0); + mpfr_set_exp (y, 0); + mpfr_set_emin (-1); + + /* the exact product is 6917529028714823679/2^64, which is rounded to + 3/8 = 0.375, which is smaller, thus tern < 0 */ + tern = mpfr_mul (z, x, y, GMP_RNDN); + MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0); + + tern = mpfr_subnormalize (z, tern, GMP_RNDN); + /* since emin = -1, and EXP(z)=-1, z should be rounded to precision + EXP(z)-emin+1 = 1, i.e., z should be a multiple of the smallest possible + positive representable value with emin=-1, which is 1/4. The two + possible values are 1/4 and 2/4, which are at equal distance of z. + But since tern < 0, we should choose the largest value, i.e., 2/4. */ + MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0); + + /* here is another test for the alternate case, where z was rounded up + first, thus we have to round down */ + mpfr_set_str_binary (x, "0.11111111111010110101011011011011"); + mpfr_set_str_binary (y, "0.01100000000001111100000000001110"); + tern = mpfr_mul (z, x, y, GMP_RNDN); + MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0); + tern = mpfr_subnormalize (z, tern, GMP_RNDN); + MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 1, -2) == 0); + + /* finally the case where z was exact, which we simulate here */ + mpfr_set_ui_2exp (z, 3, -3, GMP_RNDN); + tern = mpfr_subnormalize (z, 0, GMP_RNDN); + MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0); + + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); +} int main (int argc, char *argv[]) @@ -80,6 +178,8 @@ tests_start_mpfr (); check1 (); + check2 (); + check3 (); tests_end_mpfr (); return 0; diff -Naurd mpfr-2.3.0-a/version.c mpfr-2.3.0-b/version.c --- mpfr-2.3.0-a/version.c 2007-10-05 12:21:06.000000000 +0000 +++ mpfr-2.3.0-b/version.c 2007-10-23 03:05:51.000000000 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "2.3.0-p3"; + return "2.3.0-p4"; }