🚀 go-pugleaf

The code is straightforwardly presented in _The Science of Fractal
Images_ by Heinz-Otto Peitgen and (someone else, I forgot whom).  I
haveincluded it here, but that book is quite worthwhile for not only
fractals of all types in theory and practise, but also general random
number algorithms.

#define MAXLEVEL 200

int 		Arand, Nrand;
double 		GaussAdd, GaussFac;
double		delta[MAXLEVEL];

void InitGauss(seed)
int seed;
{
	Nrand = 4;
	Arand = INT_MAX;
	GaussAdd = sqrt(3.0*(double)Nrand);
	GaussFac = 2*GaussAdd/((double)Nrand*(double)Arand);
	srand((unsigned)seed);
}

double Gauss()
{
	double sum;
	int i;
	
	sum = 0.0;
	for(i = 1; i <= Nrand; i++)
		sum += (double)rand();
	return(GaussFac*sum - GaussAdd);
}

And viola, Gaussian random variables.

Hope this helps,

Gunther Anderson
gunther@jhunix.hcf.jhu.edu
gunther@jhuvms.hcf.jhu.edu
gunther@crabcake.cs.jhu.edu

eisen@jhunix.UUCP (Gunther Wil Anderson) presented an algorithm
that claims to produce normally-distributed random numbers.

If you read the code carefully, however, you'll observe that
it really doesn't do that.  It merely produces a crude approximation
to a normal distribution, based on the principle that the sum
of a large number of uniform distributions looks sort of like
a normal distribution.

In particular, the function will never return values with absolute
value greater than sqrt(12).

To do it right, pull out your dusty copy of Knuth's "Seminumerical
Algorithms" or just derive it yourself:  If D(x) is a probability
distribution, then let y = C(x) be a solution to the following equation:

	integral from -infinity to y   D(t) dt  = x

Then the function C(U) yields the desired distribution, where U
is a source of uniformly-distributed random numbers from 0 to 1.

(As a check, plug in D(x) = the characteristic function of [0,1], 
and observe that C(x) is the identity function on (0,1).)
--
raymond@math.berkeley.edu         mathematician by training, hacker by choice