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

[Gerson W. Barbosa]

\(\pi \approx 2\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 )\)

Write a program to evaluate the Wallis product on your HP-42S or DM-42.

Restrictions:

1. Stack-only (X, Y, Z, T & L);
2. Keep at least one copy of the previous content of the X register on the stack (no need to leave it on Y);
3. No limitation for the value of n.

Least number of steps then least number of bytes, including a 2-character label, wins.

Examples:

1 XEQ “WP” -> 2.66666666667
10 -> 3.06770380664
123 -> 3.13523960409
1E6 -> 3.141591868192120714596476681464745 (Free42 only)

Occasional diferences in the last two significant digits will be irrelevant.

Have fun!

[pinkman]

I should participate but I have a really busy weekend!
I’ll give it a try tonight instead of sleeping

[Gerson W. Barbosa]

(07-04-2020 04:57 PM)pinkman Wrote:  I should participate but I have a really busy weekend!
I’ll give it a try tonight instead of sleeping

No need to lose sleep over this one. I will wait until next Wednesday before posting my solution, when most likely it will have been surpassed more than once :-)

[Werner]

Without the END but including LBL "WP": 11 steps, 23 bytes.
Werner

[JMBaillard]

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

Happy Independence Day.

[mfleming]

Pretty off-topic, but I enjoy doing these mini-challenges in other programming languages as well

Code:


      wallis←{2××/(n-1)÷⍨n←4×(⍳⍵)*2}
      wallis 50000
3.141576946

Keep'em coming!

[Valentin Albillo]

.
Hi, Gerson:

(07-04-2020 04:42 PM)Gerson W. Barbosa Wrote:  Least number of steps then least number of bytes, including a 2-character label, wins.

A pity, should you have included a "fastest running time" win condition this "let's-think-out-of-the-box" concoction o'mine would've surely won, as its running time is about the same regardless of the input:

[19-step, 36-byte program]

01   LBL "WP"
02   STO ST Y
03   N!
04   X^2
05   RCL ST Y
06   STO+ ST X
07   N!
08   /
09   X^2
10   2
11   2
12   RCLx ST T
13   STO+ ST X
14   Y^X
15   *
16   0.5
17   RCL+ ST Z
18   /
19   END

Some Examples (run in Free42)

1 XEQ “WP” -> 2.66666666667
10 -> 3.06770380664
123 -> 3.13523960403

which virtually match your results.

As stated, the time is about the same whatever the argument, i.e., extremely fast, and it leaves the original contents of X when the program is run (i.e., the argument keyed in) in stack Y, Z and T. However, it's limited to those arguments that won't overflow the calculation, mostly the factorials, and a few steps can surely be shaved here and there, but I'm just too sleepy right now. :O

Thanks for the nice mini-challenge, I saw it at 4:00 a.m. and it was just what I needed before going to bed.

Have a nice weekend.
V.

[Gerson W. Barbosa]

Hello, Valentín!

(07-05-2020 02:07 AM)Valentin Albillo Wrote:  .
Hi, Gerson:

(07-04-2020 04:42 PM)Gerson W. Barbosa Wrote:  Least number of steps then least number of bytes, including a 2-character label, wins.

A pity, should you have included a "fastest running time" win condition this "let's-think-out-of-the-box" concoction o'mine would've surely won, as its running time is about the same regardless of the input:

...

Sorry, but the most important restriction here was evaluating the product for n as large as 1E6 or 1E7, which is impossible with that approach. Anyway, thanks for calling my attention to alternative methods as they may be useful for another application I have in mind for which a much more narrow range will suffice.

More or less in the same vein:

Code:



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

Examples:

414 -> 3.13969841676

5102 -> 3.141438733173932101937205797737162

These are the maximum values of n for the HP-42S and Free42, respectively, that will not cause overflow.

(07-05-2020 02:07 AM)Valentin Albillo Wrote:  Thanks for the nice mini-challenge, I saw it at 4:00 a.m. and it was just what I needed before going to bed.

Have a nice weekend.

You too! You are most welcome. Thanks again for you contribution.

Gerson.

[Gerson W. Barbosa]

(07-05-2020 12:56 AM)mfleming Wrote:  Pretty off-topic, but I enjoy doing these mini-challenges in other programming languages as well

Code:


      wallis←{2××/(n-1)÷⍨n←4×(⍳⍵)*2}
      wallis 50000
3.141576946

APL? It looks greek to me. But no problem, solutions in other programming languages are always welcome.

[Werner]

I came up with the same formula, Valentin.
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

Blame it on the 4 am ;-)
Cheers, Werner

[mfleming]

(07-05-2020 03:50 AM)Gerson W. Barbosa Wrote:  

(07-05-2020 12:56 AM)mfleming Wrote:  Pretty off-topic, but I enjoy doing these mini-challenges in other programming languages as well

Code:


      wallis←{2××/(n-1)÷⍨n←4×(⍳⍵)*2}
      wallis 50000
3.141576946

APL? It looks greek to me. But no problem, solutions in other programming languages are always welcome.

APL indeed, and thank you for all your posts over the years!

[Valentin Albillo]

.
Hi, Werner:

(07-05-2020 03:52 AM)Werner Wrote:  I came up with the same formula, Valentin. It can be simplified a lot:

00 { 27-Byte Prgm }
01▸LBL "WP"
02 ENTER
03 STO+ ST Y
04 COMB
[...]

Blame it on the 4 am ;-)

