Re: problem with rand ?

[David Bateman]

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