octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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

-- 
David Bateman                                address@hidden
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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]