Bullshit::Functions

# File lib/bullshit.rb, line 726
    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

# File lib/bullshit.rb, line 811
    def erf(x)
      r = sqrt(1 - exp(-x ** 2 * (4 / Math::PI + A * x ** 2) / (1 + A * x ** 2)))
      x < 0 ? -r : r
    end

# File lib/bullshit.rb, line 753
    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

# File lib/bullshit.rb, line 784
    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

# File lib/bullshit.rb, line 717
    def log_beta(a, b)
      log_gamma(a) + log_gamma(b) - log_gamma(a + b)
    rescue Errno::ERANGE, Errno::EDOM
      0 / 0.0
    end

# File lib/bullshit.rb, line 699
      def log_gamma(x)
        if x.nan? || x <= 0
          0 / 0.0
        else
          sum = 0.0
          (LANCZOS_COEFFICIENTS.size - 1).downto(1) do |i|
            sum += LANCZOS_COEFFICIENTS[i] / (x + i)
          end
          sum += LANCZOS_COEFFICIENTS[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

# File lib/bullshit.rb, line 695
      def log_gamma(x)
        lgamma(x).first
      end