Code:

for (y = 0; y < 2048; y++) {
fla = (float)y / 651.580337f; // 2047/pi
fy = cos(fla) * hfbase; fr = -sin(fla) * hfbase;
for (x = 0; x < 4096; x++) {
hf[hu][x][y] = 0;
flo = (float)x / 651.7394919613114f + adjlong; // 4095/tau
fx = cos(flo) * fr; fz = sin(flo) * fr;
// fx, fy, and fz mapped to sphere of radius 1 with center at origin if hfbase is removed from above
float sum = 0.f;
for (j = 0; j < woct; j++) {
dx = fx * pp2[j]; dy = fy * pp2[j]; dz = fz * pp2[j];
if (dx < 0.f) dx += 131072.f; ix = (int)dx; dx -= ix; x1 = ix & 0x1f;
if (dy < 0.f) dy += 131072.f; iy = (int)dy; dy -= iy; y1 = iy & 0x1f;
if (dz < 0.f) dz += 131072.f; iz = (int)dz; dz -= iz; z1 = iz & 0x1f;
x0 = x1 - 1; x2 = x1 + 1; x3 = x1 + 2;
y0 = y1 - 1; y2 = y1 + 1; y3 = y1 + 2;
z0 = z1 - 1; z2 = z1 + 1; z3 = z1 + 2;
x0 &= 0x1f; x2 &= 0x1f; x3 &= 0x1f;
y0 &= 0x1f; y2 &= 0x1f; y3 &= 0x1f;
z0 &= 0x1f; z2 &= 0x1f; z3 &= 0x1f;
p0 = tricint(dx, perlin[x0][y0][z0], perlin[x1][y0][z0], perlin[x2][y0][z0], perlin[x3][y0][z0]);
p1 = tricint(dx, perlin[x0][y1][z0], perlin[x1][y1][z0], perlin[x2][y1][z0], perlin[x3][y1][z0]);
p2 = tricint(dx, perlin[x0][y2][z0], perlin[x1][y2][z0], perlin[x2][y2][z0], perlin[x3][y2][z0]);
p3 = tricint(dx, perlin[x0][y3][z0], perlin[x1][y3][z0], perlin[x2][y3][z0], perlin[x3][y3][z0]);
pa = tricint(dy, p0, p1, p2, p3);
p0 = tricint(dx, perlin[x0][y0][z1], perlin[x1][y0][z1], perlin[x2][y0][z1], perlin[x3][y0][z1]);
p1 = tricint(dx, perlin[x0][y1][z1], perlin[x1][y1][z1], perlin[x2][y1][z1], perlin[x3][y1][z1]);
p2 = tricint(dx, perlin[x0][y2][z1], perlin[x1][y2][z1], perlin[x2][y2][z1], perlin[x3][y2][z1]);
p3 = tricint(dx, perlin[x0][y3][z1], perlin[x1][y3][z1], perlin[x2][y3][z1], perlin[x3][y3][z1]);
pb = tricint(dy, p0, p1, p2, p3);
p0 = tricint(dx, perlin[x0][y0][z2], perlin[x1][y0][z2], perlin[x2][y0][z2], perlin[x3][y0][z2]);
p1 = tricint(dx, perlin[x0][y1][z2], perlin[x1][y1][z2], perlin[x2][y1][z2], perlin[x3][y1][z2]);
p2 = tricint(dx, perlin[x0][y2][z2], perlin[x1][y2][z2], perlin[x2][y2][z2], perlin[x3][y2][z2]);
p3 = tricint(dx, perlin[x0][y3][z2], perlin[x1][y3][z2], perlin[x2][y3][z2], perlin[x3][y3][z2]);
pc = tricint(dy, p0, p1, p2, p3);
p0 = tricint(dx, perlin[x0][y0][z3], perlin[x1][y0][z3], perlin[x2][y0][z3], perlin[x3][y0][z3]);
p1 = tricint(dx, perlin[x0][y1][z3], perlin[x1][y1][z3], perlin[x2][y1][z3], perlin[x3][y1][z3]);
p2 = tricint(dx, perlin[x0][y2][z3], perlin[x1][y2][z3], perlin[x2][y2][z3], perlin[x3][y2][z3]);
p3 = tricint(dx, perlin[x0][y3][z3], perlin[x1][y3][z3], perlin[x2][y3][z3], perlin[x3][y3][z3]);
pd = tricint(dy, p0, p1, p2, p3);
o = tricint(dz, pa, pb, pc, pd);
if (j < 2) {
sum += o * pn2[j + 1];
}
else {
if (o > 32767.5f) o = 65535.f - o;
sum += o * pn2[j];
}
}
sum *= sumdiv;
i = (int)sum;
if (i < 0) i = 0;
else if (i > 65535) i = 65535;
hf[hu][x][y] = i >> 8;
} }

wow, that gradient b/g makes code a trip :p