[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: problem with rand ?
From: |
David Bateman |
Subject: |
Re: problem with rand ? |
Date: |
Wed, 28 Jul 2004 12:12:33 +0200 |
User-agent: |
Mutt/1.4.1i |
Paul,
Using rand53 makes two calls to randi to get the bits. This is going
to slow things up. With the 32-bit versions we were just about as fast
as matlab, using the 53-bit versions we will be slower.
The speed differences were
randu rand53 rand53/randu Matlab R12
tic;x=rand(1e6,1);toc 0.072 0.110 1.53 0.087
tic;x=rand(5e6,1);toc 0.360 0.550 1.53 0.434
So there is a 53% slow up by doing this, and the octave code is now slower
than Matlab. Interesting matlab claims 53 bit precision, and so we should
be capable of getting similar speed to matlab with 53 bits.
I'd thought that maybe we could regain the speed reduction by precalculating
the factor "1.0/4294967296.0" in randu, and similarly in rand53
"1.0/9007199254740992.0". However, as these are constants it seems g++
is already smart enough to know it can precalculate them and no speed
difference results in making this change.
As for randn and rande, it isn't be that difficult to convert them to
use 53 bits. Its just a matter of using rand53 in the initialization
of the table, using TWO_TO_POWER_52 and TWO_TO_POWER_53 as appropriate.
The only complication is the sign of the randn values where we previously
used a signed long to implicit get the sign. With 53 bits this is no longer
possible.
Using this code has the following effect
32-bit 53-bit 53/32 Matlab R12
tic;x=randn(1e6,1);toc 0.095 0.152 1.60 0.092
tic;x=randn(5e6,1);toc 0.470 0.760 1.62 0.463
Which is a more significant slow-up. It should be noted that the state
vector of randn in matlab is 2-element, and it appears it uses a simplier
randu generator than the Mersenne Twister, so perhaps this speed difference
is normal. The document
http://www.mathworks.com/company/newsletters/news_notes/clevescorner/spring01_cleve.html
States that the period of randn in matlab is 2^64, which will be significantly
shorter than the randn in octave using the Mersenne Twister.
The test
s1=randn(1,128*1024); ss1=sort(s1); dss1=diff(ss1); length(find(dss1==0))
averaged about 3 or 4 identical values with 32 bits, but with 53-bits I
haven't seen an identical value with a few hundred runs of the above.
I attach a patch for this to this mail to change everything to 53 bits, but
haven't applied it yet. I'm not sure whether the speed hits is worth it in
most applications, though I'm willing to be convinced if someone stand-up
and stands they really need 53 bit accuracy.
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
patch.rand
Description: Text document