diff options
Diffstat (limited to 'patches/mpfr/3.0.0/160-mpfr_sub1.patch')
| -rw-r--r-- | patches/mpfr/3.0.0/160-mpfr_sub1.patch | 628 |
1 files changed, 628 insertions, 0 deletions
diff --git a/patches/mpfr/3.0.0/160-mpfr_sub1.patch b/patches/mpfr/3.0.0/160-mpfr_sub1.patch new file mode 100644 index 00000000..2ce98e9f --- /dev/null +++ b/patches/mpfr/3.0.0/160-mpfr_sub1.patch @@ -0,0 +1,628 @@ +diff -Naurd mpfr-3.0.0-a/PATCHES mpfr-3.0.0-b/PATCHES +--- mpfr-3.0.0-a/PATCHES 2010-10-21 20:59:32.000000000 +0000 ++++ mpfr-3.0.0-b/PATCHES 2010-10-21 20:59:32.000000000 +0000 +@@ -0,0 +1 @@ ++mpfr_sub1 +diff -Naurd mpfr-3.0.0-a/VERSION mpfr-3.0.0-b/VERSION +--- mpfr-3.0.0-a/VERSION 2010-10-21 20:28:38.000000000 +0000 ++++ mpfr-3.0.0-b/VERSION 2010-10-21 20:59:32.000000000 +0000 +@@ -1 +1 @@ +-3.0.0-p5 ++3.0.0-p6 +diff -Naurd mpfr-3.0.0-a/mpfr.h mpfr-3.0.0-b/mpfr.h +--- mpfr-3.0.0-a/mpfr.h 2010-10-21 20:28:38.000000000 +0000 ++++ mpfr-3.0.0-b/mpfr.h 2010-10-21 20:59:32.000000000 +0000 +@@ -27,7 +27,7 @@ + #define MPFR_VERSION_MAJOR 3 + #define MPFR_VERSION_MINOR 0 + #define MPFR_VERSION_PATCHLEVEL 0 +-#define MPFR_VERSION_STRING "3.0.0-p5" ++#define MPFR_VERSION_STRING "3.0.0-p6" + + /* Macros dealing with MPFR VERSION */ + #define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c)) +diff -Naurd mpfr-3.0.0-a/sub1.c mpfr-3.0.0-b/sub1.c +--- mpfr-3.0.0-a/sub1.c 2010-06-10 11:00:14.000000000 +0000 ++++ mpfr-3.0.0-b/sub1.c 2010-10-21 20:59:32.000000000 +0000 +@@ -37,7 +37,9 @@ + mp_size_t cancel2, an, bn, cn, cn0; + mp_limb_t *ap, *bp, *cp; + mp_limb_t carry, bb, cc, borrow = 0; +- int inexact, shift_b, shift_c, is_exact = 1, down = 0, add_exp = 0; ++ int inexact, shift_b, shift_c, add_exp = 0; ++ int cmp_low = 0; /* used for rounding to nearest: 0 if low(b) = low(c), ++ negative if low(b) < low(c), positive if low(b)>low(c) */ + int sh, k; + MPFR_TMP_DECL(marker); + +@@ -196,7 +198,8 @@ + } + + #ifdef DEBUG +- printf ("shift_b=%d shift_c=%d diffexp=%lu\n", shift_b, shift_c, ++ printf ("rnd=%s shift_b=%d shift_c=%d diffexp=%lu\n", ++ mpfr_print_rnd_mode (rnd_mode), shift_b, shift_c, + (unsigned long) diff_exp); + #endif + +@@ -307,17 +310,18 @@ + { + if (MPFR_LIKELY(sh)) + { +- is_exact = (carry == 0); + /* can decide except when carry = 2^(sh-1) [middle] + or carry = 0 [truncate, but cannot decide inexact flag] */ +- down = (carry < (MPFR_LIMB_ONE << (sh - 1))); + if (carry > (MPFR_LIMB_ONE << (sh - 1))) + goto add_one_ulp; +- else if ((0 < carry) && down) ++ else if ((0 < carry) && (carry < (MPFR_LIMB_ONE << (sh - 1)))) + { + inexact = -1; /* result if smaller than exact value */ + goto truncate; + } ++ /* now carry = 2^(sh-1), in which case cmp_low=2, ++ or carry = 0, in which case cmp_low=0 */ ++ cmp_low = (carry == 0) ? 0 : 2; + } + } + else /* directed rounding: set rnd_mode to RNDZ iff toward zero */ +@@ -344,12 +348,32 @@ + cn -= (long int) an + cancel2; + + #ifdef DEBUG +- printf ("last %d bits from a are %lu, bn=%ld, cn=%ld\n", ++ printf ("last sh=%d bits from a are %lu, bn=%ld, cn=%ld\n", + sh, (unsigned long) carry, (long) bn, (long) cn); + #endif + ++ /* for rounding to nearest, we couldn't conclude up to here in the following ++ cases: ++ 1. sh = 0, then cmp_low=0: we can either truncate, subtract one ulp ++ or add one ulp: -1 ulp < low(b)-low(c) < 1 ulp ++ 2. sh > 0 but the low sh bits from high(b)-high(c) equal 2^(sh-1): ++ -0.5 ulp <= -1/2^sh < low(b)-low(c)-0.5 < 1/2^sh <= 0.5 ulp ++ we can't decide the rounding, in that case cmp_low=2: ++ either we truncate and flag=-1, or we add one ulp and flag=1 ++ 3. the low sh>0 bits from high(b)-high(c) equal 0: we know we have to ++ truncate but we can't decide the ternary value, here cmp_low=0: ++ -0.5 ulp <= -1/2^sh < low(b)-low(c) < 1/2^sh <= 0.5 ulp ++ we always truncate and inexact can be any of -1,0,1 ++ */ ++ ++ /* note: here cn might exceed cn0, in which case we consider a zero limb */ + for (k = 0; (bn > 0) || (cn > 0); k = 1) + { ++ /* if cmp_low < 0, we know low(b) - low(c) < 0 ++ if cmp_low > 0, we know low(b) - low(c) > 0 ++ (more precisely if cmp_low = 2, low(b) - low(c) = 0.5 ulp so far) ++ if cmp_low = 0, so far low(b) - low(c) = 0 */ ++ + /* get next limbs */ + bb = (bn > 0) ? bp[--bn] : 0; + if ((cn > 0) && (cn-- <= cn0)) +@@ -357,76 +381,115 @@ + else + cc = 0; + +- /* down is set when low(b) < low(c) */ +- if (down == 0) +- down = (bb < cc); ++ /* cmp_low compares low(b) and low(c) */ ++ if (cmp_low == 0) /* case 1 or 3 */ ++ cmp_low = (bb < cc) ? -2+k : (bb > cc) ? 1 : 0; ++ ++ /* Case 1 for k=0 splits into 7 subcases: ++ 1a: bb > cc + half ++ 1b: bb = cc + half ++ 1c: 0 < bb - cc < half ++ 1d: bb = cc ++ 1e: -half < bb - cc < 0 ++ 1f: bb - cc = -half ++ 1g: bb - cc < -half ++ ++ Case 2 splits into 3 subcases: ++ 2a: bb > cc ++ 2b: bb = cc ++ 2c: bb < cc ++ ++ Case 3 splits into 3 subcases: ++ 3a: bb > cc ++ 3b: bb = cc ++ 3c: bb < cc ++ */ + + /* the case rounding to nearest with sh=0 is special since one couldn't + subtract above 1/2 ulp in the trailing limb of the result */ +- if ((rnd_mode == MPFR_RNDN) && sh == 0 && k == 0) ++ if (rnd_mode == MPFR_RNDN && sh == 0 && k == 0) /* case 1 for k=0 */ + { + mp_limb_t half = MPFR_LIMB_HIGHBIT; + +- is_exact = (bb == cc); +- + /* add one ulp if bb > cc + half + truncate if cc - half < bb < cc + half + sub one ulp if bb < cc - half + */ + +- if (down) ++ if (cmp_low < 0) /* bb < cc: -1 ulp < low(b) - low(c) < 0, ++ cases 1e, 1f and 1g */ + { + if (cc >= half) + cc -= half; +- else ++ else /* since bb < cc < half, bb+half < 2*half */ + bb += half; ++ /* now we have bb < cc + half: ++ we have to subtract one ulp if bb < cc, ++ and truncate if bb > cc */ + } +- else /* bb >= cc */ ++ else if (cmp_low >= 0) /* bb >= cc, cases 1a to 1d */ + { + if (cc < half) + cc += half; +- else ++ else /* since bb >= cc >= half, bb - half >= 0 */ + bb -= half; ++ /* now we have bb > cc - half: we have to add one ulp if bb > cc, ++ and truncate if bb < cc */ ++ if (cmp_low > 0) ++ cmp_low = 2; + } + } + + #ifdef DEBUG +- printf (" bb=%lu cc=%lu down=%d is_exact=%d\n", +- (unsigned long) bb, (unsigned long) cc, down, is_exact); ++ printf ("k=%u bb=%lu cc=%lu cmp_low=%d\n", k, ++ (unsigned long) bb, (unsigned long) cc, cmp_low); + #endif +- if (bb < cc) ++ if (cmp_low < 0) /* low(b) - low(c) < 0: either truncate or subtract ++ one ulp */ + { + if (rnd_mode == MPFR_RNDZ) +- goto sub_one_ulp; ++ goto sub_one_ulp; /* set inexact=-1 */ + else if (rnd_mode != MPFR_RNDN) /* round away */ + { + inexact = 1; + goto truncate; + } +- else /* round to nearest: special case here since for sh=k=0 +- bb = bb0 - MPFR_LIMB_HIGHBIT */ ++ else /* round to nearest */ + { +- if (is_exact && sh == 0) +- { +- /* For k=0 we can't decide exactness since it may depend +- from low order bits. +- For k=1, the first low limbs matched: low(b)-low(c)<0. */ +- if (k) +- { +- inexact = 1; +- goto truncate; +- } +- } +- else if (down && sh == 0) +- goto sub_one_ulp; +- else +- { +- inexact = (is_exact) ? 1 : -1; ++ /* If cmp_low < 0 and bb > cc, then -0.5 ulp < low(b)-low(c) < 0, ++ whatever the value of sh. ++ If sh>0, then cmp_low < 0 implies that the initial neglected ++ sh bits were 0 (otherwise cmp_low=2 initially), thus the ++ weight of the new bits is less than 0.5 ulp too. ++ If k > 0 (and sh=0) this means that either the first neglected ++ limbs bb and cc were equal (thus cmp_low was 0 for k=0), ++ or we had bb - cc = -0.5 ulp or 0.5 ulp. ++ The last case is not possible here since we would have ++ cmp_low > 0 which is sticky. ++ In the first case (where we have cmp_low = -1), we truncate, ++ whereas in the 2nd case we have cmp_low = -2 and we subtract ++ one ulp. ++ */ ++ if (bb > cc || sh > 0 || cmp_low == -1) ++ { /* -0.5 ulp < low(b)-low(c) < 0, ++ bb > cc corresponds to cases 1e and 1f1 ++ sh > 0 corresponds to cases 3c and 3b3 ++ cmp_low = -1 corresponds to case 1d3 (also 3b3) */ ++ inexact = 1; + goto truncate; + } ++ else if (bb < cc) /* here sh = 0 and low(b)-low(c) < -0.5 ulp, ++ this corresponds to cases 1g and 1f3 */ ++ goto sub_one_ulp; ++ /* the only case where we can't conclude is sh=0 and bb=cc, ++ i.e., we have low(b) - low(c) = -0.5 ulp (up to now), thus ++ we don't know if we must truncate or subtract one ulp. ++ Note: for sh=0 we can't have low(b) - low(c) = -0.5 ulp up to ++ now, since low(b) - low(c) > 1/2^sh */ + } + } +- else if (bb > cc) ++ else if (cmp_low > 0) /* 0 < low(b) - low(c): either truncate or ++ add one ulp */ + { + if (rnd_mode == MPFR_RNDZ) + { +@@ -437,34 +500,70 @@ + goto add_one_ulp; + else /* round to nearest */ + { +- if (is_exact) ++ if (bb > cc) + { +- inexact = -1; +- goto truncate; ++ /* if sh=0, then bb>cc means that low(b)-low(c) > 0.5 ulp, ++ and similarly when cmp_low=2 */ ++ if (cmp_low == 2) /* cases 1a, 1b1, 2a and 2b1 */ ++ goto add_one_ulp; ++ /* sh > 0 and cmp_low > 0: this implies that the sh initial ++ neglected bits were 0, and the remaining low(b)-low(c)>0, ++ but its weight is less than 0.5 ulp */ ++ else /* 0 < low(b) - low(c) < 0.5 ulp, this corresponds to ++ cases 3a, 1d1 and 3b1 */ ++ { ++ inexact = -1; ++ goto truncate; ++ } + } +- else if (down) ++ else if (bb < cc) /* 0 < low(b) - low(c) < 0.5 ulp, cases 1c, ++ 1b3, 2b3 and 2c */ + { +- inexact = 1; ++ inexact = -1; + goto truncate; + } +- else +- goto add_one_ulp; ++ /* the only case where we can't conclude is bb=cc, i.e., ++ low(b) - low(c) = 0.5 ulp (up to now), thus we don't know ++ if we must truncate or add one ulp. */ + } + } ++ /* after k=0, we cannot conclude in the following cases, we split them ++ according to the values of bb and cc for k=1: ++ 1b. sh=0 and cmp_low = 1 and bb-cc = half [around 0.5 ulp] ++ 1b1. bb > cc: add one ulp, inex = 1 ++ 1b2: bb = cc: cannot conclude ++ 1b3: bb < cc: truncate, inex = -1 ++ 1d. sh=0 and cmp_low = 0 and bb-cc = 0 [around 0] ++ 1d1: bb > cc: truncate, inex = -1 ++ 1d2: bb = cc: cannot conclude ++ 1d3: bb < cc: truncate, inex = +1 ++ 1f. sh=0 and cmp_low = -1 and bb-cc = -half [around -0.5 ulp] ++ 1f1: bb > cc: truncate, inex = +1 ++ 1f2: bb = cc: cannot conclude ++ 1f3: bb < cc: sub one ulp, inex = -1 ++ 2b. sh > 0 and cmp_low = 2 and bb=cc [around 0.5 ulp] ++ 2b1. bb > cc: add one ulp, inex = 1 ++ 2b2: bb = cc: cannot conclude ++ 2b3: bb < cc: truncate, inex = -1 ++ 3b. sh > 0 and cmp_low = 0 [around 0] ++ 3b1. bb > cc: truncate, inex = -1 ++ 3b2: bb = cc: cannot conclude ++ 3b3: bb < cc: truncate, inex = +1 ++ */ + } + +- if ((rnd_mode == MPFR_RNDN) && !is_exact) ++ if ((rnd_mode == MPFR_RNDN) && cmp_low != 0) + { + /* even rounding rule */ + if ((ap[0] >> sh) & 1) + { +- if (down) ++ if (cmp_low < 0) + goto sub_one_ulp; + else + goto add_one_ulp; + } + else +- inexact = (down) ? 1 : -1; ++ inexact = (cmp_low > 0) ? -1 : 1; + } + else + inexact = 0; +diff -Naurd mpfr-3.0.0-a/tests/tfma.c mpfr-3.0.0-b/tests/tfma.c +--- mpfr-3.0.0-a/tests/tfma.c 2010-06-10 11:00:13.000000000 +0000 ++++ mpfr-3.0.0-b/tests/tfma.c 2010-10-21 20:59:32.000000000 +0000 +@@ -337,6 +337,94 @@ + mpfr_clears (x, y, z, r, (mpfr_ptr) 0); + } + ++static void ++bug20101018 (void) ++{ ++ mpfr_t x, y, z, t, u; ++ int i; ++ ++ mpfr_init2 (x, 64); ++ mpfr_init2 (y, 64); ++ mpfr_init2 (z, 64); ++ mpfr_init2 (t, 64); ++ mpfr_init2 (u, 64); ++ ++ mpfr_set_str (x, "0xf.fffffffffffffffp-14766", 16, MPFR_RNDN); ++ mpfr_set_str (y, "-0xf.fffffffffffffffp+317", 16, MPFR_RNDN); ++ mpfr_set_str (z, "0x8.3ffffffffffe3ffp-14443", 16, MPFR_RNDN); ++ mpfr_set_str (t, "0x8.7ffffffffffc7ffp-14444", 16, MPFR_RNDN); ++ i = mpfr_fma (u, x, y, z, MPFR_RNDN); ++ if (mpfr_cmp (u, t) != 0) ++ { ++ printf ("Wrong result in bug20101018 (a)\n"); ++ printf ("Expected "); ++ mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN); ++ printf ("\nGot "); ++ mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN); ++ printf ("\n"); ++ exit (1); ++ } ++ if (i <= 0) ++ { ++ printf ("Wrong ternary value in bug20101018 (a)\n"); ++ printf ("Expected > 0\n"); ++ printf ("Got %d\n", i); ++ exit (1); ++ } ++ ++ mpfr_set_str (x, "-0xf.fffffffffffffffp-11420", 16, MPFR_RNDN); ++ mpfr_set_str (y, "0xf.fffffffffffffffp+9863", 16, MPFR_RNDN); ++ mpfr_set_str (z, "0x8.fffff80ffffffffp-1551", 16, MPFR_RNDN); ++ mpfr_set_str (t, "0x9.fffff01ffffffffp-1552", 16, MPFR_RNDN); ++ i = mpfr_fma (u, x, y, z, MPFR_RNDN); ++ if (mpfr_cmp (u, t) != 0) ++ { ++ printf ("Wrong result in bug20101018 (b)\n"); ++ printf ("Expected "); ++ mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN); ++ printf ("\nGot "); ++ mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN); ++ printf ("\n"); ++ exit (1); ++ } ++ if (i <= 0) ++ { ++ printf ("Wrong ternary value in bug20101018 (b)\n"); ++ printf ("Expected > 0\n"); ++ printf ("Got %d\n", i); ++ exit (1); ++ } ++ ++ mpfr_set_str (x, "0xf.fffffffffffffffp-2125", 16, MPFR_RNDN); ++ mpfr_set_str (y, "-0xf.fffffffffffffffp-6000", 16, MPFR_RNDN); ++ mpfr_set_str (z, "0x8p-8119", 16, MPFR_RNDN); ++ mpfr_set_str (t, "0x8.000000000000001p-8120", 16, MPFR_RNDN); ++ i = mpfr_fma (u, x, y, z, MPFR_RNDN); ++ if (mpfr_cmp (u, t) != 0) ++ { ++ printf ("Wrong result in bug20101018 (c)\n"); ++ printf ("Expected "); ++ mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN); ++ printf ("\nGot "); ++ mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN); ++ printf ("\n"); ++ exit (1); ++ } ++ if (i <= 0) ++ { ++ printf ("Wrong ternary value in bug20101018 (c)\n"); ++ printf ("Expected > 0\n"); ++ printf ("Got %d\n", i); ++ exit (1); ++ } ++ ++ mpfr_clear (x); ++ mpfr_clear (y); ++ mpfr_clear (z); ++ mpfr_clear (t); ++ mpfr_clear (u); ++} ++ + int + main (int argc, char *argv[]) + { +@@ -345,6 +433,8 @@ + + tests_start_mpfr (); + ++ bug20101018 (); ++ + mpfr_init (x); + mpfr_init (s); + mpfr_init (y); +diff -Naurd mpfr-3.0.0-a/tests/tsub.c mpfr-3.0.0-b/tests/tsub.c +--- mpfr-3.0.0-a/tests/tsub.c 2010-06-10 11:00:13.000000000 +0000 ++++ mpfr-3.0.0-b/tests/tsub.c 2010-10-21 20:59:32.000000000 +0000 +@@ -201,6 +201,8 @@ + if (mpfr_cmp (z, x)) + { + printf ("Error in mpfr_sub (2)\n"); ++ printf ("Expected "); mpfr_print_binary (x); puts (""); ++ printf ("Got "); mpfr_print_binary (z); puts (""); + exit (1); + } + mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101"); +@@ -478,6 +480,156 @@ + mpfr_clear (u); + } + ++/* Bug found by Jakub Jelinek ++ * http://bugzilla.redhat.com/643657 ++ * https://gforge.inria.fr/tracker/index.php?func=detail&aid=11301 ++ * The consequence can be either an assertion failure (i = 2 in the ++ * testcase below, in debug mode) or an incorrectly rounded value. ++ */ ++static void ++bug20101017 (void) ++{ ++ mpfr_t a, b, c; ++ int inex; ++ int i; ++ ++ mpfr_init2 (a, GMP_NUMB_BITS * 2); ++ mpfr_init2 (b, GMP_NUMB_BITS); ++ mpfr_init2 (c, GMP_NUMB_BITS); ++ ++ /* a = 2^(2N) + k.2^(2N-1) + 2^N and b = 1 ++ with N = GMP_NUMB_BITS and k = 0 or 1. ++ c = a - b should round to the same value as a. */ ++ ++ for (i = 2; i <= 3; i++) ++ { ++ mpfr_set_ui_2exp (a, i, GMP_NUMB_BITS - 1, MPFR_RNDN); ++ mpfr_add_ui (a, a, 1, MPFR_RNDN); ++ mpfr_mul_2ui (a, a, GMP_NUMB_BITS, MPFR_RNDN); ++ mpfr_set_ui (b, 1, MPFR_RNDN); ++ inex = mpfr_sub (c, a, b, MPFR_RNDN); ++ mpfr_set (b, a, MPFR_RNDN); ++ if (! mpfr_equal_p (c, b)) ++ { ++ printf ("Error in bug20101017 for i = %d.\n", i); ++ printf ("Expected "); ++ mpfr_out_str (stdout, 16, 0, b, MPFR_RNDN); ++ putchar ('\n'); ++ printf ("Got "); ++ mpfr_out_str (stdout, 16, 0, c, MPFR_RNDN); ++ putchar ('\n'); ++ exit (1); ++ } ++ if (inex >= 0) ++ { ++ printf ("Error in bug20101017 for i = %d: bad inex value.\n", i); ++ printf ("Expected negative, got %d.\n", inex); ++ exit (1); ++ } ++ } ++ ++ mpfr_set_prec (a, 64); ++ mpfr_set_prec (b, 129); ++ mpfr_set_prec (c, 2); ++ mpfr_set_str_binary (b, "0.100000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000001E65"); ++ mpfr_set_str_binary (c, "0.10E1"); ++ inex = mpfr_sub (a, b, c, MPFR_RNDN); ++ if (mpfr_cmp_ui_2exp (a, 1, 64) != 0 || inex >= 0) ++ { ++ printf ("Error in mpfr_sub for b-c for b=2^64+1+2^(-64), c=1\n"); ++ printf ("Expected result 2^64 with inex < 0\n"); ++ printf ("Got "); mpfr_print_binary (a); ++ printf (" with inex=%d\n", inex); ++ exit (1); ++ } ++ ++ mpfr_clears (a, b, c, (mpfr_ptr) 0); ++} ++ ++/* hard test of rounding */ ++static void ++check_rounding (void) ++{ ++ mpfr_t a, b, c, res; ++ mpfr_prec_t p; ++ long k, l; ++ int i; ++ ++#define MAXKL (2 * GMP_NUMB_BITS) ++ for (p = MPFR_PREC_MIN; p <= GMP_NUMB_BITS; p++) ++ { ++ mpfr_init2 (a, p); ++ mpfr_init2 (res, p); ++ mpfr_init2 (b, p + 1 + MAXKL); ++ mpfr_init2 (c, MPFR_PREC_MIN); ++ ++ /* b = 2^p + 1 + 2^(-k), c = 2^(-l) */ ++ for (k = 0; k <= MAXKL; k++) ++ for (l = 0; l <= MAXKL; l++) ++ { ++ mpfr_set_ui_2exp (b, 1, p, MPFR_RNDN); ++ mpfr_add_ui (b, b, 1, MPFR_RNDN); ++ mpfr_mul_2ui (b, b, k, MPFR_RNDN); ++ mpfr_add_ui (b, b, 1, MPFR_RNDN); ++ mpfr_div_2ui (b, b, k, MPFR_RNDN); ++ mpfr_set_ui_2exp (c, 1, -l, MPFR_RNDN); ++ i = mpfr_sub (a, b, c, MPFR_RNDN); ++ /* b - c = 2^p + 1 + 2^(-k) - 2^(-l), should be rounded to ++ 2^p for l <= k, and 2^p+2 for l < k */ ++ if (l <= k) ++ { ++ if (mpfr_cmp_ui_2exp (a, 1, p) != 0) ++ { ++ printf ("Wrong result in check_rounding\n"); ++ printf ("p=%lu k=%ld l=%ld\n", p, k, l); ++ printf ("b="); mpfr_print_binary (b); puts (""); ++ printf ("c="); mpfr_print_binary (c); puts (""); ++ printf ("Expected 2^%lu\n", p); ++ printf ("Got "); mpfr_print_binary (a); puts (""); ++ exit (1); ++ } ++ if (i >= 0) ++ { ++ printf ("Wrong ternary value in check_rounding\n"); ++ printf ("p=%lu k=%ld l=%ld\n", p, k, l); ++ printf ("b="); mpfr_print_binary (b); puts (""); ++ printf ("c="); mpfr_print_binary (c); puts (""); ++ printf ("a="); mpfr_print_binary (a); puts (""); ++ printf ("Expected < 0, got %d\n", i); ++ exit (1); ++ } ++ } ++ else /* l < k */ ++ { ++ mpfr_set_ui_2exp (res, 1, p, MPFR_RNDN); ++ mpfr_add_ui (res, res, 2, MPFR_RNDN); ++ if (mpfr_cmp (a, res) != 0) ++ { ++ printf ("Wrong result in check_rounding\n"); ++ printf ("b="); mpfr_print_binary (b); puts (""); ++ printf ("c="); mpfr_print_binary (c); puts (""); ++ printf ("Expected "); mpfr_print_binary (res); puts (""); ++ printf ("Got "); mpfr_print_binary (a); puts (""); ++ exit (1); ++ } ++ if (i <= 0) ++ { ++ printf ("Wrong ternary value in check_rounding\n"); ++ printf ("b="); mpfr_print_binary (b); puts (""); ++ printf ("c="); mpfr_print_binary (c); puts (""); ++ printf ("Expected > 0, got %d\n", i); ++ exit (1); ++ } ++ } ++ } ++ ++ mpfr_clear (a); ++ mpfr_clear (res); ++ mpfr_clear (b); ++ mpfr_clear (c); ++ } ++} ++ + #define TEST_FUNCTION test_sub + #define TWO_ARGS + #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS) +@@ -491,6 +643,8 @@ + + tests_start_mpfr (); + ++ bug20101017 (); ++ check_rounding (); + check_diverse (); + check_inexact (); + bug_ddefour (); +diff -Naurd mpfr-3.0.0-a/version.c mpfr-3.0.0-b/version.c +--- mpfr-3.0.0-a/version.c 2010-10-21 20:28:38.000000000 +0000 ++++ mpfr-3.0.0-b/version.c 2010-10-21 20:59:32.000000000 +0000 +@@ -25,5 +25,5 @@ + const char * + mpfr_get_version (void) + { +- return "3.0.0-p5"; ++ return "3.0.0-p6"; + } |