Computing the pseudo-random gradient:

Precompute table of permutations P[n]
Precompute table of gradients G[n]
G = G[ i + P[ j + P[ k ] ] ]

Actually, you don't really need to call a mod function; since 256 is a power of two, I just did the addition and masked off all but the lowest 8 bits, so the actual computational cost is really more like what you now see on the left.

To precompute the table of gradients G, I did a Monte Carlo simulation, as follows:

The problem was to compute 256 gradient vectors uniformly around the surface of a sphere. I knew that I could compute points uniformly within a cubic volume, by picking uniformly random values of x,y, and z via a rand() function.

In order to make a uniformly random distribution of points on the surface of a sphere, you really just need to make a uniformly random distribution of points within a spherical volume. Then you can normalize them (ie: resize them all to unit length).

So my Monte Carlo simulation worked as follows: I generated uniformally random points within the cube which surrounds the unit sphere. I threw out all the points that lie outside the unit sphere and kept all the points that lie inside the unit sphere. These are the points for which x²+y²+z² < 1.

I stopped this process only when I had filled up a table of 256 entries. I then normalized all the entries in the table, and I was done.