diff options
| author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-08 08:08:27 +0000 |
|---|---|---|
| committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-08 08:08:27 +0000 |
| commit | 21ca2f258539b599eacc782aa19a750af299b2c2 (patch) | |
| tree | 143b648abd5be81d89bf4ff10c4e5b190fd96892 /src/jyn_asympt.c | |
| parent | 5f1a77e3de3e09714004b7a13aa2f46da756bad4 (diff) | |
| download | mpfr-21ca2f258539b599eacc782aa19a750af299b2c2.tar.gz | |
[src/jyn_asympt.c] fixed bug when sin(z)+cos(z) or sin(z)-cos(z) round to 0
[tests/tj0.c] added corresponding non-regression test
git-svn-id: https://scm.gforge.inria.fr/anonscm/svn/mpfr/trunk@14400 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src/jyn_asympt.c')
| -rw-r--r-- | src/jyn_asympt.c | 15 |
1 files changed, 13 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) { |