HP-42S (Free42, DM-42) mini-challenge – Wallis Product

[pinkman]

Significant difference between integer and real constants.
I wasn’t aware of this.

Code:


« 2. SWAP
WHILE DUP REPEAT
 DUP SQ 4. *
 ROT * LASTARG
 DROP 1. - / SWAP 1 -
END
DROP »

[Gerson W. Barbosa]

57.5 bytes (49G/50g only, though), 10000 iterations in 107.19 seconds. No need to wait until 2058... just till 2054 :-)

Code:


« 2. 1. PICK3
   START OVER SQ 4. *
DUP 1. - / * SWAP 1.
- SWAP
   NEXT NIP
»

[pinkman]

The use of START, did not remember (still doing 28S RPL!).
Nice!

[Gerson W. Barbosa]

Why not a FOR loop instead? Sure not a pure stack-only solution, but it saves six bytes and it’s faster, bringing the wait down to 31 years.

[pinkman]

Yes I wanted to use the stack only, following the philosophy of your original challenge.
And regarding the START command, it uses an anonymous local variable, it’s quite border line! (Only joking)
32 years for a few pi decimals seems acceptable, at least in comparison with the 2500 years that supposedly separate us from Pythagoras.

[Gerson W. Barbosa]

(07-06-2020 02:31 AM)Gerson W. Barbosa Wrote:  I have a 141-byte, 81-step program that gives 12 digits of \(\pi\) in less than seven seconds on the real HP-42s and 34 digits instantly on my smartphone: 3.141592653589793238462643383279504 (the last digit should be a ‘3’), using the Wallis Product as a basis.

Down to 135 bytes and 74 steps, including LBL and END. It might be optimized even more by HP-41 programmers who usually avoid using too many recall arithmetic as I do. Anyway, on the 42S I tend to value number of steps over byte count. Right now I have 5906 free bytes on the 42S and 1773076480 on Free42, way more memory than I’ll ever need :-)

Code:


00 { 135-Byte Prgm }
01▸LBL "W"
02 STO 01
03 2
04 STO 02
05 MOD
06 ENTER
07 STO+ ST X
08 1
09 -
10 4
11 STO 04
12 RCL× 01
13 1
14 RCL- ST T
15 ×
16 RCL+ ST Z
17 RCL+ ST Z
18 1
19 RCL+ 01
20 X<> 03
21 CLX
22 3
23 RCL× ST T
24 STO- 03
25 STO+ ST X
26 STO- 04
27 R↓
28 4
29 RCL× 01
30 RCL- ST Y
31 RCL+ 02
32 RCL 03
33 RCL- 04
34 RCL+ 02
35 RCL× ST T
36 RCL÷ ST Z
37 R↓
38▸LBL 00
39 STO+ ST T
40 R↓
41 R↓
42 +/-
43 STO× ST T
44 ×
45 RCL 03
46 X<>Y
47 ÷
48 STO ST T
49 R↓
50 RCL 01
51 STO+ ST X
52 X↑2
53 STO× 02
54 DSE ST X
55 STO÷ 02
56 R↓
57 RCL ST Y
58 SIGN
59 X<> ST Z
60 ABS
61 X<> 04
62 +/-
63 STO+ 03
64 NOT
65 +/-
66 NOT
67 +/-
68 X<> 04
69 DSE 01
70 GTO 00
71 R↑
72 RCL- ST T
73 RCL× 02
74 END

HP-75C program:

Code:


10 INPUT N
15 T=MOD(N,2)
20 S=2*T-1
25 L=N+1-3*T
30 M=4-6*T
35 E=4*N*(1-T)+2*T
40 D=4*N-E+2
45 C=S*(L-M+2)/E
50 W=1
55 FOR I=1 TO N
60 S=-S
65 T=4*I*I
70 W=W*T/(T-1)
75 C=S*L/(C+D)
80 L=L-M @ M=2-M
85 T=D @ D=E @ E=T
90 NEXT I
95 DISP 2*W*(C+1)

Register usage on 42S/Free42:

N, I: R01
W: R02
L: R03
M: R04
C, D, E, S, T: Stack

[Gerson W. Barbosa]

(07-05-2020 03:52 AM)Werner Wrote:  It can be simplified a lot:

Code:

00 { 27-Byte Prgm }
01▸LBL "WP"
02 ENTER
03 STO+ ST Y
04 COMB
05 4
06 LASTX
07 Y↑X
08 LASTX
09 X<> ST Z
10 ÷
11 X↑2
12 X<>Y
13 0.5
14 +
15 ÷
16 END

Yes, it can. If by simplified you don’t mean a lesser byte count, of course.

Code:



00 { 29-Byte Prgm }
01▸LBL "WP"
02 N!
03 STO× ST X
04 0.5
05 RCL+ ST L
06 GAMMA
07 1
08 RCL+ ST L
09 GAMMA
10 ×
11 ÷
12 PI
13 ×
14 END

Limits are 114 and 1166 on the 42S and Free42, respectively. These can be increased to 252 and 2123 at the cost of a couple of steps and three bytes:

Code:



00 { 32-Byte Prgm }
01▸LBL "WP"
02 N!
03 RCL ST X
04 0.5
05 RCL+ ST L
06 GAMMA
07 STO÷ ST Y
08 R↓
09 1
10 RCL+ ST L
11 GAMMA
12 ÷
13 ×
14 PI
15 ×
16 END

