Q:

Output range of Perlin noise

I'm investigating a few of the various implementations for coherent noise (I know there are libraries, but this is mostly for my own edification and curiosity) and how you can use it, and there's one problem I have with the original Perlin noise thing.
According to this frequently linked Math FAQ, the output range will be between -1 and 1, but I don't understand how the value gets to be in that range. 
As I understand it, the algorithm is basically this: each grid point has an associated random gradient vector of length 1. Then, for each point, for all four surrounding grid points, you calculate the dot product of the random gradient and the vector going from that grid-point. Then you use a fancy ease curve and linear interpolation to get that down to one value. 
But, here's my problem: these dot-products are, occasionally, going to be outside the range [-1, 1], and since you do the linear interpolation ultimately between the dot products, doesn't that mean that the final value will, on occasion, be outside the range of [-1, 1]?
Say, for instance, that one of the random vectors is (sqrt(2)/2, sqrt(2)/2) (which has a length of 1) and (0.8, 0.8) (which is in the unit square), you get a result of roughly 1.131. If that value is used in the linear interpolation, it's entirely possible that the value generated will be greater than 1. And, indeed, with my straight-forward implementation, that happens quite frequently. 
Am I missing something here?
For reference, here's my code in Java. Vec is a simple class to do simple 2d vector arithmetic, fade() is the ease curve, lerp() is linear interpolation, and gradient(x, y) gives you the gradient for that grid-point as a Vec. The gridSize variable gives you the size of the grid in pixels (it has type double):
public double getPoint(int x, int y) {
    Vec p = new Vec(x / gridSize, y / gridSize);
    Vec d = new Vec(Math.floor(p.x), Math.floor(p.y));

    int x0 = (int)d.x,
        y0 = (int)d.x;

    double d00 = gradient(x0    , y0    ).dot(p.sub(x0    , y0    )),
           d01 = gradient(x0    , y0 + 1).dot(p.sub(x0    , y0 + 1)),
           d10 = gradient(x0 + 1, y0    ).dot(p.sub(x0 + 1, y0    )),
           d11 = gradient(x0 + 1, y0 + 1).dot(p.sub(x0 + 1, y0 + 1));

    double fadeX = fade(p.x - d.x),
           fadeY = fade(p.y - d.y);

    double i1 = lerp(fadeX, d00, d10),
           i2 = lerp(fadeX, d01, d11);

    return lerp(fadeY, i1, i2);
}

Edit: here's the code for generating the random gradients: 
double theta = gen.nextDouble() * 2 * Math.PI; 
gradients[i] = new Vec(Math.cos(theta), Math.sin(theta));

Where gen is a java.util.Random. 

A:

You have y0 = (int)d.x;, but you mean d.y. This will most certainly affect your output range, and is the reason you are seeing such largely out-of-range values.

That said, the output range of Perlin noise is not actually [-1, 1]. While I'm not quite sure of the math myself (I must be getting old), this rather lengthy discussion works out that the actual range is [-sqrt(n)/2, sqrt(n)/2], where n is the dimensionality (2 in your case). So the output range of your 2D Perlin noise function should be [-0.707, 0.707]. This is somehow related to the fact that both d and the interpolation parameters are a function of p. If you read through that discussion, you may find the precise explanation you are looking for (particularly, post #7).
I am testing your implementation using the following program (which I hacked together from your example, so pardon the weird use of gridCells and gridSize):
import java.util.Random;

public class Perlin {

    static final int gridSize = 200;
    static final int gridCells = 20;
    static final Vec[][] gradients = new Vec[gridCells + 1][gridCells + 1];

    static void initializeGradient () {
        Random rand = new Random();
        for (int r = 0; r < gridCells + 1; ++ r) {
            for (int c = 0; c < gridCells + 1; ++ c) {
                double theta = rand.nextFloat() * Math.PI;
                gradients[c][r] = new Vec(Math.cos(theta), Math.sin(theta));                
            }
        }
    }

    static class Vec {
        double x;
        double y;
        Vec (double x, double y) { this.x = x; this.y = y; }
        double dot (Vec v) { return x * v.x + y * v.y; }
        Vec sub (double x, double y) { return new Vec(this.x - x, this.y - y); }
    }

    static double fade (double v) {
        // easing doesn't matter for range sample test.
        // v = 3 * v * v - 2 * v * v * v;
        return v;
    }

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

    static Vec gradient (int c, int r) {
        return gradients[c][r];
    }

    // your function, with y0 fixed. note my gridSize is not a double like yours.     
    public static double getPoint(int x, int y) {

        Vec p = new Vec(x / (double)gridSize, y / (double)gridSize);
        Vec d = new Vec(Math.floor(p.x), Math.floor(p.y));

        int x0 = (int)d.x,
            y0 = (int)d.y;

        double d00 = gradient(x0    , y0    ).dot(p.sub(x0    , y0    )),
               d01 = gradient(x0    , y0 + 1).dot(p.sub(x0    , y0 + 1)),
               d10 = gradient(x0 + 1, y0    ).dot(p.sub(x0 + 1, y0    )),
               d11 = gradient(x0 + 1, y0 + 1).dot(p.sub(x0 + 1, y0 + 1));

        double fadeX = fade(p.x - d.x),
               fadeY = fade(p.y - d.y);

        double i1 = lerp(fadeX, d00, d10),
               i2 = lerp(fadeX, d01, d11);

        return lerp(fadeY, i1, i2);

    }

    public static void main (String[] args) {

        // loop forever, regenerating gradients and resampling for range. 
        while (true) {

            initializeGradient();

            double minz = 0, maxz = 0;

            for (int x = 0; x < gridSize * gridCells; ++ x) {
                for (int y = 0; y < gridSize * gridCells; ++ y) {
                    double z = getPoint(x, y);
                    if (z < minz)
                        minz = z;
                    else if (z > maxz)
                        maxz = z;
                }
            }

            System.out.println(minz + " " + maxz);

        }

    }

}

I am seeing values within the theoretical range of [-0.707, 0.707], although I am generally seeing values between -0.6 and 0.6; which could just be a consequence of the value distribution and a low sampling rate.

