From octave-maintainers-request at bevo dot che dot wisc dot edu Fri Jan 23 06:39:55 2004 Subject: Re: randn benchmarks From: David Bateman To: Paul Kienzle Cc: octave-maintainers at bevo dot che dot wisc dot edu Date: Fri, 23 Jan 2004 13:36:03 +0100 According to Paul Kienzle (on 01/23/04): > Yeah, I saw it there, but I didn't see any indication of a > license in the ziggurat code when I looked. This is a > user contribution rather than a base package, so it > hasn't necessarily gone through the same scrutiny as > others. That said, I didn't ask the package author if he > obtained permission before using it. > > The license on the jstatsoft page says that the paper > is freely distributed. It does not say that it is freely > redistributable. > > It is probably okay for us to use it, but I prefer to get > permission first. There is also an implementation as a C++ class http://www.codeproject.com/cpp/zigurat.asp I don't think Marsaglia is going to reply, give that he was writing books in the 60's there is a strong chance he has retired... You'd probably be better of contacting teh co-author http://www.csis.hku.hk/~tsang/ and getting permission from him. In any case, the ziggurat code in this paper is for 2^32 bit seeds, which is a bit short for my liking. I'd prefer to target a 64 bit seed as does matlab, so a possible solution is just to rewrite the code from the paper for a 64 bit seed.. In any case, before that the randn in octave-forge is really terrible. I sent another mail last night, but forgot reply-all (so only you got it), which is partially quoted here > > Ok, randn in octave-forge is slow but rand is pretty good. I tried a simple > code to generate randn based on some fortran code on netlib to compare it > with what is currently in octave-forge... The code snippet is > > { > double u,v; > while(1) > { > u = randu.randExc(); > v = randu.randExc(); > v = 1.7156 * (v - 0.5); > double x = u - 0.449871; > double y = abs(v) + 0.386595; > double q = x*x + y*(0.196*y - 0.25472*x); > if (q < 0.27597) > break; > if ((q <= 0.27846) && (v*v < - 4. * log(u) *u*u)) > break; > } > X(r,c) = v / u; > } > > I've done a simple test with "a=randn(1,1e6);hist(a,[-5:0.2:5])" to see > if its normally distributed, and to my eye it appears to be. The interesting > thing is this > > before > ----- > octave:2> tic; a = randn(1500, 1500); toc > > after > ----- > octave:2> tic; a = randn(1500, 1500); toc > ans = 0.45562 > > matlab R12 > ---------- > >> tic; a = rand(1500, 1500); toc > > elapsed_time = > > 0.1762 > > So, this simplistic algorithm is not as fast as matlab, but its still a good > 4 times faster than what is currently in octave-forge. > > In fact even, the most basic way of generating a normal distribution > > { > double sq, u, v; > do { > u = randu.randExc(); > v = randu.randExc(); > sq = u*u + v*v; > } while (sq > 1.); > X(r,c) = sqrt(-2*log(sq)/sq)*u; > } > > is faster than what is in octave-forge. Is there any reason to > use the code for randn that is already in octave-forge, over one > of these? > > octave:2> tic; a = randn(1500, 1500); toc > ans = 0.86871 > Couldn't we at least replace the randn code in octave-forge with one of these to generate the distribution from the Mersenne Twister code? This would be a stop-gap measure until we can deal with the Ziggurat 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