From octave-maintainers-request at bevo dot che dot wisc dot edu Sun Jan 25 08:58:41 2004
Subject: Re: randn benchmarks
From: David Bateman <David dot Bateman at motorola dot com>
To: "Dmitri A. Sergatskov" <dmitri at unm dot edu>
Cc: Paul Kienzle <pkienzle at users dot sourceforge dot net>,    octave-maintainers@bevo.che.wisc.edu
Date: Sun, 25 Jan 2004 15:54:33 +0100

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

-- 
David Bateman                                David dot Bateman at motorola dot com
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: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary