Monte-Carlo Pi
01-11-2022, 11:37 AM (This post was last modified: 01-11-2022 11:38 AM by Ángel Martin.)
Post: #1
 Ángel Martin Senior Member Posts: 1,302 Joined: Dec 2013
Monte-Carlo Pi
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
01-11-2022, 12:19 PM
Post: #2
 Albert Chan Senior Member Posts: 1,791 Joined: Jul 2018
RE: Monte-Carlo Pi
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

01-11-2022, 12:51 PM
Post: #3
 Andres Junior Member Posts: 49 Joined: Jan 2014
RE: Monte-Carlo Pi
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

Andrés C. Rodríguez (Argentina)

My posts are mostly from old memories, not from current research.
01-11-2022, 01:11 PM
Post: #4
 J-F Garnier Senior Member Posts: 619 Joined: Dec 2013
RE: Monte-Carlo Pi
(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
01-11-2022, 02:12 PM (This post was last modified: 01-11-2022 02:25 PM by KeithB.)
Post: #5
 KeithB Member Posts: 298 Joined: Jan 2017
RE: Monte-Carlo Pi
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.
01-11-2022, 05:40 PM
Post: #6
 Csaba Tizedes Senior Member Posts: 500 Joined: May 2014
RE: Monte-Carlo Pi
(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
01-12-2022, 10:18 AM
Post: #7
 ttw Member Posts: 248 Joined: Jun 2014
RE: Monte-Carlo Pi
(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.
01-12-2022, 01:50 PM
Post: #8
 Ángel Martin Senior Member Posts: 1,302 Joined: Dec 2013
RE: Monte-Carlo Pi
Thanks to all for the good pointers and relevant information, that should give me enough to proceed with my little project.

Cheers,
ÁM
01-12-2022, 05:33 PM (This post was last modified: 01-12-2022 05:46 PM by C.Ret.)
Post: #9
 C.Ret Member Posts: 155 Joined: Dec 2013
RE: Monte-Carlo Pi
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 ?
01-13-2022, 01:04 AM
Post: #10
 ttw Member Posts: 248 Joined: Jun 2014
RE: Monte-Carlo Pi
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.
01-13-2022, 03:30 AM
Post: #11
 Dave Shaffer Member Posts: 57 Joined: Dec 2013
RE: Monte-Carlo Pi
(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?
01-13-2022, 01:58 PM
Post: #12
 KeithB Member Posts: 298 Joined: Jan 2017
RE: Monte-Carlo Pi
Sorry, I was not clear. I was a high school student at a special program that taught students like me programming on Saturdays.
01-13-2022, 06:14 PM
Post: #13
 C.Ret Member Posts: 155 Joined: Dec 2013
RE: Monte-Carlo Pi
(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.
01-13-2022, 10:10 PM
Post: #14
 Gerson W. Barbosa Senior Member Posts: 1,452 Joined: Dec 2013
RE: Monte-Carlo Pi
Just in case the thread from 2011 hasn’t been mentioned yet:

01-14-2022, 09:03 AM (This post was last modified: 01-14-2022 09:05 AM by Ángel Martin.)
Post: #15
 Ángel Martin Senior Member Posts: 1,302 Joined: Dec 2013
RE: Monte-Carlo Pi
(01-13-2022 10:10 PM)Gerson W. Barbosa Wrote:  Just in case the thread from 2011 hasn’t been mentioned yet:

... 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
01-14-2022, 05:53 PM
Post: #16
 C.Ret Member Posts: 155 Joined: Dec 2013
RE: Monte-Carlo Pi
(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 ?
01-14-2022, 08:31 PM
Post: #17
 Gerson W. Barbosa Senior Member Posts: 1,452 Joined: Dec 2013
RE: Monte-Carlo Pi
(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.
01-15-2022, 05:52 AM
Post: #18
 Joe Horn Senior Member Posts: 1,824 Joined: Dec 2013
RE: Monte-Carlo Pi
(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.

<0|ɸ|0>
-Joe-
01-19-2022, 09:07 PM (This post was last modified: 01-22-2022 10:03 AM by Gerson W. Barbosa.)
Post: #19
 Gerson W. Barbosa Senior Member Posts: 1,452 Joined: Dec 2013
RE: Monte-Carlo Pi
(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)
01-29-2022, 10:37 AM (This post was last modified: 01-29-2022 10:38 AM by Ángel Martin.)
Post: #20
 Ángel Martin Senior Member Posts: 1,302 Joined: Dec 2013
RE: Monte-Carlo Pi
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)