static public double noise(double x, double y, double z) {

      int bx, by, bz, b0, b1, b00, b10, b01, b11;
      double rx0, rx1, ry0, ry1, rz, sx, sy, sz, a, b, c, d, u, v, q[];

      bx = (int) Math.IEEEremainder(Math.floor(x), B);
      if (bx < 0)
         bx += B;
      rx0 = x - Math.floor(x);
      rx1 = rx0 - 1;

      by = (int) Math.IEEEremainder(Math.floor(y), B);
      if (by < 0)
         by += B;
      ry0 = y - Math.floor(y);
      ry1 = ry0 - 1;

      bz = (int) Math.IEEEremainder(Math.floor(z), B);
      if (bz < 0)
         bz += B;
      rz = z - Math.floor(z);

      //if (bx < 0 || bx >= B + B + 2)
      //System.out.println(bx);

      b0 = p[bx];

      bx++;

      b1 = p[bx];

      b00 = p[b0 + by];
      b10 = p[b1 + by];

      by++;

      b01 = p[b0 + by];
      b11 = p[b1 + by];

      sx = s_curve(rx0);
      sy = s_curve(ry0);
      sz = s_curve(rz);

      q = G(b00 + bz);
      u = rx0 * q[0] + ry0 * q[1] + rz * q[2];
      q = G(b10 + bz);
      v = rx1 * q[0] + ry0 * q[1] + rz * q[2];
      a = lerp(sx, u, v);
      q = G(b01 + bz);
      u = rx0 * q[0] + ry1 * q[1] + rz * q[2];
      q = G(b11 + bz);
      v = rx1 * q[0] + ry1 * q[1] + rz * q[2];
      b = lerp(sx, u, v);
      c = lerp(sy, a, b);
      bz++;
      rz--;
      q = G(b00 + bz);
      u = rx0 * q[0] + ry0 * q[1] + rz * q[2];
      q = G(b10 + bz);
      v = rx1 * q[0] + ry0 * q[1] + rz * q[2];
      a = lerp(sx, u, v);
      q = G(b01 + bz);
      u = rx0 * q[0] + ry1 * q[1] + rz * q[2];
      q = G(b11 + bz);
      v = rx1 * q[0] + ry1 * q[1] + rz * q[2];
      b = lerp(sx, u, v);
      d = lerp(sy, a, b);

      return lerp(sz, c, d);
   }