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
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
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 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
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
Code: Select all
122▸LBL 04
123 RAN
124 RCL× 00
125 IP
126 1
127 +
128 END
The computation gave me a value of π of 3.149704, i.e. an error of 0.26%, roughly.