Path: utzoo!attcan!uunet!zephyr.ens.tek.com!uw-beaver!mit-eddie!bloom-beacon!oliver
From: oliver@athena.mit.edu (James D. Oliver III)
Newsgroups: comp.lang.fortran
Subject: Re: Random Number Functions--Generate versus Input
Keywords: fast method of generating random vectors
Message-ID: <1990Nov12.022637.23612@athena.mit.edu>
Date: 12 Nov 90 02:26:37 GMT
References: <1990Nov8.170334.24579@athena.mit.edu> <1990Nov8.215752.7075@athena.mit.edu> <32565@netnews.upenn.edu>
Sender: daemon@athena.mit.edu (Mr Background)
Organization: Massachusetts Institute of Technology
Lines: 66

In article <32565@netnews.upenn.edu> efrank@truth.hep.upenn.edu (Ed Frank) writes:

>  First, it is not clear to me that this is the fastest method.  My 
>method of generating random direction cosines is as follows.  Let DIR()
>be the array of direction cosines, and let RANN() be the (uniform) random
>number generator.  PI is 3.14.  Then
>        Phi = 2.0 * PI * RANN( dum)       ! in polar coord. Phi and cos(theta)
>        Dir(3) = -1.0 + 2.0*RANN( dum)    ! follow a flat distribution.
>        Sth    = SQRT( 1.0 = Dir(3)**2)   ! convert cos(theta) to sin(theta)
>        Dir(1) = Sth * COS( phi)          ! now get X and y
>        Dir(2) = Sth * SIN( phi)          ! direction cosines
>will produce the random directions for you.  You might want to include
>this in your benchmark.  I think that this will beat your method with
>respect to speed.

>   Second, I do not completely understand your method. If I am reading
>your description correctly, you will *NOT* obtain a uniform
>distribution with your technique.  Have you histogrammed the direction
>cosines produced by your method? The distributions should be flat.
>Keep in mind that if you distribute points on the surface of a sphere
>uniformly, and then project that distribution onto a plane, you will
>not obtain a flat distrubution bounded by x**2+y**2 = 1.  Most of the
>sphere's projected area will be near x**2 + y**2 = 1.

>In short: Phi is flat from 0 to 2*pi
>          Cos(theta) is flat from -1 to 1

Well, I hadn't originally intended to go into this, but your post and
several e-mail messages I've recieved indicate that the method of
generating the uniform vectors needs to be explained a little better than
my admittedly sketchy description.  Looking back, I realized that I gave
the impression that I directly substituted the values of my two random
numbers as x and y, which of course does not give the propert distribution.

The method I used was suggested by Marsaglia (Ann. Math. Stat., 43:645-646,
1972).  I stumbled onto it entirely by accident, and it seems to be the
most elegant method out there.  Briefly, it is as follows:

1) Generate two random uniform variates u and v on (-1,1).  If S = u**2 +
v**2 > 1, then generate new values for u and v.

2) Calculate Q = sqrt(1-S), then
	x=2.*u*Q
	y=2.*v*Q
	z=1.-2.*S

The efficiency of step 1) goes as the area of a unit square divided by the
area of a unit circle, or 4/pi, so the average number of random numbers
generated per success is 8/pi, or 2.55.  Thus the method requires 2.55
random function calls plus a square root to produce a vector, as opposed to
your method, which requires 2 random function calls, one square root, and
two calls to trigonometric functions.

Marsaglia provides the proof for this method, which basically results from
the fact that the distribution of S is a uniform variate on (0,1) and
independent of (u,v).

This method is significantly faster than any method I've used requiring
trig functions.  Although I haven't tried the one you mention above, I'm
failry confident that it is faster than that one as well.

--
____________________________
	Jim Oliver   
	oliver@athena.mit.edu /	joliver@hstbme.mit.edu
	oliver%mitwccf.BITNET@MITVMA.MIT.EDU