local noise = {}
local exp,sin,cos,floor,log,acos,sqrt = math.exp,math.sin,math.cos,math.floor,math.log,math.acos,math.sqrt
local GR, PI, TAU, SQRT5, LOG_GR = (sqrt(5)+1)/2, math.pi, 2*math.pi, sqrt(5), log((sqrt(5)+1)/2)

local function cdf(x,sigma)
    return .5 + .5*(x<0 and -1 or 1)*sqrt(1.-exp(-4./TAU * x*x))
end

local function _defaultArgs(resolution,random,decayFn)
    if not resolution then resolution = 4 end
    if not random then random = math.random end
    if not decayFn then decayFn = function(x) return .1^x end end
    return resolution,random,decayFn
end

local function _amplitudesAndOffsets(decayFn,numAmplitudes, numOffsets, random)
    local sigma = 0
    local amplitudes = {}
    for i=1,numAmplitudes do
        local a = decayFn((i-1+random())/numAmplitudes)
        amplitudes[i] = a
        sigma = sigma + a^2
    end
    sigma = math.sqrt(sigma/2)
    local offsets = {}
    for i=1,numOffsets do
        offsets[i] = random()*TAU
    end
    return amplitudes,sigma,offsets
end

noise.make1d = function(resolution,random,decayFn)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,resolution,random)
    return function(x)
        local noise = 0
        for i,a in ipairs(amplitudes) do
            noise = noise + a*cos(x/a + offsets[i])
        end
        return cdf(noise/sigma);
    end
end

noise.make2d = function(resolution,random,decayFn)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,2*resolution,random)
    sigma = sigma/sqrt(2)
    return function(x,y)
        local noise = 0
        for i,a in ipairs(amplitudes) do
            local angle = ((i*GR) % 1)*TAU
            local u = x*cos(angle) - y*sin(angle)
            local v = x*cos(angle+TAU/4) - y*sin(angle+TAU/4)
            noise = noise + a/2*(cos(u/a + offsets[2*i]) + cos(v/a + offsets[2*i-1]))
        end
        return cdf(noise/sigma);
    end
end

noise.make3d = function(resolution,random,decayFn)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,3*resolution,random)
    sigma = sigma/sqrt(3)
    return function(x,y,z)
        -- Find the biggest fibonacci number F_n such that F_n < resolution
        local n = floor(log((resolution-1)*SQRT5 + .5)/LOG_GR)
        local dec = floor(.5 + (GR^n)/SQRT5) -- F_n, using closed form Fibonacci
        local inc = floor(.5 + dec/GR) -- F_(n-1)

        local noise,i,j = 0,0,0
        for i=0,resolution-1 do
            if j >= dec then
                j = j - dec
            else
                j = j + inc
                if j >= resolution then
                    j = j - dec
                end
            end
            -- Convert golden ratio sequence into polar coordinate unit vector
            local phi = ((i*GR) % 1)*TAU
            local theta = acos(-1+2*((j*GR) % 1))
            -- Make an orthonormal basis, where n1 is from polar phi/theta,
            -- n2 is roated 90 degrees along phi, and n3 is the cross product of the two
            local n1 = {sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)}
            local n2 = {sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.)}
            -- Cross product
            local n3 = {n1[2]*n2[3] - n1[3]*n2[2],
                        n1[3]*n2[1] - n1[1]*n2[3],
                        n1[1]*n2[2] - n1[2]*n2[1]}
            -- Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
            local u = n1[1]*x + n1[2]*y + n1[3]*z
            local v = n2[1]*x + n2[2]*y + n2[3]*z
            local w = n3[1]*x + n3[2]*y + n3[3]*z
            -- Pull the amplitude from the shuffled array index ("j"), not "i",
            -- otherwise neighboring unit vectors will have similar amplitudes!
            local a = amplitudes[j+1]
            -- Noise is the average of cosine of distance along each axis, shifted by offsets and scaled by amplitude.
            noise = noise + a*(cos(u/a + offsets[3*i+1]) + cos(v/a + offsets[3*i+2]) + cos(w/a + offsets[3*i+3]))/3
        end
        return cdf(noise/sigma)
    end
end

local shader1d = [[
#define TAU 6.283185307179586476925286766559005768394338798750211641949
#define MAX_RESOLUTION 64
extern int resolution;
extern float sigma, range_min, range_max;
extern float amplitudes[MAX_RESOLUTION];
extern float offsets[MAX_RESOLUTION];

vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
{
    float x = mix(range_min,range_max,texture_coords.x);
    float noise = 0.;
    for (int i=0; i < resolution; i++) {
        float a = amplitudes[i];
        noise += a*cos(x/a + offsets[i]);
    }
    noise = noise/sigma;
    noise = .5 + .5*sign(noise)*sqrt(1.-exp(-4./TAU * noise*noise));
    return vec4(noise,noise,noise,1.);
}
]]

noise.make1dShader = function(resolution,random,decayFn)
    if resolution > 64 then
        error("Resolution cannot exceed 64")
    end
    local shader = lg.newShader(shader1d)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,resolution,random)
    shader:send("sigma",sigma)
    shader:send("resolution",resolution)
    shader:send("offsets",unpack(offsets))
    shader:send("amplitudes",unpack(amplitudes))
    shader:send("range_min", 0)
    shader:send("range_max", 0)
    return shader
