Performance/bin/MonteCarlo.c

0001 #include "Random.h"
0002
0003 /**
0004  Estimate Pi by approximating the area of a circle.
0005
0006  How: generate N random numbers in the unit square, (0,0) to (1,1)
0007  and see how are within a radius of 1 or less, i.e.
0008  <pre>
0009
0010  sqrt(x^2 + y^2) < r
0011
0012  </pre>
0013   since the radius is 1.0, we can square both sides
0014   and avoid a sqrt() computation:
0015   <pre>
0016
0017     x^2 + y^2 <= 1.0
0018
0019   </pre>
0020   this area under the curve is (Pi * r^2)/ 4.0,
0021   and the area of the unit of square is 1.0,
0022   so Pi can be approximated by
0023   <pre>
0024                 # points with x^2+y^2 < 1
0025      Pi =~      --------------------------  * 4.0
0026                      total # points
0027
0028   </pre>
0029
0030 */
0031
0032 static const int SEED = 113;
0033
0034 double MonteCarlo_num_flops(int Num_samples) {
0035   /* 3 flops in x^2+y^2 and 1 flop in random routine */
0036
0037   return ((double)Num_samples) * 4.0;
0038 }
0039
0040 double MonteCarlo_integrate(int Num_samples) {
0041   Random R = new_Random_seed(SEED);
0042
0043   int under_curve = 0;
0044   int count;
0045
0046   for (count = 0; count < Num_samples; count++) {
0047     double x = Random_nextDouble(R);
0048     double y = Random_nextDouble(R);
0049
0050     if (x * x + y * y <= 1.0)
0051       under_curve++;
0052   }
0053
0054   Random_delete(R);
0055
0056   return ((double)under_curve / Num_samples) * 4.0;
0057 }