Monte-Carlo Pi

[Ángel Martin]

I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNG-related programs and MCODE functions. I'd like to include a Monte-Carlo based calculation of pi, which I *think* has already been published in the forum but can't locate anymore. Any pointers to that reference will be appreciated...

Cheers,
ÁM

[Albert Chan]

Probability of random()/random() rounded to even integer = 1 + 1/4 - PI/4 = (5-PI)/4 ≈ 0.4646

Or, we can count odd's:
Probability of random()/random() rounded to odd integer = (PI-1)/4 ≈ 0.5354

lua> s, n = 0, 5e6
lua> for i=1,n do s = s + round(random()/random())%2 end
lua> s/n
0.5354442

see https://www.hpmuseum.org/forum/thread-13...#pid129470

[Andres]

The very simple and old one I remember starts with generating two uniformly distributed numbers from the interval [0 1]. After invoking R -> P conversion, the X register contains the module of a vector which rectangular coordinates are the abovementioned random numbers.

If the module of the vector is less than 1, a counter is incremented. The ratio between this counter and the number of iterations should converge to pi/4.

From memory, it appeared on the TI-56 manual, IIRC.

I apologize in advance for any error in this post, which I am writing just from very old memories and with little time to recheck details.

(01-11-2022 11:37 AM)Ángel Martin Wrote:  I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNG-related programs and MCODE functions. I'd like to include a Monte-Carlo based calculation of pi, which I *think* has already been published in the forum but can't locate anymore. Any pointers to that reference will be appreciated...

Cheers,
ÁM

[J-F Garnier]

(01-11-2022 11:37 AM)Ángel Martin Wrote:  I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNG-related programs and MCODE functions. I'd like to include a Monte-Carlo based calculation of pi, which I *think* has already been published in the forum but can't locate anymore. Any pointers to that reference will be appreciated...

You may have a look at the work of our friend Mike (Stgt) here (copy and edit the link).

J-F

[KeithB]

It might have been mine, which I learned at Cal Tech in the seventies at a Programming for High Schoolers program they had.

1. Generate 2 uniform random numbers between 0 and 1.
2. calculate the square root of sum of the squares. This represents the distance from the origin.
3. If the sum is <= 1 you have a hit - the point is inside the unit circle
4 calculate the ratio of the number of hits / the number of total pairs
5 multiply by 4 to get Pi.

It is cool to graph the hits and see the circle gradually fill up.

Note: I see Andres beat me to it, but I had to name drop Cal Tech.

[Csaba Tizedes]

(01-11-2022 11:37 AM)Ángel Martin Wrote:  I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNG-related programs and MCODE functions. I'd like to include a Monte-Carlo based calculation of pi, which I *think* has already been published in the forum but can't locate anymore. Any pointers to that reference will be appreciated...

Cheers,
ÁM

HP-15C version

[ttw]

(01-11-2022 02:12 PM)KeithB Wrote:  It might have been mine, which I learned at Cal Tech in the seventies at a Programming for High Schoolers program they had.

1. Generate 2 uniform random numbers between 0 and 1.
2. calculate the square root of sum of the squares. This represents the distance from the origin.
3. If the sum is <= 1 you have a hit - the point is inside the unit circle
4 calculate the ratio of the number of hits / the number of total pairs
5 multiply by 4 to get Pi.

It is cool to graph the hits and see the circle gradually fill up.

Note: I see Andres beat me to it, but I had to name drop Cal Tech.

A big speedup is possible through skipping the square root.

Another trick is to run batches of size 113.

[Ángel Martin]

Thanks to all for the good pointers and relevant information, that should give me enough to proceed with my little project.

Cheers,
ÁM

[C.Ret]

Csaba Tizedes ' Wrote:  HP-15C version

Thanks for the link, this publication trigger my curiosity and inspire me a few line of code.

(01-12-2022 10:18 AM)ttw Wrote:  A big speedup is possible through skipping the square root.

Another trick is to run batches of size 113.

Yes, speedup is possible through skipping the square root and the test as well.

Making batch of any size spare a lot of time skipping all the numerous intermediate displays.
But why a size of 113 ?

Didn't the batch of hundreds or thousands spare more ?

The following codes may be of interest for HP-15C's owners volunteers to experiment that optimization.

I still being in doubt concerning the usefulness of the rectangular to polar conversion. of course, this instruction spare a few steps and bytes in the code. But, since it also compute the angle from trigonometric functions based on CORDIC algorithm...

Taken with a dreadful doubt and in order to be clear about it; I carried out a small experiment using my HP-15C:

I entered the two codes which are similar in principle. Only difference is the use of the direct rectangular to polar conversion instruction or the computation of the sum of the two squares.
In both version, the dots are drawn at random using two RAN# instructions.

Usage:
At first use:
Please clear register (pressing f CLEAR REG ) since any use of the code add a batch of random drawn and display current approximation of PI.
Registers R0 contains the count of dots in the first quadrant circle and R1 contains the total count of dot's drawn.
The register I is only use as a counter in the main loop (based on a DSE instruction).

