But using Box-Miller gives you the complication of figuring out whether
you care about being thread safety, which you have to do if you're doing
anything serious. (See the comments in the random module for the gauss
method).
If you care about thread-safety, each thread should get its own Random instance
anyways.
By the way, does anyone know why there is no Poisson random numbers in
the module? The implementation is quite simple (but not as simple as the
Linux kernel *wink*):
def poisson(lambda_=1):
L = math.exp(-lambda_)
k = 0
p = 1
while 1:
k += 1
p *= random.random()
if p<= L:
break
return k-1
although this can be improved upon for large values of lambda_.
Where large is about 10. That's where my implementation in numpy.random switches
over to a more complicated, but more efficient implementation. It does require a
log-gamma special function implementation, though.
#define LS2PI 0.91893853320467267
#define TWELFTH 0.083333333333333333333333
long rk_poisson_ptrs(rk_state *state, double lam)
{
long k;
double U, V, slam, loglam, a, b, invalpha, vr, us;
slam = sqrt(lam);
loglam = log(lam);
b = 0.931 + 2.53*slam;
a = -0.059 + 0.02483*b;
invalpha = 1.1239 + 1.1328/(b-3.4);
vr = 0.9277 - 3.6224/(b-2);
while (1)
{
U = rk_double(state) - 0.5;
V = rk_double(state);
us = 0.5 - fabs(U);
k = (long)floor((2*a/us + b)*U + lam + 0.43);
if ((us >= 0.07) && (V <= vr))
{
return k;
}
if ((k < 0) ||
((us < 0.013) && (V > us)))
{
continue;
}
if ((log(V) + log(invalpha) - log(a/(us*us)+b)) <=
(-lam + k*loglam - loggam(k+1)))
{
return k;
}
}
}
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco