module MoreMath::Functions

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