octave-maintainers
[Top][All Lists]

## Re: randn benchmarks

 From: David Bateman Subject: Re: randn benchmarks Date: Sun, 25 Jan 2004 15:54:33 +0100 User-agent: Mutt/1.4.1i

```Paul,

I'd get even simplier than what you suggest. As the most basic mapping
of [0,1) to a normal distribution can give two random variables at once
like

{
double sq, u, v;
do {
u = 2*randu.randExc() - 1;
v = 2*randu.randExc() - 1;
sq = u*u + v*v;
} while (sq > 1.);
sq = sqrt(-2*log(sq)/sq);
*x++ = sq*u;
*x++ = sq*v;
}

This will be just as fast as the other code I suggested. Implemented for
octave-forge it would be something like

{
Matrix X(nr, nc);
unsigned int n = nr * nc;
unsigned int n2 = n - 2;
double *xv = X.fortran_vec ();
double sq, u, v;
for (unsigned int i = 0; i < n2; i+=2)
{
do {
u = 2.0*randu.randExc() - 1.0;
v = 2.0*randu.randExc() - 1.0;
sq = u*u + v*v;
} while (sq > 1. || sq == 0);
sq =  sqrt(-2*log(sq)/sq);
*xv++ = sq * u;
*xv++ = sq * v;
}
// If we have a odd number of terms, do last term
if (n & 1)
{
do {
u = 2.0*randu.randExc() - 1.0;
v = 2.0*randu.randExc() - 1.0;
sq = u*u + v*v;
} while (sq > 1. || sq == 0);
sq =  sqrt(-2*log(sq)/sq);
*xv++ = sq * u;
}
}

This way there is also no question that the the algorithm is correct.

However, and a big however, I've reimplemented the Ziggurat technique
from Marsaglia and Tsang's article without reference to their code,
using 256 level ziggurat and the Mersenne Twister as the base generator
of [0,1). I attach the code here for your perusal.

Thus I'd say there is no reason for me not to commit this version of
the Ziggurat rather than any other code... This is 2.5 times faster
than the above. I intend to commit the patch for this, but being aware
that rand.cc is your code at the moment, you of course have the right
to turn it down (ie. reverse my patch), though I hope you won't. I
plan to do this today, with my sort code.

Cheers
David

--
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: