Path: utzoo!utgpu!news-server.csri.toronto.edu!mailrus!tut.cis.ohio-state.edu!purdue!haven!decuac!shlump.nac.dec.com!ryn.esg.dec.com!allvax!jroth
From: jroth@allvax.dec.com (Jim Roth)
Newsgroups: comp.graphics
Subject: Re: random points on surface of sphere
Summary: a random walk is very fast
Keywords: uniform spherical distribution, random walk
Message-ID: <1523@ryn.esg.dec.com>
Date: 3 May 90 06:52:30 GMT
Sender: guest@ryn.esg.dec.com
Organization: Digital Equipment Corporation
Lines: 50


In article <49791@lanl.gov>, matthias@mpx1.lampf.lanl.gov (Matthias, Bjorn E.) writes...
>This topic has certainly been beaten to death recently.  Still, I believe
>that there are some more useful things to say.
.. stuff deleted ...
>in Monte Carlo simulations to generate quantities with certain distributions.
>Often it is not possible to invert the distribution in question and then one
>resorts to rejection methods (e.g. von Neumann rejection) or to random walk
>traversal algorithms (e.g. Metropolis algorithm).  But when it is possible,
>it is usually faster to do things directly as outlined above for this case.

In this case, it is considerably faster to do a random walk on the sphere,
because this avoids relatively expensive transcendental function calls; in
fact the only library call besides the random number generator is a square root
which could be done faster as an inverse square root if such a routine were
available (or better if a normalize to unit vector routine were available.)

Try the following code, which cyclically permutes the coordinates and
performs a random plane rotation; the distribution of theta is not too
important, so rejection is not needed here.  The limit distribution is the
same as long as it is not pathologically bad.

	x = sqrt(1.0/3.0)	! initial unit vector
	y = sqrt(1.0/3.0)
	z = sqrt(1.0/3.0)
	do i=1,succ
		cs = 1.0-2.0*ran(ix)
		sn = 1.0-2.0*ran(ix)
		t = 1.0/sqrt(cs**2+sn**2)
		cs = cs*t
		sn = sn*t
		t = x
		x = cs*y-sn*z
		y = sn*y+cs*z
		z = t				
		ax(i)=x
		ay(i)=y
		az(i)=z
	enddo

Needless to say, the rejection method could be tweaked by randomly selecting
a sub-cube not completely outside the unit sphere and uniformly choosing a
point within it, or some other table driven technique to improve the succes
rate.  It would approach the speed of the random walk then.

However, for high dimensions the random walk is very attractive, since a
plane rotation in one of n choose 2 planes can be done just as cheaply as in
2-D (as long as successive vectors need not be completely independant.)

- Jim