+- 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
Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Gerson W. Barbosa - 10-23-2021 02:49 PM
FWIW,
A new method for the fast evaluation of ζ(2) using the definition as a basis
The series, 1 + 1/4 + 1/9 + 1/25 + 1/36 + 1/49…, converges very slowly to the exact result, π^2/6. In order to obtain n correct digits the series should be evaluated up to the (10^n)th term. However, the addition of n+1 terms of a simple continued fraction after the evaluation of the first n terms of the series will significantly speed up the rate of convergence, yielding slightly more than 2n correct digits.
For example,
for n = 3,
1+1/4+1/9+1/((3+1/2)+1/(12*(3+1/2)+16/(5*(3+1/2)+81/(28*(3+1/2))))) = 55783/33912 = 1.6449339
The coefficients of the denominators of the continued fraction, 12, 5, 28, 9, 44, 13…, obey the formula k(i) = (5 - 3*(-1)^i)*(i + 1/2). The numerators, 1, 16, 81, 256, 625, 1296…, are quite obvious.
HP-42S/Free42 program:
Code:
00 { 65-Byte Prgm }
01▸LBL "z"
02 0.5
03 +
04 STO 01
05 IP
06 0
07 STO 02
08▸LBL 00
09 RCL ST Y
10 X↑2
11 1/X
12 +
13 -1
14 RCL ST Z
15 Y↑X
16 +/-
17 3
18 ×
19 5
20 +
21 0.5
22 RCL+ ST T
23 ×
24 RCL× 01
25 RCL+ 02
26 R↑
27 X↑2
28 X↑2
29 X<>Y
30 ÷
31 STO 02
32 R↓
33 DSE ST Y
34 GTO 00
35 RCL 02
36 RCL+ 01
37 1/X
38 +
39 ENDn = 12 on the HP-42S and n = 16 on Free42 will suffice for 12 and 34 correct digits, respectively.
6 XEQ “z” → 1.64493406685
16 XEQ “z” →
1.644934066848226436472415166646025
HP-71B BASIC program;
Code:
10 S=0
15 C=0
20 INPUT N
25 K=N+.5
30 A=3-6*MOD(N,2)
35 FOR I=N TO 1 STEP -1
40 S=S+1/(I*I)
45 C=I^4/(C+K*(5-A)*(I+.5))
50 A=-A
55 NEXT I
60 DISP S+1/(K+C)RUN
? 6
1.64493406685
Interested readers are invited to provide
- optimized versions of the given programs;
- versions for other calculator, such as the HP-41;
- a proof (I don’t have any – this is the result of a Friday afternoon work only, which until minutes ago I thought to be a Saturday afternoon. Still looking like Sunday morning to me).
Pointing out typos and mistakes, either math or grammar related, are welcome.
Gerson.
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 10-25-2021 01:29 PM
Code:
function zeta2(n)
local a, b, c, N, t = 0, 0, (n%2)*6+2, n+0.5
local N2 = N*N
for i = n, 2, -1 do
t = i*i
a = a + 1/t
b = t*t/(b + c*N2)
c, N2 = 10-c, N2-N
end
return 1/(N+1/(b+c*N2)) + a + 1
endThis lua code loop from n down to 2 instead of n to 1.
Note: for n = 0, this return 1.4 instead of 2.
lua> for n=0,8 do print(n, zeta2(n)) end
0 1.4
1 1.6428571428571428
2 1.644949494949495
3 1.644933946685539
4 1.6449340678010402
5 1.6449340668406054
6 1.6449340668482877
7 1.644934066848226
8 1.6449340668482264
Convergence speed very impressive, n=8 converged to float(pi^2/6)
How did you derive the continued fraction form ?
I was reading Why pi^2 so close to 10 ? With pi^2 ≈ 10, pi^2/6 ≈ 5/3
5/3 = 1 + sum(4/(4*k^2-1), k=2..inf) // terms telescopically cancel, leaving only 5/3
ζ(2) = 1 + sum(1/k^2, k=2..inf)
5/3 - ζ(2) = sum(4/(4*k^2-1) - 1/k^2, k=2..inf)
ζ(2) = 5/3 - sum(1/(k^2*(4*k^2-1)), k=2..inf)
We might as well pull 1 term out. (pi^2/6 ≈ 9.9/6 = 1.65)
ζ(2) = 1.65 - sum(1/(k^2*(4*k^2-1)), k=3..inf)
With denominator of O(k^4), this converge much faster than ζ(2) "definition"
How can we accelerate this with continued fraction ?
---
This converge even faster, A cute sum of Ramanujan
zeta2(n) := (10 - sum(1/(k*(k+1))^3, k=1..n)) / 6
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Gerson W. Barbosa - 10-26-2021 02:12 AM
(10-25-2021 01:29 PM)Albert Chan Wrote:Code:
function zeta2(n)
local a, b, c, N, t = 0, 0, (n%2)*6+2, n+0.5
local N2 = N*N
for i = n, 2, -1 do
t = i*i
a = a + 1/t
b = t*t/(b + c*N2)
c, N2 = 10-c, N2-N
end
return 1/(N+1/(b+c*N2)) + a + 1
end
This lua code loop from n down to 2 instead of n to 1.
Note: for n = 0, this return 1.4 instead of 2.
lua> for n=0,8 do print(n, zeta2(n)) end
0 1.4
1 1.6428571428571428
2 1.644949494949495
3 1.644933946685539
4 1.6449340678010402
5 1.6449340668406054
6 1.6449340668482877
7 1.644934066848226
8 1.6449340668482264
Convergence speed very impressive, n=8 converged to float(pi^2/6)
That gives linear convergence ( 25/12 digits per iteration ), as you can see by the Decimal Basic code and output:
Code:
OPTION ARITHMETIC DECIMAL_HIGH
INPUT PROMPT "number of decimal digits: ":nd
LET tm = TIME
LET n = CEIL(12*nd/25)
PRINT "n = ";n
LET p = n*n
LET q = p*p
LET a = 0
LET b = 0
LET c = 6*MOD(n,2) + 2
LET d = q - (n - 1)^4
LET e = 2*n - 1
LET f = 12*p
LET g = 24*n - 12
LET m = n + 0.5
LET m2 = m*m
FOR i = 2 TO n
LET a = a + 1/p
LET b = q/(b + c*m2)
LET c = 10 - c
LET m2 = m2 - m
LET p = p - e
LET e = e - 2
LET q = q - d
LET f = f - g
LET d = d - f - 2
LET g = g - 24
NEXT i
LET z2 = 1/(m + 1/(b + c*m2)) + a + 1
LET r = TIME - tm
LET r$ = STR$(z2 - INT(z2))
PRINT
PRINT STR$(INT(z2));
PRINT r$(0:1);
FOR i = 2 TO nd + 1
PRINT r$(i:i);
IF MOD((i - 1),10) = 0 THEN PRINT " ";
IF MOD((i - 1),50) = 0 THEN
PRINT USING " (%%%%)": i - 1
PRINT REPEAT$(" ",1 + LEN(STR$(INT(z2))));
END IF
NEXT i
IF MOD (i - 2,50) <> 0 OR nd = 0 THEN PRINT
PRINT
PRINT "Runtime: ";
PRINT USING "0.##": r;
PRINT " seconds"
END
Output:
number of decimal digits: 999
n = 480
1.6449340668 4822643647 2415166646 0251892189 4990120679 (0050)
8437735558 2293700074 7040320087 3833628900 6197587053 (0100)
0400431896 2337190679 6287246870 0500778793 5102946330 (0150)
8662768317 3330936776 2605095251 0068721400 5479681155 (0200)
8794890360 8232777619 1984075645 5876963235 6367097100 (0250)
9694890208 5932008051 6364788783 3884604444 5184059825 (0300)
1452506833 8763142276 5879392958 8063204472 1979084773 (0350)
4091059020 8378289549 2782638903 7976358334 3942045159 (0400)
1208180995 9345444877 4587965008 8088940870 1111634710 (0450)
6931614618 4288798154 8624483590 9183448757 3874283940 (0500)
8276028756 3214346010 0135766209 8204872069 0400073826 (0550)
6356030240 2284462963 0324566097 1719514277 2131595125 (0600)
5679986190 8719315439 5352410638 0440721421 3396547505 (0650)
8015872316 5839947624 3491422433 4836290488 7009665059 (0700)
8622630341 0959673655 2811371670 3269114987 8403435716 (0750)
1605776676 3330672527 3689423841 6640889536 2275954007 (0800)
7279474812 7102520498 3784332300 1716574481 0302860434 (0850)
9668847942 1672843359 7281997793 8100084665 6078053778 (0900)
2885947278 6259316186 6458829216 0658193859 2324153258 (0950)
0646178120 1884649777 6259849775 6060938460 605146769
Runtime: 0.33 secondsI used your optimizations (thanks!) and avoided all exponentiation-related multiplications inside the loop.
(10-25-2021 01:29 PM)Albert Chan Wrote: How did you derive the continued fraction form ?
Numerically. I computed the difference of the exact result and the partial sums for a few n and examined their inverses. For example:
10 -> 10.507917
20 -> 20.504062
90 -> 90.500921
That gives the first term of the continued fraction, 1/(n + 1/2)
By doing the calculation with 34 digits, it was possible to get the next three terms:
1/((n + 1/2) + 1/(12*(n + 1/2) + 16/(5*(n + 1/2) + 81/(28*(n+1/2) + ... ))))
The numerators are obvious, but the coefficients of the denominators don't appear to follow a clear pattern. Then I tried the solver:
solve 1 + 1/((1 + 1/2) + 1/(12*(1 + 1/2) + 16/(5*(1 + 1/2) + 81/(28*(1 + 1/2) + 256/(x*(1 + 1/2)))))) = pi^2/6 --> x ~ 13
But neither 13 nor 12 appeared to be correct. Then another try:
solve 134225/90000 + 1/((6 + 1/2) + 1/(12*(6 + 1/2) + 16/(5*(6 + 1/2) + 81/(28*(6 + 1/2) + 256/(x*(6 + 1/2)))))) = pi^2/6 --> x ~ 9
This appeared to be right. Next step was searching OEIS for a sequence starting with 1, 12, 5, 28, 9 -- and there it was:
https://oeis.org/A110185
0, 1, 12, 5, 28, 9, 44, 13, 60, 17, 76, 21, 92, 25, 108, 29...
"a(n) = (5+3(-1)^n)(2n-1)/2, with a(0)=0. Sum(a(i), i=0..n) = A085787(A042948(n)). - Bruno Berselli, Jan 20 2012"
Apparently no connection with the original problem, but it did work.
(10-25-2021 01:29 PM)Albert Chan Wrote: This converge even faster, A cute sum of Ramanujan
zeta2(n) := (10 - sum(1/(k*(k+1))^3, k=1..n)) / 6
Is that ok? 50 terms of the above series give only 10 correct digits:
1.644934066(94481120)
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 10-26-2021 09:47 AM
(10-26-2021 02:12 AM)Gerson W. Barbosa Wrote: Numerically. I computed the difference of the exact result and the partial sums for a few n and examined their inverses.
For example:
10 -> 10.507917
20 -> 20.504062
90 -> 90.500921
That gives the first term of the continued fraction, 1/(n + 1/2)
Interesting many CF correction formula involve the tag along +1/2
Quote:By doing the calculation with 34 digits, it was possible to get the next three terms:
1/((n + 1/2) + 1/(12*(n + 1/2) + 16/(5*(n + 1/2) + 81/(28*(n+1/2) + ... ))))
The numerators are obvious, but the coefficients of the denominators don't appear to follow a clear pattern ...
Impressive you can deduce denominator pattern with only 4 numbers: 1, 12, 5, 28, ...
Anther way to check units (dimensional analysis).
Numerator: 1, 16, 81, ... seems to involve unit of U^4
Denominator should have unit U^2, to have everything consistent.
But we only have unit U from (n+1/2), so coefficients somehow also have unit of U
To pull units off coefficients, we interpolate:
(1,1), (3,5) ⇒ y = 2*x-1 = 2*(x-0.5)
(2,12), (4,28) ⇒ y= 8*x-4 = 8*(x-0.5)
It would be better if we had more coefficients to confirm the pattern.
But, this suggested "unit free" coefficients are alternating 2, 8, 2, 8 ...
1/(2*0.5*(n+0.5) + 1/(8*1.5*(n+0.5) + 16/(2*2.5*(n+0.5) + 81/(8*3.5*(n+0.5) + ... ))))
Quote:(10-25-2021 01:29 PM)Albert Chan Wrote: This converge even faster, A cute sum of Ramanujan
zeta2(n) := (10 - sum(1/(k*(k+1))^3, k=1..n)) / 6
Is that ok? 50 terms of the above series give only 10 correct digits:
1.644934066(94481120)
Considering number of terms needed to sum 1/k^2, this is a huge improvement.
Maybe we can accelerate convergence with CF correction ?
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Ren - 10-26-2021 02:13 PM
I beg your kindness and pardon,
but to a Mathoklutz like me,
What does ζ mean?
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - floppy - 10-26-2021 03:04 PM
(10-26-2021 02:13 PM)Ren Wrote: I beg your kindness and pardon,an explanation here https://mathworld.wolfram.com/RiemannZetaFunctionZeta2.html
but to a Mathoklutz like me,
What does ζ mean?
This bring me to the idea, I must have a look in the relation to GAGM (quite new thing in math, from a russian; because there is i in it, too) https://indico-hlit.jinr.ru/event/187/contributions/1769/attachments/543/931/SAdlaj.pdf
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 10-26-2021 03:24 PM
(10-26-2021 02:13 PM)Ren Wrote: What does ζ mean?
(10-25-2021 01:29 PM)Albert Chan Wrote: ζ(2) = 1.65 - sum(1/(k^2*(4*k^2-1)), k=3..inf)
With denominator of O(k^4), this converge much faster than ζ(2) "definition"
How can we accelerate this with continued fraction ?
Only 1 term CF correction. It does help.
CAS> zeta2(n) := 1.65 - sum(1/(k^2*(4*k^2-1)), k=3..n)
CAS> corr2(n) := -1/horner([12,0,4.2,0],n+0.5)
CAS> zeta2(10) → 1.6450058264
CAS> Ans+corr2(10) → 1.64493406774
CAS> pi*pi/6. → 1.64493406685
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Gerson W. Barbosa - 10-26-2021 03:58 PM
(10-26-2021 02:13 PM)Ren Wrote: What does ζ mean?
Basically it is a Greek letter, zeta, used do denote the mathematical function zeta(x), defined as
ζ(x) = 1 + 1/2ˣ + 1/3ˣ + 1/4ˣ + 1/5ˣ + …
where x is a real number greater than zero
When x = 1, for instance,it is equivalent to the harmonic function, ζ(x) = H(x) = 1 + 1/2 + 1/3+ 1/4 + 1/5 + …
Later the domain was extended to all the real axis and eventually to all the complex plane, when things started to become really interesting.
https://www.britannica.com/science/Riemann-zeta-function
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 10-26-2021 08:28 PM
(10-26-2021 09:47 AM)Albert Chan Wrote: Interesting many CF correction formula involve the tag along +1/2
I get it now ! It is because correction is bounded by integral test.
∫(1/x^2, x=n .. inf) > Σ(1/k^2, k=n+1 .. inf) > ∫(1/x^2, x=n+1 .. inf)
Σ(1/k^2, k=n+1 .. inf) ≈ ∫(1/x^2, x=n+0.5 .. inf) = 1/(n+0.5)
Example, 2 terms + correction: ζ(2) ≈ 1 + 1/4 + 1/2.5 = 1.65
This explained why CF formula have form 1/((n+0.5) + ...)
(10-26-2021 03:24 PM)Albert Chan Wrote: CAS> zeta2(n) := 1.65 - sum(1/(k^2*(4*k^2-1)), k=3..n)
CAS> corr2(n) := -1/horner([12,0,4.2,0],n+0.5)
CAS> zeta2(10) → 1.6450058264
CAS> Ans+corr2(10) → 1.64493406774
CAS> pi*pi/6. → 1.64493406685
Above example, f(x) = 1/(x^2*(4*x^2-1) = -1/x^2 - 1/(x+1/2) + 1/(x-1/2)
∫(f(x), x)
= 1/x - (ln(x+1/2) - ln(x-1/2))
= 1/x - 2*atanh(1/(2x))
= 1/x - 2*(1/(2x) + 1/(2x)^3/3 + ...)
≈ -1/(12*x^3)
This explained corr2(n) denominator big term: 12*(n+0.5)^3
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Ren - 10-27-2021 01:32 PM
Thank you Albert Chan, Gerson W. Barbosa, and floppy for enlightening me!
Mr. Chan, I found the video most helpful (and overwhelming!) B^)
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 10-27-2021 05:12 PM
(10-25-2021 01:29 PM)Albert Chan Wrote: Convergence speed very impressive, n=8 converged to float(pi^2/6)
2nd version, by putting some reasonable number to b, before loop start (*)
Now, n=6 already converged to float(pi^2/6)
Code:
function zeta2(n)
local N, c, t = n+0.5, (n%2)*6+2
local N2 = N*N
local a, b = 0, (N2+N) * 0.93/(10-c)
for i = n, 2, -1 do
t = i*i
a = a + 1/t
b = t*t/(b + c*N2)
c, N2 = 10-c, N2-N
end
return 1/(N+1/(b+c*N2)) + a + 1
endlua> for n=0,8 do print(n, zeta2(n)) end
0 1.4539195555018722
1 1.6448912932530366
2 1.6449340922532778
3 1.6449340670471453
4 1.6449340668467034
5 1.6449340668482317
6 1.6449340668482264
7 1.6449340668482264
8 1.6449340668482264
(*) 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
2312/4913*2 ≈ 0.9412
If N is large, this should be the constant to use instead of 0.93
Update: Looking ahead a few more, 0.944 is about optimum
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 10-29-2021 02:16 PM
(10-26-2021 09:47 AM)Albert Chan Wrote: It would be better if we had more coefficients to confirm the pattern.
But, this suggested "unit free" coefficients are alternating 2, 8, 2, 8 ...
1/(2*0.5*(n+0.5) + 1/(8*1.5*(n+0.5) + 16/(2*2.5*(n+0.5) + 81/(8*3.5*(n+0.5) + ... ))))
We can also sum odd terms: sum(1/(2k-1)^2, k=1..inf) = pi^2/8 = 3/4*ζ(2)
ζ(2) = 4/3 * sum(1/(2k-1)^2, k=1..inf)
If we sum only n terms, correction term ≈ ∫(1/(2x-1)^2, x=n+1/2 .. inf) = 1/(4n)
Doing numerically for the continued fraction form, we get a very similar form as above.
1/(4*n + 1/(2*1.5*n + 16/(8*2.5*n + 81/(2*3.5*n + 256/(8*4.5*n + ... )))))
Example, estimate ζ(2) with n=3
lua> n = 3
lua> 1 + 1/3^2 + 1/5^2
1.1511111111111112
lua> 1/(4*n + 1/(2*1.5*n + 16/(8*2.5*n + 81/(2*3.5*n))))
0.08258933029740148
lua> 4/3 * (1.1511111111111112 + 0.08258933029740148)
1.6449339218780168
lua> pi*pi/6
1.6449340668482264
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 10-31-2021 03:40 PM
Another approach, getting ζ(2) with alternating series.
Σ(all terms) = Σ(even terms) * 4
Σ(odd terms) = Σ(all terms) - Σ(even terms) = Σ(even terms) * 3
Σ(odd terms) - Σ(even terms) = Σ(even terms) * 2
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 - ...
Note: Tk - Tk-1 = k*(k+1)/2 - (k-1)*k/2 = k
zeta2_alt(n) is computationally less expensive than zeta2(n)
Code:
function zeta2_alt(n)
local x, t = 3*n*(n+1)+1
local a, b = 0, 0
for i = n, 2, -1 do
t = i*i
a = 2/t - a
b = t*t / (x-b)
x = x - 4*i
end
return (1-n%2*2)/(x-4-1/(x-b)) - a + 2
endUpdate:
revert zeta2_alt(n) with b=0, to fairly compare against zeta2(n) (also, b=0)
zeta2_alt(n) is not as good, less accurate by about 2/5 of a decimal digit.
lua> err = function(z2) return abs(pi*pi/6-z2) end
lua> for n=1,6 do print(n, log10(err(zeta2_alt(n))/err(zeta2(n)))) end
1 0.3872388516915953
2 0.40466274755687004
3 0.4105758432666672
4 0.4132661850460007
5 0.41470949038470145
6 0.41521536217708266
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-01-2021 12:56 AM
We can get the CF correction formula, using Euler–Maclaurin formula (see #9)
XCas> C(k) := bernoulli(k)/k!; // Euler-Maclaurin formula coefs
XCas> f, a, b := 1/x^2, N, inf; // N = n+1
XCas> corr := int(f,x,a,b) → 1/N
XCas> corr += preval(f,a,b) * (-1/2) → 1/N+1/(N^2*2)
XCas> f := f' :; corr += preval(f,a,b) * C(k:=2)
XCas> f := f'':; corr += preval(f,a,b) * C(k+=2)
Run last line a few times, we have this corr, as polynomial of 1/N:
XCas> e2r(corr(N=1/x))
[-691/2730, 0, 5/66, 0, -1/30, 0, 1/42, 0, -1/30, 0, 1/6, 1/2, 1, 0]
Confirm numerically:
XCas> n:=5; sum(1./k^2, k=1..n) + corr(N=n+1), pi*pi/6.
(5, 1.64493406685, 1.64493406685)
Now, we are ready to convert corr into CF formula.
corr assumed N=n+1, but we wanted N=n+1/2, so we shift, and flip it.
Note, we replace N by N+1/2 in 1 step, instead of N by n+1 then n by N-1/2
XCas> [top,bot] := f2nd(1/corr(N=N+1/2)) :;
XCas> q:=quo(top,bot,N); [top,bot] := [bot,top-q*bot]:; → N
Run last line a few more times, we have the other quotients:
12*N, 5/16*N, 448/81*N, 729/4096*N, 180224/50625*N, 8125/65536*N
Convert simple CF to generalized CF, we have (note: now N = n+0.5)
corr = 1/(N+ 1/(12N + 16/(5N + 81/(28*N + 256/(9*N + 5^4/(44*N + 6^4/(13*N + ...
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Gerson W. Barbosa - 11-01-2021 05:04 AM
(11-01-2021 12:56 AM)Albert Chan Wrote: We can get the CF correction formula, using Euler–Maclaurin formula (see #9)
.
.
.
Convert simple CF to generalized CF, we have (note: now N = n+0.5)
corr = 1/(N+ 1/(12N + 16/(5N + 81/(28*N + 256/(9*N + 5^4/(44*N + 6^4/(13*N + ...
Thank you very much! I am not an expert and cannot fully grasp what you have done, but it looks like you have a proof.
Almost three hundred years ago, before tackling the Basel problem, Euler evaluated the series to a few decimal digits using the Euler-MacLaurin method, incidentally the same method used in these RPL programs. The continued fraction correction appears to require less computational effort to obtain the same number of digits.
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-01-2021 10:42 PM
(10-26-2021 02:12 AM)Gerson W. Barbosa Wrote: That gives linear convergence ( 25/12 digits per iteration ), as you can see by the Decimal Basic code and output ...
This may have explained the mysterious 25/12 digits per iteration.
According to Loch's theorem, for "typical" real number between 0 and 1
Asymptote when accuracy is extremely good:
Gain in decimal accuracy (per iteration) ≈ pi^2/(6*ln(2)*ln(10)) ≈ 1.0306 digit
With Decimal Basic code, I checked for for digits accuracy, difference to ζ(2) = pi^2/6
Note: digits = 1 - log10(abs(pi^2/6 - x)). So, if x=1, it is 1.1905 digits accurate.
Anyway, we are only interested in differences.
n=100: 210.6882
n=101: 212.7781 → gained 2.0899 digit
n=102: 214.8680 → gained 2.0899 digit
As experiment, I added an extra CF term to b, and compare against original b=0
n=100: 211.9306 → gained 1.2424 digit
n=101: 214.0206 → gained 1.2425 digit
n=102: 216.1107 → gained 1.2427 digit
This is better than 1.0306 digits. Perhaps CF's not "typical".
Simple CF coefs is decreasing (both odd and even, showing side-by-side)
Code:
[1.0, 12.0,
0.3125, 5.53086419753,
0.177978515625, 3.55998024691,
0.123977661133, 2.62034818825,
0.0950344838202, 2.07210260024,
0.0770232130308, 1.71322636883,
0.0647422966949, 1.46016607257 ] .* NThe difference (~ 0.85 digit) may be due to changes to CF constants.
Going from n to n+1, all the simple CF coefs scaled up, by factor (N+1)/N
We can estimate the gain in accuracy I learn from giac source
When CF coefs increases, estimate is better, roughly by its square. (and, we have n+1 of them)
log10(((N+1)/N)^(2*(n+1))) = 2*(n+1) * log10(1+1/N) = (2N+1)/ln(10) * log1p(1/N)
limit((2N+1)/ln(10) * log1p(1/N), N=inf) = 2/ln(10) ≈ 0.8686 digit
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Gerson W. Barbosa - 11-02-2021 12:28 PM
(11-01-2021 10:42 PM)Albert Chan Wrote:(10-26-2021 02:12 AM)Gerson W. Barbosa Wrote: That gives linear convergence ( 25/12 digits per iteration ), as you can see by the Decimal Basic code and output ...
...
With Decimal Basic code, I checked for for digits accuracy, difference to ζ(2) = pi^2/6
Note: digits = 1 - log10(abs(pi^2/6 - x)). So, if x=1, it is 1.1905 digits accurate.
Anyway, we are only interested in differences.
n=100: 210.6882
n=101: 212.7781 → gained 2.0899 digit
n=102: 214.8680 → gained 2.0899 digit
I had noticed n = 478 was enough for 999 decimal digits. Your calculations suggests 1883/901 is a better estimate for the convergence rate. I've replaced line LET n = CEIL(12*nd/25) with LET n = CEIL(901*nd/1883). On another computer I tried I get 999 decimal digits in 0.27 seconds (same processor, same clock speed, but a cleaner Windows installation).
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-03-2021 12:38 AM
(10-31-2021 03:40 PM)Albert Chan Wrote: Another approach, getting ζ(2) with alternating series.
I translated the code to Decimal Basic, to check its convergence rate.
Amazingly, it is almost the same as original! 2.0899 digit per iteration
However, accuracy is less by 0.4180 digit, compared with non alternating sum.
Cost per loop is a lot smaller though.
Code:
OPTION ARITHMETIC DECIMAL_HIGH
10 INPUT PROMPT "n = ":n
LET x = 3*n*(n+1)+1
LET a = 0
LET b = 0
FOR i = n TO 2 step -1
LET t = i*i
LET a = 2/t - a
LET b = t*t / (x-b)
LET x = x - 4*i
NEXT i
LET z2 = (1-MOD(n,2)*2)/(x-4-1/(x-b)) - a + 2
PRINT "Accurate digits ="; 1-LOG10(ABS(PI*PI/6-z2))
GOTO 10
ENDNote: I did not add the b initial correction, for fair comparison.
n = 100
Accurate digits = 210.27025201561072
n = 101
Accurate digits = 212.36014715889508
n = 102
Accurate digits = 214.45004193935936
n = 400
Accurate digits = 837.2346060430264
n = 401
Accurate digits = 839.32448364856168
n = 402
Accurate digits = 841.41436124813102
For n=478, it almost reached 1000 digits full precision (1 ULP error, last digit = 8 instead of 9)
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-03-2021 01:14 AM
(11-03-2021 12:38 AM)Albert Chan Wrote: However, accuracy is less by 0.4180 digit, compared with non alternating sum.
Accuracy gap of 2/5 of a digit start from the beginning, from n=1.
For n=400, gap = 0.417975 digit, which is very close to log10(1+ϕ) ≈ 0.4179752805
I don't know if this is related, but it is the same ratio if we add 1 more CF term to estimate ϕ
ϕ = [1; 1, 1, 1, 1, ...] = [1;ϕ] = [1;1,ϕ] = [1;1,1,ϕ] = ...
Adding 1 CF term will make estimate better, roughly by denominator square, ϕ^2 = 1+ϕ
XCas> phi := (1+sqrt(5))/2. → 1.61803398875
XCas> e1 := dfc2f(makelist(1, 1, 20)) - phi → 9.77190839357e-09
XCas> e2 := dfc2f(makelist(1, 1, 21)) - phi → -3.73253694619e-09
XCas> abs(e1/e2) → 2.61803393629
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Albert Chan - 11-03-2021 11:28 PM
(11-03-2021 12:38 AM)Albert Chan Wrote: Note: I did not add the b initial correction, for fair comparison.
...
For n=478, it almost reached 1000 digits full precision (1 ULP error, last digit = 8 instead of 9)
Fibonacci numbers pops out again, when I tried to estimate initial b
XCas> X := n*n+n+1
XCas> T(k) := k*(k+1)/2
XCas> b2 := (n+2)^4/(X+4*T(n+2)) :;
XCas> b1 := (n+1)^4/(X+4*T(n+1) - b2) :;
XCas> expand(partfrac(b1)[1])
3/8*n^2 + 23/32*n + 29/128
Note that 3/8 = F2 / F4
Rinse and repeat, add all the way to b10, we get this initial b:
2584/6765*n^2 + 159050/203401*n + 7976465719/20640116475
Again, 2584/6765 = F18 / F20
We might as well replace top coef with 1/ϕ^2 = 2-ϕ
The others are not as important. Patch with initial b:
LET b = ((n+2.0472)*n+1.0118) * (3-SQR(5))/2
n = 100
Accurate digits = 216.76041181744226
n = 101
Accurate digits = 218.85708591701648
n = 102
Accurate digits = 220.95367741927092
n = 400
Accurate digits = 844.54363777628652
n = 401
Accurate digits = 846.63481219743616
n = 402
Accurate digits = 848.72598292186942
For n=475, it reached 1000 digits full precision.