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: 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