end

local shader2d = [[
#define TAU 6.283185307179586476925286766559005768394338798750211641949
#define PHI 1.618033988749894848204586834365638117720309179805762862135
#define MAX_RESOLUTION 64
extern int resolution;
extern float sigma;
extern float amplitudes[MAX_RESOLUTION];
extern vec2 offsets[MAX_RESOLUTION];
extern vec2 range_min, range_max;

vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
{
    vec2 pos = mix(range_min,range_max,texture_coords);
    float noise = 0.;
    for (int i=0; i < resolution; i++) {
        float angle = mod(float(i)*PHI, 1.)*TAU;
        float u = pos.x*cos(angle) - pos.y*sin(angle);
        float v = pos.x*cos(angle+TAU/4.) - pos.y*sin(angle+TAU/4.);
        float a = amplitudes[i];
        noise += a*mix(cos(u/a + offsets[i].x), cos(v/a + offsets[i].y), .5);
    }
    noise = noise/sigma;
    noise = .5 + .5*sign(noise)*sqrt(1.-exp(-4./TAU * noise*noise));
    return vec4(noise,noise,noise,1.);
}
]]

noise.make2dShader = function(resolution,random,decayFn)
    if resolution > 64 then
        error("Resolution cannot exceed 64")
    end
    local shader = lg.newShader(shader2d)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,2*resolution,random)
    sigma = sigma/sqrt(2)
    local offsets2 = {}
    for i=1,#offsets-1,2 do
        table.insert(offsets2, {offsets[i],offsets[i+1]})
    end
    shader:send("sigma",sigma)
    shader:send("resolution",resolution)
    shader:send("offsets",unpack(offsets2))
    shader:send("amplitudes",unpack(amplitudes))
    shader:send("range_min", {0,0})
    shader:send("range_max", {1,1})
    return shader
end

local shader3d = [[
#define TAU 6.283185307179586476925286766559005768394338798750211641949
#define PHI 1.618033988749894848204586834365638117720309179805762862135
#define LOG_PHI 0.481211825059603447497758913424368423135184334385660519660
#define SQRT5 2.236067977499789696409173668731276235440618359611525724270
#define MAX_RESOLUTION 64
extern int resolution;
extern float sigma, z;
extern float amplitudes[MAX_RESOLUTION];
extern vec3 offsets[MAX_RESOLUTION];
extern vec2 range_min, range_max;

// https://www.graphics.rwth-aachen.de/media/papers/jgt.pdf
vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
{
    vec3 pos = vec3(mix(range_min,range_max,texture_coords), z);
    // Find the biggest fibonacci number F_n such that F_n < RESOLUTION
    int n = int(log((float(resolution)-1.)*SQRT5 + .5)/LOG_PHI);
    int dec = int(.5 + pow(PHI,n)/SQRT5); // F_n, using closed form Fibonacci
    int inc = int(.5 + dec/PHI); // F_(n-1)

    float noise = 0.;
    for (int i=0, j=0; i<resolution; ++i) {
        if (j >= dec) {
            j -= dec;
        } else {
            j += inc;
            if (j >= resolution)
                j -= dec;
        }
        // Convert golden ratio sequence into polar coordinate unit vector
        float phi = mod(float(i)*PHI,1.)*TAU;
        float theta = acos(mix(-1.,1.,mod(float(j)*PHI,1.)));
        // Make an orthonormal basis, where n1 is from polar phi/theta,
        // n2 is roated 90 degrees along phi, and n3 is the cross product of the two
        vec3 n1 = vec3(sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi));
        vec3 n2 = vec3(sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.));
        vec3 n3 = cross(n1,n2);
        // Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
        float u = dot(n1, pos);
        float v = dot(n2, pos);
        float w = dot(n3, pos);
        // Pull the amplitude from the shuffled array index ("j"), not "i",
        // otherwise neighboring unit vectors will have similar amplitudes!
        float a = amplitudes[j];
        //float a = pow(mod(float(i+1)*(PHI-1.), 1.), .3);
        // Noise is the average of cosine of distance along each axis, shifted by offsets and scaled by amplitude.
        noise += a*(cos(u/a + offsets[i].x) + cos(v/a + offsets[i].y) + cos(w/a + offsets[i].z))/3.;
    }
    noise = noise/sigma;
    noise = .5 + .5*sign(noise)*sqrt(1.-exp(-4./TAU * noise*noise));
    return vec4(noise,noise,noise,1.);
}
]]

noise.make3dShader = function(resolution,random,decayFn)
    if resolution > 64 then
        error("Resolution cannot exceed 64")
    end
    local shader = lg.newShader(shader3d)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,3*resolution,random)
    sigma = sigma/sqrt(3)
    local offsets2 = {}
    for i=1,#offsets-1,3 do
        table.insert(offsets2, {offsets[i],offsets[i+1],offsets[i+2]})
    end
    shader:send("sigma",sigma)
    shader:send("resolution",resolution)
    shader:send("offsets",unpack(offsets2))
    shader:send("amplitudes",unpack(amplitudes))
    shader:send("range_min", {0,0})
    shader:send("range_max", {1,1})
    return shader
end

return noise