1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
|
#include <iostream>
#include <iomanip>
#include <chrono>
#include <random>
using namespace std;
unsigned seed = chrono::system_clock::now().time_since_epoch().count();
default_random_engine generator(seed);
uniform_real_distribution<double> distribution( 0.0, 1.0 );
double myRandom() // needs to return a uniformly-distributed value on [0,1]
{
double x;
x = distribution( generator );
return x;
}
double getPi( long N ) // from ratio of area of a circle to the unit square
{
long ntrue = 0;
double x, y;
for ( long i = 0; i < N; i++ )
{
x = myRandom();
y = myRandom();
if ( x * x + y * y <= 1.0 ) ntrue++;
}
return 4.0 * ntrue / N;
}
double getPiByVolume( long N ) // from ratio of volume of a sphere to a unit cube
{
long ntrue = 0;
double x, y, z;
for ( long i = 0; i < N; i++ )
{
x = myRandom();
y = myRandom();
z = myRandom();
if ( x * x + y * y + z * z <= 1.0 ) ntrue++;
}
return 6.0 * ntrue / N;
}
double getLog2( long N ) // from area under curve y=1/x (i.e. xy=1) from x = 1 to 2
{
long ntrue = 0;
double x, y;
for ( long i = 0; i < N; i++ )
{
x = myRandom() + 1;
y = myRandom();
if ( x * y <= 1.0 ) ntrue++;
}
return ( ntrue + 0.5 ) / N;
}
double getE( long N ) // see Wikipedia: I don't get it though!
{
double sum;
long long total = 0;
for ( long i = 0; i < N; i++ )
{
sum = 0;
while ( sum <= 1 )
{
sum += myRandom();
total++;
}
}
return ( total + 0.5 ) / N;
}
double getRoot2( long N ) // from pdf; sqrt(A) = A * Prob( X * X < 1/A )
{
long ntrue = 0;
double x;
for ( long i = 0; i < N; i++ )
{
x = myRandom();
if ( x * x <= 0.5 ) ntrue++;
}
return 2.0 * ntrue / N;
}
double getPhi( long N ) // from sqrt(5) as above
{
long ntrue = 0;
double x;
for ( long i = 0; i < N; i++ )
{
x = myRandom();
if ( x * x <= 0.2 ) ntrue++;
}
return 0.5 + 2.5 * ntrue / N;
}
void printIt( string target, double (*f)( long ), long N, string actual )
{
cout << setw(11) << target << setw(12) << setprecision(5) << f(N) << " ( " << actual << " )" << endl;
}
int main()
{
long N;
cout << "Input number of random points: "; cin >> N;
cout << "\nMonte-Carlo estimates for: " << endl;
printIt( "pi:" , getPi , N, "3.14159..." );
printIt( "pi(again):", getPiByVolume, N, "3.14159..." );
printIt( "ln(2):" , getLog2 , N, "0.693147.." );
printIt( "e:" , getE , N, "2.71828..." );
printIt( "sqrt(2):" , getRoot2 , N, "1.41421..." );
printIt( "phi:" , getPhi , N, "1.61803..." );
}
|