Wallis' product exploration

[pinkman]

Last week my cousin told me Wallis' product was quite his favorite formula in mathematics.
   
I remembered the thread opened on this forum about PI approximations (like 355 / 113), and we decided to program Wallis product on Free42 while having a coffee...

We built it this way:

Code:


01 LBL "WALLIS"
02 1
03 STO+ 00     // n++ (n is register 00)
04 RCL 00      // recl n
05 X^2         // this is n^2
06 4
07 x           // 4n^2 (numerator)
08 ENTER       // duplicate
09 X<>Y        // enable stack lift
10 1
11 -           // 4n^2 - 1 (denominator)
12 /           // 4n^2 / (4n^2 - 1)
13 STOx 01     // pi/2 approximation is in register 01
14 RCL 01
15 2
16 x           // (pi/2 approximation) x 2
17 RTN

The usage is simple:
- initiate n with value 0 in register 00: 0 STO 00
- then initiate product with value 1 in register 01: 1 STO 01
- execute the product as many times as you want: XEQ WALLIS

We were really disapointed!
First iteration: 2.6666...
2nd: 2.844444...
3rd: 2.9...
10th: 3.067...
50th: 3.126...

Yes the product seems to approach PI/2, but soooooo slowly!

My other program included a loop until the precision of the calculator has been reached:

Code:


01 LBL "PICALC"
02 RCL 01      // take last PI/2 approximation...
03 STO 02      // ... and save it
04 1
05 STO+ 00     // n++ (n is register 00)
06 RCL 00      // recl n
07 X^2         // this is n^2
08 4
09 x           // 4n^2 (numerator)
10 ENTER       // duplicate
11 X<>Y        // enable stack lift
12 1
13 -           // 4n^2 - 1 (denominator)
14 /           // 4n^2 / (4n^2 - 1)
15 STOx 01     // pi/2 approximation is in register 01
16 RCL 01      // recl pi/2 current approximation
17 RCL 02      //  recl pi/2 previous approximation
18 -
19 X≠0?        // if maximum calculator precision has not been reached then...
20 GTO "PICALC"   // ... loop and calculate another one
21 RCL 01      // we're done, prove it
22 2
23 x           // by calculating 2 x pi/2 approximation
24 RTN

I tried it a few times and breaked the program after several minutes to get 7 correct significant digits of pi after... 100E6 iterations.

Despite the fact that Wallis' product is really slow (and it's very deceptive), I have one question:
What do you think about calculations precisions?
- I mean, at first iteration the number is rounded, because pi/2 is approximated by 1.3333.... so errors will appear, and will be multiplied at each loop. I even wondered if those approximations could be a reason for the convergence being so slow.
- I also mean, each time we iterate, we calculate a square of big numbers, approaching the limit of the calculator, leading to (4n^2) = (4n^2 - 1) => (4n^2)/(4n^2 - 1) = 1, but, unless my algorithm is wrong, I did not reach the limit.

[Albert Chan]

You can check the result with a direct formula, for n terms, here

Or, using lgamma(), we have

\(\log \left(2 \prod _{k=1} ^{k=n} {4k^2 \over 4k^2-1}\right) =
4\log(\Gamma(n)) - 2\log(\Gamma(2n)) + n \log(16) + 2 \log(n) - \log(4n+2) \)

Example, using log1p and sum the logarithm of the terms backwards, this matched above formula

lua> require'mathx'
lua> log1p = mathx.log1p
lua> s=0; for i=100e6, 1, -1 do s = s - log1p(-1/(4*i*i)) end
lua> = 2*math.exp(s)
3.1415926457358117

[pinkman]

(02-08-2020 02:13 AM)Albert Chan Wrote:  You can check the result with a direct formula, for n terms, here

Really interesting document, thanks! I'll send it to my cousin.
(By the way I won't say using factorial or double factorial is a "direct formula", as it is an iteration also. But... of course if there were a direct formula for PI, nobody would have tried to calculate its decimals with years and years of computing! )
I used the double factorial function of Wolfram Alpha to evaluate pi for n=10, then 100, ... and compare with my algorithm accuracy.

Code:


2/(2*n+1)*((2*n)!!/(2*n-1)!!)^2 where n=10

Wolfram & Free42 agree with the following values:
n=10: 3.06770380664[...more digits with Wolfram but I fixed accuracy on 12]
n=100: 3.13378749063
n=1E3: 3.14080774603
n=1E6: 3.14159186819
n=1E8: 3.14159264574

I can conclude I have no problems of accuracy with.
With Free42 on iPhone XR it takes approximately 6 minutes to compute for 1E8.

(02-08-2020 02:13 AM)Albert Chan Wrote:  Or, using lgamma(), we have

\(\log \left(2 \prod _{k=1} ^{k=n} {4k^2 \over 4k^2-1}\right) =
4\log(\Gamma(n)) - 2\log(\Gamma(2n)) + n \log(16) + 2 \log(n) - \log(4n+2) \)

Example, using log1p and sum the logarithm of the terms backwards, this matched above formula

lua> require'mathx'
lua> log1p = mathx.log1p
lua> s=0; for i=100e6, 1, -1 do s = s - log1p(-1/(4*i*i)) end
lua> = 2*math.exp(s)
3.1415926457358117

I found quite the same for Free42 and Wolfram ! n=1E8: 3.14159264574

Now with the help of your link to calculus7.org and Wolfram, I computed approx. values of PI with the faster converging formula:

Code:


(4*n+3)*((2*n)!!/(2*n+1)!!)^2 where n=10

n=10: 3.14074437347[... 12 significant numbers]
n=100: 3.14158298190
n=1E3: 3.14159255556
n=1E4: 3.14159265261 -> better accuracy than 1st algorithm and n=1E8
n=1E5: 3.14159265358
n=1E6: 3.14159265359 -> 12 significant digits, it will not change anymore

Yes it is a good version of the algorithm.

Next step (tonight?) I'll try to create an optimized algorithm to compute the fractions of double factorials (!!) on Free42.

Thanks again for your reply, Albert.
Regards.

[Allen]

I think the "slowness" of the clever asymptotic behavior of \( \prod_{k=1} ^{k=n} {4k^2 \over 4k^2-1} \) can best be seen by factoring the numerator and denominator for the first few \( n\) to see what is going on:

n=2
\( 2^{6}\over 3^{2} 5^{1}\)

n=3
\( 2^{8}\over 5^{2} 7^{1}\)

n=5
\( 2^{16}\over 3^{4} 7^{2} 11^{1}\)

n=7
\( 2^{22}\over 3^{3} 5^{1} 11^{2} 13^{2}\)

n=11
\( 2^{38}\over 3^{2} 7^{2} 13^{2} 17^{2} 19^{2} 23^{1}\)

n=13
\( 2^{46}\over 3^{3} 5^{4} 7^{2} 17^{2} 19^{2} 23^{2}\)

n=17
\( 2^{64}\over 3^{6} 5^{3} 7^{1} 11^{2} 19^{2} 23^{2} 29^{2} 31^{2}\)

n=19
\( 2^{70}\over 3^{3} 5^{4} 7^{2} 11^{2} 13^{1} 23^{2} 29^{2} 31^{2} 37^{2}\)

n=23
\( 2^{84}\over 3^{6} 5^{4} 13^{2} 29^{2} 31^{2} 37^{2} 41^{2} 43^{2} 47^{1}\)

n=29
\( 2^{108}\over 3^{2} 5^{2} 7^{2} 11^{2} 17^{2} 19^{2} 31^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{1}\)

Do you notice anything interesting about the relationship between \( n \) and the denominator when \( n \) is prime?

What about when it's not prime?

n=28
\( 2^{106}\over 3^{1} 5^{2} 7^{2} 11^{2} 17^{2} 19^{1} 29^{2} 31^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2}\)

n=29
\( 2^{108}\over 3^{2} 5^{2} 7^{2} 11^{2} 17^{2} 19^{2} 31^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{1}\)

n=30
\( 2^{112}\over 7^{2} 11^{2} 17^{2} 19^{2} 31^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{2} 61^{1}\)

n=31
\( 2^{114}\over 3^{2} 7^{3} 11^{2} 17^{2} 19^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{2} 61^{2}\)

n=32
\( 2^{126}\over 3^{4} 5^{1} 7^{4} 11^{2} 13^{1} 17^{2} 19^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{2} 61^{2}\)

n=33
\( 2^{128}\over 3^{2} 5^{2} 7^{4} 13^{2} 17^{2} 19^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{2} 61^{2} 67^{1}\)

n=34
\( 2^{132}\over 3^{3} 5^{2} 7^{4} 13^{2} 19^{2} 23^{1} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{2} 61^{2} 67^{2}\)


Reminds me the highly composite binomial coefficients \( \binom{2n}{n} \).

To represent as \( \pi \) rather than \( \pi \over 2 \), one can just add 1 to the numerator like so:

\( \pi_{29} \approx \frac {2^{108+1}} {3^{2} 5^{2} 7^{2} 11^{2} 17^{2} 19^{2} 31^{2} 37^{2} 41^{2} 43^{2} 47^{2} 53^{2} 59^{1} }\)

[Albert Chan]

(02-08-2020 07:13 PM)Allen Wrote:  Reminds me the highly composite binomial coefficients \( \binom{2n}{n} \).

That is because the formula had this denominator (squared!)

\(\Large 2\prod _{k=1} ^{k=n} {4k^2 \over 4k^2-1} =
{2^{4n+2} \over 4n+2} ÷ \binom{2n}{n}^2 \)

[Allen]

maybe a more compact way of looking at the prime powers in the denominator:

here \( 08002222022000022222222221 \rightarrow 2^{0} 3^{8} 5^{0} 7^{0} 11^{2} 13^{2} 17^{2} 19^{2} 23^{0} 29^{2} 31^{2} 37^{0} 41^{0} 43^{0} 47^{0} 53^{2} 59^{2} 61^{2} 67^{2} 71^{2} 73^{2} 79^{2} 83^{2} 89^{2} 97^{2} 101^{1} \)

The carats under each line indicate the prime factors of n... and a * next to n represents prime n.

Code:

 
  n   prime exponents in denominator from p=3..101]
 50 08002222022000022222222221
    ^ ^                       
 49 06401222022000022222222220
       ^                      
 48 04440222022000022222222210
    ^^                        
*47 06340221022000022222222200
                  ^           
 46 05240220021000222222222200
    ^       ^                 
 45 04230120220000222222222200
     ^^                       
 44 08420020220000222222222100
    ^   ^                     
*43 07422020210000222222222000
                 ^            
 42 06322010200002222222222000
    ^^ ^                      
*41 08242000200002222222221000
                ^             
 40 04242000200022222222220000
    ^ ^                       
 39 00442000200022222222210000
     ^   ^                    
 38 02431200200022222222200000
    ^      ^                  
*37 01220202200022222222200000
               ^              
 36 00020202200222222222100000
    ^^                        
 35 04020202200222222221000000
      ^^                      
 34 03240202100222222220000000
    ^     ^                   
 33 02240222000222222210000000
     ^  ^                     
 32 04142122000222222200000000
    ^                         
*31 02032022000222222200000000
              ^               
 30 00022022002222222100000000
    ^^^                       
*29 02222022002222221000000000
             ^                
 28 01222021022222220000000000
    ^  ^                      
 27 00141020022222220000000000
     ^                        
 26 06040020022222210000000000
    ^    ^                    
 25 05040210022222200000000000
      ^                      
 24 04420200022222200000000000
    ^^                      
*23 06400200022222100000000000
            ^               
 22 04300200222222000000000000
    ^   ^                  
 21 02202200222221000000000000
     ^ ^                   
 20 04222200222210000000000000
    ^ ^                   
*19 03422100222200000000000000
           ^             
 18 02422002222100000000000000
    ^^                  
