Use Chudnovsky rational approximation of pi^2 for zetasEven

[Bodigrim]

To evaluate Riemann zeta function on even arguments (zetasEven) we use formula: The main issue here is that for reasonably big n we divide one huge number by another one, meaning a catastrophic loss of precision. For instance,

> approximateValue (zetasEven !! 25) :: Double
0.9999999999999996
> approximateValue (zetasEven !! 999) :: Double
NaN

which does not make any sense, because the zeta function is always finite and > 1. Currently we fix this behaviour by flooring values of zetasEven by 1.

This discrepancy hurts evaluation of Riemann zeta function on odd arguments and propagates further to numeric errors in Dirichlet beta function.

---

This is a rough idea, I have not given it much thought.

Chudnovsky algorithm gives us a rapidly converging series of rational approximations to pi^2. We can use it to generate a rational approximation of zetasEven !! n and consequently calculate not only zetaEven !! n itself (which is any way close to 1), but also zetasEven !! n - 1 using usual Double with a fairly good relative precision.

[b-mehta]

Surely the issue here is in approximateValue from ExactPi, rather than in zetasEven? Specifically, zetasEven returns an exact value, and the loss of precision occurs when converting the exact value into a Double.

That said, I've implemented Chudnovsky's algorithm in my branch.

[rockbmb]

:+1: for loss of precision in Dirichlet beta. The code introduced by #120 works, but not very well. Error propagation from zetasOdd and zetasEven is quite strong.

[Bodigrim]

Cool!

@b-mehta do you know any rapidly converging rational series for pi? If yes, then we could possibly combine it with Chudnovsky algorithm for pi^2 and squeeze them together into exact-pi instead of Data.ExactPi. rationalApproximations.

[b-mehta]

@Bodigrim I don't know of any off the top of my head, but I'll research. In the meantime, I've just commited to b-mehta/arithmoi/chudnovsky for a solution which should work if the answer turns out to be no: I've made a datatype ExactPiSq as a simple analogue to ExactPi along with a new approximateValue which uses Chudnovsky's algorithm to get the best approximation for a given Fractional instance.

The test value in the docs for zeta(50) used to differ from the true value after only around 10 non-zero places, but now agrees to a Fixed Prec50 at least. In addition, we can now get values for zetasEven !! n-1 as you mentioned (cf WolframAlpha):

> import Data.Number.Fixed
> let Exact z q = zetasEven !! 50
> approximateValue (Exact z q) :: Fixed Prec50
1.00000000000000000000000000000078886090522101180736
> getRationalLimit [q * p^z - 1 | p <- chudnovsky]
7.888609052210118e-31

[b-mehta]

One possible solution would be to use Chudnovsky's algorithm directly, and use a rational approximation to the square root using continued fraction convergents which is easy to calculate and should be fairly easy to accelerate also.

[Bodigrim]

Looks awesome!

The approach with rational approximation to the square root sounds reasonable.

I suggest passing comparison function a -> a -> Bool to getRationalLimit as an argument and do not rely on Eq instance. This is not only more flexible in terms of precision, but is also a cure against NaN /= NaN thing, which is a common reason for fix-point evaluations cycling forever.

[b-mehta]

Sounds good, I'd like to do some numerical testing myself before making changes in exact-pi, but the modifications to getRationalLimit sound good. Would it be better to move the Chudnovsky work here to that repo instead?

[Bodigrim]

Yes, IMHO Chudnovsky algorithm itself fits into exact-pi better.

[b-mehta]

I'm hesitant to replace approximateValue since it uses the target type's pi value which may be preferable in some cases.

But, for Chudnovsky to have an impact on the issues in zetasEven' and Dirichlet beta, it seems preferable to have something of the type ExactPi -> a for Fractional/Floating a. This was my original aim for getRationalLimit but as you point out, NaN would create issues there. I'm unsure how to resolve this, since it is desirable to find the limit. Any suggestions?

[Bodigrim]

Sure, we'd rather replace rationalApproximations than approximateValue. Let us return to an impact on zetas after that.

[rockbmb]

Unfortunate that a relatively simple issue is taking so long to get merged 💤💤. The package's maintainer seems to be very busy with other things. We cannot fault him, open source takes time and effort that isn't always available.

@Bodigrim in a few weeks (2019), if nothing happens, will you release your fork of exact-pi? Should the maintainer return, you can just merge your work into his repository and close yours.

[Bodigrim]

Yes, if nothing happens I will proceed with my fork.

[Bodigrim]

I managed to upload exact-pi-0.5.0.0 with @b-mehta's improvements to rationalApproximations.

[rockbmb]

The issue has been merged, I'm guessing we can proceed with a PR to improve zetas/betas.

[Bodigrim]

I patched arithmoi to benefit from rationalApproximations. It has an immediate effect even for Double. For example, betas 1e-15 !! 31 gives a correct value of 0.9999999999999983 instead of 0.9999999999999976 previously.

@b-mehta thanks for valuable contribution!

[rockbmb]

@Bodigrim That's great. Thanks for doing the change yourself.