I know that if I want to generate directions uniformly within an n-dimensional hypersphere, I need to:
1) generate points $d_i, \, i\in{1......n}$ on a standard normal distribution $~N(0,1)$
2) and then get the unit vector with all of the $d_i$'s as the directions:
$D=(d_1...d_i...d_n)(\sum\limits_{i=1}^{n}d_i^2)^{-1/2}$
Somehow this results in a uniform distribution, even though I picked the $d_i$ values from a normal distribution. How?
https://mathoverflow.net/questions/26098/how-to-pick-a-random-direction-in-n-dimensional-space
You don't need to use a normal distribution to get a uniformly distributed direction. There are other methods.
For example, you could generate a point within the unit ball around the origin by using $n$ uniformly-distributed variables to get coordinates within an $n$-dimensional cube circumscribed around that ball and then rejecting any points that fall outside the ball. This gets progressively less efficient as $n$ increases, because the ratio of the volume of the ball to the volume of the cube tends toward zero.
Additionally, as noted in Is there something special about the normal distribution in the context of the random vector being equally likely to point in any direction?, it is possible to make a spherically-symmetric $n$-dimensional distribution that is not a multivariate Gaussian, and then you can use this to get a vector in a uniformly-distributed direction without using rejection.
But let's just look at why the multivariate Gaussian works, and not worry about other ways you could have generated a random direction.
The probability density of a single standard normal variable is $$ f_{N(0,1)}(x) = \frac{e^{-x^2/2}}{\sqrt{2\pi}} . $$
But we can take $n$ independently distributed variables with this same density and using them to determine the $n$ coordinates of a vector in an $n$-dimensional Euclidean vector space. We get a vector $\mathbf x$ with coordinates $(x_1, x_2, x_3, \ldots, x_n).$ Since the variables are all independent, we can get the density of their joint distribution at that point just by multiplying together their univariate densities at $x_1,$ $x_2,$, $x_3,$ etc. That is, the density of the joint distribution at $\mathbf x$ is \begin{align} f(x_1, x_2, x_3, \ldots, x_n) &= f_{N(0,1)}(x_1) f_{N(0,1)}(x_2) f_{N(0,1)}(x_3) \cdots f_{N(0,1)}(x_n) \\ &= \left(\frac{e^{-x_1^2/2}}{\sqrt{2\pi}} \right) \left(\frac{e^{-x_2^2/2}}{\sqrt{2\pi}} \right) \left(\frac{e^{-x_3^2/2}}{\sqrt{2\pi}} \right)\cdots \left(\frac{e^{-x_n^2/2}}{\sqrt{2\pi}} \right) \\ &= \frac{e^{-(x_1^2 + x_2^2 + x_3^2 + \cdots + x_n^2)/2}}{(2\pi)^{n/2}} \end{align}
But in $n$ dimensions, the length of a vector with coordinates $(x_1, x_2, x_3, \ldots, x_n)$ is $\sqrt{x_1^2 + x_2^2 + x_3^2 + \cdots + x_n^2}.$ So if we let $r = \sqrt{x_1^2 + x_2^2 + x_3^2 + \cdots + x_n^2}$ (that is, $r$ is the length of the random vector $\mathbf x$), then $$ f(x_1, x_2, x_3, \ldots, x_n) = \frac{e^{-r^2/2}}{(2\pi)^{n/2}}. $$
That is, the distribution $f$ is spherically symmetric--it is a function purely of the length of the vector, without regard for its direction. The density at distance $r$ from the origin is the same in every direction; every direction is equally likely.
OK, so now we have a random vector of some random length, equally likely to be in any direction. But you want a unit vector pointing in a random direction. So we take our random vector of random length, and find the unit vector in the same direction (which is still equally likely to be any direction). We do this by dividing by the length of the vector: $$ \mathbf d = \frac{\mathbf x}{r} = \frac{(x_1, x_2, x_3, \ldots, x_n)} {\sqrt{x_1^2 + x_2^2 + x_3^2 + \cdots + x_n^2}} = \left(\sum_{i=0}^n x_i^2\right)^{-1/2} (x_1, x_2, x_3, \ldots, x_n). $$
And there's the formula given in the mathoverflow answer.