*17 06312002222000000000000000
          ^            
 16 05201022222000000000000000
    ^                 
 15 04200022221000000000000000
     ^^              
 14 06400022210000000000000000
    ^  ^            
*13 03420022200000000000000000
         ^         
 12 00220222200000000000000000
    ^^             
*11 02020222100000000000000000
        ^         
 10 01012222000000000000000000
    ^ ^         
  9 00202221000000000000000000
     ^         
  8 04202210000000000000000000
    ^          
 *7 03102200000000000000000000
       ^     
  6 02022100000000000000000000
    ^^       
 *5 04021000000000000000000000
      ^     
  4 02220000000000000000000000
    ^     
 *3 00210000000000000000000000
     ^   
 *2 02100000000000000000000000
    ^   
  1 01000000000000000000000000

[Allen]

(02-08-2020 07:53 PM)Albert Chan Wrote:  That is because the formula had this denominator (squared!)
\(\Large 2\prod _{k=1} ^{k=n} {4k^2 \over 4k^2-1} =
{2^{4n+2} \over 4n+2} ÷ \binom{2n}{n}^2 \)

Interesting!!

Since the central binomial coefficients are related to Catalan numbers, Wallis' product \( \pi \) approximation for each \( n \) is inversely proportional to the number of ways a regular \( n \)-gon can be divided into \( n-2 \) triangles.

As \( n \rightarrow \infty \) then the \(n\)-gon shape approaches a circle.

Wait, \( \pi \) is related to circles because of triangles?

[pinkman]

(02-08-2020 08:58 PM)Allen Wrote:  

(02-08-2020 07:53 PM)Albert Chan Wrote:  That is because the formula had this denominator (squared!)
\(\Large 2\prod _{k=1} ^{k=n} {4k^2 \over 4k^2-1} =
{2^{4n+2} \over 4n+2} ÷ \binom{2n}{n}^2 \)

Interesting!!

Since the central binomial coefficients are related to Catalan numbers, Wallis' product \( \pi \) approximation for each \( n \) is inversely proportional to the number of ways a regular \( n \)-gon can be divided into \( n-2 \) triangles.

As \( n \rightarrow \infty \) then the \(n\)-gon shape approaches a circle.

Wait, \( \pi \) is related to circles because of triangles?

Really interesting. And funny
So, this apparently abstract formula of Wallis is the same approach to find PI than any other method: converge to a circle. Well, what did I expect? No PI without circles!
It also leaded me to wake up my memory about binomial coefficients, that I had a bit (completely) forgotten...

[Allen]

What I find mind blowing about it is that it's a combinatoric function, not of just triangles, but how many different ways triangles can tile a nearly-circular polygon.

I tried yesterday to find a printed/published formula for \( \pi \) in terms of catalan numbers, but could not. combining terms with Albert's observation above, we get

\(\pi = \lim\limits_{n\to\infty} \frac {2^{4n+2} } { (2n+1) (n+1)^2 C_n^2} \)

interesting scaling polynomial in the bottom

\( 2 n^3 + 5 n^2 + 4 n + 1 \)

Out of curiosity, I think we can replace the \( (n+1)^2\) with a Triangluar number \( T_n\) and further drop a factor of 4 from the numerator:

\( \pi = \lim\limits_{n\to\infty} \frac {n^2 2^{4n} } { (2n+1) T_n^2 C_n^2} \)

in effort to reduce some polynomial terms, but I can't tell if that's a step forward or a step back.. Ideally, one could represent \( \pi \) as an integer divided by a hand full of (discrete) combinatoric functions.

[Allen]

Using identity 5.36 ( p. 186 of Concrete Mathematics )

\( \binom{n-\frac{1}{2}} {n} = \binom{2n} {n} ÷ 2^{2n} \)

I think we can also reduce the limit to:
\( \pi = \lim\limits_{n\to\infty} \frac {1} {n+\frac {1}{2}} ÷ \binom{n-\frac{1}{2}} {n}^{2} \)

