diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-01-20 15:29:28 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-01-20 15:29:28 +0000 |
commit | 6c94f4e1e21d6c5d2eee5082c53fecc774c18740 (patch) | |
tree | 3c62332d413e91cf9dced3583a4d385b4e4edb09 | |
parent | f010fe07b165160a455217ce143244a38eb72704 (diff) | |
download | mpfr-6c94f4e1e21d6c5d2eee5082c53fecc774c18740.tar.gz |
Fixed bug for zeta(s) with s near an even negative integer.
(merged changesets r9852-9854 from the trunk)
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/3.1@9855 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | doc/algorithms.tex | 4 | ||||
-rw-r--r-- | src/zeta.c | 31 | ||||
-rw-r--r-- | tests/tzeta.c | 22 |
3 files changed, 47 insertions, 10 deletions
diff --git a/doc/algorithms.tex b/doc/algorithms.tex index 3b987b0ad..208dec1e1 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -3097,7 +3097,9 @@ The algorithm for the Riemann Zeta function is due to Jean-Luc R\'emy and Sapphorain P\'etermann \cite{PeRe06,PeRe07}. For $s < 1/2$ we use the functional equation \[ \zeta(s) = 2^s \pi^{s-1} \sin\left(\frac{\pi s}{2}\right) \Gamma(1-s) - \zeta(1-s). \] + \zeta(1-s) \] +(in that case, one should take care of cancellation when $\sin(\pi s/2)$ is +small, i.e., when $s$ is near an even negative integer). For $s \geq 1/2$ we use the Euler-MacLaurin summation formula, applied to the real function $f(x) = x^{-s}$ for $s > 1$: \[ \zeta(s) = \sum_{k=1}^{N-1} \frac{1}{k^s} + \frac{1}{2N^s} diff --git a/src/zeta.c b/src/zeta.c index f29e5d0f5..e7042d0a6 100644 --- a/src/zeta.c +++ b/src/zeta.c @@ -377,8 +377,8 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) } } - /* Check for case s= 1 before changing the exponent range */ - if (mpfr_cmp (s, __gmpfr_one) ==0) + /* Check for case s=1 before changing the exponent range */ + if (mpfr_cmp (s, __gmpfr_one) == 0) { MPFR_SET_INF (z); MPFR_SET_POS (z); @@ -420,7 +420,7 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) MPFR_ZIV_INIT (loop, prec1); for (;;) { - mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN);/* s1 = 1-s */ + mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN); /* s1 = 1-s */ mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */ mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */ if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k, @@ -432,17 +432,32 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) break; } mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */ - mpfr_const_pi (p, MPFR_RNDD); - mpfr_mul (y, s, p, MPFR_RNDN); - mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* s*Pi/2 */ - mpfr_sin (y, y, MPFR_RNDN); /* sin(Pi*s/2) */ - mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); + + mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */ + + /* multiply z_pre by 2^s*Pi^(s-1) where p=Pi, s1=1-s */ mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */ mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */ mpfr_pow (y, y, s1, MPFR_RNDN); /* (2*Pi)^(s-1) */ mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN); + /* multiply z_pre by sin(Pi*s/2) */ + mpfr_mul (y, s, p, MPFR_RNDN); + mpfr_div_2ui (p, y, 1, MPFR_RNDN); /* p = s*Pi/2 */ + mpfr_sin (y, p, MPFR_RNDN); /* y = sin(Pi*s/2) */ + if (MPFR_GET_EXP(y) < 0) /* take account of cancellation in sin(p) */ + { + mpfr_t t; + mpfr_init2 (t, prec1 - MPFR_GET_EXP(y)); + mpfr_const_pi (t, MPFR_RNDD); + mpfr_mul (t, s, t, MPFR_RNDN); + mpfr_div_2ui (t, t, 1, MPFR_RNDN); + mpfr_sin (y, t, MPFR_RNDN); + mpfr_clear (t); + } + mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); + if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz, rnd_mode))) break; diff --git a/tests/tzeta.c b/tests/tzeta.c index 7bdf65cdd..ab0665983 100644 --- a/tests/tzeta.c +++ b/tests/tzeta.c @@ -243,7 +243,6 @@ main (int argc, char *argv[]) mpfr_set_str_binary (s, "1.10010e4"); mpfr_zeta (z, s, MPFR_RNDZ); - mpfr_set_prec (s, 53); mpfr_set_prec (y, 53); mpfr_set_prec (z, 53); @@ -394,6 +393,27 @@ main (int argc, char *argv[]) mpfr_nextabove (s); MPFR_ASSERTN (mpfr_equal_p (z, s) && inex > 0); + /* bug reported by Fredrik Johansson on 19 Jan 2016 */ + mpfr_set_prec (s, 536); + mpfr_set_ui_2exp (s, 1, -424, MPFR_RNDN); + mpfr_sub_ui (s, s, 128, MPFR_RNDN); /* -128 + 2^(-424) */ + for (prec = 6; prec <= 536; prec += 8) /* should go through 318 */ + { + mpfr_set_prec (z, prec); + mpfr_zeta (z, s, MPFR_RNDD); + mpfr_set_prec (y, prec + 10); + mpfr_zeta (y, s, MPFR_RNDD); + mpfr_prec_round (y, prec, MPFR_RNDD); + if (! mpfr_equal_p (z, y)) + { + printf ("mpfr_zeta fails near -128 for inprec=%lu outprec=%lu\n", + (unsigned long) mpfr_get_prec (s), (unsigned long) prec); + printf ("expected "); mpfr_dump (y); + printf ("got "); mpfr_dump (z); + exit (1); + } + } + mpfr_clear (s); mpfr_clear (y); mpfr_clear (z); |