[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
- Re: randn benchmarks, (continued)
- Re: randn benchmarks, Paul Kienzle, 2004/01/23
- Re: randn benchmarks, David Bateman, 2004/01/23
- Re: randn benchmarks, Paul Kienzle, 2004/01/23
- Re: randn benchmarks, David Bateman, 2004/01/24
- Re: randn benchmarks, Dirk Eddelbuettel, 2004/01/23
- Re: randn benchmarks, Dmitri A. Sergatskov, 2004/01/23
- Re: randn benchmarks, Dirk Eddelbuettel, 2004/01/24
- Re: randn benchmarks, Paul Kienzle, 2004/01/24
- Re: randn benchmarks, Paul Kienzle, 2004/01/24
- Re: randn benchmarks, Dmitri A. Sergatskov, 2004/01/24
- Re: randn benchmarks,
David Bateman <=
- eigenvalues 3 times speedup patch [Was: benchmarks], David Bateman, 2004/01/23
- Re: benchmarks, David Bateman, 2004/01/26