combination function in python

[Jussi Piitulainen]

Steven D'Aprano writes:
> bearophileHUGS wrote:
...
>>     return factorial(n) // (factorial(k) * factorial(n-k))
> 
> That's a naive and slow implementation. For even quite small values
> of n and k, you end up generating some seriously big long ints, and
> then have to (slowly!) divide them.

A better _definition_ of the binomial coefficient with upper index r
and lower index k is (r * (r - 1) * ...) / (k * (k - 1) * ...) with k
factors in both products. These are called falling factorial powers by
Graham, Knuth and Patashnik. Their notation is to write n^k and k^k
but with the exponent underlined; the latter is just k!, when k > 0. A
straightforward implementation below.

> A better implementation would be something like this:
> 
> def binomial(n, k):
>     if not 0 <= k <= n:
>         return 0
>     if k == 0 or k == n:
>         return 1
>     # calculate n!/k! as one product, avoiding factors that 
>     # just get canceled
>     P = k+1
>     for i in xrange(k+2, n+1):
>         P *= i
>     # if you are paranoid:
>     # C, rem = divmod(P, factorial(n-k))
>     # assert rem == 0
>     # return C
>     return P//factorial(n-k)
> 
> There's probably even a really clever way to avoid that final
> division, but I suspect that would cost more in time and memory than
> it would save.

Here's one non-clever one for integers n, k that uses n^k / k^k
(falling powers) with the smaller of k and n - k as lower index:

def choose(n, k):
   if 0 <= k <= n:
       ntok = 1
       ktok = 1
       for t in xrange(1, min(k, n - k) + 1):
          ntok *= n
          ktok *= t
          n -= 1
       return ntok // ktok
   else:
       return 0