[Biblioteca] SciLua Math

--Creator: Bolodefchoco
--Made in: 17/07/2016
--Last update: 17/07/2016
--[[ Notes:
    math.round
        Does:
            Arredonda pro número inteiro mais próximo
        Args:
            x --> Número
    math.step
        Does:
            {1x>=0/0x<0}
        Args:
            x --> Número
    math.sign
        Does:
            {1x>=0/-1x<0}
        Args:
            x --> Número
    math.phi
        Does:
            Φ(x)=(∫x−∞)(1/√2π)exp{−(t^2/2)}dt
        Args:
            x --> Número
    math.iphi
        Does:
            (IΦ)(x)=Φ−1(x)
        Args:
            x --> Número
    math.gamma
        Does:
            Γ(z)=(∫∞0)e^−t t^(z−1)dt
        Args:
            z --> Número
    math.loggamma
        Does:
            Valor absoluto do valor da função math.gamma
        Args:
            z --> Número
    math.logbeta
        Does:
            Logaritmo da função math.beta
        Args:
            a --> Número
            b --> Número
    math.beta
        Does:
        B(a,b)=(Γ(a)Γ(b))/Γ(a+b)

        Args:
            a --> Número
            b --> Número
]]--


math.round = function(x)
    return x<0 and math.ceil(x-.5) or math.floor(x+.5)
end

math.step = function(x)
    return math.max(0,math.min(math.floor(x)-1,1))
end

math.sign = function(x)
    return -1 + math.step(x)*2
end

math.phi = function(x)
    if x<=-8 then return 0 elseif x>=8 then return 1 else
        local s,b,q = x,x,x^2
        for i = 3,1/0,2 do
            b = b*q/i
            local t = s
            s = t + b
            if s == t then break end
        end
        return .5+s*math.exp(-.5*q-.91893853320467274178)
    end
end

local iphifast,iphi
do
    local a = {
        [1] = 0,
        [2] = -3.969683028665376e+01,
        [3] = 2.209460984245205e+02,
        [4] = -2.759285104469687e+02,
        [5] = 1.383577518672690e+02,
        [6] = -3.066479806614716e+01,
        [7] = 2.506628277459239e+00
    }
 local b = {
        [1] = 0,
        [2] = -5.447609879822406e+01,
        [3] = 1.615858368580409e+02,
        [4] = -1.556989798598866e+02,
        [5] = 6.680131188771972e+01,
        [6] = -1.328068155288572e+01
    }
  local c = {
        [1] = 0,
        [2] = -7.784894002430293e-03,
        [3] = -3.223964580411365e-01,
        [4] = -2.400758277161838e+00,
        [5] = -2.549732539343734e+00,
        [6] = 4.374664141464968e+00,
        [7] = 2.938163982698783e+00
    }
  local d = {
        [1] = 0,
        [2] = 7.784695709041462e-03,
        [3] = 3.224671290700398e-01,
        [4] = 2.445134137142996e+00,
        [5] = 3.754408661907416e+00
    }
    math.iphifast = function(p)
        if math.abs(p-.5) < .47575 then
            local q = p-.5
            local r = q^2
            return (((((a[1]*r+a[2])*r+a[3])*r+a[4])*r+a[5])*r+a[6])*q/(((((b[1]*r+b[2])*r+b[3])*r+b[4])*r+b[5])*r+1)
        else
            local iu = math.ceil(p-.97575)
            local z = (1-iu)*p+iu*(1-p)
            local sign = 1-2*iu
            local q = math.sqrt(-2*math.log(z))
            return sign*(((((c[1]*q+c[2])*q+c[3])*q+c[4])*q+c[5])*q+c[6])/((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1)
        end
    end
    math.iphi = function(p)
        if p<=0 then
            return -1/0
        elseif p>=1 then
            return 1/0
        else
            local x = math.iphifast(p)
            local e = math.phi(x)-p
            local u = e*math.sqrt(2*math.pi)*math.exp(x^2/2)
            return x-u/(1+x*u/2)
        end
    end
    local gamma,loggamma
end

do
    local gamma_r10 = 10.900511
    local gamma_dk = {
        [1] = 2.48574089138753565546e-5,
        [2] = 1.05142378581721974210,
        [3] = -3.45687097222016235469,
        [4] = 4.51227709466894823700,
        [5] = -2.98285225323576655721,
        [6] = 1.05639711577126713077,
        [7] = -1.95428773191645869583e-1,
        [8] = 1.70970543404441224307e-2,
        [9] = -5.71926117404305781283e-4,
        [10] = 4.63399473359905636708e-6,
        [11] = -2.71994908488607703910e-9
  }
    local gamma_c = 2*math.sqrt(math.exp(1)/math.pi)
    math.gamma = function(z)
        if z<0 then
            return math.pi/(math.sin(math.pi*z)*math.gamma(1-z))
        end
        local sum = gamma_dk[0];sum = sum + gamma_dk[1]/(z+0);sum = sum + gamma_dk[2]/(z+1);sum = sum + gamma_dk[3]/(z+2);sum = sum + gamma_dk[4]/(z+3);sum = sum + gamma_dk[5]/(z+4);sum = sum + gamma_dk[6]/(z+5);sum = sum + gamma_dk[7]/(z+6);sum = sum + gamma_dk[8]/(z+7);sum = sum + gamma_dk[9]/(z+8);sum = sum + gamma_dk[10]/(z+9)
        return gamma_c*((z+gamma_r10-.5)/math.exp(1))^(z-.5)*sum
    end
    math.loggamma = function(z)
        if z<0 then
            return math.log(math.pi)-math.log(math.abs(math.sin(math.pi*z)))-math.loggamma(1-z)
        end
        local sum = gamma_dk[0];sum = sum + gamma_dk[1]/(z+0);sum = sum + gamma_dk[2]/(z+1);sum = sum + gamma_dk[3]/(z+2);sum = sum + gamma_dk[4]/(z+3);sum = sum + gamma_dk[5]/(z+4);sum = sum + gamma_dk[6]/(z+5);sum = sum + gamma_dk[7]/(z+6);sum = sum + gamma_dk[8]/(z+7);sum = sum + gamma_dk[9]/(z+8);sum = sum + gamma_dk[10]/(z+9)
        return math.log(gamma_c)+(z-.5)*math.log(z+gamma_r10-.5)-(z-.5)+math.log(sum)
    end
end

local logbeta = function(a,b)
  if a<=0 or b<=0 then return 0/0 end
  return math.loggamma(a)+math.loggamma(b)-math.loggamma(a+b)
end

local beta = function(a,b)
  return math.exp(logbeta(a,b))
end