Post Reply 
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
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
Find all posts by this user
Quote this message in a reply
01-11-2022, 12:19 PM
Post: #2
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

see https://www.hpmuseum.org/forum/thread-13...#pid129470
Find all posts by this user
Quote this message in a reply
01-11-2022, 12:51 PM
Post: #3
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)

Please disregard idiomatic mistakes.
My posts are mostly from old memories, not from current research.
Find all posts by this user
Quote this message in a reply
01-11-2022, 01:11 PM
Post: #4
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
Visit this user's website Find all posts by this user
Quote this message in a reply
01-11-2022, 02:12 PM (This post was last modified: 01-11-2022 02:25 PM by KeithB.)
Post: #5
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.
Find all posts by this user
Quote this message in a reply
01-11-2022, 05:40 PM
Post: #6
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
Find all posts by this user
Quote this message in a reply
01-12-2022, 10:18 AM
Post: #7
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.
Find all posts by this user
Quote this message in a reply
01-12-2022, 01:50 PM
Post: #8
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
Find all posts by this user
Quote this message in a reply
01-12-2022, 05:33 PM (This post was last modified: 01-12-2022 05:46 PM by C.Ret.)
Post: #9
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 ?
Find all posts by this user
Quote this message in a reply
01-13-2022, 01:04 AM
Post: #10
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.
Find all posts by this user
Quote this message in a reply
01-13-2022, 03:30 AM
Post: #11
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?
Find all posts by this user
Quote this message in a reply
01-13-2022, 01:58 PM
Post: #12
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.
Find all posts by this user
Quote this message in a reply
01-13-2022, 06:14 PM
Post: #13
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.
Find all posts by this user
Quote this message in a reply
01-13-2022, 10:10 PM
Post: #14
RE: Monte-Carlo Pi
Just in case the thread from 2011 hasn’t been mentioned yet:

https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=180167
Find all posts by this user
Quote this message in a reply
01-14-2022, 09:03 AM (This post was last modified: 01-14-2022 09:05 AM by Ángel Martin.)
Post: #15
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:

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
Find all posts by this user
Quote this message in a reply
01-14-2022, 05:53 PM
Post: #16
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 ?
Find all posts by this user
Quote this message in a reply
01-14-2022, 08:31 PM
Post: #17
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.
Find all posts by this user
Quote this message in a reply
01-15-2022, 05:52 AM
Post: #18
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-
Visit this user's website Find all posts by this user
Quote this message in a reply
01-19-2022, 09:07 PM (This post was last modified: 01-22-2022 10:03 AM by Gerson W. Barbosa.)
Post: #19
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)
Find all posts by this user
Quote this message in a reply
01-29-2022, 10:37 AM (This post was last modified: 01-29-2022 10:38 AM by Ángel Martin.)
Post: #20
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
Find all posts by this user
Quote this message in a reply
Post Reply 




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