math.special = math.special or {}
do
    local special = math.special

    function special.factorial(n)
        if n % 1 ~= 0 then return nil end
        if n == 0 or n == 1 then return 1 end
        return n*factorial(n-1)  
    end
    
    -- используем приближение Ланцоша
    function special.gamma(x)
        local g = 7
        local p = {
            0.99999999999980993,
            676.5203681218851,
           -1259.1392167224028,
            771.32342877765313,
           -176.61502916214059,
            12.507343278686905,
           -0.13857109526572012,
            9.9843695780195716e-6,
            1.5056327351493116e-7
        }
    
        if x < 0.5 then
            -- Формула дополнения
            return math.pi / (math.sin(math.pi * x) * special.gamma(1 - x))
        else
            x = x - 1
            local a = p[1]
            for i = 2, #p do
                a = a + p[i] / (x + i - 1)
            end
            local t = x + g + 0.5
            return math.sqrt(2 * math.pi) * t^(x + 0.5) * math.exp(-t) * a
        end
    end
    
        
   -- Функция для вычисления логарифма гамма-функции
    function special.loggamma(k)
        -- Для целых значений k используем приближение Стирлинга
        if k < 8 then
            -- Прямое вычисление для небольших k
            local res = 0
            for i = 2, k do
                res = res + math.log(i)
            end
            return res
        else
            -- Приближение Стирлинга
            local k_float = k
            local term1 = (k_float + 0.5) * math.log(k_float)
            local term2 = -k_float + 0.5 * math.log(2 * math.pi)
            local term3 = 1 / (12 * k_float)
            return term1 + term2 + term3
        end
    end
    
        -- Функция erf(x)
    function special.erf(x)
        -- Константы
        local a1 =  0.254829592
        local a2 = -0.284496736
        local a3 =  1.421413741
        local a4 = -1.453152027
        local a5 =  1.061405429
        local p  =  0.3275911
    
        -- Сохранение знака x
        local sign = 1
        if x < 0 then
            sign = -1
        end
        x = math.abs(x)
    
        -- Вычисление по формуле
        local t = 1.0 / (1.0 + p * x)
        local y = 1.0 - (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t) * math.exp(-x * x)
    
        return sign * y
    end

    
    -- Функция beta(a, b)
    function special.beta(a, b)
        return math.exp(special.loggamma(a) + special.loggamma(b) - special.loggamma(a + b))
    end
    
end

do -- reqire "math.special"
    local uniform = {}
    uniform.rvs = (type(math.random) == "function") and math.random or math.uniform.rvs
    
    local random = {
      uniform = uniform,
      poisson = {},
      exponential = {},
      normal = {}
    }
    setmetatable(random, {__call = function(tbl,...) return random.uniform.rvs(...) end})
    
    local meta_distr = {
      __call = function(distr,...)
          return distr.rvs(...)
      end
    }
  
    for key, val in pairs(random) do
        setmetatable(val, meta_distr)
    end
    
    math.random = random
end


do -- poisson distribution
    local poisson = math.random.poisson
    local loggamma = math.special.loggamma
  
    function poisson.ptrs(lam)
      local slam = math.sqrt(lam)
      local loglam = math.log(lam)
      local b = 0.931 + 2.53 * slam
      local a = -0.059 + 0.02483 * b
      local invalpha = 1.1239 + 1.1328 / (b - 3.4)
      local vr = 0.9277 - 3.6224 / (b - 2)
      
      while true do
          local U = math.random() - 0.5
          local V = math.random()
          local us = 0.5 - math.abs(U)
          local k = math.floor((2 * a / us + b) * U + lam + 0.43)
          if (us >= 0.07) and (V <= vr) then
              return k
          end
          if (k < 0) or ((us < 0.013) and (V > us)) then
              -- Переходим к следующей итерации цикла
          else
              local lhs = math.log(V) + math.log(invalpha) - math.log(a / (us * us) + b)
              local rhs = -lam + k * loglam - loggamma(k + 1)
              if lhs <= rhs then
                  return k
              end
          end
        end
    end
    
    poisson.rvs = function(lambda) 
        return poisson.ptrs(lambda)
    end 
end

do -- exponential distribution
    local exponential = math.random.exponential
    
    function exponential.standart(lambda)
        local u = math.random()
        return - math.log(u) * lambda
    end
    
    exponential.rvs = function(lambda) 
        return exponential.standart(lambda)
    end
end


do -- normal distrebution
    local normal = math.random.normal
    
    function normal.box_muller()
        local u1 = math.random()
        local u2 = math.random()
        return math.sqrt(-2 * math.log(u1)) * math.cos(2 * math.pi * u2)
    end
    
    function normal.rvs(mean, std)
        mean = mean or 0
        std = std or 1
        return normal.box_muller()*std + mean 
    end
    
end