initial commit

1 files changed, 157 insertions, 0 deletions

diff --git a/reference/sphere.rand b/reference/sphere.rand
new file mode 100644
index 0000000..bf56705
--- /dev/null
+++ b/reference/sphere.rand
@@ -0,0 +1,157 @@
+From: ags@seaman.cc.purdue.edu (Dave Seaman)
+Newsgroups: sci.math.num-analysis
+Subject: Re: N-dim spherical random number drawing
+Date: 20 Sep 1996 13:04:58 -0500
+
+In article <wmEgVDG00VUtMEhUwg@andrew.cmu.edu>,
+Shing-Te Li  <sl2x+@andrew.cmu.edu> wrote:
+>Does anyone know a good algorithm that can draw a vector of random
+>numbers uniformly from a N dimensional sphere?
+
+Here are four methods for generating uniformly-distributed points on a unit
+sphere:
+
+  (1) The normal-deviate method.  Choose x, y, and z, each
+      normally-distributed with mean 0 and variance 1.  (Actually, the
+      variance does not matter, provided it remains fixed.)  Normalize
+      the result to obtain a uniform density on the unit sphere.
+
+      This method also generalizes nicely to n dimensions.  It works
+      because the vector chosen (before normalization) has a density
+      that depends only on the distance from the origin.  The method is
+      mentioned in Knuth, _Seminumerical Algorithms_, 2nd. ed., pp.
+      130-131.
+
+  (2) The hypercube rejection method.  Choose x, y, and z, uniformly
+      distributed on [-1,1].  Compute s = x^2 + y^2 + z^2.  If s > 1,
+      reject the triplet and start over.  Once you have an acceptable
+      vector, normalize it to obtain a point on the unit sphere.
+
+      This method also generalizes to n dimensions, but as the dimension
+      increases the probability of rejection becomes high and many random
+      numbers are wasted.
+
+  (3) The trig method.  This method works only in 3-space, but it is
+      very fast.  It depends on the slightly counterintuitive fact (see
+      proof below) that each of the three coordinates of a uniformly
+      distributed point on S^2 is uniformly distributed on [-1,1] (but
+      the three are not independent, obviously).  Therefore, it
+      suffices to choose one axis (Z, say) and generate a uniformly
+      distributed value on that axis.  This constrains the chosen point
+      to lie on a circle parallel to the X-Y plane, and the obvious
+      trig method may be used to obtain the remaining coordinates.
+
+	(a) Choose z uniformly distributed in [-1,1].
+	(b) Choose t uniformly distributed on [0, 2*pi).
+	(c) Let r = sqrt(1-z^2).
+	(d) Let x = r * cos(t).
+	(e) Let y = r * sin(t).
+
+  (4) A variation of method (3) that doesn't use trig may be found in
+      Knuth, _loc. cit._.  The method is equivalent to the following:
+
+	(a) Choose u and v uniformly distributed on [-1,1].
+	(b) Let s = u^2 + v^2.
+	(c) If s > 1, return to step (a).
+	(d) Let a = 2 * sqrt(1-s).
+	(e) The desired point is (a*u, a*v, 2*s-1).
+      
+      This method uses two-dimensional rejection, which gives a higher
+      probability of acceptance compared to algorithm (2) in three
+      dimensions.  I group this with the trig method because both
+      depend on the fact that the projection on any axis is uniformly
+      distributed on [-1,1], and therefore both methods are limited to
+      use on S^2.
+
+I have found the trig method to be fastest in tests on an IBM RS/6000
+model 590, using the XLF 3.2 Fortran compiler, with the IMSL
+subroutines DRNUN and DRNNOA generating the random numbers.  The trig
+routines were accelerated by using the MASS library, and sqrt
+computations were performed by the hardware sqrt instruction available
+on the 590.  Timings for 1 million random points on S^2:
+
+	(1) normal-deviate method                   9.60 sec.  
+	(2) 3-D rejection                           8.84 sec.  
+	(3) trig method                             2.71 sec.  
+	(4) polar with 2-D rejection                5.08 sec.
+
+The fact that algorithm (3) does not involve rejection turns out to be
+a huge advantage because it becomes possible to generate all the points
+in a single loop with no conditional statements, which optimizes very
+effectively and more than offsets the cost of using trig functions.
+The alternative with algorithms (2) and (4) is to generate the points
+one at a time and test each one for acceptance before proceeding to the
+next, which requires nested loops and many more subroutine calls to
+generate the random numbers.  Algorithm (1) does not explicitly use
+rejection, but the IMSL subroutine DRNNOA uses rejection in generating
+the normally-distributed numbers that are required.
+
+In higher dimensions, the normal-deviate method is preferred
+because the probability of rejection in the hypercube method
+grows very rapidly with the dimension.
+
+--------------------------------------------------------------------
+
+Here is a proof of correctness for the fact that the z-coordinate is
+uniformly distributed.  The proof also works for the x- and y-
+coordinates, but note that this works only in 3-space.
+
+The area of a surface of revolution in 3-space is given by
+
+	A = 2 * pi * int_a^b f(x) * sqrt(1 + [f'(x}]^2) dx
+
+Consider a zone of a sphere of radius R.  Since we are integrating in
+the z direction, we have
+
+	f(z) = sqrt(R^2 - z^2)
+	f'(z) = -z / sqrt(R^2-z^2)
+
+	1 + [f'(z)]^2 = r^2 / (R^2-z^2)
+
+	A = 2 * pi * int_a^b sqrt(R^2-z^2) * R/sqrt(R^2-z^2) dz
+
+	  = 2 * pi * R int_a^b dz
+
+	  = 2 * pi * R * (b-a)
+
+	  = 2 * pi * R * h, 
+	  
+where h = b-a is the vertical height of the zone.  Notice how the integrand
+reduces to a constant.  The density is therefore uniform.
+
+-- 
+Dave Seaman			dseaman@purdue.edu
+      ++++ stop the execution of Mumia Abu-Jamal ++++
+    ++++ if you agree copy these lines to your sig ++++
+++++ see http://www.xs4all.nl/~tank/spg-l/sigaction.htm ++++
+==============================================================================
+
+From: George Marsaglia <geo@stat.fsu.edu>
+Subject: Random points on a sphere.
+Date: Tue, 15 Jun 1999 18:03:07 -0400
+Newsgroups: sci.math.num-analysis,sci.math,sci.stat.math
+
+Questions on methods for random sampling from the surface
+of a sphere keep coming up, and seemingly get lost as
+new generations of computer users encounter the problem.
+
+I described the methods that I had used for years in
+"Sampling from the surface of a sphere", Ann Math Stat, v43, 1972,
+recommending the two that seemed the fastest:
+
+1) Choose standard normal variates x1,x2,x3,
+put R=1/sqrt(x1^2+x2^2+x3^2) and return the point
+(x1*R,x2*R,x3*R).
+
+2) Generate u and v, uniform in [-1,1]  until S=u^2+v^2 < 1,
+then return (1-2*S,u*r,v*r), with r=2*sqrt(1-S).
+
+The second method is fast and easiest to program---unless
+one already has a fast normal generator in his library.
+
+Method 1) seems by far the fastest (and easily applies to higher
+dimensions) if my latest normal generator is used.
+It produces normal variates at the rate of seven million
+per second in a 300MHz Pc.
+
+George Marsaglia