local noise = {}
local exp,sin,cos,floor,log,acos,sqrt,abs = math.exp,math.sin,math.cos,math.floor,math.log,math.acos,math.sqrt,math.abs
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)
    return .5 + .5*(x<0 and -1 or 1)*sqrt(1.-exp(-4./TAU * x*x))
end

local function cdf_prime(x, dx)
    return (0.31831 * exp(-2/PI * x*x) * abs(x)*dx)/sqrt(1-exp(-2/PI*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
        assert(a > 0.0001)
        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

local function sign(x)
    if x == 0 then return 0
    elseif x < 0 then return -1 else return 1 end
end
noise.make1d = function(resolution,random,decayFn)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,resolution,random)
    local function noise(x)
        local n = 0
        for i,a in ipairs(amplitudes) do
            n = n + a*cos(x/a + offsets[i])
        end
        n = n/sigma
        return cdf(n)
    end

    local function gradient(x)
        local n, dndx = 0, 0
        for i,a in ipairs(amplitudes) do
            n = n + a*cos(x/a + offsets[i])
            dndx = dndx - sin(x/a + offsets[i])
        end
        dndx = dndx/sigma
        n = n/sigma
        return cdf_prime(n,dndx)
    end
    return noise, gradient
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)
    local function noise(x,y)
        local n = 0
        for i,a in ipairs(amplitudes) do
            local angle = ((i*GR) % 1)*TAU
            local sina, cosa = sin(angle), cos(angle)
            local u = x*cosa - y*sina
            local v = -x*sina - y*cosa
            local k, w = offsets[2*i], offsets[2*i-1]
            n = n + a*(cos(u/a + k) + cos(v/a + w))
        end
        n = n/(2*sigma)
        return cdf(n)
    end
    local function gradient(x,y)
        local n, dx, dy = 0,0,0
        for i,a in ipairs(amplitudes) do
            local angle = ((i*GR) % 1)*TAU
            local sina, cosa = sin(angle), cos(angle)
            local u = x*cosa - y*sina
            local v = -x*sina - y*cosa
            local k, w = offsets[2*i], offsets[2*i-1]
            n = n + a*(cos(u/a + k) + cos(v/a + w))
            dx = dx + (sina*sin(v/a + w) - cosa*sin(u/a + k))
            dy = dy + (sina*sin(u/a + k) + cosa*sin(v/a + w))
        end
        n = n/(2*sigma)
        dx = dx/(2*sigma)
        dy = dy/(2*sigma)
        return cdf_prime(n,dx), cdf_prime(n,dy)
    end
    return noise, gradient
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)
    local function noise(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 n,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 n1x,n1y,n1z = sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)
            local n2x,n2y,n2z = sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.)
            -- Cross product
            local n3x,n3y,n3z = n1y*n2z - n1z*n2y,
                        n1z*n2x - n1x*n2z,
                        n1x*n2y - n1y*n2x
            -- Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
            local u = n1x*x + n1y*y + n1z*z
            local v = n2x*x + n2y*y + n2z*z
            local w = n3x*x + n3y*y + n3z*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.
            n = n + a*(cos(u/a + offsets[3*i+1]) + cos(v/a + offsets[3*i+2]) + cos(w/a + offsets[3*i+3]))
        end
        return cdf(n/(3*sigma))
    end
    local function gradient(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 n,i,j = 0,0,0
        local dndx,dndy,dndz = 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 n1x,n1y,n1z = sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)
            local n2x,n2y,n2z = sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.)
            -- Cross product
            local n3x,n3y,n3z = n1y*n2z - n1z*n2y,
                        n1z*n2x - n1x*n2z,
                        n1x*n2y - n1y*n2x
            -- Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
            local u = n1x*x + n1y*y + n1z*z
            local v = n2x*x + n2y*y + n2z*z
            local w = n3x*x + n3y*y + n3z*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.
            n = n + a*(cos(u/a + offsets[3*i+1]) + cos(v/a + offsets[3*i+2]) + cos(w/a + offsets[3*i+3]))

            local kx,ky,kz = -sin(u/a+offsets[3*i+1]),-sin(v/a+offsets[3*i+2]),-sin(w/a + offsets[3*i+3])
            dndx = dndx + (n1x*kx + n2x*ky + n3x*kz)
            dndy = dndy + (n1y*kx + n2y*ky + n3y*kz)
            dndz = dndz + (n1z*kx + n2z*ky + n3z*kz)
        end
        n = n / (3*sigma)
        return cdf_prime(n, dndx/(3*sigma)), cdf_prime(n, dndy/(3*sigma)), cdf_prime(n, dndz/(3*sigma))
    end
    return noise, gradient
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.;
#ifdef GRADIENT
    float dndx = 0.;
#endif
    for (int i=0; i < resolution; i++) {
        float a = amplitudes[i];
        noise += a*cos(x/a + offsets[i]);
#ifdef GRADIENT
        dndx -= sin(x/a + offsets[i]);
#endif
    }
    noise /= sigma;
#ifdef GRADIENT
    dndx /= sigma;
    dndx = (0.31831 * exp(-4./TAU * noise*noise) * abs(noise)*dndx)/sqrt(1.0-exp(-4./TAU * noise*noise));
    // TODO: normalize properly
    dndx = .5 + .5*dndx;
    return vec4(dndx,dndx,dndx, 1.);
#else
    noise = .5 + .5*sign(noise)*sqrt(1.-exp(-4./TAU * noise*noise));
    return vec4(noise,noise,noise, 1.);
#endif
}
]]

noise.make1dShader = function(resolution,random,decayFn)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    if type(resolution) ~= 'number' or (resolution % 1 > 0) or (resolution <= 0) then
        error("Resolution must be a positive integer.")
    end
    local shaderCode = "#define RESOLUTION "..tostring(resolution).."\n"..shader1d
    local noiseShader = lg.newShader(shaderCode)
    local gradShader = lg.newShader("#define GRADIENT\n"..shaderCode)
    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,resolution,random)
    do -- Dumb hack to work around a bug in Love 0.10.2 not sending the last value
        table.insert(amplitudes, 1.)
        table.insert(offsets, 1.)
    end
    for _,shader in ipairs{noiseShader, gradShader} do
        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)
    end
    return noiseShader,gradShader
end

local shader2d = [[
#define TAU 6.283185307179586476925286766559005768394338798750211641949
#define PHI 1.618033988749894848204586834365638117720309179805762862135
extern float sigma;
extern float amplitudes[RESOLUTION];
extern vec2 offsets[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.;
#ifdef GRADIENT
    float dndx = 0., dndy = 0.;
#endif
    for (int i=0; i < RESOLUTION; i++) {
        float angle = mod(float(i)*PHI, 1.)*TAU;
        float cosa = cos(angle), sina = sin(angle);
        float u = pos.x*cosa - pos.y*sina;
        float v = -pos.x*sina - pos.y*cosa;
        float a = amplitudes[i];
        float k = offsets[i].x, w = offsets[i].y;
        noise += a*(cos(u/a + k) + cos(v/a + w));
#ifdef GRADIENT
        dndx += -cosa*sin(u/a + k) + sina*sin(v/a + w);
        dndy += sina*sin(u/a + k) + cosa*sin(v/a + w);
#endif
    }
    noise /= 2.*sigma;

#ifdef GRADIENT
    dndx /= 2.*sigma;
    dndx = (0.31831 * exp(-4./TAU * noise*noise) * abs(noise)*dndx)/sqrt(1.0-exp(-4./TAU * noise*noise));
    dndx = .5 + .5*dndx;

    dndy /= 2.*sigma;
    dndy = (0.31831 * exp(-4./TAU * noise*noise) * abs(noise)*dndy)/sqrt(1.0-exp(-4./TAU * noise*noise));
    dndy = .5 + .5*dndy;

    return vec4(dndx,dndy,0.,1.);
#else
    noise = .5 + .5*sign(noise)*sqrt(1.-exp(-4./TAU * noise*noise));
    return vec4(noise,noise,noise,1.);
#endif
}
]]

