459,272 Members | 1,442 Online Need help? Post your question and get tips & solutions from a community of 459,272 IT Pros & Developers. It's quick & easy.

# combination function in python

 P: n/a how to use the combination function in python ? For example 9 choose 2 (written as 9C2) = 9!/7!*2!=36 Please help, I couldnt find the function through help. Apr 15 '07 #1
17 Replies

 P: n/a DanielJohnson >import gmpygmpy.comb(9,2) mpz(36) However, there's no equivalent function built into Python, if that's what you're asking (gmpy's a third-party extension). It's of course easy to write one for yourself, if you want the functionality but don't want to download gmpy and don't care for gmpy's superb speed. Alex Apr 15 '07 #2

 P: n/a DanielJohnson: Please help, I couldnt find the function through help. You can't find it because it's not there: def factorial(n): """factorial(n): return the factorial of the integer n. factorial(0) = 1 factorial(n) with n<0 is -factorial(abs(n)) """ result = 1 for i in xrange(1, abs(n)+1): result *= i if n >= 0: return result else: return -result def binomial(n, k): """binomial(n, k): return the binomial coefficient (n k).""" assert n>0 and isinstance(n, (int, long)) and isinstance(k, (int, long)) if k < 0 or k n: return 0 if k == 0 or k == n: return 1 return factorial(n) // (factorial(k) * factorial(n-k)) Bye, bearophile Apr 15 '07 #3

 P: n/a On Sun, 15 Apr 2007 02:38:31 -0700, bearophileHUGS wrote: DanielJohnson: >Please help, I couldnt find the function through help. You can't find it because it's not there: def factorial(n): """factorial(n): return the factorial of the integer n. factorial(0) = 1 factorial(n) with n<0 is -factorial(abs(n)) """ result = 1 for i in xrange(1, abs(n)+1): result *= i if n >= 0: return result else: return -result def binomial(n, k): """binomial(n, k): return the binomial coefficient (n k).""" assert n>0 and isinstance(n, (int, long)) and isinstance(k, (int, long)) if k < 0 or k n: return 0 if k == 0 or k == n: return 1 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 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. -- Steven. Apr 15 '07 #4

 P: n/a Steven D'Aprano: 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 implementation would be something like this: You are right, thank you for the improvement (the advantage of the older implementation is that it's naive, so it's a bit more probably correct compared to more complex code I may write. For Python code I often tend to write a naive version first, create many testcases, slowly fixing all the corner cases (like factorial(-5)), and only later find a faster/better implementation if I have some time to do it or if I need it. If I need to do lot of binomials the gmpy by Alex helps). Bye, bearophile Apr 15 '07 #5

 P: n/a 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 Apr 15 '07 #6

 P: n/a On Apr 15, 8:37 am, Jussi Piitulainen wrote: 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 It might be even better to do the divisions as you go, rather than leaving them all to the end. That way the intermediate results stay smaller. So (leaving out the bounds checking) one just does: def choose(n, k): ntok = 1 for t in xrange(min(k, n-k)): ntok = ntok*(n-t)//(t+1) return ntok Mark Apr 15 '07 #7

 P: n/a Jussi Piitulainen wrote: >There's probably even a really clever way to avoid that finaldivision, but I suspect that would cost more in time and memory thanit would save. We're getting closer and closer to something I already posted a few times here. This implementation was unfortunate because I consistently used an uncommon name for it so people couldn't easily find it (mea culpa), and maybe also because it uses the despised reduce builtin. def noverk(n,k): return reduce(lambda a,b: a*(n-b)/(b+1),xrange(k),1) BTW I hereby give up the strange name for this and request a better name that everybody will use so there will be no confusion anymore. Maybe comb(n,k) ? 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 Ha, my function uses smaller subproducts :-) A. Apr 15 '07 #8

 P: n/a Mark Dickinson writes: Jussi Piitulainen wrote: 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 It might be even better to do the divisions as you go, rather than leaving them all to the end. That way the intermediate results stay smaller. So (leaving out the bounds checking) one just does: def choose(n, k): ntok = 1 for t in xrange(min(k, n-k)): ntok = ntok*(n-t)//(t+1) return ntok Yes, that's what I once saw done. Thanks. Its correctness is more subtle, so I prefer to add the parentheses below to emphasise that the product has to be computed before the division. I also renamed the variable to p: it's no longer n^k (falling). And I put the bounds back in. def choose(n, k): if 0 <= k <= n: p = 1 for t in xrange(min(k, n - k)): p = (p * (n - t)) // (t + 1) return p else: return 0 Apr 15 '07 #9

 P: n/a Mark Dickinson: def choose(n, k): ntok = 1 for t in xrange(min(k, n-k)): ntok = ntok*(n-t)//(t+1) return ntok With few tests, it seems this is faster than the version by Jussi only with quite big n,k. (Another test to be added, is the assert n>0 of my original version (testing that n,k seem useful too)). Bye, bearophile Apr 15 '07 #10

 P: n/a On Apr 15, 11:34ï¿½am, Anton Vredegoor wrote: Jussi Piitulainen wrote: 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. We're getting closer and closer to something I already posted a few times here. This implementation was unfortunate because I consistently used an uncommon name for it so people couldn't easily find it But then, who's looking for it? (mea culpa), and maybe also because it uses the despised reduce builtin. def noverk(n,k): * * *return reduce(lambda a,b: a*(n-b)/(b+1),xrange(k),1) BTW I hereby give up the strange name for this and request a better name that everybody will use so there will be no confusion anymore. Maybe comb(n,k) ? No, that name's already used by gmpy. And a single function doesn't replace all of gmpy's other functionality, many of which are needed in the same applications where comb() is used, so there's really no need for your function. Your time is better spent applying the tools provided by gmpy rather than trying to re-invent the wheel. > 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 Ha, my function uses smaller subproducts :-) A. Apr 15 '07 #11

 P: n/a On Apr 15, 3:33 pm, bearophileH...@lycos.com wrote: With few tests, it seems this is faster than the version by Jussi only with quite big n,k. True---and for large n and k it's difficult to imagine any real-world applications that wouldn't be better served by using the lngamma function instead. That is, the natural log of n choose k can be computed much more quickly as lngamma(n+1) - lngamma(k+1) - lngamma(n-k+1) I assume that lngamma is implemented somewhere in either numpy or scipy, but right now I can't find it... Mark Apr 15 '07 #12

 P: n/a me********@aol.com wrote: >We're getting closer and closer to something I already posted a fewtimes here. This implementation was unfortunate because I consistentlyused an uncommon name for it so people couldn't easily find it But then, who's looking for it? The OP was trying to find it in the docs, assuming it was some kind of builtin function. I can see the logic of that and I can also see the logic of including some other smallish functions like for example fac. Failing that, the next recourse would be the Internet, or a Usenet search but that would imply well named Usenet posts and function names. >(mea culpa), and maybe also because it uses the despised reduce builtin.def noverk(n,k):ï¿½ ï¿½ ï¿½return reduce(lambda a,b: a*(n-b)/(b+1),xrange(k),1) This is a rather concise function which has the added advantage that it returns 0 when k>n. >BTW I hereby give up the strange name for this and request a better namethat everybody will use so there will be no confusion anymore. Maybecomb(n,k) ? No, that name's already used by gmpy. And a single function doesn't replace all of gmpy's other functionality, many of which are needed in the same applications where comb() is used, so there's really no need for your function. Actually, by proposing comb I am *honoring* gmpy and I'm also complying with it. In Python we use namespaces to differentiate between such things. You might want to read the documentation about it some time. Your time is better spent applying the tools provided by gmpy rather than trying to re-invent the wheel. Please let me be the judge of how to spend my time. In this case it seems rather unjustified to introduce dependencies on external modules when only a few functions are needed. Since I'm also interested in the functions themselves -in this case- I'd rather have a few lines of code in my source than importing an optimized code library. There *are* situations where such things are completely justified, but I don't think this is one of them. You could take it up with the gmpy author and induce him to get gmpy included in the standard distro if you are so inclined. A. Apr 15 '07 #13

 P: n/a On Sun, 15 Apr 2007 13:58:38 -0700, Mark Dickinson wrote: ... for large n and k it's difficult to imagine any real-world applications that wouldn't be better served by using the lngamma function instead. That is, the natural log of n choose k can be computed much more quickly as lngamma(n+1) - lngamma(k+1) - lngamma(n-k+1) Yes, but if you actually want the number of combinations, not the log of the number of combinations, you then have to raise e to that power to get the answer you want -- and that takes time. Is it still faster? How about round-off error due to floating point? Should you be truncating or rounding? What about overflow error for large enough n, k? What exactly are we doing with a five hundred digit long integer anyway? Exact (long) integer maths still is useful. I assume that lngamma is implemented somewhere in either numpy or scipy, but right now I can't find it... You know what they say about what happens when you ASS_U_ME :-) -- Steven. Apr 15 '07 #14 