module MoreMath::Functions

def beta_regularized(x, a, b, epsilon = 1E-16, max_iterations = 1 << 16)
  x, a, b = x.to_f, a.to_f, b.to_f
  case
  when a.nan? || b.nan? || x.nan? || a <= 0 || b <= 0 || x < 0 || x > 1
    0 / 0.0
  when x > (a + 1) / (a + b + 2)
    1 - beta_regularized(1 - x, b, a, epsilon, max_iterations)
  else
    fraction = ContinuedFraction.for_b do |n, x|
      if n % 2 == 0
        m = n / 2.0
        (m * (b - m) * x) / ((a + (2 * m) - 1) * (a + (2 * m)))
      else
        m = (n - 1) / 2.0
        -((a + m) * (a + b + m) * x) / ((a + 2 * m) * (a + 2 * m + 1))
      end
    end
    exp(a * log(x) + b * log(1.0 - x) - log(a) - log_beta(a, b)) / 
      fraction[x, epsilon, max_iterations]
  end
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

def cantor_pairing(*xs)
  CantorPairingFunction.cantor_pairing(*xs)
end

def cantor_pairing_inv(c, n = 2)
  CantorPairingFunction.cantor_pairing_inv(c, n)
end

def erf(x)
  erf_a = MoreMath::Constants::FunctionsConstants::ERF_A
  r = sqrt(1 - exp(-x ** 2 * (4 / Math::PI + erf_a * x ** 2) / (1 + erf_a * x ** 2)))
  x < 0 ? -r : r
end

def erf(x)
  Math.erf(x)
end

def gammaP_regularized(x, a, epsilon = 1E-16, max_iterations = 1 << 16)
  x, a = x.to_f, a.to_f
  case
  when a.nan? || x.nan? || a <= 0 || x < 0
    0 / 0.0
  when x == 0
    0.0
  when 1 <= a && a < x
    1 - gammaQ_regularized(x, a, epsilon, max_iterations)
  else
    n = 0
    an = 1 / a
    sum = an
    while an.abs > epsilon && n < max_iterations
      n += 1
      an *= x / (a + n)
      sum += an
    end
    if n >= max_iterations
      raise Errno::ERANGE
    else
      exp(-x + a * log(x) - log_gamma(a)) * sum
    end
  end
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

def gammaQ_regularized(x, a, epsilon = 1E-16, max_iterations = 1 << 16)
  x, a = x.to_f, a.to_f
  case
  when a.nan? || x.nan? || a <= 0 || x < 0
    0 / 0.0
  when x == 0
    1.0
  when a > x || a < 1
    1 - gammaP_regularized(x, a, epsilon, max_iterations)
  else
    fraction = ContinuedFraction.for_a do |n, x|
      (2 * n + 1) - a + x
    end.for_b do |n, x|
      n * (a - n)
    end
    exp(-x + a * log(x) - log_gamma(a)) *
      fraction[x, epsilon, max_iterations] ** -1
  end
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

def log_beta(a, b)
  log_gamma(a) + log_gamma(b) - log_gamma(a + b)
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

def log_gamma(x)
  lgamma(x).first
end

def log_gamma(x)
  if x.nan? || x <= 0
    0 / 0.0
  else
    sum = 0.0
    lc = Constants::FunctionsConstants::LANCZOS_COEFFICIENTS
    half_log_2_pi = Constants::FunctionsConstants::HALF_LOG_2_PI
    (lc.size - 1).downto(1) do |i|
      sum += lc[i] / (x + i)
    end
    sum += lc[0]
    tmp = x + 607.0 / 128 + 0.5
    (x + 0.5) * log(tmp) - tmp + half_log_2_pi + log(sum / x)
  end
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end