MonteCarlo Pi

01112022, 11:37 AM
(This post was last modified: 01112022 11:38 AM by Ángel Martin.)
Post: #1




MonteCarlo Pi
I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNGrelated programs and MCODE functions. I'd like to include a MonteCarlo 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 

01112022, 12:19 PM
Post: #2




RE: MonteCarlo Pi
Probability of random()/random() rounded to even integer = 1 + 1/4  PI/4 = (5PI)/4 ≈ 0.4646
Or, we can count odd's: Probability of random()/random() rounded to odd integer = (PI1)/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/thread13...#pid129470 

01112022, 12:51 PM
Post: #3




RE: MonteCarlo 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 TI56 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. (01112022 11:37 AM)Ángel Martin Wrote: I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNGrelated programs and MCODE functions. I'd like to include a MonteCarlo 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... Andrés C. Rodríguez (Argentina) Please disregard idiomatic mistakes. My posts are mostly from old memories, not from current research. 

01112022, 01:11 PM
Post: #4




RE: MonteCarlo Pi
(01112022 11:37 AM)Ángel Martin Wrote: I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNGrelated programs and MCODE functions. I'd like to include a MonteCarlo 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). JF 

01112022, 02:12 PM
(This post was last modified: 01112022 02:25 PM by KeithB.)
Post: #5




RE: MonteCarlo 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. 

01112022, 05:40 PM
Post: #6




RE: MonteCarlo Pi
(01112022 11:37 AM)Ángel Martin Wrote: I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNGrelated programs and MCODE functions. I'd like to include a MonteCarlo 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... HP15C version 

01122022, 10:18 AM
Post: #7




RE: MonteCarlo Pi
(01112022 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. A big speedup is possible through skipping the square root. Another trick is to run batches of size 113. 

01122022, 01:50 PM
Post: #8




RE: MonteCarlo Pi
Thanks to all for the good pointers and relevant information, that should give me enough to proceed with my little project.
Cheers, ÁM 

01122022, 05:33 PM
(This post was last modified: 01122022 05:46 PM by C.Ret.)
Post: #9




RE: MonteCarlo Pi
Csaba Tizedes ' Wrote: HP15C version Thanks for the link, this publication trigger my curiosity and inspire me a few line of code. (01122022 10:18 AM)ttw Wrote: A big speedup is possible through skipping the square root. 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 HP15C'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 HP15C: 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 userkeylabel to the code). VERSION#1 Code: 001 STO+0 STO+1 STO I VESION#2 Code: 001 STO+0 STO+1 STO I 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 HP15C, 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 ? 

01132022, 01:04 AM
Post: #10




RE: MonteCarlo 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=XP/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. 

01132022, 03:30 AM
Post: #11




RE: MonteCarlo Pi  
01132022, 01:58 PM
Post: #12




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


01132022, 06:14 PM
Post: #13




RE: MonteCarlo Pi
(01132022 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=XP/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,... 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 HP15C enthusiasts have already found how to do it. 

01132022, 10:10 PM
Post: #14




RE: MonteCarlo Pi
Just in case the thread from 2011 hasn’t been mentioned yet:
https://www.hpmuseum.org/cgisys/cgiwrap...ead=180167 

01142022, 09:03 AM
(This post was last modified: 01142022 09:05 AM by Ángel Martin.)
Post: #15




RE: MonteCarlo Pi
(01132022 10:10 PM)Gerson W. Barbosa Wrote: Just in case the thread from 2011 hasn’t been mentioned yet: ... and here's the HP41 version of Valentín's brilliant HP71 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 HP41, but very reasonable using V41 in Turbo mode. Code: 01 LBL "MCPI" Thanks to all again for your inputs. Cheers, ÁM 

01142022, 05:53 PM
Post: #16




RE: MonteCarlo Pi
(01142022 09:03 AM)Ángel Martin Wrote: ... and here's the HP41 version of Valentín's brilliant HP71 code. It uses SEED and RNDM, both functions from the AMC_OS/X module Hello, Where did you put the raw file ? What will be the equivalent program on a HP41C in the base configuration without any module ? 

01142022, 08:31 PM
Post: #17




RE: MonteCarlo Pi
(01142022 09:03 AM)Ángel Martin Wrote: ... and here's the HP41 version of Valentín's brilliant HP71 code. It uses SEED and RNDM, both functions from the AMC_OS/X module The AMCOSX module is missing in i41CX+ here (probably accidentally deleted), but it runs on Free42: Code:
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(((40011/(8002–1/(4001^2+(240065/80021)/4)))/5)^624)/sqrt(163) where all the constants greater than 163 are 4001 or multiples thereof. 

01152022, 05:52 AM
Post: #18




RE: MonteCarlo Pi
(01142022 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: According to the 50g's LongFloat Library, this is greater than pi by approximately 2.299E34. <0ɸ0> Joe 

01192022, 09:07 PM
(This post was last modified: 01222022 10:03 AM by Gerson W. Barbosa.)
Post: #19




RE: MonteCarlo Pi
(01112022 01:11 PM)JF Garnier Wrote:(01112022 11:37 AM)Ángel Martin Wrote: I'm putting together a (very humble) RANDOM ROM with a (biased) selection of RNGrelated programs and MCODE functions. I'd like to include a MonteCarlo 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... Mike (Stgt) has written another HP41CX program that computes only one random number per iteration, but it requires his VORANOGE (Voyagers' Pseudorandom Number Generator) module. His equivalent HP15C program might illustrate the method, though: Code:
1356 f A > 3.141592920 ( ~10s, HP15C 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 HP41CX plus VORANOGE module and the HP15C: Code:
Code:
This should run also on the HP11C as RCL RAN# is not needed. We can do some cheating to get an excellent result on the original HP15C in a reasonable time by choosing proper seed and number of points: 14 STO RAN# 452 f A > 3.141592920 (7m34s, HP15C) 

01292022, 10:37 AM
(This post was last modified: 01292022 10:38 AM by Ángel Martin.)
Post: #20




RE: MonteCarlo Pi
(01192022 09:07 PM)Gerson W. Barbosa Wrote: We can do some cheating to get an excellent result on the original HP15C in a reasonable time by choosing proper seed and number of points: 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)