module MoreMath::Functions

def beta(a, b)
  if a > 0 && b > 0
    exp(log_beta(a, b))
  else
    0.0 / 0
  end
end

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: epsilon, max_iterations: max_iterations)
  else
    fraction = ContinuedFraction.for_b do |n, y|
      if n % 2 == 0
        m = n / 2.0
        (m * (b - m) * y) / ((a + (2 * m) - 1) * (a + (2 * m)))
      else
        m = (n - 1) / 2.0
        -((a + m) * (a + b + m) * y) / ((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: epsilon, max_iterations: 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 erfc(x)
  1.0 - erf(x)
end

def gamma(x)
  if x < 0.0
    return PI / (sin(PI * x) * exp(log_gamma(1 - x)))
  else
    exp(log_gamma(x))
  end
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: epsilon, max_iterations: 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: epsilon, max_iterations: max_iterations)
  else
    fraction = ContinuedFraction.for_a do |n, y|
      (2 * n + 1) - a + y
    end.for_b do |n, y|
      n * (a - n)
    end
    exp(-x + a * log(x) - log_gamma(a)) *
      fraction[x, epsilon: epsilon, max_iterations: 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_ceil(n, b = 2)
  raise ArgumentError, "n is required to be > 0" unless n > 0
  raise ArgumentError, "b is required to be > 1" unless b > 1
  e, result = 1, 0
  until e >= n
    e *= b
    result += 1
  end
  result
end

def log_floor(n, b = 2)
  raise ArgumentError, "n is required to be > 0" unless n > 0
  raise ArgumentError, "b is required to be > 1" unless b > 1
  e, result = 1, 0
  until e * b > n
    e *= b
    result += 1
  end
  result
end

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

def log_gamma(x)
  x = x.to_f
  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

def logb(x, b = 2)
  Math.log(x) / Math.log(b)
end

def numberify_string(string, alphabet = 'a'..'z')
  NumberifyStringFunction.numberify_string(string, alphabet)
end

def stringify_number(number, alphabet = 'a'..'z')
  NumberifyStringFunction.stringify_number(number, alphabet)
end