Threaded Mode | Linear Mode

Albert Chan · (This post was last modified: 11-07-2021 12:45 PM by Albert Chan.)

Wallis product part is *very slowly* converging; probably not a good idea to cut down CF terms.

Another issue is to limit rounding erors of Wallis product.
To improve accuracy, we keep track of partial product, less 1.

(1+a)*(1+1/t) - 1 = a + (1+a)/t

To recover pi, we just add back 1 (after loops finished)

pi ≈ (1+a)*(1+b)*8/3 = (a*b + b + (a-1/8)) * 8/3 + 3 (*)

Code:

OPTION ARITHMETIC DECIMAL_HIGH

10 INPUT  PROMPT "n = ":n

   LET x = 8*n+4

   LET a = 0

   LET b = 0

   FOR i = n TO 2 STEP -1

      LET t = 4*i*i-1

      LET a = a + (1+a)/t

      LET b = t / (x+b)

   NEXT i

   LET b = 2/(x-1+3/(x+b)) ! CF correction

   LET r = (a*b+b+(a-0.125)) * 8/3 + 3

   PRINT "Accurate digits ="; 1-LOG10(MAX(ABS(PI-r),1E-1000))

   GOTO 10

END

n = 100
Accurate digits = 211.18296908289588
n = 101
Accurate digits = 213.27288787760364
n = 300
Accurate digits = 629.16112325948984
n = 301
Accurate digits = 631.25100446831984

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

Previous version take more total terms (3p=3*326=978) vs here (2n=2*478=956)
But, pairwise sum can be greatly simplified → Previous version is faster, by about 30%

(*) By Sterbenz lemma, term (a-1/8) is exact

n<2: it never enter the loop, thus a = 0
n=2: a = 16/15 - 1 = 1/15, within [1/16,1/4]
With more products, a is increasing. limit(a, n=inf) = 3/8*pi - 1 ≈ 1/5.6 ≤ 1/4