(I could have messed something up there..)

[pinkman]

You’re right!
   

[Allen]

Nice, Thank you for checking in CAS!

of course at infinite limits the extra +1/2 is irrelevant, and since \( \binom{n-\frac{1}{2}} {n}^{2} = \binom{-1/2} {n}^{2} \), we can remove some more symbols.


as we can see on a discrete limit Wolfram alpha (or with continuous n)

\( \pi = \lim\limits_{n\to\infty} \frac {1} {n \binom{-\frac{1}{2}} {n}^{2}} \)

or circumference \(c \) is

\( c = \lim\limits_{n\to\infty} 2 \binom{-\frac{1}{2}} {n}^{-2} \)

[Gerson W. Barbosa]

(02-07-2020 10:55 PM)pinkman Wrote:  First iteration: 2.6666...
2nd: 2.844444...
3rd: 2.9...
10th: 3.067...
50th: 3.126...

Yes the product seems to approach PI/2, but soooooo slowly!

.
.
.

My other program included a loop until the precision of the calculator has been reached:

I tried it a few times and breaked the program after several minutes to get 7 correct significant digits of pi after... 100E6 iterations.
.

Try this:

Code:


00 { 73-Byte Prgm }
01▸LBL "WALLIS"
02 0
03 STO 00
04 1
05 STO 01
06▸LBL 00
07 1
08 STO+ 00
09 RCL 00
10 X↑2
11 4
12 ×
13 ENTER
14 X<>Y
15 1
16 -
17 ÷
18 STO× 01
19 RCL 01
20 STO+ ST X
21 64
22 RCL× 00
23 RCL ST X
24 72
25 +
26 RCL× 00
27 23
28 +
29 X<>Y
30 56
31 +
32 RCL× 00
33 15
34 +
35 ÷
36 ×
37 STOP
38 GTO 00
39 END

Usage: XEQ WALLIS R/S R/S...

or that:

Code:


00 { 76-Byte Prgm }
01▸LBL "WALLIS"
02 STO 03
03 STO 04
04 0
05 STO 00
06 1
07 STO 01
08▸LBL 00
09 1
10 STO+ 00
11 RCL 00
12 X↑2
13 4
14 ×
15 ENTER
16 X<>Y
17 1
18 -
19 ÷
20 STO× 01
21 RCL 01
22 DSE 03
23 GTO 00
24 64
25 RCL× 04
26 RCL ST X
27 72
28 +
29 RCL× 04
30 23
31 +
32 X<>Y
33 56
34 +
35 RCL× 04
36 15
37 +
38 ÷
39 ×
40 STO+ ST X
41 END

Usage: n XEQ WALLIS

Convergence gets worse as n increases, but it won’t take long before you can recognize the constant.

[EdS2]

Interesting thread!

(02-07-2020 10:55 PM)pinkman Wrote:  Despite the fact that Wallis' product is really slow (and it's very deceptive), I have one question:
What do you think about calculations precisions?

I was a bit surprised that we can multiply millions of numbers and still get recognisable results. But then I had a think: a good multiplication or division algorithm will be accurate to within one unit in the last place (an ULP) or indeed, I think, to within half an ULP, because of rounding.

The average error then will be a quarter ULP per multiplication. I'm going to assume the error can be treated as random.

So after a very long series of multiplications, we'll get the exact result, plus an error term which is a random walk of millions of quarter ULPs - which will only add up to an average error on the order of thousands of ULPs, because there's a square root law in such cases.

And thousands of ULPs is only three decimal digits of error. Which is why, I think, we still recognise pi coming out as the result.

[Gerson W. Barbosa]

(02-09-2020 10:58 PM)Gerson W. Barbosa Wrote:  Convergence gets worse as n increases, but it won’t take long before you can recognize the constant.

This is better:

Code:



00 { 97-Byte Prgm }
01▸LBL "WALLIS"
02 STO 03
03 STO 04
04 0
05 STO 00
06 1
07 STO 01
08▸LBL 00
09 1
10 STO+ 00
11 RCL 00
12 X↑2
13 4
14 ×
15 ENTER
16 X<>Y
17 1
18 -
19 ÷
20 STO× 01
21 RCL 01
22 DSE 03
23 GTO 00
24 STO+ ST X
25 512
26 RCL× 04
27 RCL ST X
28 832
29 +
30 RCL× 04
31 592
32 +
33 RCL× 04
34 167
35 +
36 X<>Y
37 704
38 +
39 RCL× 04
40 464
41 +
42 RCL× 04
43 105
44 +
45 ÷
46 ×
47 END

19 XEQ WALLIS ->

3.1415926535(94170710602783580375138)

1000 XEQ WALLIS ->

3.14159265358979323846264(8085575516)

The first program should be changed accordingly. After 15 iterations a physical 42S will return 3.14159265353.

The correction factor is

(n (n (512 n + 832) + 592) + 167)/(n (n (512 n + 704) + 464) + 105)

That’s five terms of a continued fraction in Horner form. I have yet to check whether it generalizes for infinite terms or not.

[pinkman]

A few more readings about convergence speed here: https://pdfs.semanticscholar.org/c767/92...0feefd.pdf

[Gerson W. Barbosa]

Try the following for even n and see what you get.


\(\frac{\pi }{2}\approx \left ( \frac{4}{3} \cdot \frac{16}{15}\cdot \frac{36}{35}\cdot\frac{64}{63} \cdots \frac{ 4n ^{2}}{ 4n ^{2}-1}\right )\left ( 1+\frac{1}{4n+\frac{3}{2-\frac{1}{4n+\frac{5}{2-\frac{3}{4n+\frac{7}{2-\frac{5}{4n+\frac{9}{2-\frac{7}{\dots \frac{ \ddots }{2-\frac{n-3}{4n+\frac{n+1}{2}}}}}}}}}}}} \right )\)



---

P.S.:

This will handle both parities:

\(\frac{\pi }{2}\approx \left ( \frac{4}{3} \times \frac{16}{15}\times \frac{36}{35}\times\frac{64}{63} \times \cdots \times \frac{ 4n ^{2}}{ 4n ^{2}-1}\right )\left ( 1+\frac{1}{4n+\frac{3}{2-\frac{1}{4n+\frac{5}{2-\frac{3}{4n+\frac{7}{2-\frac{5}{4n+\frac{9}{2-\frac{7}{\dots \frac{ \ddots }{4n \left ( n \bmod 2\right) + 2\left ( n+1 \bmod 2 \right )+\frac{1-n+\left ( 2n+1 \right )\left ( n \bmod 2 \right )}{4n \left (n+1 \bmod 2\right)
+ 2\left ( n \bmod 2 \right )}}}}}}}}}}} \right )\)

This approximation gives \(\frac{4}{3}n\) correct significant digits.

[Gerson W. Barbosa]

(02-11-2020 04:11 PM)Gerson W. Barbosa Wrote:  \(\frac{\pi }{2}\approx \left ( \frac{4}{3} \cdot \frac{16}{15}\cdot \frac{36}{35}\cdot\frac{64}{63} \cdots \frac{ 4n ^{2}}{ 4n ^{2}-1}\right )\left ( 1+\frac{1}{4n+\frac{3}{2-\frac{1}{4n+\frac{5}{2-\frac{3}{4n+\frac{7}{2-\frac{5}{4n+\frac{9}{2-\frac{7}{\dots \frac{ \ddots }{2-\frac{n-3}{4n+\frac{n+1}{2}}}}}}}}}}}} \right )\)

That's slightly better:

\[\frac{\pi }{2}\approx \left ( \frac{4}{3} \cdot \frac{16}{15}\cdot \frac{36}{35}\cdot\frac{64}{63} \cdots \frac{ 4n ^{2}}{ 4n ^{2}-1}\right )\left ( 1+\frac{1}{4n+\frac{3}{2-\frac{1}{4n+\frac{5}{2-\frac{3}{4n+\frac{7}{2-\frac{5}{4n+\frac{9}{2-\frac{7}{\dots \frac{ \ddots }{4n+\frac{n+1}{2-\frac{n-1}{4n}}}}}}}}}}}} \right )\]


