# 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±2^-B)^2 => 1±2^(1-b)+2^(-2b) ≈ 1±2^(1-b)$?

— Sander Huisman 2009-12-24