diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-10 23:27:03 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-10 23:27:03 +0000 |
commit | 30a7e3ccc8525d6a1300ce03fbcbeb8b59143ae3 (patch) | |
tree | 3a932455fc65194be45607ccdf410322e0d7f38a | |
parent | 32a9b16a85bc52289be725edb06ae6b95cafa9c8 (diff) | |
download | mpfr-30a7e3ccc8525d6a1300ce03fbcbeb8b59143ae3.tar.gz |
[src/jyn_asympt.c] Fixed bug when sin(z)±cos(z) rounds to 0.
[tests/tj0.c] Testcase.
(merged changesets r14400,14405 from the trunk)
git-svn-id: https://scm.gforge.inria.fr/anonscm/svn/mpfr/branches/4.1@14426 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | src/jyn_asympt.c | 15 | ||||
-rw-r--r-- | tests/tj0.c | 21 |
2 files changed, 34 insertions, 2 deletions
diff --git a/src/jyn_asympt.c b/src/jyn_asympt.c index 86fefe4c4..6e5c7c37d 100644 --- a/src/jyn_asympt.c +++ b/src/jyn_asympt.c @@ -69,6 +69,8 @@ FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r) MPFR_ZIV_INIT (loop, w); for (;;) { + int ok = 1; + mpfr_set_prec (c, w); mpfr_init2 (s, w); mpfr_init2 (P, w); @@ -92,6 +94,13 @@ FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r) /* now s approximates sin(z)+cos(z), and c approximates sin(z)-cos(z), with total absolute error bounded by 2^(1-w). */ + /* if s or c is zero, MPFR_GET_EXP will fail below */ + if (MPFR_IS_ZERO(s) || MPFR_IS_ZERO(c)) + { + ok = 0; + goto clear; + } + /* precompute 1/(8|z|) */ mpfr_si_div (iz, MPFR_IS_POS(z) ? 1 : -1, z, MPFR_RNDN); /* err <= 1 */ mpfr_div_2ui (iz, iz, 3, MPFR_RNDN); @@ -257,6 +266,9 @@ FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r) err = (err >= err2) ? err + 1 : err2 + 1; /* the absolute error on c is bounded by 2^(err - w) */ + err -= MPFR_GET_EXP (c); + + clear: mpfr_clear (s); mpfr_clear (P); mpfr_clear (Q); @@ -266,8 +278,7 @@ FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r) mpfr_clear (err_s); mpfr_clear (err_u); - err -= MPFR_GET_EXP (c); - if (MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r))) + if (ok && MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r))) break; if (diverge != 0) { diff --git a/tests/tj0.c b/tests/tj0.c index 14e30bf6d..6f6f02124 100644 --- a/tests/tj0.c +++ b/tests/tj0.c @@ -27,6 +27,25 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define REDUCE_EMAX 262143 /* otherwise arg. reduction is too expensive */ #include "tgeneric.c" +/* bug found in revision 14399 with GMP_CHECK_RANDOMIZE=1612721106588971 */ +static void +bug20210208 (void) +{ + mpfr_t x, y; + int inex; + + mpfr_init2 (x, 79); + mpfr_init2 (y, 1); + mpfr_set_str (x, "2.552495117262005805960565e+02", 10, MPFR_RNDN); + mpfr_clear_flags (); + inex = mpfr_j0 (y, x, MPFR_RNDZ); + MPFR_ASSERTN (mpfr_cmp_si_2exp (y, -1, -5) == 0); + MPFR_ASSERTN (inex > 0); + MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_INEXACT); + mpfr_clear (x); + mpfr_clear (y); +} + int main (int argc, char *argv[]) { @@ -35,6 +54,8 @@ main (int argc, char *argv[]) tests_start_mpfr (); + bug20210208 (); + mpfr_init (x); mpfr_init (y); |