#include
#include
#include
#include
#include //Random number stuff
using namespace std;
unsigned int MonteCarlo(double (*)(double), double (*)(double), double (*)(double), double, double, double, double&, double&, double&); //Monte Carlo method
inline double Integrand(double); //The function to be integrated
inline double Inverse_CDFunc(double); //Returns the Inverse CDF of the random number distribution
inline double PDFunc(double); //Returns the PDF of the random number distribution
inline double RandF(double (*)(double));//Returns a random number according to the distribution provided
inline double Rand(); //Returns a uniformally distributed random number between 0 and 1
int main()
{
srand(time(NULL)); //Seed the random number generator
double Value;
double RMS;
double Error;
unsigned int Number;
Number = MonteCarlo(Integrand, PDFunc, Inverse_CDFunc, 0., 1., .0001, Value, RMS, Error);
cout << "Value = " << Value << endl; //Output information for use by the user.
cout << "RMS = " << RMS << endl;
cout << "Error = " << Error << endl;
cout << "Number = " << Number << endl;
return(0);
}
unsigned int MonteCarlo(double (*Func)(double), double (*PDF)(double), double (*CDF)(double), double a, double b, double Tolerance, double& Value, double& RMS, double& Error) //Monte Carlo method
{
unsigned long int i; //Counter
double Entry; //Sample
double First = 0; //Expectation of the first moment
double Second = 0; //Expectation of the second moment
double Random; //Random number holder
double* FirstPtr; //Pointer to the array where each thread will place it's own answer to First
double* SecondPtr; //Pointer to the array where each thread will place it's own answer to Second
unsigned int* iPtr; //Pointer to the array where each thread will place it's own number samples, smaller than i so that i can hold it
#pragma omp parallel private(i, Entry, First, Second, Random)
{
#pragma omp master
{
Tolerance *= sqrt(double(omp_get_num_threads()));//Decrease the tolerance so each thread doesn't go too far so that the combined effort is tolereance
FirstPtr = new double[omp_get_num_threads()]; //Create new First array
SecondPtr = new double[omp_get_num_threads()]; //Create new Second array
iPtr = new unsigned int[omp_get_num_threads()]; //Create new i array
}
#pragma omp barrier //To prevent the threads other than the master from continuing on without the master finishing it's task of varible prep
for(i = 0; i < UINT_MAX; i++) //Take up to UINT_MAX (Maximum value of unsigned interger) samples
{
Random = RandF(CDF); //Create a random number distributed according to the distribution function
if(Random >= a && Random <= b) //Normal order
{
Entry = Func(Random)/PDF(Random); //Find the value of F(x) and normalize for the distribution of random numbers
First = (double(i)*First+Entry)/(i+1.); //Calculate the running average of the first moment
Second = (double(i)*Second+pow(Entry,2))/(i+1.);//Calculate the running average of the second moment
}
else if(Random >= b && Random <= a) //Reverse order
{
Entry = -Func(Random)/PDF(Random); //Find the value of F(x) and normalize for the distribution of random numbers
First = (double(i)*First+Entry)/(i+1.); //Calculate the running average of the first moment
Second = (double(i)*Second+pow(Entry,2))/(i+1.);//Calculate the running average of the second moment
}
else //Average in a zero for out of range
{
First *= double(i)/(i+1.); //Calculate the running average of the first moment
Second *= double(i)/(i+1.); //Calculate the running average of the second moment
}
if((Second-pow(First,2)) > 0 && sqrt((Second-pow(First,2))/double(i))/abs(First) < Tolerance) //Early break conditions
break;
}
FirstPtr[omp_get_thread_num()] = First; //Store First, Second, and i of the respective threads in an array where the master thread can get to it
SecondPtr[omp_get_thread_num()] = Second;
iPtr[omp_get_thread_num()] = i;
#pragma omp barrier //Allows all threads to complete their repective tasks
#pragma omp master //Ensure that the master thread alone averages results over the threads
{
cout << "Thread\tFirst\t\tSecond\t\tSamples" << endl;
cout << "0\t" << First << "\t" << Second << "\t" << i << endl;
for(int j = 1; j < omp_get_num_threads(); j++)
{
cout << j << "\t" << FirstPtr[j] << "\t" << SecondPtr[j] << "\t" << iPtr[j] << endl;
First = (double(j)*First+FirstPtr[j])/(j+1.); //Average all Firsts
Second = (double(j)*Second+SecondPtr[j])/(j+1.); //Average all Seconds
i += iPtr[j]; //Sum all i's
}
FirstPtr[0] = First;
SecondPtr[0] = Second;
iPtr[0] = i;
}
}
Value = FirstPtr[0]; //To make sure that the last value has been placed
RMS = sqrt((SecondPtr[0]-pow(FirstPtr[0],2)));
Error = sqrt((SecondPtr[0]-pow(FirstPtr[0],2))/double(iPtr[0]))/abs(FirstPtr[0]);
return(iPtr[0]);
}
inline double Integrand(double Arg)
{
return(.01/(M_PI*(pow(Arg-.5,2)+.0001)));
}
inline double PDFunc(double Arg)
{
return(.01/(M_PI*(pow(Arg-.5,2)+pow(.01,2))));
}
inline double Inverse_CDFunc(double x)
{
return(.5+.01*tan(M_PI*(x-.5)));
}
inline double RandF(double (*CDF)(double))
{
return(CDF(Rand()));
}
inline double Rand()
{
return(rand()/double(RAND_MAX));
}