🚀 go-pugleaf

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