 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))) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → (2312*N^2+1904*N+480)/4913 XCAS> simplify(partfrac(b1(c=2))) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → (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, 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)) [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 ...