Proving optimized function is still correct

[Fernando Pérez]

Huaiyu Zhu wrote:

> On Tue, 13 Nov 2001 16:55:48 GMT, Terry Reedy <tjreedy at home.com> wrote:
>>A recursive implementation like
 [snip]
> Since m and n are symmetrical, one could also compare them and choose
> between space optimization vs recursion depth optimization.

The problem with all these forms is that they grow in execution time with the 
parameters (though they are exact, if they use longs). Here's a constant-time 
version, which is however limited to the range where the answer fits in a 
float (I fixed some things from yesterday's version). For low values of m 
this is expensive, but once m gets big it's immensely faster.

If one really needs speed, the ideal solution is probably a combination of a 
good recursive solution for low m with this version for high m, and caching 
of results for low (m,n) values if those are commonly needed.

from Numeric import *
#-----------------------------------------------------------------------------
def gammln(x):
    """Return the natural log of gamma(x)"""

    cof = array([76.18009172947146,-86.50532032941677,
                 24.01409824083091,-1.231739572450155,
                 0.1208650973866179e-2,-0.5395239384953e-5],Float)

    y = arange(x+1,x+7,typecode = Float)
    tmp = x+5.5
    tmp -= (x+0.5)*log(tmp)
    ser = 1.000000000190015 + sum(cof/y)
    return -tmp+log(2.5066282746310005*ser/x)

#-----------------------------------------------------------------------------
def bico(n,k):
    """Return binomial coefficient (n k)"""
    return floor(0.5+exp(gammln(n+1.0)-gammln(k+1.0)-gammln(n-k+1.0)))

#-----------------------------------------------------------------------------
def combi_dist(m,n):
    """Return the number of combinations of m mutually distinguishable
    elements in which each element may occur 0,1,2,...,n times in any
    combination.

    Constant time, uses floats."""

    # this is just bico(m+n-1,n). Write it out explicitly for speed
    nn = m+n-1
    kk = n
    return floor(0.5 + 
           exp(gammln(nn+1.0)-gammln(kk+1.0)-gammln(nn-kk+1.0)))

This bypasses the issue of explicitly building the factorials by computing 
them through a gamma function (the implementation is from Numerical Recipes, 
so not necessarily the best).

Cheers,

f