Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B]

[Albert Chan]

(10-31-2021 03:40 PM)Albert Chan Wrote:  Let s = (-1)^n, x=n*n+n+1, T's = triangular number

ζ(2) ≈ (2 - 2/2^2 + 2/3^2 - ... + s*2/n^2)

       - s/(x - 1/(x+4*T1 - 2^4/(x+4*T2 - 3^4/(x+4*T3 - 4^4/(x+4*T4 - ...

Alternating sum CF correction, compared with non-alternating sum version, is lousy.
The only reason RHS also have 2.0899 digits per iteration is because error to fix is shrinking fast.

lua> z2 = pi^2/6
lua> z2-(1+1/4) , z2-(1-1/4)*2
0.3949340668482264      0.1449340668482264
lua> z2-(1+1/4+1/9) , z2-(1-1/4+1/9)*2
0.28382295573711525   -0.07728815537399591
lua> z2-(1+1/4+1/9+1/16) , z2-(1-1/4+1/9-1/16)*2
0.22132295573711525    0.047711844626004085

But, this may actually be useful !

We can spend more time summing the alternating series, less time doing CF corrections.
In the end, we may come out ahead.

As an experiment, we cut down CF terms in half. (p pairs alternating series, do p CF terms)
Bonus: summing alternating series pairwise also made final sum more accurate.

Code:

OPTION ARITHMETIC DECIMAL_HIGH
10 INPUT  PROMPT "p (pairs) = ":p
   LET x = (6*p+4)*p+1 ! half CF terms
   LET a = 0
   LET b = 0
   FOR i = p TO 2 STEP -1
      let t = i*i
      LET k = i - 0.25
      LET a = a + k / ((t-k)*t)
      let b = t*t / (x-b)
      LET x = x - 4*i
   NEXT i   
   LET z2 = 1/(x-4-1/(x-b)) + a*0.5 + 1.5
   PRINT "Accurate digits ="; 1-LOG10(ABS(PI*PI/6-z2))
   GOTO 10
END

p (pairs) = 66
Accurate digits = 204.17629078832496
p (pairs) = 68
Accurate digits = 210.30106011838216
p (pairs) = 266
Accurate digits = 816.6266359578656
p (pairs) = 268
Accurate digits = 822.75106992654514

For p=326, it almost reached 1000 digits full precision (1 ULP error, last digit = 8 instead of 9)

Because of asymmetry of terms, 2p (pairs) here, we expect cost about n = 3p, for previous version.
2*(2p) + (2p) = 6p = (3p) + (3p)

Using previous standard, for gain in digits per "iteration"

p = 66 to 68:     (210.3011 - 204.1763) / 3 = 2.0416
p = 266 to 268: (822.7511 - 816.6266) / 3 = 2.0415

It seems convergence is slightly worse. But, timings suggested otherwise.
For almost full precision case (p = 326) vs previous version (n=478), this version is faster, by quite a bit:

time(old version, n=478) / time(new version, p=326) = (0.22 sec) / (0.14 sec) ≈ 1.57

I cannot explain this ...