#!/usr/bin/ruby

func P(n, k) { n! / ((n-k)!) }
func C(n, k) { binomial(n, k) }
 
class Logarithm(value) {
    method to_s {
        var e = int(value/10.log)
        "%.8fE%+d" % (exp(value - e*10.log), e)
    }
}
 
func lstirling(n) {
    n < 10 ? (lstirling(n+1) - log(n+1))
           : (0.5*log(2*Num.pi*n) + n*log(n/Num.e + 1/(12*Num.e*n)))
}
 
func P_approx(n, k) {
    Logarithm((lstirling(n) - lstirling(n -k)))
}
 
func C_approx(n, k) {
    Logarithm((lstirling(n) - lstirling(n -k) - lstirling(k)))
}
 
say "=> Exact results:"
for n (1..12) {
    var p = n//3
    say "P(#{n}, #{p}) = #{P(n, p)}"
}

for n (10..60 `by` 10) {
    var p = n//3
    say "C(#{n}, #{p}) = #{C(n, p)}"
}

say '';
say "=> Floating point approximations:"
for n ([5, 50, 500, 1000, 5000, 15000]) {
    var p = n//3
    say "P(#{n}, #{p}) = #{P_approx(n, p)}"
}

for n (100..1000 `by` 100) {
    var p = n//3
    say "C(#{n}, #{p}) = #{C_approx(n, p)}"
}
