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"; }