[VA] SRC #010 - Pi Day 2022 Special

[Albert Chan]

(03-22-2022 12:14 AM)Valentin Albillo Wrote:  

(03-18-2022 04:22 PM)Albert Chan Wrote:  Lets recover true PN, and compare errors of products vs exp(sum of logs) [...] Note that ln(C) is odd function. Rewrite ln(C) as polynomial of 1/N, we have:

\(\displaystyle
\ln(C) = \frac{1}{N}
+ \frac{5/9}{N^3}
+ \frac{13/45}{N^5}
+ \frac{127/315}{N^7}
- \frac{89/135}{N^9}
\;+\; ... \) [...]

I must point out that this formal series of correction factors is asymptotic and divergent, i.e., its coefficients might be small and even decreasing for a while but eventually they grow bigger and bigger, both numerators and denominators, and thus can't be used to obtain arbitrary precision, as I explained in another case in post #27 of my Short & Sweet Math Challenge #24.
...

To settle down the question, if someone with access to Mathematica or some other arbitrary-precision software can compute the product for N=100,000 using 100 digits, say, or as many as necessary to ensure full 34 correct digits or more, and post here the resulting value I'd appreciate it. Thanks in advance.

>>> from mpmath import *
>>> mp.dps = 100
>>> pn = lambda n: exp(nsum(lambda x: 1+log1p(-1/(x*x))*x*x,[2,n]) + 1.5)
>>> n = mpf(100000)
>>> N = 2*n+1
>>> x = pn(n)
>>> print x
3.141608361513791562872866895754895278060325823725833279147116393910631517290786764227775828378244404

It does matched my 34-digits "true" PN.

ln(C) correction (terms upto 1/N^9) seems safe to use.

>>> err = lambda c: float(pi - x * exp(-c))
>>> err(13/(45*N**5) + 5/(9*N**3) + 1/N)
-9.8950471946808673e-38
>>> err(127/(315*N**7) + 13/(45*N**5) + 5/(9*N**3) + 1/N)
4.0449821226917704e-48
>>> err(-89/(135*N**9) + 127/(315*N**7) + 13/(45*N**5) + 5/(9*N**3) + 1/N)
-1.229817502771026e-57