diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2017-12-23 01:39:00 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2017-12-23 01:39:00 +0000 |
commit | 7733b3e1a28763816d3df3d4c763a72b5951b0f4 (patch) | |
tree | 8afb3142f3cf684571cb343af900e34e9775f4de /src | |
parent | 7f6fb5048e81630e682b1784403a0fea77dc53af (diff) | |
download | mpfr-7733b3e1a28763816d3df3d4c763a72b5951b0f4.tar.gz |
Merged changesets r12026-12045 from the trunk (bug fixes and tests).
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/4.0@12046 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src')
-rw-r--r-- | src/hypot.c | 8 | ||||
-rw-r--r-- | src/lngamma.c | 13 | ||||
-rw-r--r-- | src/sin.c | 5 | ||||
-rw-r--r-- | src/subnormal.c | 2 |
4 files changed, 20 insertions, 8 deletions
diff --git a/src/hypot.c b/src/hypot.c index 07ceccfb7..246f55d0b 100644 --- a/src/hypot.c +++ b/src/hypot.c @@ -98,7 +98,13 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode) /* If z > abs(x), then it was already rounded up; otherwise z = abs(x), and we need to add one ulp due to y. */ if (mpfr_abs (z, x, rnd_mode) == 0) - mpfr_nexttoinf (z); + { + mpfr_nexttoinf (z); + /* since mpfr_nexttoinf does not set the overflow flag, + we have to check manually for overflow */ + if (MPFR_UNLIKELY (MPFR_IS_INF (z))) + MPFR_SET_OVERFLOW (); + } MPFR_RET (1); } else /* MPFR_RNDZ, MPFR_RNDD, MPFR_RNDN */ diff --git a/src/lngamma.c b/src/lngamma.c index 147ef6bfc..5511fd1dc 100644 --- a/src/lngamma.c +++ b/src/lngamma.c @@ -350,7 +350,10 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) ulp(u)/2 + (2-z0)*max(1,log(2-z0))*2^(1-w) = (1/2 + (2-z0)*max(1,log(2-z0))*2^(1-E(u))) ulp(u) */ d = (double) MPFR_GET_EXP(s) * 0.694; /* upper bound for log(2-z0) */ - err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u); + if (MPFR_IS_ZERO(u)) /* in that case the error on u is zero */ + err_u = 0; + else + err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u); err_u = (err_u >= 0) ? err_u + 1 : 0; /* now the error on u is bounded by 2^err_u ulps */ @@ -397,11 +400,15 @@ GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd) } else { - err_s += 1 - MPFR_GET_EXP(v); + /* if v = 0 here, it was 1 before the call to mpfr_log, + thus the error on v was zero */ + if (!MPFR_IS_ZERO(v)) + err_s += 1 - MPFR_GET_EXP(v); err_s = (err_s >= 0) ? err_s + 1 : 0; /* the error on v is bounded by 2^err_s ulps */ err_u += MPFR_GET_EXP(u); /* absolute error on u */ - err_s += MPFR_GET_EXP(v); /* absolute error on v */ + if (!MPFR_IS_ZERO(v)) /* same as above */ + err_s += MPFR_GET_EXP(v); /* absolute error on v */ mpfr_sub (s, v, u, MPFR_RNDN); /* the total error on s is bounded by ulp(s)/2 + 2^(err_u-w) + 2^(err_s-w) <= ulp(s)/2 + 2^(max(err_u,err_s)+1-w) */ @@ -148,9 +148,8 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) sign = MPFR_SIGN(xx); /* now that the argument is reduced, precision m is enough */ mpfr_set_prec (c, m); - mpfr_cos (c, xx, MPFR_RNDZ); /* can't be exact */ - mpfr_nexttoinf (c); /* now c = cos(x) rounded away */ - mpfr_mul (c, c, c, MPFR_RNDU); /* away */ + mpfr_cos (c, xx, MPFR_RNDA); /* c = cos(x) rounded away */ + mpfr_mul (c, c, c, MPFR_RNDU); /* away */ mpfr_ui_sub (c, 1, c, MPFR_RNDZ); mpfr_sqrt (c, c, MPFR_RNDZ); if (MPFR_IS_NEG_SIGN(sign)) diff --git a/src/subnormal.c b/src/subnormal.c index 11a6af3ca..0a135aa9f 100644 --- a/src/subnormal.c +++ b/src/subnormal.c @@ -144,7 +144,7 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd) { if (SAME_SIGN (inexact, MPFR_INT_SIGN (y))) mpfr_nexttozero (dest); - else + else /* subnormal range, thus no overflow */ mpfr_nexttoinf (dest); inexact = -inexact; } |