For each further use:
Input size of the batch and start the program by pressing the R/S key.
Eventually for the really first use or in case of an interrupted previous use, be sure to restart at the beginning by pressing g RTN or any appropriate instruction (GTO CH 000 , f CLEAR PGRM in RUN mode, or add a user-key-label to the code).


VERSION#1

Code:

001- STO+0  STO+1  STO I
004- f LBL 0
005-    f RAN#  f RAN#  g →P  g INT  STO-0  f DSE I  GTO 0
012-  4  RCL*0  RCL/1

This version is the shortest code, but a batch of size 100 need 3'04" to complete on my HP-15C

VESION#2

Code:

001- STO+0  STO+1  STO I
004- f LBL 0
005-    f RAN#  ENTER  *  f RAN#  ENTER  *  +  g INT  STO-0  f DSE I  GTO 0
016-  4  RCL*0  RCL/1

This longer version is about 31.5% faster than the shortest, since a batch of size 100 need only 2'06" to complete.

So I regret to announce to the whole community, that despite its elegance and sobriety, the idea of using the conversion →P instruction is not the most efficient and that nothing beats a good calculation conventionally placed in the operational stack.
If its true on a HP-15C, it is certainly the same on other models or brands ?

However, I still have some doubts about the effectiveness of my method and the possible existence of a bias related to my somewhat brutal way of counting ... Maybe? Peut-être ?

[ttw]

The use of multiplies 113 (and some other numbers) comes from the theory of continued fractions. Any number X may be approximated by a fraction P/Q (possibly improper). Trivially the error e=|X-P/Q| is less than 1/2Q. For each X, there exists a sequence of Qs with an error of less than 1/2Q^2. The sequence for Pi is 7, 106, 113, 33102, 33215,...

If one takes samples of size 33000, one gets at best errors of size 1/33000; with size 33215, there may be some errors of size 1/1103236225. Some samples are extraordinarily good.

Of course, in "real -ife" computation, one wouldn't know what Q to use.

[Dave Shaffer]

(01-11-2022 02:12 PM)KeithB Wrote:  Note: I see Andres beat me to it, but I had to name drop Cal Tech.

Here's another Techer (Ph.D. Astronomy 1974 - we must have been there together at some occasion!)

Anybody else out there?

[KeithB]

Sorry, I was not clear. I was a high school student at a special program that taught students like me programming on Saturdays.

[C.Ret]

(01-13-2022 01:04 AM)ttw Wrote:  The use of multiplies 113 (and some other numbers) comes from the theory of continued fractions. Any number X may be approximated by a fraction P/Q (possibly improper). Trivially the error e=|X-P/Q| is less than 1/2Q. For each X, there exists a sequence of Qs with an error of less than 1/2Q^2. The sequence for Pi is 7, 106, 113, 33102, 33215,...

If one takes samples of size 33000, one gets at best errors of size 1/33000; with size 33215, there may be some errors of size 1/1103236225. Some samples are extraordinarily good.

Of course, in "real -ife" computation, one wouldn't know what Q to use.

Thank you for this explanation, I will now have to deepen this somewhat and perfect my modest knowledge on this subject ...

In the meantime, I have found a way to shorten VERSION#2 while making it even faster.
But I won't say anything about it here, I hope that other HP-15C enthusiasts have already found how to do it.

[Gerson W. Barbosa]

Just in case the thread from 2011 hasn’t been mentioned yet:

https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=180167

[Ángel Martin]

(01-13-2022 10:10 PM)Gerson W. Barbosa Wrote:  Just in case the thread from 2011 hasn’t been mentioned yet:

https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=180167

... and here's the HP-41 version of Valentín's brilliant HP-71 code. It uses SEED and RNDM, both functions from the AMC_OS/X module

Just enter the number of points and press R/S.
Obviously very slow on a real HP-41, but very reasonable using V41 in Turbo mode.

Code:

01  LBL "MCPI"
02  STO 00
03  0
04  SEED
05  FIX 0
06  LBL 00
07  RNDM
08  RNDM
09  /
10  RND
11  2
12  MOD
13  -
14  DSE  Y
15  GTO 00
16  FIX 9
17  RCL 00
18  /
19  -4
20  *
21  1
22  +
23  END

Thanks to all again for your inputs.

Cheers,
ÁM

[C.Ret]

(01-14-2022 09:03 AM)Ángel Martin Wrote:  ... and here's the HP-41 version of Valentín's brilliant HP-71 code. It uses SEED and RNDM, both functions from the AMC_OS/X module

Just enter the number of points and press R/S.
Obviously very slow on a real HP-41, but very reasonable using V41 in Turbo mode.

Hello,

Where did you put the raw file ?

What will be the equivalent program on a HP-41C in the base configuration without any module ?

[Gerson W. Barbosa]

(01-14-2022 09:03 AM)Ángel Martin Wrote:  ... and here's the HP-41 version of Valentín's brilliant HP-71 code. It uses SEED and RNDM, both functions from the AMC_OS/X module

Just enter the number of points and press R/S.
Obviously very slow on a real HP-41, but very reasonable using V41 in Turbo mode.

Code:

