Post Reply 
Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B]
10-26-2021, 02:12 AM (This post was last modified: 10-26-2021 03:58 AM by Gerson W. Barbosa.)
Post: #3
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B]
(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) Smile

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 seconds

I 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)
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: Evaluation of ζ(2) by the definition (sort of) [HP-42S & HP-71B] - Gerson W. Barbosa - 10-26-2021 02:12 AM



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