/*     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.              
	   
	   
	A revised version that uses the STL vector class template
	
	*/
    

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

using namespace std;

int randn(                      // The return value will be 0 if it worked, 1 if it did not.
           vector<double> & Z,  // The array of numbers to return, filled with independent standard normals
           int    nZ ) {        // The number of normals requested, 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.
   const double pi = 3.1415926535;     // Duh
   double       X, R, T;
   int         iU;                     // Hold the integer random number from rand() to test for zero.
   
   if ( nZ <= 0 ) {
      cout << "randn quitting because nZ = " << nZ << endl;
	  return 1;  }
   if ( nZ > Z.size() ) {          // 
      cout << "randn quitting because the vector Z is too small" << endl;
	  cout << "   nZ = " << nZ << ", Z size = " << Z.size() << 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.
      iU     = rand();
      if ( iU == 0 ) iU = 1;                       // rand() should never give zero, but tell that to Microsoft.
      R      = ( (double) iU )/RAND_MAX;
	     R      = sqrt( -2*log( R ));                 // rand()+1 avoids the possibility that rand() returns a zero.
	     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 ) + 1 );
	  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.
}