On Sat, Sep 22, 2007 at 10:35:16AM +0200, Gael Varoquaux wrote:
> I would go for the "generate_fourplets" solution if I had a way to
> calculate the binomial coefficient without overflows.
Sorry, premature optimisation is the root of all evil, but turning ones
brain on early is good.
"""
##############################################################################
# Some routines for calculation of binomial coefficients
def gcd(m,n):
while n:
m,n=n,m%n
return m
def binom_(n,k):
if k==0:
return 1
else:
g = gcd(n,k)
return binomial(n-1, k-1)/(k/g)*(n/g)
def binomial(n,k):
if k > n/2: # Limit recursion depth
k=n-k
return binom_(n,k)
"""
This is surprisingly fast (surprising for me, at least).
Using this and the C code I have, I can generate the quadruplets of 200
integers quite quickly:
In [5]: %timeit b = [1 for i in generate_quadruplets(200)]
10 loops, best of 3: 1.61 s per loop
With generate_quadruplets given by:
"""
##############################################################################
def generate_quadruplets(size):
""" Returns an iterator on tables listing all the possible unique
combinations of four integers below size. """
C_code = """
int index = 0;
for (int j=0; j<i+1; j++) {
for (int k=0; k<j+1; k++) {
for (int l=0; l<k+1; l++) {
quadruplets(index, 0) = i;
quadruplets(index, 1) = j;
quadruplets(index, 2) = k;
quadruplets(index, 3) = l;
index++ ;
}
}
}
"""
for i in xrange(size):
multiset_coef = binomial(i+3, 3)
quadruplets = empty((multiset_coef, 4), dtype=uint32)
inline(C_code, ['quadruplets', 'i'],
type_converters=converters.blitz)
yield quadruplets
"""
This fits my needs.
Maybe I should add this on the cookbook. It is a bit specific, but it
shows a serie of interesting problems, I think.
I also think the binomial and gcd should go in scipy (I could not find
them, but maybe they are already there). Maybe in a new module, as I
don't really see where this fits in.
Gaël
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
