module MoreMath::Functions

def gammaQ_regularized(x, a, epsilon = 1E-16, max_iterations = 1 << 16)

+max_iterations+-times.
+a+ with an error of <= +epsilon+, but only iterate
Return an approximation of the regularized gammaQ function for +x+ and 
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