noise.make2dShader = function(resolution,random,decayFn)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    if type(resolution) ~= 'number' or (resolution % 1 > 0) or (resolution <= 0) then
        error("Resolution must be a positive integer.")
    end
    local shaderCode = "#define RESOLUTION "..tostring(resolution).."\n"..shader2d
    local noiseShader = lg.newShader(shaderCode)
    local gradShader = lg.newShader("#define GRADIENT\n"..shaderCode)
    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
    do -- Dumb hack to work around a bug in Love 0.10.2 not sending the last value
        table.insert(amplitudes, 1.)
        table.insert(offsets2, {1,1})
    end
    for _,shader in ipairs{noiseShader, gradShader} do
        shader:send("sigma",sigma)
        shader:send("offsets",unpack(offsets2))
        shader:send("amplitudes",unpack(amplitudes))
        shader:send("range_min", {0,0})
        shader:send("range_max", {1,1})
    end
    return noiseShader,gradShader
end

local shader3d = [[
#define TAU 6.283185307179586476925286766559005768394338798750211641949
#define PHI 1.618033988749894848204586834365638117720309179805762862135
#define LOG_PHI 0.481211825059603447497758913424368423135184334385660519660
#define SQRT5 2.236067977499789696409173668731276235440618359611525724270
extern float sigma, z;
extern float amplitudes[RESOLUTION];
extern vec3 offsets[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 fib_n = int(log((float(RESOLUTION)-1.)*SQRT5 + .5)/LOG_PHI);
    int dec = int(.5 + pow(PHI,fib_n)/SQRT5); // F_n, using closed form Fibonacci
    int inc = int(.5 + dec/PHI); // F_(fib_n-1)

    float n = 0.;
#ifdef GRADIENT
    float dndx = 0., dndy = 0., dndz = 0.;
#endif
    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.
        n += a*(cos(u/a + offsets[i].x) + cos(v/a + offsets[i].y) + cos(w/a + offsets[i].z));
#ifdef GRADIENT
        vec3 k = vec3(-sin(u/a+offsets[i].x),-sin(v/a+offsets[i].y),-sin(w/a + offsets[i].z));
        dndx += (n1.x*k.x + n2.x*k.y + n3.x*k.z)/3.;
        dndy += (n1.y*k.x + n2.y*k.y + n3.y*k.z)/3.;
        dndz += (n1.z*k.x + n2.z*k.y + n3.z*k.z)/3.;
#endif
    }
    n /= 3.*sigma;
#ifdef GRADIENT
    dndx /= sigma;
    dndx = (0.31831 * exp(-4./TAU * n*n) * abs(n)*dndx)/sqrt(1.0-exp(-4./TAU * n*n));
    dndx = .5 + .5*dndx;

    dndy /= sigma;
    dndy = (0.31831 * exp(-4./TAU * n*n) * abs(n)*dndy)/sqrt(1.0-exp(-4./TAU * n*n));
    dndy = .5 + .5*dndy;

    dndz /= sigma;
    dndz = (0.31831 * exp(-4./TAU * n*n) * abs(n)*dndz)/sqrt(1.0-exp(-4./TAU * n*n));
    dndz = .5 + .5*dndz;

    return vec4(dndx,dndy,dndz, 1.);
#else
    n = .5 + .5*sign(n)*sqrt(1.-exp(-4./TAU * n*n));
    return vec4(n,n,n,1.);
#endif
}
]]

noise.make3dShader = function(resolution,random,decayFn)
    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
    if type(resolution) ~= 'number' or (resolution % 1 > 0) or (resolution <= 0) then
        error("Resolution must be a positive integer.")
    end
    local shaderCode = "#define RESOLUTION "..tostring(resolution).."\n"..shader3d
    local noiseShader = lg.newShader(shaderCode)
    local gradShader = lg.newShader("#define GRADIENT\n"..shaderCode)
    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
    do -- Dumb hack to work around a bug in Love 0.10.2 not sending the last value
        table.insert(amplitudes, 1.)
        table.insert(offsets, {1,1,1})
    end
    for _,shader in ipairs{noiseShader, gradShader} do
        shader:send("sigma",sigma)
        shader:send("offsets",unpack(offsets2))
        shader:send("amplitudes",unpack(amplitudes))
        shader:send("range_min", {0,0})
        shader:send("range_max", {1,1})
    end
    return noiseShader, gradShader
end

return noise