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

noise.make1dShader = function(resolution,random,decayFn)
    local shader = lg.newShader("noise1d.glsl")
    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

noise.make2dShader = function(resolution,random,decayFn)
    local shader = lg.newShader("noise2d.glsl")
    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

noise.make3dShader = function(resolution,random,decayFn)
    local shader = lg.newShader("noise3d.glsl")
    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