Various Monte Carlo integration results due to exit

Question

I just stumbled upon a C ++ error / function that I cannot fully understand, and was hoping that someone who knows C ++ better can point me in the right direction.

Below you will find my attempt to find out the area under the Gaussian curve using Monte Carlo integration. Recipe:

• Create a large number of normally distributed random variables (with an average value of zero and a standard deviation of one).
• Dial these numbers.
• Take the average of all squares. The average will be a very close estimate of the area under the curve (1.0 in the case of the Gaussian).

The code below consists of two simple functions: rand_uni , which returns a random variable evenly distributed between zero and one, and rand_norm , which is (rather poor, but "sufficient for the government to work") an approximation of a normally distributed random variable.

main runs through the loop one billion times, each time calling rand_norm , squaring it with pow and adding it to the cumulative variable. After this cycle, the accumulated result is simply divided by the number of runs and printed on the terminal as Result=<SOME NUMBER> .

The problem is the very bizarre behavior of the code below: when each generated random variable is printed on cout (yes, one billion times), the final result is correct regardless of the compiler used (1.0015, which is pretty close to what I want). If I do not print a random variable in each iteration of the loop, I get inf under gcc and 448314 under clang .

Honestly, it just amazes the mind, and since this is my first meeting with C ++, I really don't know what the problem is: is it something with pow ? cout acting weird?

Any hints would be greatly appreciated!

A source

 // Monte Carlo integration of the Gaussian curve #include <iostream> #include <cstdlib> #include <cmath> using namespace std; enum { no_of_runs = 1000000 }; // uniform random variable double rand_uni() { return ((double) rand() / (RAND_MAX)); }; // approximation of a normaly distributed random variable double rand_norm() { double result; for(int i=12; i > 0; --i) { result += rand_uni(); } return result - 6; }; int main(const int argc, const char** argv) { double result = 0; double x; for (long i=no_of_runs; i > 0; --i) { x = pow(rand_norm(), 2); #ifdef DO_THE_WEIRD_THING cout << x << endl; // MAGIC?! #endif result += x; } // Prints the end result cout << "Result=" << result / no_of_runs << endl << endl; } 

Makefile

 CLANG=clang++ GCC=g++ OUT=normal_mc default: *.cpp $(CLANG) -o $(OUT).clang.a *.cpp $(CLANG) -o $(OUT).clang.b -DDO_THE_WEIRD_THING *.cpp $(GCC) -o $(OUT).gcc.a *.cpp $(GCC) -o $(OUT).gcc.b -DDO_THE_WEIRD_THING *.cpp 

Answer 1

result in double rand_norm() not initialized 0 to the beginning of + = random values.

All cout use and reset some of the stack memory, which is later used by calling rand_norm, "fixing" your problem when you execute cout.

Change the first line to double result = 0.0; to fix it.

 double rand_norm() { double result = 0.0; for(int i=12; i > 0; --i) { result += rand_uni(); } return result - 6; }; 

In addition, you have to sow your random number generator, you will always get exactly the same sequence of pseudo random numbers without it, add

 srand(time(NULL)); 

before the first call to rand (), or check out some of the best random number generators, for example. Boost.random

Answer 2

In your rand_norm () function, the result is not initialized, this is your error. Just in case, you are wondering why this worked when you printed the values, because the unified variable will have the value that it has in the memory area, it is a stack in your code, so calling cout <<changed your stack values. Here is a fixed version of your function:

 double rand_norm() { double result = 0.0; for(int i=12; i > 0; --i) { result += rand_uni(); } return result - 6; };