# HP Forums

Full Version: Wallis' product exploration
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Last week my cousin told me Wallis' product was quite his favorite formula in mathematics.
[attachment=8075]
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.
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
(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.

Regards.
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} }$$
(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$$
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
(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?

(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...
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.
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..)
You’re right!
[attachment=8076]
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}$$
(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.

(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.
(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.
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.
(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
(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.
(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.
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :