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

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

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] (/thread-17613.html)

Pages: 1 2

RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-04-2021 08:35 PM

(10-27-2021 05:12 PM)Albert Chan Wrote: (*) this is how initial b is estimated, by looking ahead.

XCAS> b2 := (N+3/2)^4 / (c*N*(N+2))
XCAS> b1 := (N+1/2)^4 / ((10-c)*N*(N+1) + b2)

XCAS> simplify(partfrac(b1(c=8))[1]) → (2312*N^2+1904*N+480)/4913
XCAS> simplify(partfrac(b1(c=2))[1]) → (578*N^2+476*N+120)/4913

Let's get generalized CF convergents from the top down, instead of bottom up.

I set c=8.0, to reduce expensive symbolic calculations (setting c=2.0 will get ¼th as big)
Also, I use x = n+0.5 instead of N = n+0.5, to match Decimal Basic code.
Warning: Decimal Basic is case-insensitive, N and n are the same variable.

XCas> nextv(v, cf) := [v[1], normal(cf*v)] // 2nd row = next convergent

XCas> v := identity(2) // initial b convergents, start with 1/0, b0/1 = 0/1
XCas> c, n, x2 := 8.0, x-0.5, x*x

XCas> v := nextv(v, [(n+=1)^4, (c:=10-c)*(x2+=x)]):; e2r(quo(v[1]))

[0.5, 0.5, 0.25] // b = horner([0.5, 0.5, 0.25], n+0.5)

Repeat the last command, we have:

[0.470588235294, 0.387543252595, 0.0976999796458]
[0.472222222222, 0.401234567901, 0.141117969822]
[0.472131147541, 0.399892502016, 0.133170617805]
...
[0.472135955, 0.4, 0.13416407865] // converged 12 digits

Guessing that coefficients are somehow related to ϕ, above is likely this:

lua> r = sqrt(5) -- = 2*phi - 1
lua> 2*(r-2), 0.4, 0.06*r
0.4721359549995796 0.4 0.1341640786499874

Let's test this, for Decimal Basic version of zeta2(n)

Code:

OPTION ARITHMETIC DECIMAL_HIGH

10 INPUT  PROMPT "n = ":n

   LET c = MOD(n,2)*6+2   

   LET x = n+0.5

   LET x2 = x*x

   LET a = 0

   LET b = SQR(5)

   LET b = (4*(b-2)*x2+0.8*x+0.12*b)/(10-c)

   FOR i = n TO 2 STEP -1

      LET t = i*i

      LET a = 1/t + a

      LET b = t*t / (c*x2+b)

      LET c = 10-c

      LET x2 = x2-x

   NEXT i   

   LET z2 = 1/(x+1/(c*x2+b)) + a + 1

   PRINT "Accurate digits ="; 1-LOG10(ABS(PI*PI/6-z2))

   GOTO 10

END

n = 100
Accurate digits = 217.89547073326986
n = 101
Accurate digits = 219.99818460987848
n = 102
Accurate digits = 222.10077322982628
n = 400
Accurate digits = 846.6549411985644
n = 401
Accurate digits = 848.74806270861004
n = 402
Accurate digits = 850.84117615610204

For n=474, it reached 1000 digits full precision.

RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-04-2021 10:42 PM

If we add a 1 CF correction term to initial b, result is even better.

LET b = SQR(5)
LET b = (4*(b-2)*x2+0.8*x+0.12*b - 3.2/(50*x+23*b))/(10-c)

n = 100
Accurate digits = 223.10474793591716
n = 101
Accurate digits = 225.21595297592808
n = 102
Accurate digits = 227.3269505519926
n = 400
Accurate digits = 853.05677365480504
n = 401
Accurate digits = 855.15205424416744
n = 402
Accurate digits = 857.24732141737914

For n=471, it reached 1000 digits full precision.

RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-05-2021 03:55 PM

(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 ...