As part of a more complex algorithm, I need a fast method to find random points of the nsphere, $S^n$, starting with a RNG (random number generator). A simple way to do this (in low dimensions at least) is to select a random point of the (n+1)ball and normalize it. And to get a random point of the (n+1)ball select a random point of the (n+1)cube $[1,1]^{n+1}$ (by selecting (n+1) points of $[0,1)$ with the RNG and scaling using $x \mapsto 2x1$) and then use "rejection", i.e., just ignore a point if it is not in the (n+1)ball. This works fine if n is reasonably small, however for large n the volume of the ball is such a tiny fraction of the volume of the cube that rejection is enormously inefficient. So what is a good alternative approach.

17$\begingroup$ I don't know how to do this, but one abstract way to choose a random point on the sphere is to choose the coordinates according to a Gaussian distribution. The resulting random point will be chosen rotationallysymmetrically, so one can then divide by the length to get a point on the sphere chosen randomly and rotationally symmetrically. So if one had a way to approximate a Gaussian distribution, then this should be possible (and I suppose the law of large numbers says that one should be able to do this starting with any probability distribution...). $\endgroup$ – Ian Agol Jul 10 '13 at 19:11

$\begingroup$ (I meant the central limit theorem, not the law of large numbers, but it looks like you've received some answers addressing approximating a normal distribution) $\endgroup$ – Ian Agol Jul 10 '13 at 19:26

2$\begingroup$ Just an aside: This is treated, e.g., on the Wikipedia page for $n$sphere and (very closely related) versions of it have been asked previously on stats.SE and math.SE. Similar algorithms exist for generating uniform points on the simplex. $\endgroup$ – cardinal Jul 11 '13 at 1:01
The usual approach is to generate $n+1$ i.i.d. mean zero Gaussian random variables $X_1, \dotsc, X_{n+1}$ to get a random point $X$ in $(n+1)$space with rotationally invariant distribution and normalize.
Incidentally, if you ever actually need to generate a random point in an $n$ball, the best way is probably to generate a random point in the $(n+1)$sphere as above and drop the last two coordinates.

14$\begingroup$ To complete Mark's answer, to get two Gaussian variates from two uniform variates, use the BoxMuller transform: en.wikipedia.org/wiki/Box%E2%80%93Muller_transformation $\endgroup$ – Yoav Kallus Jul 10 '13 at 19:21

1$\begingroup$ @Mark Meckes Many thanks, Mark. In retrospect, if I had been more adept with my Googling I might have discovered that the Gaussian RV approach was a "wellknown" answer to my question, and would not have needed to ask it on MO. But as usual, I have come away very impressed with the utility of this great website and learned a lot from the answers I received. $\endgroup$ – Dick Palais Jul 11 '13 at 16:02
Here some algorithm that produce normal distributed random numbers :
1) Polar algorithm
2) Ziggurat algorithm
3) BoxMuller algorithm
By dividing the random vector of n+1 random numbers by the norm of the vector you get random numbers on the nsphere.
I can't resist linking to this gem which describes a simple method for a generalisation of your problem.
Nice question! I ran into a similar problem a few years ago  even for dimension $10$, the rejection method was annoyingly slow. One of the problems is that such questions straddle at least three huge fields (discrete mathematics, statistics and CS), so often one doesn't know where to start looking.
Here's what we finally used: a nice, explicit solution for how to generate uniformly random points on the $n$sphere can be found in Section III of the unpublished paper:
Cumbus, Damien, Walker, Uniform sampling in the hypersphere via latent variables and the Gibbs sampler (1996).
You can find a copy of their work here. If you want to sample from some measurable subset of the $n$sphere instead of the whole thing, try the (much less explicit) technical report
Shao, Badler, Spherical Sampling by Archimedes' Theorem (also 1996)
which can be found here. Good luck with your algorithm.
If you convert to an ndimensional spherical coordinate system, then you would need n1 angles and a radius (in [0,1]) to uniquely determine a point within (or on) the nsphere. (If you want points strictly on the nsphere, force the radius to be 1).
To get the angles, utilize a uniform distribution over [0,2*pi). The radius is another uniform variate over [0,1]. Then you get a uniform distribution over the nsphere by conversion (from the wikipedia page for the nsphere):
x_1 = r*cos(phi_1)
x_2 = r*sin(phi_1)*cos(phi_2)
. . .
x_{n1} = r*sin(phi_1)*sin(phi_2)*...*sin(phi_{n2})*cos(phi_{n1})
x_n = r*sin(phi_1)*sin(phi_2)*...*sin(phi_{n2})*sin(phi_{n1})
The cost of this is n1 uniform variates (and another if you want r in [0,1]) and the subsequent cosine and sine operations (most of which can be reused...)

3$\begingroup$ This will not generate uniform distribution. $\endgroup$ – Emil Jeřábek Jul 11 '13 at 16:52