Using random numbers to compute π

Contributions to this software library are always welcome. Please ensure that you post program listings rather than .raw files. They give a reasonable idea of what your program does without having to load them into a DM42 and you can also include comments in your code. Check out the following link for a decoder/encoder: https://technical.swissmicros.com/decoders/dm42/

You can then copy/paste the listing and post it in "code" tags.
Post Reply
grsbanks
Posts: 1122
Joined: Tue Apr 25, 2017 11:23 am
Location: Preston, Lancs, UK
Contact:

Using random numbers to compute π

Post by grsbanks »

Here's a little bit of maths trivia that can be fun to explore.

As n→∞, the probability that two natural numbers between 1 and n are coprime (i.e. that they have no factors in common other than 1) approaches 6/π². I didn't just make this up, it was established by Euler in 1735 and proven in 1741: https://en.wikipedia.org/wiki/Coprime_i ... oprimality

So, what's to stop us from sampling pairs of natural numbers, counting how many of them are coprime and examining the proportion of coprime pairs to the total number of pairs? I'll call this proportion 'p'.

In theory, p should be equal to 6/π², so by reorganising that a bit, if we take sqrt(6/p) we should end up with π or something close to it.

A few caveats first. Firstly, this only gets close to reality for values of n (the highest random natural number that we work with) that approach infinity. However, we get an error of only about 0.3% if we stick to n=200 or about 0.06% for n=1000 (see the "Basel" problem, compute the sum of the reciprocals of the squares of x for x ranging from 1 to n and compare that to π²/6). It's not going to take any longer to generate random numbers up to 1000 than it is random numbers up to 200, so I'll stick with 1000 for the purposes of this experiment.

Secondly, for the sample of data points to be significant, it must be as large as possible. The DM42 has plenty of computing power under the hood, especially if hooked up to USB power, so sample sizes of 1000, 10000 and even 100000 are not a problem. The higher this number, the better, but the longer it'll take to run.

Thirdly, the random number generator in Free42 is already pretty darn good but it'll never be perfect, so that will introduce a few errors as well.

Given the three sources of errors above, I'll be very happy if I get results within one percent, i.e. if it is able to give me a value of π somewhere between, say, 3.1 and 3.2.

Now, that really great screen on the DM42 should not be allowed to go to waste while this is running, so I added instructions to plot the value p computed after each iteration of the loop.

The program is attached as a .raw file, but here's the code with a few comments on its operation.

First, here's the initial setup. It displays its current setup and provides a menu for the user to specify the parameters of the experiment.

Code: Select all

00 { 267-Byte Prgm }
01▸LBL "PI2"
02 ALL
03 6
04 1
05 NEWMAT
06 LSTO "REGS"
07 6
08 PI
09 X↑2
10 ÷
11 STO 02
12 CLST
13 CF 29
14▸LBL 10
15 "HIGHEST: "
16 ARCL 00
17 ├"[LF]SAMP SZ: "
18 ARCL 01
19 AVIEW
20 CLMENU
21 "HIGH"
22 KEY 1 GTO 21
23 "SAMP"
24 KEY 2 GTO 22
25 KEY 9 GTO 12
26 RCL 00
27 100
28 X>Y?
29 GTO 11
30 RCL 01
31 1ᴇ3
32 X>Y?
33 GTO 11
34 "GO→"
35 KEY 6 GTO 30
36▸LBL 11
37 MENU
38 STOP
39 GTO 10
40▸LBL 12
41 EXITALL
42 RTN
43▸LBL 21
44 ABS
45 IP
46 STO 00
47 GTO 10
48▸LBL 22
49 ABS
50 IP
51 STO 01
52 GTO 10
No doubt you'll notice that the "GO→" menu entry is only displayed if the highest random number to use ('n') is 100 or greater and if the sample size is at least 1000. But if it's there and you hit it, we go off to LBL 30 and start the run.

Code: Select all

53▸LBL 30
54 RCL 01
55 STO 03
56 3
57 STO "GrMod"
58 CLLCD
59 -120
60 1
61 PIXEL
This initialises the counter for the number of iterations to perform, sets the display to the highest resolution of 400x240 and draws a horizontal line halfway down the screen. This line shows the level that we want the proportion 'p' to converge on during the course of the test run.

The main loop runs here:

Code: Select all

62▸LBL 06
63 XEQ 04
64 XEQ 04
65 XEQ 01
66 1
67 X≠Y?
68 GTO 07
69 STO+ 04
70▸LBL 07
71 RCL 01
72 RCL- 03
73 1
74 +
75 STO ST Y
76 RCL 04
77 X<>Y
78 ÷
79 STO 05
80 RCL- 02
81 20
82 ×
83 120
84 ×
85 RCL+ ST L
86 X<0?
87 GTO 08
88 X<>Y
89 RCL÷ 01
90 400
91 ×
92 PIXEL
93▸LBL 08
94 DSE 03
95 GTO 06
The subroutine at LBL 04 just generates a random number from 1 to n.

The subroutine at LBL 01 computes the greatest common divisor of the numbers in X and Y using the Euclidean algorithm. If that number is 1 then we know that the two random numbers generated were coprime and we increment the number of coprime pairs found. We then compute the new 'p', apply some scaling, plot the point on the graph and go back for the next iteration.

When everything is finished, we compute π from the proportion 'p' that we have and then calculate the error between that and the real value of π as a percentage. Our π computed from 'p' is left in Y and the error percentage in X.

Code: Select all

96 RCL 05
97 1/X
98 6
99 ×
100 SQRT
101 PI
102 RCL- ST Y
103 RCL÷ ST L
104 ABS
105 100
106 ×
107 FIX 04
108 GTO 12
This is the subroutine that calculates the gcd of X and Y:

Code: Select all

109▸LBL 01
110 X>Y?
111 X<>Y
112▸LBL 02
113 MOD
114 X=0?
115 GTO 03
116 LASTX
117 X<>Y
118 GTO 02
119▸LBL 03
120 LASTX
121 RTN
And finally, this generates the random number from 1 to 'n':

Code: Select all

122▸LBL 04
123 RAN
124 RCL× 00
125 IP
126 1
127 +
128 END
This graph was generated using 5000 pairs of random numbers from 1 to 1000:

Image

The computation gave me a value of π of 3.149704, i.e. an error of 0.26%, roughly.

Image
Attachments
pi_2.raw
(270 Bytes) Downloaded 159 times
There are only 10 kinds of people in the world: those who understand binary and those who do not.
cdmackay
Posts: 281
Joined: Fri Oct 05, 2018 8:33 pm
Location: Cambridge, UK
Contact:

Re: Using random numbers to compute π

Post by cdmackay »

nice! love the plot.
Cambridge, UK
41CL/DM41X 12/15C/16C DM15/16 17B/II/II+ 28S 42S/DM42 32SII 48GX 50g 35s WP34S PrimeG2 WP43S/pilot
Casio, Rockwell 18R
Post Reply