diff options
| -rw-r--r-- | src/fma.c | 8 | ||||
| -rw-r--r-- | tests/tfma.c | 24 |
2 files changed, 29 insertions, 3 deletions
@@ -234,7 +234,7 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, /* Let's eliminate the obvious case where x*y and z have the same sign. No possible cancellation -> real overflow. Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3, - then |x*y| >= 2^(emax+1), and |x*y + z| >= 2^emax. This case + then |x*y| >= 2^(emax+1), and |x*y + z| > 2^emax. This case is also an overflow. */ if (MPFR_SIGN (u) == MPFR_SIGN (z) || e >= __gmpfr_emax + 3) { @@ -253,7 +253,8 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, inexact = mpfr_mul (u, u, y, MPFR_RNDN); MPFR_ASSERTN (inexact == 0); - /* Now, we need to add z/4... But it may underflow! */ + /* Now, we need to add z/4, which is in fact a subtraction since + x*y and z have different signs. But it may underflow! */ { mpfr_t zo4; mpfr_srcptr zz; @@ -262,7 +263,8 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, if (MPFR_GET_EXP (u) > MPFR_GET_EXP (z) && MPFR_GET_EXP (u) - MPFR_GET_EXP (z) > MPFR_PREC (u)) { - /* |z| < ulp(u)/2, therefore one can use z instead of z/4. */ + /* |z| < ulp(u)/2, therefore one can use z instead of z/4, + except maybe when u is a power of 2! */ zz = z; } else diff --git a/tests/tfma.c b/tests/tfma.c index f3a3b30e7..b8918a192 100644 --- a/tests/tfma.c +++ b/tests/tfma.c @@ -196,6 +196,29 @@ test_overflow2 (void) } static void +test_overflow3 (void) +{ + mpfr_t x, y, z, r; + int inex; + mpfr_prec_t p = 8; + + mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0); + mpfr_init2 (r, 2 * p); + mpfr_set_ui_2exp (x, 1, mpfr_get_emax () - 1, MPFR_RNDN); + mpfr_set_ui (y, 2, MPFR_RNDN); /* y = 2 */ + mpfr_set_si_2exp (z, -1, mpfr_get_emax () - 2*p - 2, MPFR_RNDN); + mpfr_nextabove (z); + mpfr_clear_flags (); + /* We have x*y = 2^emax and z = -2^(emax-2p-2)*(1-2^(-p)) thus + x*y+z = 2^emax - 2^(emax-2p-2) + 2^(emax-3p-2) should overflow, + since it is closest from 2^emax than from 2^emax - 2^(emax-2p). */ + inex = mpfr_fma (r, x, y, z, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_inf_p (r) && mpfr_sgn (r) > 0); + mpfr_clears (x, y, z, r, (mpfr_ptr) 0); +} + +static void test_underflow1 (void) { mpfr_t x, y, z, r; @@ -833,6 +856,7 @@ main (int argc, char *argv[]) set_emax (MPFR_EMAX_MAX); test_overflow1 (); test_overflow2 (); + test_overflow3 (); test_underflow1 (); test_underflow2 (); test_underflow3 (2); |