2009-12-13 (updated 2014-02-27)

Source of error in K₀ calculation found

The Staudt-Clausen theorem requires a conversion factor between the ζ-function and the Bernoulli numbers such that

$\zeta(2n) = f_n \left\vert B_{2n} \right\vert$

where

$f_n = \frac{\left( 2 \pi \right)^{2n}}{2 (2n)!}$

Since the exponentiation and factorial calculation is pretty demanding at high precision and high n it is done by iteratively updating the previous value:

$f_n = \frac{4 \pi^2}{2n (2n - 1)} f_{n-1}$

Or, in the more recent versions where n is decreasing:

$f_n = \frac{(2n +1) (2n + 2)}{4 \pi^2} f_{n+1}$

The implementation was as follows:

// factor *= (2n + 1) * (2n + 2) / (2π)²
mpfr_div(factor, factor, twopi2, GMP_RNDN);
mpfr_mul_ui(factor, factor, (2*n+1)*(2*n+2), GMP_RNDN);

Where twopi2 is set to $4 \pi^2$ .

Now, both factor and twopi2 are accurate up to B bits, that is:

factor $= f_n (1 \pm 2^{-B})$ twopi2 $= 4\pi^2 (1 \pm 2^{-B})$

Such that the product is:

$4 \pi^2 f_n (1 \pm 2^{-B})^2 = 4 \pi^2 f_n (1 \pm 2^{-B+1})$

So we lost one bit of accuracy! To compensate for this, I add N extra bits to the accuracy of factor and twopi2, where N is the number of iterations required. This is in fact a bit of overkill, since twopi2 stays B bits accurate and the above derivation assumed factor and twopi2 to have the same error.

The program was then tested with several different parameters and run at 1, 10, 100, 300 and 500 thousand digits of accuracy to verify that the bug is now completely gone.

The impact of this change is about a 3% increase in the number of bits required in factor and twopi2, there was a very slight performance penalty, but this was easily offset by discovering that I could apply this improvement in a few more places.

In conclusion, my Khinchin calculator is back in the game! At this moment I’m running two calculations of one million and one point one million digits on the Qubis server. If all goes well I should have a verified set of one million K₀ digits before the year is over!

Check! Minor detail: The last term after the expansion is always positive ;)

— Sander Huisman 2009-12-24

@SH: I started out with the general case of a number A with a bits accuracy and a number B with b bits accuracy:

$A(1\pm2^{-a}) \times B(1\pm2^{-b})$

$=AB(1\pm2^{-a}\pm2^{-b}\pm2^{-a-b})$

The $2^{-a-b}$ term is insignificant so can be dropped.

$=AB(1\pm2^{-a}\pm2^{-b})$

In case both number are equally accurate (worst case), then $a=b$ and:

$=AB(1\pm2^{-a}\pm2^{-a})$

$=AB(1\pm2\cdot2^{-a})$

$=AB(1\pm2^{-a+1})$

So we lost one bit!

PS: With $(1\pm2^{-n})$ I of course mean a number in the range $[1-2^{-n} {,} 1+2^{-n})$

— Remco Bloemen 2009-12-24

To comment on my own post, this $(2*n+1)*(2*n+2)$ factor is exactly the reason why the program can’t be used in 32 bits, this number gets to about 40 bits when you calculate a million digits.

— Remco Bloemen 2009-12-24

Looking good! Was this on purpose: $(1\pm2^-B)^2 => 1\pm2^(1-b)+2^(-2b) \approx 1\pm2^(1-b)$ ?

— Sander Huisman 2009-12-24