HP-42S (Free42, DM-42) mini-challenge – Wallis Product
07-13-2020, 08:33 PM (This post was last modified: 07-16-2020 11:07 AM by Gerson W. Barbosa.)
Post: #29
 Gerson W. Barbosa Senior Member Posts: 1,487 Joined: Dec 2013
RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product
(07-13-2020 07:04 AM)Werner Wrote:  The limit of my routine on Free42 is 10203, which is quite a lot larger though.
It could've been slightly larger still, but COMB errors out on C(20408,10204) even though the result is less than 1e6145. COMB should be able to compute up to C(20420,10210) using a little trick, as demonstrated in this RPN version:

Code:
00 { 71-Byte Prgm }01▸LBL "CMB" 02 LSTO "L" 03 R↓ 04 R↓ 05 LSTO "Z" 06 R↓ 07 LSTO "T" 08 R↓ 09 STO- ST Y 10 X>Y? 11 X<>Y 12 STO+ ST Y 13 1 14 X>Y? 15 GTO 00 16 1ᴇ-4 17▸LBL 02 18 RCL× ST T 19 RCL÷ ST Y 20 DSE ST T 21 ISG ST Y 22 X<>Y 23 DSE ST Z 24 GTO 02 25 1ᴇ4 26 × 27▸LBL 00 28 RCL "L" 29 ABS 30 RCL "T" 31 STO ST Y 32 RCL "Z" 33 R↑ 34 END

Yes, I had noticed that. The farther you can go, the better, no problem a few extra bytes. I thought I would need more, but n = 1000 on Free42 was enough for me to find out that

$\frac{2\pi }{\pi -2W_{n}}= 8n+5+\frac{3}{8n+4+\frac{15}{8n+4+\frac{35}{8n+4+\frac{63}{8n+4+\frac{99}{\ddots​ } }}}}$

The denominators of the continued fraction in the right-hand side of the expression are easily recognizable. If $$\pi$$ is isolated and the continued fraction, taken with only the first three terms, is factored we get

$\pi \approx W_{n}\cdot \frac{4096n^{4}+8704n^{3}+10304n^{2}+6224n+1473}{4096n^{4}+7680n^{3}+8768n^{2}+4​656n+945 }$

or

$\pi \approx W_{n}\cdot \frac{n^{4}+\frac{17}{8}n^{3}+\frac{161}{64}n^{2}+\frac{389}{256}n+\frac{1473}{4​096}}{n^{4}+\frac{15}{8}n^{3}+\frac{137}{64}n^{2}+\frac{291}{256}n+\frac{945}{40​96}}$

The above polynomial exactly matches the best approximation polynomial presented in the paper linked by Thibault here, which is kind of a loose proof that the formula holds. Notice that the method used by the author is not suitable for the obtention of higher order approximation polynomials. Given the regularity of the continued fraction such polynomials are easily obtainable, but they are not necessary as the continued fraction allows for easier programming.

We are now ready to present the new, the astounding, the amazing... and the unknown Wallis-Wasicki formula :-) Well, since it was easy to get with virtually no math, I don't think it might be that unknown. If any of you find it in a book or over the internet, please let me know so I can name it properly.

$\pi = W_{n}\cdot \left ( 2+\frac{4}{8n+3+\frac{3}{8n+4+\frac{15}{8n+4+\frac{35}{8n+4+\frac{63}{\ddots }}}} } \right )$

where $$W_{n}$$ is the Wallis Product evaluated to n terms.

Or, written as an approximation,

$\pi \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 ( 2+\frac{4}{8n+3+\frac{3}{8n+4+\frac{15}{8n+4+ \frac{35}{8n+4 + \frac{63}{\dots\frac{\ddots }{8n+4+\frac{4n^{2}-1}{8n+4}}} }}} } \right )$

This is of course exact when n → ∞

This appears to produce slightly more than two decimal digits per term, I don't know the exact rate of convergence but the computation of 1000 decimal digits of pi required 478 iterations. Because that has rendered the older formula obsolete, I decided to rewrite my just optimized HP-42S/Free42 program above in this thread. This is the first, still unoptimized version:

Code:
 00 { 57-Byte Prgm } 01▸LBL "W" 02 STO 01 03 FP 04 8 05 RCL× ST L 06 4 07 + 08 2 09▸LBL 00 10 RCL 01 11 STO+ ST X 12 X↑2 13 STO× ST Y 14 DSE ST X 15 STO÷ ST Y 16 X<> ST Z 17 RCL+ ST T 18 STO÷ ST Z 19 X<> ST L 20 X<>Y 21 DSE 01 22 GTO 00 23 R↓ 24 + 25 1 26 + 27 RCL× ST Z 28 -2 29 RCL+ ST L 30 ÷ 31 END

01: 3.140740740740740740740740740740741
02: 3.141598489616110761485210824417872
03: 3.141592609971230454956940475005568
04: 3.141592653927057641418818559928142
05: 3.141592653587141242124980978994893
06: 3.141592653589814292336411237641248
07: 3.141592653589793070317559127943467
08: 3.141592653589793239810871784159143
09: 3.141592653589793238451802835169985
10: 3.141592653589793238462730724196084
11: 3.141592653589793238462642678516210
12: 3.141592653589793238462643388972955
13: 3.141592653589793238462643383233466
14: 3.141592653589793238462643383279871
15: 3.141592653589793238462643383279500
16: 3.141592653589793238462643383279502

