Perlin Noise function - values way out of range

Question

So, I've been trying to use Perlin Noise in my project, and bizarrely the noise function outputs mostly values way out of the [-1, 1] bounds. I've implemented it in C through Perlin's code for his improved noise in Java, and another useful implementation in C++ by sol. I've ran through it at least 20 times and made sure it was exactly the same, but couldn't find any causes for the issue. Here is how it looks:

//   DEFAULT_PERMUTATION is defined through the preprocessor.
// It is basically an initializer of the same array used in Perlin's
// implementation, except it has been duplicated into a 512 element array


double simpleNoise(double x, double y, double z){

double res;

int p[512] = DEFAULT_PERMUTATION;

int X = (int) floor(x) & 255;
int Y = (int) floor(y) & 255;
int Z = (int) floor(z) & 255;

double u = fade(x);
double v = fade(y);
double w = fade(z);

int A = p[X] + Y;
int AA = p[A] + Z;
int AB = p[A + 1] + Z;
int B = p[X + 1] + Y;
int BA = p[B] + Z;
int BB = p[B + 1] + Z;

x -= floor(x);
y -= floor(y);
z -= floor(z);


res = lerp(w, lerp(v, lerp(u, grad(p[AA  ], x  , y  , z   ),  
                                 grad(p[BA  ], x-1, y  , z   )), 
                         lerp(u, grad(p[AB  ], x  , y-1, z   ),  
                                 grad(p[BB  ], x-1, y-1, z   ))),
                 lerp(v, lerp(u, grad(p[AA+1], x  , y  , z-1 ),  
                                 grad(p[BA+1], x-1, y  , z-1 )), 
                         lerp(u, grad(p[AB+1], x  , y-1, z-1 ),
                                 grad(p[BB+1], x-1, y-1, z-1 ))));
return res;
}

The other private functions are pretty straightfoward too:

double fade(double t){
return t * t * t * (t * (t * 6 - 15) + 10);
}

double lerp(double t, double a, double b){
return a + t * (b - a); 
}

double grad(int hash, double x, double y, double z){
int h = hash & 15;                     
  double u = h<8 ? x : y,                 
         v = h<4 ? y : h==12||h==14 ? x : z;
  return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
}

In my project, the noise looks like this, due to the fact that the color only cares about values between 0 and 1, but tests with some positions give absurd values like -15, 32, some bigger than thousands. Feeding the values (1,-4,1), (2,5,1), (2,0,0) are interesting to say the least.

What could the problem be?

EDIT: For investigation purposes, this program lists inconsistencies found:

#define DEFAULT_PERMUTATION { \
    151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142, \
    8,99,37,240,21,10,23,190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117, \
    35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71, \
    134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41, \
    55,46,245,40,244,102,143,54, 65,25,63,161,1,216,80,73,209,76,132,187,208, 89, \
    18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226, \
    250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182, \
    189,28,42,223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167,  \
    43,172,9,129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246, \
    97,228,251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239, \
    107,49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254, \
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,\
    151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142, \
    8,99,37,240,21,10,23,190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117, \
    35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71, \
    134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41, \
    55,46,245,40,244,102,143,54, 65,25,63,161,1,216,80,73,209,76,132,187,208, 89, \
    18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226, \
    250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182, \
    189,28,42,223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167,  \
    43,172,9,129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246, \
    97,228,251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239, \
    107,49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254, \
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180 }
#define SIZE 10
double fade(double t){
return t * t * t * (t * (t * 6 - 15) + 10);
}

double lerp(double t, double a, double b){
return a + t * (b - a); 
}

double grad(int hash, double x, double y, double z){
int h = hash & 15;                     
double u = h<8 ? x : y,                 
     v = h<4 ? y : h==12||h==14 ? x : z;
return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
}
double simpleNoise(double x, double y, double z){
double res;
int p[512] = DEFAULT_PERMUTATION;
int X = (int) floor(x) & 255;
int Y = (int) floor(y) & 255;
int Z = (int) floor(z) & 255;
double u = fade(x);
double v = fade(y);
double w = fade(z);
int A = p[X] + Y;
int AA = p[A] + Z;
int AB = p[A + 1] + Z;
int B = p[X + 1] + Y;
int BA = p[B] + Z;
int BB = p[B + 1] + Z;
x -= floor(x);
y -= floor(y);
z -= floor(z);
res = lerp(w, lerp(v, lerp(u, grad(p[AA  ], x  , y  , z   ),  
                             grad(p[BA  ], x-1, y  , z   )), 
                     lerp(u, grad(p[AB  ], x  , y-1, z   ),  
                             grad(p[BB  ], x-1, y-1, z   ))),
             lerp(v, lerp(u, grad(p[AA+1], x  , y  , z-1 ),  
                             grad(p[BA+1], x-1, y  , z-1 )), 
                     lerp(u, grad(p[AB+1], x  , y-1, z-1 ),
                             grad(p[BB+1], x-1, y-1, z-1 ))));
return res;
}
main(){
int i, j, k, n=0;
double e;
for(i=-SIZE/2; i<SIZE/2; i++){
    for(j=-SIZE/2; j<SIZE/2; j++){
        for(k=-SIZE/2; k<SIZE/2; k++){
            e = simpleNoise(i, j, k);
            if(fabs(e) > 1){
               printf("P(%d,%d,%d): %f\n", i, j, k, e);
               n++;
            }
        }
    }
}
printf("%d inconsistencies found", n);
}

Answer 1

The reference implementation here does

x -= floor(x);
y -= floor(y);
z -= floor(z);

before such steps as

double u = fade(x);
double v = fade(y);
double w = fade(z);

Your implementation moved that step much further down the code. Why?

If you move it back to its original location, in your "investigation" code the noise function will always expectedly return zero, since Perlin noise function is supposed to return zero at integer points.