Implementation of BTRD binomial random varirate generation algo from http://epub.wu.ac.at/1242/1/document.pdf . THIS IS EXPLOSIVE CODE IN PUBLIC DOMAIN. USE WITH CAUTION.

BtrdAlgorithm.rb

# When the mean is still small, but the probability of success
# is not miniscule. Covers n > 20
# Taken from "The Generation of Binomial Random Variates"
# Wolfgang Hormann (http://epub.wu.ac.at/1242/1/document.pdf)
# This is algorithm BTRD from the text

# This is a precalculated table for correction terms in
# stirling's approximation to log(k!) each index corresponds
# to a seperate value of k
fc_table = [0.08106146679532726, 0.04134069595540929,
0.02767792568499834, 0.02079067210376509, 0.01664469118982119,
0.01387612882307075, 0.01189670994589177, 0.01041126526197209,
0.009255462182712733, 0.008330563433362871]

# lambda for larger values
fc = lambda do |d|
  next fc_table[d] if d < 10
  dp1 = d + 1
  dp1_sq = dp1 ** 2
  # (1/12 - (1/360 - 1/1260x**2)/x**2)/x; where x is d + 1
  # Boy does this expression raise a lot of cool warnings...
  ((1.0 /12 - (1.0 / 360 - (1.0 / 1260) / dp1_sq)) / dp1_sq) / dp1
end
m = ((n + 1) * prob).floor
r = prob / (1 - prob)
nr = (n + 1) * r
npq = n * prob * (1 - prob)
b = 1.15 + 2.53 * Math.sqrt(npq)
a = -0.0873 + 0.0248 * b + 0.01 * prob
c = n * prob + 0.5
alpha = (2.83 + 5.1 / b) * Math.sqrt(npq)
vr = 0.92 - 4.2 / b
urvr = 0.86 * vr
u = 0
v = 0

loop do
  # Steps 1, 2 and first half of 3.0
  loop do
    v = prng.rand
    if v <= urvr
      u = v / vr - 0.43
      return ((((2 * a) / (0.5 - u.abs)) + b) * u + c).floor
    elsif v >= vr
      u = prng.rand - 0.5
    else
      u = v / vr - 0.93
      u = (u <=> 0) * 0.5 - u
      v = vr * prng.rand
    end
    
    us = 0.5 - u.abs
    k = ((((2 * a) / us) + b) * u + c).floor
    
    # Equivalent to the goto directive talked about in the paper
    # This restarts the proess if k is out of bounds
    break unless k < 0 || k > n
  end
  
  # Second half of step 3.0
  v = v * alpha / (b + a / (us ** 2))
  km = (k - m).abs
  
  # Step 3.1
  f = 1.0
  if m < k
    m.upto(k) {|i| f = f * (nr / i - r)}
  elsif m > k
    k.upto(m) {|i| v = v * (nr / i - r)}
  end
  
  # This is yet another goto-1
  # Return value if condition fulfilled, else goto-1
  # Since the next part is a 
  if v <= f
    return k
  
  # Step 3.2
  if km > 15
    v = log(v)
    rho = (km / npq) * (((km / 3.0 + 0.625) * km + 1.0 / 6) / npq + 0.5)
    t = - (km ** 2) / (2 * npq)
    return k if v < t - rho
    
    # If no value was returned, then we need to break loop
    # Paper says goto-1 if v > t + rho
    break if v <= (t + rho)
  end
  
  # Step 3.3 of the paper
  nm = n - m + 1
  h = (m + 0.5) * log((m + 1) / (r * nm)) + fc.call(m) + fc.call(n - m)
  
  # Step 3.4
  nk = n - k + 1
  if v <= h + (n + 1) * log(nm / nk) + (k + 0.5) * log(nk * r / (k + 1)) -fc.call(k) - fc.call(n - k)
    return k
  else
    # GO TO 1 NOW
    continue
  end
end