For readability, here is an HP-75C version:

Code:
 10 INPUT N 15 C=0 20 D=8*N+4 25 W=1 30 FOR I=N TO 1 STEP -1 35 T=4*I*I 40 W=W*T/(T-1) 45 C=(T-1)/(C+D) 50 NEXT I 55 C=C+D+1 60 DISP 2*C*W/(C-2)

According to Pi: A Source Book, by Lennart Berggren, Jonathan Borwein and Peter Borwein, "not even 100 years of computing on a supercomputer programmed to add or multiply the terms of either sequence [Wallis' Product (1665) and Gregory's Series (1671)] would yield 100 digits of pi" (page 113). And yet, with the simple correction factor above, 1000 digits of $$\pi$$ can be computed in less than half a second on an ordinary desktop computer running Decimal BASIC:

Code:
 OPTION ARITHMETIC DECIMAL_HIGH INPUT  PROMPT "number of decimal digits: ":nd LET tm = TIME  LET n = INT(nd*239/500)  PRINT "n = ";n LET c = 0 LET d = 8*n + 4  LET w = 1 FOR i = n TO 1 STEP -1     LET t = 4*i*i    LET w = w*t/(t - 1)    LET c = (t - 1)/(c + d) next i LET c = c + d + 1 LET p = 2*c*w/(c - 2) LET r = TIME - tm LET r$= STR$(2*c*w/(c - 2) - INT(2*c*w/(c - 2))) 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 "0.##": r; PRINT " seconds" END

number of decimal digits: 1000
n =  478

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 2164201989

Runtime:  0.43 seconds

Best regards,

Gerson.

———

P.S.: Stack-only HP-42S/Free42 version:

Code:
 00 { 59-Byte Prgm } 01▸LBL "W" 02 RCL ST X 03 8 04 × 05 4 06 + 07 0 08 2 09 X<> ST Z 10▸LBL 00 11 R↑ 12 RCL+ ST X 13 STO× ST X 14 STO× ST T 15 DSE ST X 16 STO÷ ST T 17 X<>Y 18 STO+ ST Z 19 X<> ST Z 20 STO÷ ST Y 21 X<> ST L 22 R↓ 23 X<>Y 24 DSE ST T 25 GTO 00 26 1 27 + 28 + 29 × 30 -2 31 RCL+ ST L 32 ÷ 33 END

5 XEQ W →

Y: 3.00217595456
X: 3.14159265359

~ 1.7 seconds on the real HP-42S.

———

P.P.S.: We can save four bytes and four steps, if all we want is the best $$\pi$$ approximation:

Code:
 00 { 55-Byte Prgm } 01▸LBL "W" 02 4 03 0 04 8 05 RCL× ST T 06 RCL+ ST Z 07▸LBL 00 08 R↑ 09 RCL+ ST X 10 STO× ST X 11 STO× ST T 12 DSE ST X 13 STO÷ ST T 14 X<>Y 15 STO+ ST Z 16 X<> ST Z 17 STO÷ ST Y 18 X<> ST L 19 R↓ 20 X<>Y 21 DSE ST T 22 GTO 00 23 DSE ST X 24 + 25 1/X 26 0.5 27 + 28 × 29 END

5 XEQ W →

3.14159265358
( ~ 1.6 seconds )

The Wallis Product approximation is still accessible after program completion:

RCL ÷ . ST L 2 ÷ →

3.00217595456
 « Next Oldest | Next Newest »

 Messages In This Thread HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-04-2020, 04:42 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - pinkman - 07-04-2020, 04:57 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-04-2020, 05:07 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Werner - 07-04-2020, 09:25 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - JMBaillard - 07-04-2020, 09:41 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 04-22-2021, 11:27 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - mfleming - 07-05-2020, 12:56 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-05-2020, 03:50 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - mfleming - 07-05-2020, 05:32 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Valentin Albillo - 07-05-2020, 02:07 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-05-2020, 03:43 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Werner - 07-05-2020, 03:52 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Valentin Albillo - 07-05-2020, 11:20 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-13-2020, 03:32 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-06-2020, 02:31 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - EdS2 - 07-06-2020, 11:06 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-06-2020, 12:40 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-11-2020, 06:43 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - pinkman - 07-06-2020, 06:32 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-06-2020, 10:47 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - pinkman - 07-06-2020, 06:04 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-06-2020, 06:53 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - pinkman - 07-06-2020, 08:39 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - pinkman - 07-06-2020, 09:36 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-06-2020, 10:57 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - pinkman - 07-06-2020, 11:12 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-07-2020, 03:53 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - pinkman - 07-07-2020, 06:46 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Werner - 07-13-2020, 07:04 AM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-13-2020 08:33 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-20-2020, 02:49 PM RE: HP-42S (Free42, DM-42) mini-challenge – Wallis Product - Gerson W. Barbosa - 07-20-2020, 03:47 PM