01  LBL "MCPI"
02  STO 00
03  0
04  SEED
05  FIX 0
06  LBL 00
07  RNDM
08  RNDM
09  /
10  RND
11  2
12  MOD
13  -
14  DSE  Y
15  GTO 00
16  FIX 9
17  RCL 00
18  /
19  -4
20  *
21  1
22  +
23  END

The AMCOSX module is missing in i41CX+ here (probably accidentally deleted), but it runs on Free42:

Code:



00 { 41-Byte Prgm }
01▸LBL "MCPI"
02 STO 00
03 1
04 SEED
05 FIX 00
06▸LBL 00
07 RAN
08 RAN
09 ÷
10 RND
11 2
12 MOD
13 -
14 DSE ST Y
15 GTO 00
16 FIX 09
17 RCL 00
18 ÷
19 -4
20 ×
21 1
22 +
23 END

145000 XEQ "MCPI" -> 3.141627586 ( ~ 2 seconds)

BTW, here’s an alternative to Valentin’s “faster, very precise, approximate way” (to get pi). Something more accurate than Free42 might be necessary to verify it is not exactly pi, though:

ln(((4001-1/(8002–1/(4001^2+(24006-5/8002-1)/4)))/5)^6-24)/sqrt(163)

where all the constants greater than 163 are 4001 or multiples thereof.

[Joe Horn]

(01-14-2022 08:31 PM)Gerson W. Barbosa Wrote:  BTW, here’s an alternative to Valentin’s “faster, very precise, approximate way” (to get pi). Something more accurate than Free42 might be necessary to verify it is not exactly pi, though:

ln(((4001-1/(8002–1/(4001^2+(24006-5/8002-1)/4)))/5)^6-24)/sqrt(163)

where all the constants greater than 163 are 4001 or multiples thereof.

According to the 50g's LongFloat Library, this is greater than pi by approximately 2.299E-34.

[Gerson W. Barbosa]

(01-11-2022 01:11 PM)J-F Garnier Wrote:  

(01-11-2022 11:37 AM)Ángel Martin Wrote:  I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNG-related programs and MCODE functions. I'd like to include a Monte-Carlo based calculation of pi, which I *think* has already been published in the forum but can't locate anymore. Any pointers to that reference will be appreciated...

You may have a look at the work of our friend Mike (Stgt) here (copy and edit the link).

Mike (Stgt) has written another HP-41CX program that computes only one random number per iteration, but it requires his VORANOGE (Voyagers' Pseudorandom Number Generator) module. His equivalent HP-15C program might illustrate the method, though:

Code:


001- LBL A
002- STO I
003- 0
004- STO RAN#
005- FIX 0
006- LBL 1 
007- RCL RAN#
008- RAN#
009- /
010- RND
011- 2
012- /
013- FRAC
014- +
015- DSE I
016- GTO 1
017- FIX 9
018- x<>y
019- /
020- 8
021- *
022- 1
023- +
024- RTN

1356 f A -> 3.141592920 ( ~10s, HP-15C LE)

Update:

Mike is interested neither in Pi nor in the Monte Carlo method. His motivation is finding the fastest possible implementations. I think he has succeeded, judging by his solutions for both the HP-41CX plus VORANOGE module and the HP-15C:

Code:



01•LBL "PIMC"
02 STO 01
03 CLX
04 STO 00
05 RCL 01
06 ENTER^
07•LBL 01
08 RCL 00
09 RN
10 X^2
11 STO 00
12 +
13 INT
14 -
15 DSE Y
16 GTO 01
17 RCL 01
18 /
19 4
20 *
21 END

2187000 CPU cycles for input 1E3 if PR detached and VORANOGE on page 8.

Code:



  1-42,21,11  LBL A
  2-   44 25  STO I
  3-       0  0      ; no need to square this one
  4-   44  0  STO 0
  5-      34  x<>y
  6-      36  ENTER^
  7-42,21, 0  LBL 0
  8-   45  0  RCL 0  ; previous RN^2
  9-   42 36  RAN#   ; get a new RN
 10-   43 11  x^2
 11-   44  0  STO 0  ; save for future use
 12-      40  +
 13-   43 44  INT
 14-      30  -      ; hit-=1 if outside
 15-42, 5,25  DSE I
 16-   22  0  GTO 0
 17-      34  x<>y
 18-      10  /
 19-       4  4
 20-      20  *

3842552 CPU counts for input 1E3

This should run also on the HP-11C as RCL RAN# is not needed.

We can do some cheating to get an excellent result on the original HP-15C in a reasonable time by choosing proper seed and number of points:

14 STO RAN#

452 f A -> 3.141592920 (7m34s, HP-15C)

[Ángel Martin]

(01-19-2022 09:07 PM)Gerson W. Barbosa Wrote:  We can do some cheating to get an excellent result on the original HP-15C in a reasonable time by choosing proper seed and number of points:

14 STO RAN#

452 f A -> 3.141592920 (7m34s, HP-15C)

This program doesn't seem to produce more accurate results as the number of iterations gets larger (always starting with the same seed of course), which is counterintuitive to me. Do you see the same ?

ÁM