Here are 1000 digits of \(\pi\), in a reasonable time, courtesy of John Wallis and yours truly :-)

------------------------------------------------------------------------------

n: 842

3.1415926535 8979323846 2643383279 5028841971 6939937510 
  5820974944 5923078164 0628620899 8628034825 3421170679 
  8214808651 3282306647 0938446095 5058223172 5359408128 
  4811174502 8410270193 8521105559 6446229489 5493038196 
  4428810975 6659334461 2847564823 3786783165 2712019091 
  4564856692 3460348610 4543266482 1339360726 0249141273 
  7245870066 0631558817 4881520920 9628292540 9171536436 
  7892590360 0113305305 4882046652 1384146951 9415116094 
  3305727036 5759591953 0921861173 8193261179 3105118548 
  0744623799 6274956735 1885752724 8912279381 8301194912 
  9833673362 4406566430 8602139494 6395224737 1907021798 
  6094370277 0539217176 2931767523 8467481846 7669405132 
  0005681271 4526356082 7785771342 7577896091 7363717872 
  1468440901 2249534301 4654958537 1050792279 6892589235 
  4201995611 2129021960 8640344181 5981362977 4771309960 
  5187072113 4999999837 2978049951 0597317328 1609631859 
  5024459455 3469083026 4252230825 3344685035 2619311881 
  7101000313 7838752886 5875332083 8142061717 7669147303 
  5982534904 2875546873 1159562863 8823537875 9375195778 
  1857780532 1712268066 1300192787 6611195909 216420198 
  
Runtime:   .96 seconds

------------------------------------------------------------------------------

Only linear convergence, but better than waiting forever for a just few digits.

Decimal BASIC program:

Code:



OPTION ARITHMETIC DECIMAL_HIGH
INPUT  PROMPT "n: ":n
LET tm = TIME 
LET nd = 1000
LET pr = 1
LET n1 = n - 1
LET n2 = n1 + 2
LET d = 4*n
LET c = 0
FOR i = 1 TO n/2
   LET np = 4*(2*i - 1)*(2*i - 1)
   LET dp = np - 1
   LET pr = pr*np/dp
   LET np = 16*i*i
   LET dp = np - 1
   LET pr = pr*np/dp
   LET c = n2/(2 - n1/(d + c))
   LET n2 = n2 - 2
   LET n1 = n1 - 2
next i
LET c = 2 + 2/(c + d)
LET p = pr*c
LET r = TIME - tm
LET r$ = STR$(p - INT(p))                           
PRINT
PRINT STR$(INT(p));
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
      PRINT REPEAT$(" ",1 + LEN(STR$(INT(p))));
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "##.##": r;
PRINT " seconds"
END

[pinkman]

(02-11-2020 10:48 PM)Gerson W. Barbosa Wrote:  Only linear convergence, but better than waiting forever for a just few digits.

Thank you for all, this is beautiful.

[Gerson W. Barbosa]

(02-12-2020 10:01 PM)pinkman Wrote:  

(02-11-2020 10:48 PM)Gerson W. Barbosa Wrote:  Only linear convergence, but better than waiting forever for a just few digits.

Thank you for all, this is beautiful.

Thank you and your cousin for having started that conversation about the Wallis Product and posting here, otherwise I probably wouldn’t have tried. Unlike approximation polynomials, the equivalent continued fractions usually follow beautiful patterns. Once you have a continued fraction, those polynomials are just their successive approximants. For instance,

Simplify [1 + 1/(4 n + 3/(2 - 1/(4 n + 5/2)))]

-> (64 n^2 + 72 n + 23)/(64 n^2 + 56 n + 15)

= (n^2 + 9n/8 + 23/64)/(n^2 + 7n/8 + 15/64)

This technique had worked previously for other slowly convergent series, so I guessed it would work for Wallis as well.