🚀 go-pugleaf

Gnu g++ v2.3.2 for Sparc-Sun-SunOS 4.1.1.  libg++ v2.3.  

In Normal.cc, I am not sure that the given algorithm gives an exactly
normally distributed pair of results.  I have included a modification
implementing the Box-Muller algorithm, which is exact, at the cost of
a sin/cos evaluation.  

James F. Carter        (213) 825-2897
UCLA-Mathnet;  6221 MSA; 405 Hilgard Ave.; Los Angeles, CA, USA  90024-1555
Internet: jimc@math.ucla.edu            BITNET: jimc%math.ucla.edu@INTERBIT
UUCP:...!{ucsd,ames,ncar,gatech,purdue,rutgers,decvax,uunet}!math.ucla.edu!jimc

----------- From libg++.2.3/libg++/src/Normal.cc

double Normal::operator()()
{
    
    if (haveCachedNormal == 1) {
	haveCachedNormal = 0;
	return(cachedNormal * pStdDev + pMean );
    } else {
	
	for(;;) {
	    double u1 = pGenerator -> asDouble();
	    double u2 = pGenerator -> asDouble();
	    double v1 = 2 * u1 - 1;
	    double v2 = 2 * u2 - 1;
	    double w = (v1 * v1) + (v2 * v2);
	    
//
//	We actually generate two IID normal distribution variables.
//	We cache the one & return the other.
// 
	    if (w <= 1) {
		double y = sqrt( (-2 * log(w)) / w);
		double x1 = v1 * y;
		double x2 = v2 * y;
		
		haveCachedNormal = 1;
		cachedNormal = x2;
		return(x1 * pStdDev + pMean);
	    }
	}
    }
}

------------ Suggested change

// Reference: Ripley, Brian D.  Stochastic Simulation.  Wiley (NY) 1987.
// p. 54, algorithm 3.1, Box-Muller (not Mu\"ller).
double Normal::operator()()
{
    
    if (haveCachedNormal == 1) {
	haveCachedNormal = 0;
	return(cachedNormal * pStdDev + pMean );
    } else {

//	Get 2 flat random variables in 0..1.  Generator gives 0 <= u1 < 1.
	double u1 = pGenerator -> asDouble();
//	This one must never be zero.
	double u2;
	while (0.0 == (u2 = pGenerator -> asDouble())) /*null stmt*/;
	    
//
//	We actually generate two IID normal distribution variables.
//	We cache the one & return the other.
// 
	double th = 6.28318530717958648 * u1
	double r = sqrt(-2 * log(u2));
		
	haveCachedNormal = 1;
	cachedNormal = r * sin(th);
	return(pStdDev*r*cos(th) + pMean);
    }
}