Add a function for computing binomial coefficients to the math module #79612
Comments
|
A recuring pain point, for me and for many others who use Python for mathematical computations, is that the standard library does not provide a function for computing binomial coefficients. I would like to suggest adding a function, in the math module, with the following signature: binomial(n: Integral, k: Integral) -> Integral A simple native Python implementation would be: from functools import reduce
from math import factorial
from numbers import Integral
def binomial(n: Integral, k: Integral) -> Integral:
if k < 0 or n < 0:
raise ValueError("math.binomial takes non-negative parameters.")
k = min(k, n-k)
num, den = 1, 1
for i in range(k):
num = num * (n - i)
den = den * (i + 1)As far as I can tell, all of the math module is implemented in C, so this should be done in C too, but the implemented behaviour should be equivalent. I will submit a Github pull request once I have a ready-to-review patch. Not starting a PEP, per PEP-1:
|
KellerFuchs
mannequin
added
3.8
|
You import reduce but never use it :-) +1 for this, I certainly miss it too. |
|
Yes, that was a copypasta mistake (and I also import factorial needlessly) as the file I prototyped it in had some other code for testing my proposed implementation. :) |
|
+1 I have wanted this a number of times. FWIW, most recently I wrote it like this: def comb(n, r):
'Equivalent to factorial(n) // (factorial(r) * factorial(n-r))'
c = 1
r = min(r, n-r)
for i in range(1, r+1):
c = c * (n - r + i) // i
return cI'm not sure is this is better than a single divide, but it kept the intermediate values from exploding in size, taking advantage of cancellation at each step. Also, I'm not sure what the predominant choice for variable names should be, "n things taken r at a time" or "n things taken k at time". Also, it's worth considering whether the function should be called "binomial", "choose", "combinations", or "comb". The word "binomial" seems too application specific but would make sense if we ever introduced a "multinomial" counterpart. The word "choose" is how we usually pronounce it. The word "combinations" fits nicely with the related counting functions, "combinations, permutations, and factorial". The word "comb" is short, works well with "perm" and "fact", and nicely differentiates itself as the integer counterparts of the combinatoric functions in the itertools module. Wolfram uses both choose and Binomial[n,m] |
|
For some ranges of inputs, it may make sense to use the apparently naive implementation Here are some timings (Python 3.7.1, macOS 10.13), using Raymond's In [5]: %timeit comb(1000, 600) In [6]: %timeit factorial(1000) // factorial(600) // factorial(400) In [7]: %timeit factorial(1000) // (factorial(600) * factorial(400)) But that's just one set of inputs, on one system; your results may vary. |
|
There's also the question of what inputs should be considered valid: Note that this needs fixing with both of the code snippets shown so far: they both return 1 for k > n. |
|
@Raymond Hettinger:
Here goes the bike shed! :) I'd rather not call it comb, because it collides with a completely-unrelated English word, and it's not obvious what it stands for unless one already knows.
That's a good point; math.permutations doesn't exist, but itertools.permutations does, so I would expect (by analogy) a “combinations” functions to produce all possible k-sized subsets (rather than counting them), and that's exactly what itertools.combinations does. On the other hand, combinations and permutations are names that make it perhaps more obvious what is being counted, so perhaps math.{combinations,permutations} should be aliases for math.{binomial,factorial} ? Is the name collision with itertools a problem? TL;DR: “binomial” is more consistent with the current naming in math and itertools, but perhaps it makes sense to introduce math.{combination,permutations} as aliases? (Note that I used math.binomial as the name in the PR so far, but that's only because I needed a name, not because I consider the discussion ended.) @mark Dickinson:
Yes, I avoided pushing specifically for a given algorithm (rather than getting upfront agreement on the functional behaviour) because the performance characteristics will likely be quite different once implemented in C in the math module. (Unless I'm mistaken and there is a way to add pure-Python functions to the math module?)
From a mathematical standpoint, (n choose k) is defined for all non-negative k, n, with (n chooze k) = 0 when k>n or k=0. It's necessary behaviour for the usual equations to hold (like (n + 1 choose k + 1) = (n choose k) + (n choose k + 1)). Negative k and n, on the other hand, should clearly be a ValueError, and so should non-integers inputs; this is consistent with factorial's behaviour. I started a pull request and (for now) only added tests which document that (proposed) behaviour, so we can more easily discuss whether that's what we want.
Yes :) |
|
I think that it is better to add this function in a new module imath, which could contain other integer functions: factorial, gcs, as_integer_ration, isqrt, isprime, primes. See https://mail.python.org/pipermail/python-ideas/2018-July/051917.html |
|
Mathematically, |
It's not so clear cut. You can find different definitions out there. Knuth et. al., for example, in the book "Concrete Mathematics", extend the definition not just to negative k, but to negative n as well. Mathematicians aren't very good at agreeing on things. :-) But that doesn't really matter: what we need to decide is what behaviour is useful for the users of the function. |
You don't mean the "k=0" part of that, right? |
|
One more decision that needs to be made: should the new function accept integer-valued floats? Or should any I'd personally prefer that floats not be accepted; I think this was a misfeature of |
|
On Fri, Dec 07, 2018 at 12:04:44AM +0000, Raymond Hettinger wrote:
I've done a quick survey of some of the most common/popular scientific TI Nspire all call this nCr, and nPr for the permutation version. This matches For what its worth, the colour I prefer for this bikeshed are "comb" and |
I think the key word there is *extend*. To the degree that any However, I think that it is too easy to get the order of n and r (n and |
|
On Fri, Dec 07, 2018 at 01:37:36PM +0000, Mark Dickinson wrote:
Agreed. We can always add support for floats later, but its hard to We ought to require n, r (or n, k) to be non-negative ints with 0 <= r |
|
@serhiy Storchaka:
imath is a fine idea, and you already started a discussion in python-ideas@, but it's a much bigger undertaking than just adding this one function, and you can move it there once imath happens. |
|
@mark Dickinson:
Indeed not; the tests in the PR actually assert binomial(n, n) == binomial(n, 0) == 1. |
OK; I only went with that because I assumed there were Good Reasons© that factorial did it, but if rejecting integral floats isn't going to be a controversial move I'm also in favor of it. Updating the PR momentarily to check that binomial rejects floats. |
|
@steven D'Aprano:
It's also not universal; in my experience, most calculators are localized for a given market, and may use different notations (in particular, the notation for combinations/binomial numbers changes across countries). |
brettcannon
changed the title
|
Brett, what was the purpose of the title change?
|
+1 These would be my preferences as well :-) |
|
My two cents:
Not concerned about speed. It's possible to do this with no divisions involving integers bigger than (*) Sketch: make lists of the k numerators and k-1 denominators (skip 1). For each prime P <= k, a modulus operation can determine the index of the first multiple of P in each list. For that, and for each P'th later list element, divide out the multiples of P, adding 1 to a running count for numerators, subtracting 1 for denominators, and reducing each list element by the Ps divided out. Then if the final P count isn't 0 (it will always be >= 0), append pow(P, count) to a list of factors. pow() is efficient. After that, all the denominators will be reduced to 1, so can be ignored thereafter. It just remains to multiply all the reduced numerators and prime-power factors. Catenate them all in a single list, heapify it (min heap), and so long as at least 2 factors remain pop the two smallest and push their product. This attempts to balance bit lengths of incoming factors, allowing close-to-best-case-for-it Karatsuba multiplication to kick in. But that's nuts ;-) To get even nutsier, special-case P=2 to use shifts instead, skip adding pow(2, count) to the list of factors, and just shift left by the count at the end. That said, even the "dumbest possible loop" may benefit in C by shifting out all trailing zeroes, multiplying/dividing only odd integers, and shifting left at the end. |
|
[Tim]
+1 to all of this. |
|
Just for fun, here's a gonzo implementation (without arg-checking) using ideas from the sketch. All factors of 2 are shifted out first, and all divisions are done before any multiplies. For large arguments, this can run much faster than a dumb loop. For example, combp(10**100, 400) takes about a quarter the time of a dumb-loop divide-each-time-thru implementation. # Return number of trailing zeroes in `n`.
def tzc(n):
result = 0
if n:
mask = 1
while n & mask == 0:
result += 1
mask <<= 1
return result
# Return exponent of prime `p` in prime factorization of
# factorial(k).
def facexp(k, p):
result = 0
k //= p
while k:
result += k
k //= p
return result
def combp(n, k):
from heapq import heappop, heapify, heapreplace
if n-k < k:
k = n-k
if k == 0:
return 1
if k == 1:
return n
firstnum = n - k + 1
nums = list(range(firstnum, n+1))
assert len(nums) == k
# Shift all factors of 2 out of numerators.
shift2 = 0
for i in range(firstnum & 1, k, 2):
val = nums[i]
c = tzc(val)
assert c
nums[i] = val >> c
shift2 += c
shift2 -= facexp(k, 2) # cancel all 2's in factorial(k)
assert shift2 >= 0
# Any prime generator is fine here. `k` can't be
# huge, and we only want the primes through `k`.
pgen = psieve()
p = next(pgen)
assert p == 2
for p in pgen:
if p > k:
break
pcount = facexp(k, p)
assert pcount
# Divide that many p's out of numerators.
i = firstnum % p
if i:
i = p - i
for i in range(i, k, p):
val, r = divmod(nums[i], p)
assert r == 0
pcount -= 1
while pcount:
val2, r = divmod(val, p)
if r:
break
else:
val = val2
pcount -= 1
nums[i] = val
if pcount == 0:
break
assert pcount == 0
heapify(nums)
while len(nums) > 1:
a = heappop(nums)
heapreplace(nums, a * nums[0])
return nums[0] << shift2I'm NOT suggesting to adopt this. Just for history in the unlikely case there's worldwide demand for faster |
|
Can I work on C implementation if no-one else is doing it right now? |
Sounds fine to me. You might want to coordinate with @KellerFuchs to see what the status of their PR is; maybe the two of you can collaborate? @KellerFuchs: are you still planning to work on this? |
|
Kellar and Yash, my suggestion is to separate the work into two phases. Start with an initial patch that implements this simplest possible implementation, accompanied by clean documentation and thorough testing. Once everyone has agreed on the API (i.e. calling it "comb()", how to handle various input datatypes, and handling of corner cases), and the patch is applied, only then turn to a second pass for optimizations (special casing various types, minimizing how big intermediate values can get by doing early cancellation, exploiting even/odd patterns etc). |
|
Why do we want __index__ support? This seems like an unnecessary extension without relevant use cases. |
|
NumPy integer types are not subclasses of int. |
|
I'd expect any math module function accepting integers to be sufficiently duck-typed to accept integer-like things (i.e., objects that implement |
|
Mark, please see: https://twitter.com/daveinstpaul/status/1134919179361034240 What are your thoughts? |
|
One other idea: it may be worth putting in an alias for "binomial" to improve findability and readability in contexts where a person really does what binomial coefficients and is not otherwise thinking about counting contexts. |
Sigh. I don't object to extending to I'd say leave it as-is for 3.8, see what the reaction is, and maybe relax constraints in 3.9 if that seems appropriate. But if a majority of others really want to make the change now, that's okay with me. |
I concur. That said, the referenced tweet was a reaction :-) FWIW, itertools.combinations(seq, r) returns 0 values when r > len(seq). Tim, what do you think? |
|
I'm particularly sceptical about real-world use-cases for k < 0. Maybe we should consider removing just the k > n requirement. |
|
Strongly prefer requiring 0 <= k <= n at first. This is a programming language: it will be applied to real problems, not to abstract proofs where some slop can be helpful in reducing the number of cases that need to be considered. The Twitter fellow is certainly right that "0" is the clear answer to "how many 5-element subsets does have a 4-element set have?", but in the context of real problems it's more likely (to my eyes) that a programmer asking that question of Raymond, I'm not clear on what you mean by "alias". You mean supplying the same function under more than one name? I'd be -1 on that (where would it end? the answer to every future name bikeshedding issue then would be "all of the above"). But +1 on, e.g., adding a doc index entry for "binomial" pointing to "comb". |
|
I think that binomial(n, k) should return 0 when k > n or k < 0. This is a practical consideration. I'm concerned about evaluating sums involving binomial coefficients. Mathematicians are often rather loose about specifying the upper and lower bounds of summation, because the unwanted terms are zero anyway. These formulas should not result in value errors when they are implemented directly. To give a simplistic example, the sum of the first n positive integers is binomial(n+1, 2), and the formula should still work if n is zero. |
|
@david Ratcliffe
I agree that's a convincing use-case. Can you think of any practical uses for allowing negative |
Radcliffe! Apologies for the misspelling. It's still early in this timezone and I haven't had any coffee yet. |
|
I'm not convinced, although I agree relaxing k <= n is less damaging than relaxing k >= 0. Python isn't aimed at mathematicians (although some 3rd-party packages certainly are, and they're free to define things however they like). We have to trade off convenience for experts in edge cases against the likelihood that an ordinary user is making a mistake. For example, that's why, despite that Python supports complex numbers, math.sqrt(-1) raises ValueError. For _most_ users, trying that was probably an error in their logic that led to the argument being negative. Experts can use cmath.sqrt() instead. Ordinary users think comb(n, k) is the value of n!//(k!*(n-k)!), and as far as they're concerned factorial is only defined for non-negative integers. 0 <= k <= n follows from that. (Indeed, the current docs for And ordinary users think of the sum of the first n integers as "number of terms times the average", or memorize directly that the answer is n*(n+1)/2. That works fine for n=0. Only an expert thinks of it as So I would still rather enforce 0 <= k <= n at the start, and relax it later only if there's significant demand. In going on three decades of using Python and having written my own |
|
Here are a few other data points. There is no consensus but either returning 0 or erroring out seem reasonable. MS Excel's COMBIN: If number < 0, number_chosen < 0, or number < number_chosen, COMBIN returns the #NUM! error value. Scipy.misc.comb: If k > N, N < 0, or k < 0, then a 0 is returned. Matlib.nchoose: Error: K must an integer between 0 and N. Maple.numbcomb: Unclear from the docs what this does Wolfram alpha: Returns 0 |
|
I'm going to repeat part of an earlier comment :-) """ Why can't we take universal agreement for an initial answer? ;-) Looking at what other packages do doesn't really interest me, unless they differ on the cases everyone agreed would be useful - they have different audiences than Python has. Guess what Wolfram returns for (-4) choose (-5) Go ahead. That's right: -4. Sheesh - I don't care who their audience is, any code calling |
|
I understand that pure mathematics is not the primary use case. But I would On Sun, Jun 2, 2019 at 8:59 PM Raymond Hettinger <report@bugs.python.org>
|
|
I'm not fatally opposed to relaxing k <= n. David makes some good points about it, and - as Raymond already noted - "0" is consistent with the behavior of itertools.combinations(). The docs would need to change, though, because the factorial "definition" would no longer make sense. Defining it in terms of the falling factorial would, but that's too obscure for Python's audience. Probably a combinatorial definition would be best: `comb(n, k) is the number of k-element subsets of an n-element set". Then it's clear that n and k are cardinals (integers >= 0), and that the result is 0 when k > n. |
|
For what its worth, there are concrete, practical applications for which in turn has applications in physics, chemistry and other sciences: https://en.wikipedia.org/wiki/Fractional_calculus#Applications Take that however you see fit :-) My own preference is for a comb() function that returns 0 for out of And yes, Mathematica does accept fractional and complex arguments to https://www.wolframalpha.com/input/?i=choose%28sqr%28-3.5%29,+sqrt%28-4%2Bi%29%29 |
|
Given that (a) we're right up against feature freeze, and (b) the discussions about API already happened in plenty of good time some months ago, leading to agreement on the initial API, I think we should avoid any last-minute changes and stick with what we have for 3.8. Then we can have a more relaxed discussion about what, if anything, should be changed for 3.9. |
|
I don't have a preference one way or the other (both ways have plausible arguments, error checking vs utility beyond simple combinatorics, and both ways have ample precedents). That said, if we're going to ultimately relax the constraints, it would be better to do it now before the API is released. Even if we were already in feature freeze, it would be okay to change it (the purpose of a beta is to elicit end-user feedback which is what David has just given). |
|
Python needs a benevolent dictator ;-) I'd be happy for Mark, Raymond, or me to play that role here. If it were me, I'd drop the k <= n requirement: both arguments must be non-negative integers. Because a combinatorial (not factorial-, and certainly not gamma-, based) definition of |
+1 from me. |
|
Fine by me. The same change should be applied to |
Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.
Show more details
GitHub fields:
bugs.python.org fields:
The text was updated successfully, but these errors were encountered: