/*     A procedure that fills an array Z with nZ standard normals using the Box Muller method.
       Written by Jonathan Goodman for assignment 8 in the class 
	   http://www.math.nyu.edu/faculty/goodman/teaching/DerivSec09/index.html
	   
	   It uses the native C random number generator rand(), which may not be that good.
	   Feel free to replace this with a better one, but remember to set the seed.                */
    

#include <iostream>
#include <math.h>

using namespace std;

int randn(                // The return value will be 0 if it worked, 1 if it did not.
           double Z[],    // The array of numbers to return, filled with independent standard normals
           int    nZ ) {  // The size of the array Z, must be allocated in the calling procedure.

   int    i;                     // The array index, not the loop index -- see below.
   int    nZBy2;                 // nZ/2 if nZ is even, (nZ-1)/2 if nZ is odd.
   double pi = 3.1415926535;     // Duh
   double X, R, T;
   
   if ( nZ <= 0 ) {
      cout << "randn quitting because nZ = " << nZ << endl;
	  return 1;  }
   
   nZBy2 = nZ/2;
   
   for ( int l = 0; l < 10; l++) rand();         //  Discard the first few random numbers, which can be bad.
   
   i = 0;
   for ( int k = 0; k < nZBy2; k++) {              // Box Muller makes normals two at a time.
      R      = ( (double) rand() )/RAND_MAX;
	  R      = sqrt( -2*log( R ) );
	  T      = 2*pi*( (double) rand() )/RAND_MAX;
	  Z[i++] = R * cos(T);
	  Z[i++] = R * sin(T);
     }
   if ( 2*nZBy2 < nZ ) {                       // Fill in the last value if nZ is odd.
	  R      = ( (double) rand() )/RAND_MAX;
	  R      = sqrt( -2*log( R ) );
	  T      = 2*pi*( (double) rand() )/RAND_MAX;
	  Z[i]   = R * cos(T);                        // Don't need i++ because this is the last time.
	}
	  
	     
   return 0;       // It worked, hopefully.
}