circular random number distribution

[Robert G. Brown]

Jochen Küpper writes:

> Heiko Bauke <Heiko.Bauke@Physik.Uni-Magdeburg.DE> writes:
>
>> | do {
>> |     x = gsl_ran_flat(rng, -1, 1);
>> |     y = gsl_ran_flat(rng, -1, 1);
>> | } while(sqrt(x*x + y*y) > 1.);
>
> Oops, of course! (Stupid me.)
>
> Ok, I wrote a "real" benchmark, see attached code, generating 1e8
> random pairs in a 1/1-circle using gsl_rng_taus and the different
> methods discussed (see attached code for details).
>
> Intel(R) Xeon(TM) CPU 2.80GHz machine with gcc-3.3.3:
> ,----
> | reject: 15.59 s
> | polar:  27.25 s
> | dir_2d: 24.38 s
> | trig2d: 27.68 s
> `----
> AMD Opteron(tm) Processor 850 2.4GHz with gcc-3.3.3:
> ,----
> | reject: 4.48 s
> | polar:  16.25 s
> | dir_2d: 8.29 s
> | trig2d: 16.87 s

These numbers make nearly perfect sense.  taus in my benchmarking
program is one of the fastest at 20-35 nsec/rand.  This means that the
resulting times are heavily dominated by the relative efficiency of
trancendentals.  The opteron times are truly impressive, as well,
suggesting that taus is able to take advantage of the wide data path and
pipeline of the opteron CPU.

Note also that the efficiency of accept-reject over alternative
hyperspherical methods goes like the fraction of the volume outside of
the hypersphere:

   2d:  ((2r)^2 - \pi r^2)/(2r)^2 = (4 - \pi)/4 \approx 0.22

   3d:  ((2r)^3 - 4\pi/3 r^3)/(2r)^3 = (8 - 4\pi/3)/8 \approx 0.47

   n-sphere ((2r)^n - S_n r^2/n)/(2r)^n
              = (2^n - 2^{(n+1)/2} \pi^{(n-1)/2}/((n-2)!! n))/2^n (n odd)
              = (2^n - 2 \pi^{n/2}/((n/2-1)! n))/2^2 (n odd)/2^n (n even)

which unfortunately goes to 1 as n-> infinity.  If you like, the volume
of any n-volume for large n is concentrated on its surface.  This means
that there will generally be a crossover where accept-reject ON AN
N-CUBE becomes inefficient and loses out to direct methods.  I'd expect
them to pretty much break even in 3d and go downhill from there.

However, accept reject methods can often be "rescued" by means of Clever
Tricks -- usually by means of generating samples to be tested from a
volume that is a better approximation of the desired volume than the
N-cube is of the N-sphere.

When I run taus through dieharder (which doesn't yet implement all of
the diehard tests but does implement a few others) it seems to do OK
(makes it through 5-bit random, which weak ones generally do not).  For
most purposes it should be a reasonable choice and it is fast.

    rgb
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20050809/5ab0b06f/attachment.sig>