Hehe, I certainly could, I was more or less a zombie at those wee hours but I actually did notice the COMB possibility.

However I preferred to use N! instead because I thought that it added some extra value to the computation, thinking that N! behaved just like X! does in many HP models, i.e., if the argument is integer it would provide the factorial as expected, but if x wasn't integer it would return the equivalent Gamma function.

Alas, in the morning I checked it and no, in the HP42S N! doesn't work for noninteger arguments, so you have to code the equivalent Gamma evaluation in its place, which for my code above merely consist in putting the two steps "1 +" before each of the two factorials and changing both factorials to GAMMA, regardless of further optimizations (which surely are possible).

After both insertions my code runs fine and now allows for these interesting possibilites:

a) What's the value of the Wallis Product if adding up 10.5 terms ?

            10.5 XEQ "WP" -> 3.07102204549

b) What's the value of the Wallis Product if adding up Pi terms ?

            PI XEQ "WP" -> 2.93377409474

c) How many terms do we need to add up to get 3 ?

            4.9132932375 (terms) XEQ "WP" -> 3

Best regards.
V.

[Gerson W. Barbosa]

The Wallis Product converges very slowly towards \(\frac{\pi}{2}\). J-M Baillard’s program above, which is about 30% faster than mine, takes 173.2 seconds on my HP-42S to run 1000 iterations, giving only two correct decimal digits of \(\pi\): 3.14080774637. Quite disappointing, isn’t it? How long would it take to compute 12 significant digits? As you have figured out by now, 10^d iterations are required for d correct significant digits. How long would it take to return 3.14159265359 on the HP-42S, or to return \(\pi\) to 34 digits +/- a couple of ULPs on your computer or smartphone running Free42? Just do the math and be astounded by the results.

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. These are achieved with only 8 and 24 iterations, respectively. I’ve managed to save a few bytes and steps by shamelessly borrowing part of J-M Baillard’s code. I will publish it here later, if you are interested.

Thank you all for your participation, contributions and great insights so far.

Gerson.

Edited do fix a few typos.

[pinkman]

After a first try of 19 steps, I finally found a 12 steps solutions, the same than JMBaillard :/

[Gerson W. Barbosa]

Here is mine, 30 bytes, 15 steps:

Code:


00 { 30-Byte Prgm }
01▸LBL "WP"
02 2
03▸LBL 00
04 RCL ST Y
05 X↑2
06 STO+ ST X
07 STO+ ST X
08 SIGN
09 X<> ST L
10 RCL- ST L
11 RCL÷ ST L
12 ÷
13 DSE ST Y
14 GTO 00
15 END

Thank you all for the great three-step saving!

Gerson.

[EdS2]

(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. These are achieved with only 8 and 24 iterations, respectively. I’ve managed to save a few bytes and steps by shamelessly borrowing part of J-M Baillard’s code. I will publish it here later, if you are interested.

Interested? But of course!

[Gerson W. Barbosa]

(07-06-2020 11:06 AM)EdS2 Wrote:  

(07-06-2020 02:31 AM)Gerson W. Barbosa Wrote:  ...

I’ve managed to save a few bytes and steps by shamelessly borrowing part of J-M Baillard’s code. I will publish it here later, if you are interested.

Interested? But of course!

So, here it is, based on the amazing, on the astounding Wallis-Wasicki formula :-)

Code:


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

01: 3.151515151515151515151515151515152
02: 3.140740740740740740740740740740741
03: 3.141577263316393751176359872012044
04: 3.141593864074484229522989212911695
05: 3.141592682610236909552589409774771
06: 3.141592651452292998395352098655921
07: 3.141592653532218907530988374852819
08: 3.141592653593856174058282094810482
09: 3.141592653589909875031699539089613
10: 3.141592653589785247473875899289923
11: 3.141592653589792999646253245914736
12: 3.141592653589793254475588473668716
13: 3.141592653589793238954647785071681
14: 3.141592653589793238430189382757839
15: 3.141592653589793238461625819757986
16: 3.141592653589793238462709648277746
17: 3.141592653589793238462645493302460
18: 3.141592653589793238462643247285411
19: 3.141592653589793238462643378896107
20: 3.141592653589793238462643383559629
21: 3.141592653589793238462643383288622
22: 3.141592653589793238462643383278927
23: 3.141592653589793238462643383279480
24: 3.141592653589793238462643383279504

The HP-71B program might be more readable:

Code:



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

RUN

? 8

3.14159265359

For more efficiency it is better to process two terms at a time, but then the results for odd n are not accessible:

Code:



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

RUN

? 8

3.14159265361

[pinkman]

Amazing!

I know it is off topic, but here is my attempt in RPL: only stack, no local vars, preserve stack, only 4 levels (just for fun!):

Code:


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

[Gerson W. Barbosa]

(07-06-2020 06:04 PM)pinkman Wrote:  I know it is off topic, but here is my attempt in RPL: only stack, no local vars, preserve stack, only 4 levels (just for fun!):

Code:


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

Not off-topic at all, since it’s still on Wallis.

67.5 bytes and (de)counting! Not that I intend to do an attempt, at least for now. On the 50g it should run faster if you use real constants: 119.17 seconds for ten thousand iterations. Should we decide to try 1E11 iterations on the real 50g to possibly get 10 correct decimal digits, we would have to wait until the Summer of 2058 before eventually getting an answer.

[pinkman]

2058... but we will all be dead because of Y2038 bug!