[VA] SRC #010 - Pi Day 2022 Special

[Albert Chan]

It is not hard to see why for log sums, we have error O(√(n)

1 + log1p(-1/k^2) * k^2
= 1 - (k^-2 + k^-4/3 + k^-6/4 + ...) * k^2
= 1 - (1 + k^-2/3 + k^-4/4 + ...)
= -(k^-2/3 + k^-4/4 + ...)

Because of 1 in front, term errors are in orders of machine epsilon.
Worst case, we have errors of O(n).

But, because errors spread-out somewhat randomly, we have O(√n)

You might try sum terms from index of 2 to n, instead of in reverse.
I would guess you would produce similar sized error for PN

--

For products of factors, (1-1/k^2)^(k^2):

We expected base have errors, also in order of machine epsilon.
However, errors are not random, but clustered when k is huge.
Example, for 10-digits calculator, this is the rounded base.

b(k) = 1-1/k^2

b(99999) = 0.99999 99998 99998, rounded up
b(99998) = 0.99999 99998 99996, rounded up
...
b(82000) = 0.99999 99998 51279, *still* rounded up

(1+ε)^(n^2) = 1 + n^2 ε

Product of n-1 terms, we expected worst case errors of O(n^3)
Of course, errors are not totally skewed, we expected O(n^2+)

From previous post:
PN(n=1e5) errors = 15,684,238,090 ULP ≈ 1e5 ^ 2.04

(03-14-2022 08:03 PM)Valentin Albillo Wrote:  
PN(N) = 3.14159 42243 85727 33446 22511 05879 403 ( 7 correct digits save 2 ulp )

Using ln(C) correction, "true" PN = 3.14159 42243 85727 33456 11796 83910 689

PN(n=1e6) errors = 98,928,578,031,286 ULP ≈ 1e6 ^ 2.33