diff options
| -rw-r--r-- | src/mpfr-impl.h | 3 | ||||
| -rw-r--r-- | src/sum.c | 37 | ||||
| -rw-r--r-- | tests/tsum.c | 39 |
3 files changed, 64 insertions, 15 deletions
diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index 7a4d00569..00a38bad2 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -1876,7 +1876,8 @@ __MPFR_DECLSPEC int mpfr_round_raw_4 _MPFR_PROTO ((mp_limb_t *, __MPFR_DECLSPEC int mpfr_check _MPFR_PROTO ((mpfr_srcptr)); __MPFR_DECLSPEC int mpfr_sum_sort _MPFR_PROTO ((mpfr_srcptr *const, - unsigned long, mpfr_srcptr *)); + unsigned long, mpfr_srcptr *, + mpfr_prec_t *)); __MPFR_DECLSPEC int mpfr_get_cputime _MPFR_PROTO ((void)); @@ -45,9 +45,13 @@ static void count_sort (mpfr_srcptr *const, unsigned long, mpfr_srcptr *, mpfr_exp_t, mpfr_uexp_t); /* Either sort the tab in perm and returns 0 - Or returns 1 for +INF, -1 for -INF and 2 for NAN */ + Or returns 1 for +INF, -1 for -INF and 2 for NAN. + Also set *maxprec to the maximal precision of tab[0..n-1] and of the + initial value of *maxprec. +*/ int -mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm) +mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm, + mpfr_prec_t *maxprec) { mpfr_exp_t min, max; mpfr_uexp_t exp_num; @@ -79,6 +83,8 @@ mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm) if (MPFR_GET_EXP (tab[i]) > max) max = MPFR_GET_EXP(tab[i]); } + if (MPFR_PREC (tab[i]) > *maxprec) + *maxprec = MPFR_PREC (tab[i]); } if (MPFR_UNLIKELY (sign_inf != 0)) return sign_inf; @@ -213,7 +219,8 @@ heap_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm) /* Sum a list of float with order given by permutation perm, - * intermediate size set to F. + * intermediate size set to F. Return non-zero if at least one of + * the operations is inexact (thus 0 implies that the sum is exact). * Internal use function. */ static int @@ -230,16 +237,19 @@ sum_once (mpfr_ptr ret, mpfr_srcptr *const tab, unsigned long n, mpfr_prec_t F) for (i = 1; i < n - 1; i++) { MPFR_ASSERTD (!MPFR_IS_NAN (sum) && !MPFR_IS_INF (sum)); - error_trap |= mpfr_add (sum, sum, tab[i], MPFR_RNDN); + if (mpfr_add (sum, sum, tab[i], MPFR_RNDN)) + error_trap = 1; } - error_trap |= mpfr_add (ret, sum, tab[n - 1], MPFR_RNDN); + if (mpfr_add (ret, sum, tab[n - 1], MPFR_RNDN)) + error_trap = 1; mpfr_clear (sum); return error_trap; } /* Sum a list of floating-point numbers. + * If the return value is 0, then the sum is exact. + * Otherwise the return value gives no information. */ - int mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mpfr_rnd_t rnd) { @@ -266,7 +276,8 @@ mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mpfr_rnd_t rnd) /* Sort and treat special cases */ MPFR_TMP_MARK (marker); perm = (mpfr_srcptr *) MPFR_TMP_ALLOC (n * sizeof *perm); - error_trap = mpfr_sum_sort (tab, n, perm); + prec = MPFR_PREC (ret); + error_trap = mpfr_sum_sort (tab, n, perm, &prec); /* Check if there was a NAN or a INF */ if (MPFR_UNLIKELY (error_trap != 0)) { @@ -281,8 +292,7 @@ mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mpfr_rnd_t rnd) MPFR_RET (0); } - /* Initial precision */ - prec = MAX (MPFR_PREC (tab[0]), MPFR_PREC (ret)); + /* Initial precision is max(prec(ret),prec(tab[0]),...,prec(tab[n-1])) */ k = MPFR_INT_CEIL_LOG2 (n) + 1; prec += k + 2; mpfr_init2 (cur_sum, prec); @@ -295,8 +305,7 @@ mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mpfr_rnd_t rnd) error_trap = sum_once (cur_sum, perm, n, prec + k); if (MPFR_LIKELY (error_trap == 0 || (!MPFR_IS_ZERO (cur_sum) && - mpfr_can_round (cur_sum, - MPFR_GET_EXP (cur_sum) - prec + 2, + mpfr_can_round (cur_sum, prec - 2, MPFR_RNDN, rnd, MPFR_PREC (ret))))) break; MPFR_ZIV_NEXT (loop, prec); @@ -305,11 +314,13 @@ mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mpfr_rnd_t rnd) MPFR_ZIV_FREE (loop); MPFR_TMP_FREE (marker); - error_trap |= mpfr_set (ret, cur_sum, rnd); + if (mpfr_set (ret, cur_sum, rnd)) + error_trap = 1; mpfr_clear (cur_sum); MPFR_SAVE_EXPO_FREE (expo); - error_trap |= mpfr_check_range (ret, 0, rnd); + if (mpfr_check_range (ret, 0, rnd)) + error_trap = 1; return error_trap; /* It doesn't return the ternary value */ } diff --git a/tests/tsum.c b/tests/tsum.c index 5d9c8b598..38ebf52e3 100644 --- a/tests/tsum.c +++ b/tests/tsum.c @@ -126,6 +126,7 @@ test_sort (mpfr_prec_t f, unsigned long n) mpfr_ptr *tabtmp; mpfr_srcptr *perm; unsigned long i; + mpfr_prec_t prec = MPFR_PREC_MIN; /* Init stuff */ tab = (mpfr_t *) (*__gmp_allocate_func) (n * sizeof (mpfr_t)); @@ -140,7 +141,7 @@ test_sort (mpfr_prec_t f, unsigned long n) tabtmp[i] = tab[i]; } - mpfr_sum_sort ((mpfr_srcptr *)tabtmp, n, perm); + mpfr_sum_sort ((mpfr_srcptr *)tabtmp, n, perm, &prec); if (check_is_sorted (n, perm) == 0) { @@ -300,6 +301,41 @@ void check_special (void) mpfr_clears (tab[0], tab[1], tab[2], r, (mpfr_ptr) 0); } +/* bug reported by Joseph S. Myers on 2013-10-27 + https://sympa.inria.fr/sympa/arc/mpfr/2013-10/msg00015.html */ +static void +bug20131027 (void) +{ + mpfr_t r, t[4]; + mpfr_ptr p[4]; + char *s[4] = { + "0x1p1000", + "-0x0.fffffffffffff80000000000000001p1000", + "-0x1p947", + "0x1p880" + }; + int i; + + mpfr_init2 (r, 53); + for (i = 0; i < 4; i++) + { + mpfr_init2 (t[i], i == 0 ? 53 : 1000); + mpfr_set_str (t[i], s[i], 0, MPFR_RNDN); + p[i] = t[i]; + } + mpfr_sum (r, p, 4, MPFR_RNDN); + + if (MPFR_NOTZERO (r)) + { + printf ("mpfr_sum incorrect in bug20131027: expected 0, got\n"); + mpfr_dump (r); + exit (1); + } + + for (i = 0; i < 4; i++) + mpfr_clear (t[i]); + mpfr_clear (r); +} int main (void) @@ -310,6 +346,7 @@ main (void) tests_start_mpfr (); check_special (); + bug20131027 (); test_sort (1764, 1026); for (p = 2 ; p < 444 ; p += 17) for (n = 2 ; n < 1026 ; n += 42 + p) |