Regards,

Gerson.

[Werner]

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

Cheers, Werner

[Gerson W. Barbosa]

(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

[Gerson W. Barbosa]

Lest you think I am an expert in continued fractions, which definitely I AM NOT,
I will disclose my reply to someone who wanted to know how I had found the above formula.

Please do not be too disappointed at my "method".

-------

If you compute the Wallis product to 10000 digits,

3.140807746030394560485868201553507

and take its difference to pi,

0.000078534907871191484085332565629

you’ll notice the first significant digits are the same of those of pi/4,

0.7853981633

If you divide pi/4 by that difference you’ll get

10000.62500468726562573348383657523

Multiply this by 8 (because the fractional part is .625, whose multiplication by 8 will lead to a near-integer) and get

80005.00003749812500586787069260184

That’s the ‘8*n + 5’ part in the formula.

Take the fractional part of that,

0.00003749812500586787069260184

Invert it:

26668.00006249687481715340569058372

Multiply by 3:

80004.00018749062445146021707175116

That’s the first term of the continued fraction:

3/(8*n + 4)

Again, take the fractional part of that and invert it:

5333.600028938469868026787967169474

If you multiply this by 3 you’ll notice you’ll have to multiply it again by 5 so that the result is a near integer. So, let’s multiply it by 15:

80004.0004340770480204018195075421

Hence the next term of the continued fraction:

15/(8*n + 4)

So far we have

(pi/4)/(pi - Wn) = 8n + 5 + 3/(8n + 4 + 15/(8n + 4 + ...))

If you keep doing this for the next terms you’ll soon run out of precision to get meaningful results. So, now that you know the first
two terms, start over with W(1000) instead of W(10000). Thus, you can go up to 6 terms, which suffices for unveiling the pattern:

The first part before the CF starts is 8*n + 5, then the numerators are 3, 15, 35, 63, 99... and the denominators are 8*n + 4 + ...

Isolate pi and there you are.

-------

HP 50g program

(263 bytes, CHKSUM = # 9159h)

Code:


« PUSH RAD -105 CF -3 CF DUP 10 * 21 / CEIL 8 OVER * 4 + SWAP 0 1 ROT OVER
  FOR i i DUP + SQ * LASTARG NIP 1 - DUP UNROT / 4 ROLLD PICK3 ROT + / ROT -1
  STEP UNROT + 1 - 2 4 ROT / + * EXPAND FXND DUP SIZE R→I ALOG OVER - PICK3 *
SWAP IQUOT + →STR DUP HEAD 0 I→R →STR TAIL + SWAP TAIL + 1 ROT 2 + SUB POP
»

Time to compute 314 decimal digits of \(\pi\) on the real 50g: 2979.4707 seconds (14.9503 seconds on the emulator).

The formula is linearly convergent (about 2.1 digits per iteration).
Blame it on the slowness of EXPAND to process one algebraic expression containing 150 product and 151 continued fraction terms.

3.141592653589793238462643383279502884197169399375105820974944592307816406286208​99862803482534211706798214808651328230664709384460955058223172535940812848111745​02841027019385211055596446229489549303819644288109756659334461284756482337867831​6527120190914564856692346034861045432664821339360726024914127372458700660631

[Gerson W. Barbosa]

Back on topic, I've received two solutions from former forum member Mike (Stgt). J-M Baillard's and Werner's solution, which I assume to be exactly the same, are unbeatable. Please consider these as hors-concours submissions. The first Mike's solution is one byte less, but one step longer, hence his attempt to overcome that limitation by producing his own plug-in ROM (user-code only) with some more trig funcions, which I haven't tested yet. Again, please don't take that as some kind of lack of sportmanship (Mike used the German Fußball term "Nachtreten"), rather as another valuable contribution to this interesting topic.

Code:

 
00 { 22-Byte Prgm }
01▸LBL "WP"
02 2
03▸LBL 01
04 RCL ST Y
05 STO+ ST X
06 1/X
07 ASIN
08 COS
09 X↑2
10 ÷
11 DSE ST Y
12 GTO 01
13 END

(24 bytes on the HP-41)

Plug-in ROM HP-41 version, 12-step, 24 bytes:

Code:


01▶LBL "WP"
02 2
03▶LBL 01
04 RCL Y
05 ST+ X
06 XROM "ACSC"
07 COS
08 X^2
09 /
10 DSE Y
11 GTO 01
12 END

"ACSC" atands for arccosecant, the inverse cosecant function.

[Gerson W. Barbosa]

(07-04-2020 09:41 PM)JMBaillard Wrote:  Hi everyone,
here is one solution:

01 LBL "WP"
02 2
03 LBL 01
04 RCL Y
05 ST+ X
06 X^2
07 ST* Y
08 DSE X
09 /
10 DSE Y
11 GTO 01
12 END

23 bytes on free42
25 bytes on HP41

Same number of bytes, same number of steps, HP-42 only, inspired in a near-solution by former forum member Mike (Stgt):

Code:



00 { 23-Byte Prgm }
01▸LBL "WP"
02 2
03▸LBL 01
04 RCL ST Y
05 STO+ ST X
06 ACOSH
07 TANH
08 X↑2
09 ÷
10 DSE ST Y
11 GTO 01
12 END