summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/zeta.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/zeta.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/zeta.c31
1 files changed, 24 insertions, 7 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/zeta.c b/Build/source/libs/mpfr/mpfr-src/src/zeta.c
index dc7e12b2501..ec2d88de7ec 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/zeta.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/zeta.c
@@ -1,6 +1,6 @@
/* mpfr_zeta -- compute the Riemann Zeta function
-Copyright 2003-2019 Free Software Foundation, Inc.
+Copyright 2003-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -82,16 +82,33 @@ mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
/* Input: p - an integer
Output: fills tc[1..p], tc[i] = bernoulli(2i)/(2i)!
tc[1]=1/12, tc[2]=-1/720, tc[3]=1/30240, ...
+ Assumes all the tc[i] have the same precision.
+
+ Uses the recurrence (4.60) from the book "Modern Computer Arithmetic"
+ by Brent and Zimmermann for C_k = bernoulli(2k)/(2k)!:
+ sum(C_k/(2k+1-2j)!/4^(k-j), j=0..k) = 1/(2k)!/4^k
+ If we put together the terms involving C_0 and C_1 we get:
+ sum(D_k/(2k+1-2j)!/4^(k-j), j=1..k) = 0
+ with D_1 = C_0/4/(2k+1)/(2k)+C_1-1/(2k)/4=(k-1)/(12k+6),
+ and D_k = C_k for k >= 2.
+
+ FIXME: we have C_k = (-1)^(k-1) 2/(2pi)^(2k) * zeta(2k),
+ see for example formula (4.65) from the above book,
+ thus since |zeta(2k)-1| < 2^(1-2k) for k >= 2, we have:
+ |C_k - E_k| < E_k * 2^(1-2k) for k >= 2 and E_k := (-1)^(k-1) 2/(2pi)^(2k).
+ Then if 2k-1 >= prec we can evaluate E_k instead, which only requires one
+ multiplication per term, instead of O(k) small divisions.
*/
static void
mpfr_zeta_c (int p, mpfr_t *tc)
{
- mpfr_t d;
- int k, l;
-
if (p > 0)
{
- mpfr_init2 (d, MPFR_PREC (tc[1]));
+ mpfr_t d;
+ int k, l;
+ mpfr_prec_t prec = MPFR_PREC (tc[1]);
+
+ mpfr_init2 (d, prec);
mpfr_div_ui (tc[1], __gmpfr_one, 12, MPFR_RNDN);
for (k = 2; k <= p; k++)
{
@@ -306,7 +323,7 @@ compute_add (mpfr_srcptr s, mpfr_prec_t precz)
/* since 1/eps = 2^(precz+14), if EXP(sd) >= precz+14, then
sd >= 1/2*2^(precz+14) thus 2*sd >= 2^(precz+14) >= 1/eps */
if (mpfr_get_exp (t) >= precz + 14)
- mpfr_mul_2exp (t, t, 1, MPFR_RNDU);
+ mpfr_mul_2ui (t, t, 1, MPFR_RNDU);
else
mpfr_set_ui_2exp (t, 1, precz + 14, MPFR_RNDU);
/* now t = max(1/eps,2*sd) */
@@ -318,7 +335,7 @@ compute_add (mpfr_srcptr s, mpfr_prec_t precz)
else
mpfr_set (t, m1, MPFR_RNDU);
/* now t = max(8,m1) */
- mpfr_div_2exp (t, t, precz + 14, MPFR_RNDU); /* eps*max(8,m1) */
+ mpfr_div_2ui (t, t, precz + 14, MPFR_RNDU); /* eps*max(8,m1) */
mpfr_add_ui (t, t, 1, MPFR_RNDU); /* 1+eps*max(8,m1) */
mpfr_mul (t, t, u, MPFR_RNDU); /* t = c */
mpfr_add_ui (u, m1, 13, MPFR_RNDU); /* 13+m1 */