How does this hash function work in this perlin noise implementation?

In Perlin's implementation of his noise function, he uses an array to hash the input coordinates and later find the gradient vectors. But in Stefan Gustavson's implementation, a strange "permute" function is used:

vec4 permute(vec4 x)
{
  return mod289(((x*34.0)+1.0)*x);
}

vec4 mod289(vec4 x)
{
  return x - floor(x * (1.0 / 289.0)) * 289.0;

}

Then, it is used like this:

vec4 i = permute(permute(ix) + iy);

vec4 gx = fract(i * (1.0 / 41.0)) * 2.0 - 1.0 ;
vec4 gy = abs(gx) - 0.5 ;
vec4 tx = floor(gx + 0.5);
gx = gx - tx;

I understand that the "*2.0 - 1.0" is for getting gx to the range -1 to +1, but I do not understand the rest of the hashing. How does this "permute" function work and why does it divided by 41.0?

Thanks.