\(\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!

I should participate but I have a really busy weekend!

I’ll give it a try tonight instead of sleeping

(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 :-)

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

Werner

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.

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!

.

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.

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.

(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.

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

(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!

.

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.

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.
After a first try of 19 steps, I finally found a 12 steps solutions, the same than JMBaillard :/

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.

(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!

(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

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 »

(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.

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