diff options
author | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2004-04-20 09:29:48 +0000 |
---|---|---|
committer | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2004-04-20 09:29:48 +0000 |
commit | 87e63696c2bfac2de0ff5978d608d6dfb8ea05db (patch) | |
tree | d2d151575bc2e66465604ff8b27518c6ab6a6a94 | |
parent | c2a44912370e67df5da882ae044c29989d5511af (diff) | |
download | mpfr-87e63696c2bfac2de0ff5978d608d6dfb8ea05db.tar.gz |
Fix a bug when both op are < 0 (It seems it was my fault).
Add a test to check it.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2869 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | agm.c | 127 | ||||
-rw-r--r-- | tests/tagm.c | 6 |
2 files changed, 74 insertions, 59 deletions
@@ -43,21 +43,21 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) MPFR_RET_NAN; } /* now one of a or b is Inf or 0 */ - /* If a or b is -Inf or +Inf, the result is +Inf if both operands are - stricly positive, and NaN otherwise */ + /* If a and b is +Inf, the result is +Inf. + Otherwise if a or b is -Inf, the result is NaN */ else if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2)) { - if (!(MPFR_IS_STRICTPOS(op1) && MPFR_IS_STRICTPOS(op2))) - { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - else /* one of a, b is +Inf, the other is a "normal" number */ + if (MPFR_IS_STRICTPOS(op1) && MPFR_IS_STRICTPOS(op2)) { MPFR_SET_INF(r); MPFR_SET_SAME_SIGN(r, op1); MPFR_RET(0); /* exact */ } + else + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } } else /* a and b are neither NaN nor Inf, and one is zero */ { /* If a or b is 0, the result is +0 since a sqrt is positive */ @@ -69,7 +69,14 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) } MPFR_CLEAR_FLAGS(r); - /* precision of the following calculus */ + /* If a or b is negative (excluding -Infinity), the result is NaN */ + if (MPFR_UNLIKELY(MPFR_IS_NEG(op1) || MPFR_IS_NEG(op2))) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + + /* precision of the following calculus */ q = MPFR_PREC(r); p = q + 15; @@ -114,57 +121,59 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) /* Main loop */ - while (go_on) { - int err, can_round; - mp_prec_t eq; - double erraux; - + while (go_on) + { + int err, can_round; + mp_prec_t eq; + double erraux; + #if 0 - erraux = (vo == uo) ? 0.0 - : __gmpfr_ceil_exp2 (-2.0 * (double) p * uo / (vo - uo)); - err = 1 + (int) ((3.0 / 2.0 * (double) __gmpfr_ceil_log2 ((double) p) - + 1.0) * __gmpfr_ceil_exp2 (- (double) p) - + 3.0 * erraux); + erraux = (vo == uo) ? 0.0 + : __gmpfr_ceil_exp2 (-2.0 * (double) p * uo / (vo - uo)); + err = 1 + (int) ((3.0 / 2.0 * (double) __gmpfr_ceil_log2 ((double) p) + + 1.0) * __gmpfr_ceil_exp2 (- (double) p) + + 3.0 * erraux); #else - /* since the argument of __gmpfr_ceil_exp2() is negative, - erraux is always bounded by 1. - We also have floor((3/2*ceil(log2(p))+1) * 2^(-p)) = 0, - since it the value f(p) inside floor() is 5/8 for p=2, 1/2 for p=3, - and f(2p) <= f(p)/2 + 1/2^p. So err=4. - */ - erraux = 1.0; - err = 4; + /* since the argument of __gmpfr_ceil_exp2() is negative, + erraux is always bounded by 1. + We also have floor((3/2*ceil(log2(p))+1) * 2^(-p)) = 0, + since it the value f(p) inside floor() is 5/8 for p=2, 1/2 for p=3, + and f(2p) <= f(p)/2 + 1/2^p. So err=4. + */ + erraux = 1.0; + err = 4; #endif #if 0 /* this can not happen since p = q+15, and err=4 */ - if (p - err - 3 <= q) - { - p = q + err + 4; - err = 1 + (int) ((3.0 / 2.0 * __gmpfr_ceil_log2 ((double) p) + 1.0) - * __gmpfr_ceil_exp2 (- (double) p) + 3.0 * erraux); - } + if (p - err - 3 <= q) + { + p = q + err + 4; + err = 1 + (int) ((3.0 / 2.0 * __gmpfr_ceil_log2 ((double) p) + 1.0) + * __gmpfr_ceil_exp2 (- (double) p) + 3.0 * erraux); + } #endif - - /* Calculus of un and vn */ - do - { - mpfr_mul (tmp, u, v, GMP_RNDN); - mpfr_sqrt (tmpu, tmp, GMP_RNDN); - mpfr_add (tmp, u, v, GMP_RNDN); - mpfr_div_2ui (tmpv, tmp, 1, GMP_RNDN); - mpfr_set (u, tmpu, GMP_RNDN); - mpfr_set (v, tmpv, GMP_RNDN); - } - while (mpfr_cmp2 (u, v, &eq) != 0 && eq <= p - 2); - - /* Roundability of the result */ + + /* Calculus of un and vn */ + do + { + mpfr_mul (tmp, u, v, GMP_RNDN); + mpfr_sqrt (tmpu, tmp, GMP_RNDN); + mpfr_add (tmp, u, v, GMP_RNDN); + mpfr_div_2ui (tmpv, tmp, 1, GMP_RNDN); + mpfr_set (u, tmpu, GMP_RNDN); + mpfr_set (v, tmpv, GMP_RNDN); + } + while (mpfr_cmp2 (u, v, &eq) != 0 && eq <= p - 2); + + /* Roundability of the result */ can_round = mpfr_can_round (v, p - err - 3, GMP_RNDN, GMP_RNDZ, - q + (rnd_mode == GMP_RNDN)); - + q + (rnd_mode == GMP_RNDN)); + if (can_round) go_on = 0; - - else { + + else + { go_on = 1; p += 5; s = (p - 1) / BITS_PER_MP_LIMB + 1; @@ -175,8 +184,8 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) MPFR_TMP_INIT(tmpp, tmp, p, s); mpfr_set (u, a, GMP_RNDN); mpfr_set (v, b, GMP_RNDN); - } - } + } + } /* End of while */ /* Setting of the result */ @@ -187,10 +196,10 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) TMP_FREE(marker); return inexact; /* agm(u,v) can be exact for u, v rational only for u=v. - Proof (due to Nicolas Brisebarre): it suffices to consider - u=1 and v<1. Then 1/AGM(1,v) = 2F1(1/2,1/2,1;1-v^2), - and a theorem due to G.V. Chudnovsky states that for x a - non-zero algebraic number with |x|<1, then - 2F1(1/2,1/2,1;x) and 2F1(-1/2,1/2,1;x) are algebraically - independent over Q. */ + Proof (due to Nicolas Brisebarre): it suffices to consider + u=1 and v<1. Then 1/AGM(1,v) = 2F1(1/2,1/2,1;1-v^2), + and a theorem due to G.V. Chudnovsky states that for x a + non-zero algebraic number with |x|<1, then + 2F1(1/2,1/2,1;x) and 2F1(-1/2,1/2,1;x) are algebraically + independent over Q. */ } diff --git a/tests/tagm.c b/tests/tagm.c index af6c5aecd..73066af79 100644 --- a/tests/tagm.c +++ b/tests/tagm.c @@ -176,6 +176,12 @@ check_nans (void) mpfr_agm (m, x, y, GMP_RNDN); MPFR_ASSERTN (mpfr_cmp_ui (m ,1) == 0); + /* agm(-1,-2) == NaN */ + mpfr_set_si (x, -1, GMP_RNDN); + mpfr_set_si (y, -2, GMP_RNDN); + mpfr_agm (m, x, y, GMP_RNDN); + MPFR_ASSERTN (mpfr_nan_p (m)); + mpfr_clear (x); mpfr_clear (y); mpfr_clear (m); |