Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]

10232021, 02:49 PM
Post: #1




Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
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. HP42S/Free42 program: Code:
n = 12 on the HP42S and n = 16 on Free42 will suffice for 12 and 34 correct digits, respectively. 6 XEQ “z” → 1.64493406685 16 XEQ “z” → 1.644934066848226436472415166646025 HP71B BASIC program; Code:
RUN ? 6 1.64493406685 Interested readers are invited to provide  optimized versions of the given programs;  versions for other calculator, such as the HP41;  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. 

10252021, 01:29 PM
(This post was last modified: 10252021 07:16 PM by Albert Chan.)
Post: #2




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
Code: function zeta2(n) 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) 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^21), 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^21)  1/k^2, k=2..inf) ζ(2) = 5/3  sum(1/(k^2*(4*k^21)), 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^21)), 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 

10262021, 02:12 AM
(This post was last modified: 10262021 03:58 AM by Gerson W. Barbosa.)
Post: #3




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10252021 01:29 PM)Albert Chan Wrote: That gives linear convergence ( 25/12 digits per iteration ), as you can see by the Decimal Basic code and output: Code:
I used your optimizations (thanks!) and avoided all exponentiationrelated multiplications inside the loop. (10252021 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)(2n1)/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. (10252021 01:29 PM)Albert Chan Wrote: This converge even faster, A cute sum of Ramanujan Is that ok? 50 terms of the above series give only 10 correct digits: 1.644934066(94481120) 

10262021, 09:47 AM
Post: #4




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10262021 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. 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: 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*x1 = 2*(x0.5) (2,12), (4,28) ⇒ y= 8*x4 = 8*(x0.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:(10252021 01:29 PM)Albert Chan Wrote: This converge even faster, A cute sum of Ramanujan Considering number of terms needed to sum 1/k^2, this is a huge improvement. Maybe we can accelerate convergence with CF correction ? 

10262021, 02:13 PM
Post: #5




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
I beg your kindness and pardon,
but to a Mathoklutz like me, What does ζ mean? 10B, 10BII, 12C, 14B, 15C, 16C, 17B, 18C, 19BII, 20b, 22, 29C, 35, 38G, 39G, 41CV, 48G, 97 

10262021, 03:04 PM
(This post was last modified: 10262021 03:08 PM by floppy.)
Post: #6




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10262021 02:13 PM)Ren Wrote: I beg your kindness and pardon,an explanation here https://mathworld.wolfram.com/RiemannZet...Zeta2.html 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://indicohlit.jinr.ru/event/187/co...SAdlaj.pdf HP71 & Multimod, HP41CV/CX & Nov64d, PILBOX, HPILPRINTER 82162A, RS232 HPIL interface, Debian11 64bitsPC w/ pyILPER, V41R9F & EMU71 

10262021, 03:24 PM
Post: #7




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10262021 02:13 PM)Ren Wrote: What does ζ mean? (10252021 01:29 PM)Albert Chan Wrote: ζ(2) = 1.65  sum(1/(k^2*(4*k^21)), k=3..inf) Only 1 term CF correction. It does help. CAS> zeta2(n) := 1.65  sum(1/(k^2*(4*k^21)), 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 

10262021, 03:58 PM
Post: #8




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10262021 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/Riema...afunction 

10262021, 08:28 PM
(This post was last modified: 10272021 01:21 PM by Albert Chan.)
Post: #9




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10262021 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) + ...) (10262021 03:24 PM)Albert Chan Wrote: CAS> zeta2(n) := 1.65  sum(1/(k^2*(4*k^21)), k=3..n) Above example, f(x) = 1/(x^2*(4*x^21) = 1/x^2  1/(x+1/2) + 1/(x1/2) ∫(f(x), x) = 1/x  (ln(x+1/2)  ln(x1/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 

10272021, 01:32 PM
Post: #10




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
Thank you Albert Chan, Gerson W. Barbosa, and floppy for enlightening me!
Mr. Chan, I found the video most helpful (and overwhelming!) B^) 10B, 10BII, 12C, 14B, 15C, 16C, 17B, 18C, 19BII, 20b, 22, 29C, 35, 38G, 39G, 41CV, 48G, 97 

10272021, 05:12 PM
(This post was last modified: 11022021 10:49 PM by Albert Chan.)
Post: #11




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10252021 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) lua> 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 / ((10c)*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 

10292021, 02:16 PM
Post: #12




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10262021 09:47 AM)Albert Chan Wrote: It would be better if we had more coefficients to confirm the pattern. We can also sum odd terms: sum(1/(2k1)^2, k=1..inf) = pi^2/8 = 3/4*ζ(2) ζ(2) = 4/3 * sum(1/(2k1)^2, k=1..inf) If we sum only n terms, correction term ≈ ∫(1/(2x1)^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 

10312021, 03:40 PM
(This post was last modified: 11062021 03:00 PM by Albert Chan.)
Post: #13




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
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: T_{k}  T_{k1} = k*(k+1)/2  (k1)*k/2 = k zeta2_alt(n) is computationally less expensive than zeta2(n) Code: function zeta2_alt(n) Update: 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/6z2) 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 

11012021, 12:56 AM
Post: #14




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
We can get the CF correction formula, using Euler–Maclaurin formula (see #9)
XCas> C(k) := bernoulli(k)/k!; // EulerMaclaurin 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 N1/2 XCas> [top,bot] := f2nd(1/corr(N=N+1/2)) :; XCas> q:=quo(top,bot,N); [top,bot] := [bot,topq*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 + ... 

11012021, 05:04 AM
Post: #15




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(11012021 12:56 AM)Albert Chan Wrote: We can get the CF correction formula, using Euler–Maclaurin formula (see #9) 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 EulerMacLaurin 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. 

11012021, 10:42 PM
Post: #16




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10262021 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 sidebyside) Code: [1.0, 12.0, The 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 

11022021, 12:28 PM
Post: #17




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(11012021 10:42 PM)Albert Chan Wrote:(10262021 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 ... 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). 

11032021, 12:38 AM
Post: #18




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(10312021 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 Note: 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) 

11032021, 01:14 AM
(This post was last modified: 11032021 07:13 AM by Albert Chan.)
Post: #19




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(11032021 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.77190839357e09 XCas> e2 := dfc2f(makelist(1, 1, 21))  phi → 3.73253694619e09 XCas> abs(e1/e2) → 2.61803393629 

11032021, 11:28 PM
Post: #20




RE: Evaluation of ζ(2) by the definition (sort of) [HP42S & HP71B]
(11032021 12:38 AM)Albert Chan Wrote: Note: I did not add the b initial correction, for fair comparison. 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 = F_{2} / F_{4} 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 = F_{18} / F_{20} 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) * (3SQR(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. 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)