Source of error in K₀ calculation found

Remco Bloemen

2009-12-23, last updated 2014-02-27

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

ζ(2n)=fn|B2n| \zeta(2n) = f_n \left\vert B_{2n} \right\vert

where

fn=(2π)2n2(2n)! 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:

fn=4π22n(2n1)fn1 f_n = \frac{4 \pi^2}{2n (2n - 1)} f_{n-1}

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

fn=(2n+1)(2n+2)4π2fn+1 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π24 \pi^2.

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

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

Such that the product is:

4π2fn(1±2B)2=4π2fn(1±2B+1) 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±2a)×B(1±2b)A(1\pm2^{-a}) \times B(1\pm2^{-b})

=AB(1±2a±2b±2ab)=AB(1\pm2^{-a}\pm2^{-b}\pm2^{-a-b})

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

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

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

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

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

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

So we lost one bit!

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

— Remco Bloemen 2009-12-24

To comment on my own post, this (2*n+1)*(2*n+2)(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±2B)2=>1±2(1b)+2(2b)1±2(1b)(1±2^-B)^2 => 1±2^(1-b)+2^(-2b) ≈ 1±2^(1-b)?

— Sander Huisman 2009-12-24