Utility for calculating the cumulative binomial distribution in ruby

cumulative_binomial_distribution.rb

module CumulativeBinomialDistribution
  extend self

  # s = successes
  # p = probability of success per trial
  # n = number of trials
  def normal_estimate(s, p, n)
    u = n * p
    o = (u * (BigDecimal.new('1')-p)) ** BigDecimal.new('0.5')
    # continuity correction
    s = s + BigDecimal.new('0.5')

    return BigDecimal.new('0.5') * (BigDecimal.new('1') + erf((s-u)/(o * BigDecimal.new('2') ** BigDecimal.new('0.5'))))
  end

  def erf(z)
    t = BigDecimal.new('1') / (BigDecimal.new('1') + BigDecimal.new('0.5') * z.abs)
    # use Horner's method
    ans = BigDecimal.new('1') - t * BigMath.exp( -z*z - BigDecimal.new('1.26551223') +
                            t * ( BigDecimal.new('1.00002368') +
                            t * ( BigDecimal.new('0.37409196') +
                            t * ( BigDecimal.new('0.09678418') +
                            t * ( BigDecimal.new('-0.18628806') +
                            t * ( BigDecimal.new('0.27886807') +
                            t * ( BigDecimal.new('-1.13520398') +
                            t * ( BigDecimal.new('1.48851587') +
                            t * ( BigDecimal.new('-0.82215223') +
                            t * ( BigDecimal.new('0.17087277')))))))))), 30)
    if z >= BigDecimal.new('0')
      return ans
    else
      return -ans
    end
  end
end