I'm working on a game for which I want deterministic demo playback that is portable between architectures that treat floating point numbers differently. I'm using the Racket language, which conveniently has, as a primitive data type, non-floating-point representations of rational-number fractions. I want to use these to implement an approximately normally-distributed random function that accepts parameters for mean and standard deviation (skewness would be gold-plating).
Because of the limitations I've mentioned, any operation that takes in rational numbers and puts out irrational ones will need to be reimplemented from scratch in a way that produces approximations based on Racket's native fractions, not based on floating points. I've looked around at various algorithms for normal random functions, but of these, even many of the "simplest" ones like the Box-Muller transform involve things like square roots, logarithms, and trig functions. Iterated averaging is easy, so square roots aren't a problem, but I don't want to reinvent any more wheels than I need to here.
What are some algorithms I can use for generating approximately normal random numbers without invoking irrational operations like roots, logarithms, and trig functions?
I settled on a solution after typing up this question but before sending it, so I'll Share My Knowledge Q&A-Style.
After poring over several different SO posts on normally-distributed random numbers, I found that the best solution for my purposes was actually the most naive one: abuse the Central Limit Theorem. Random variables of any distribution, when added up, approximate a normal distribution just fine. In Racket, my solution turned out to be the delightfully concise
(define (random/normal μ σ)
(+ (* (- (for/sum ([i 12])
(random/uniform 0 1))
6)
σ)
μ))
where random/uniform is my function for generating uniformly random rational numbers.
In infix, imperative pseudocode, this means:
Function random_normal(μ, σ):
iterations := 12
sum := 0
for i from 1 to iterations:
sum += random_uniform(0, 1)
sum -= iterations / 2 # center the distribution on 0
return σ * sum + μ
A few SO answers mention this solution, but don't explain why 12 is the magic number here. When we add up those numbers, we want the standard deviation of that random sum to equal 1 so that we can stretch out or squish down the bell curve by the desired amount in a single multiplicative step.
If you sum a sample of random variables X1, X2, … XN, the standard deviation of the approximately normal distribution this creates is equal to
σXN/√N
where σx is the standard deviation of the variables themselves.* The standard deviation of a uniform random distribution from 0 to 1 is equal to 1/√12† so by substituting this in for σX we see that what we want is just
N/√(12N)=1,
which works out easily to N=12.
* See "Central Limit Theorem" on Wolfram MathWorld. Equation is given under identity (2), here multiplied by N to give the standard deviation of the sum rather than of the average.
† See "Continuous uniform distribution" on Wikipedia. Table on the right, "variance" (σ2), square-rooted.
It does, but the range of your distribution has to be truncated somewhere unless you have infinite memory, and ±6σ is A) almost as good as Box-Muller on a 32-bit machine and B) already huge.