Issue35431

Created on **2018-12-06 21:18** by **kellerfuchs**, last changed **2019-01-04 18:43** by **FR4NKESTI3N**.

Pull Requests | |||
---|---|---|---|

URL | Status | Linked | Edit |

PR 11018 | open | kellerfuchs, 2018-12-07 11:20 | |

PR 11414 | open | FR4NKESTI3N, 2019-01-02 19:31 |

Messages (34) | |||
---|---|---|---|

msg331251 - (view) | Author: (kellerfuchs) | Date: 2018-12-06 21:18 | |

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) return num//den 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]: > Small enhancements or patches often don't need a PEP and can be injected into the Python development workflow with a patch submission to the Python issue tracker. [PEP 1]: https://www.python.org/dev/peps/pep-0001/#id36 |
|||

msg331255 - (view) | Author: Steven D'Aprano (steven.daprano) * | Date: 2018-12-06 22:31 | |

You import reduce but never use it :-) +1 for this, I certainly miss it too. |
|||

msg331256 - (view) | Author: (kellerfuchs) | Date: 2018-12-06 23:22 | |

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. :) |
|||

msg331257 - (view) | Author: Raymond Hettinger (rhettinger) * | Date: 2018-12-07 00:04 | |

+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 c I'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] SciPy uses comb(n, k). Maple uses both numbcomb(n,m) and binomial(n,m). TeX uses {n \choose k} |
|||

msg331280 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2018-12-07 08:42 | |

For some ranges of inputs, it may make sense to use the apparently naive implementation `factorial(n) // factorial(k) // factorial(n - k)`. The current factorial implementation is significantly optimised, and using it directly may be faster than using an iterative solution. Here are some timings (Python 3.7.1, macOS 10.13), using Raymond's `comb` function from msg331257: In [5]: %timeit comb(1000, 600) 186 µs ± 442 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) In [6]: %timeit factorial(1000) // factorial(600) // factorial(400) 97.8 µs ± 256 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) In [7]: %timeit factorial(1000) // (factorial(600) * factorial(400)) 91.1 µs ± 789 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) But that's just one set of inputs, on one system; your results may vary. |
|||

msg331281 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2018-12-07 09:03 | |

There's also the question of what inputs should be considered valid: `binomial(n, k)` for `k > n` should either return 0 or be a `ValueError`, but which? Same question for `k < 0`. There's a symmetry argument for allowing `k < 0` if you allow `k > n`, but I can't think of pragmatic reasons to allow `k < 0`, while allowing `k > n` _does_ seem potentially useful. Note that this needs fixing with both of the code snippets shown so far: they both return 1 for k > n. |
|||

msg331293 - (view) | Author: (kellerfuchs) | Date: 2018-12-07 11:25 | |

@Raymond Hettinger: > it's worth considering whether the function should be called "binomial", "choose", "combinations", or "comb" Here goes the bike shed! :) Kidding, this is a pretty important ergonomics/discoverability concern, and I hadn't given it much thought yet. 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. > The word "combinations" fits nicely with the related counting functions, "combinations, permutations, and factorial". 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: > The current factorial implementation is significantly optimised, and using it directly may be faster than using an iterative solution. 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?) > `binomial(n, k)` for `k > n` should either return 0 or be a `ValueError`, but which? 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)). As such, I'd argue that returning 0 is both more likely to be the thing the user wants (as in, it's necessary behaviour for combinatorics) and easier to reason about. 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. > Note that this needs fixing with both of the code snippets shown so far: they both return 1 for k > n. Yes :) I noticed last night, as I wrote Hypothesis tests for the snippet, but didn't think it was super important to send an updated version. |
|||

msg331295 - (view) | Author: Serhiy Storchaka (serhiy.storchaka) * | Date: 2018-12-07 11:30 | |

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 |
|||

msg331296 - (view) | Author: Serhiy Storchaka (serhiy.storchaka) * | Date: 2018-12-07 11:32 | |

Mathematically, `binomial(n, k)` for `k > n` is defined as 0. |
|||

msg331308 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2018-12-07 13:08 | |

> Mathematically, `binomial(n, k)` for `k > n` is defined as 0. 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. |
|||

msg331309 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2018-12-07 13:10 | |

@kellerfuchs: > 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. You don't mean the "k=0" part of that, right? |
|||

msg331312 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2018-12-07 13:37 | |

One more decision that needs to be made: should the new function accept integer-valued floats? Or should any `float` input give a TypeError. I'd personally prefer that floats not be accepted; I think this was a misfeature of `factorial` that we shouldn't compound. |
|||

msg331318 - (view) | Author: Steven D'Aprano (steven.daprano) * | Date: 2018-12-07 14:44 | |

On Fri, Dec 07, 2018 at 12:04:44AM +0000, Raymond Hettinger wrote: > 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". I've done a quick survey of some of the most common/popular scientific calculators: TI Nspire TI-84 Plus Casio Classpad Casio FX-82AU Plus II all call this nCr, and nPr for the permutation version. This matches the notation taught in secondary school maths classes in Australia. That's common and familiar notation for secondary school students, but personally I'm not super-keen on it. For what its worth, the colour I prefer for this bikeshed are "comb" and "perm", which are the names used by the HP 48GX calculator. Second choice would be to spell the names out in full, "combinations" and "permutations". |
|||

msg331323 - (view) | Author: Steven D'Aprano (steven.daprano) * | Date: 2018-12-07 14:52 | |

> > Mathematically, `binomial(n, k)` for `k > n` is defined as 0. > > 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. :-) I think the key word there is *extend*. To the degree that any mathematical definition is "obvious", the obvious definition for number of combinations ("n choose r") is going to be 1 for r == 0 and 0 for r > n. However, I think that it is too easy to get the order of n and r (n and k) mixed up, and write combinations(5, 10) when you wanted to choose 5 from 10. I know I make that mistake on my calculator *all the time*, and so do my students, even with the nPr notation. So I recommend we raise ValueError for r > n. |
|||

msg331325 - (view) | Author: Steven D'Aprano (steven.daprano) * | Date: 2018-12-07 14:56 | |

On Fri, Dec 07, 2018 at 01:37:36PM +0000, Mark Dickinson wrote: > I'd personally prefer that floats not be accepted; Agreed. We can always add support for floats later, but its hard to remove it if it turns out to be problematic. We ought to require n, r (or n, k) to be non-negative ints with 0 <= r <= n. Extending this to negative ints or floats is probably YAGNI, but if somebody does need it, they can request an enhancement in the future. |
|||

msg331335 - (view) | Author: (kellerfuchs) | Date: 2018-12-07 16:47 | |

@Serhiy Storchaka: > I think that it is better to add this function in a new module imath, which could contain other integer functions 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. As such, I think it's out-of-scope in this issue. |
|||

msg331336 - (view) | Author: (kellerfuchs) | Date: 2018-12-07 16:49 | |

@Mark Dickinson: > You don't mean the "k=0" part of that, right? Indeed not; the tests in the PR actually assert binomial(n, n) == binomial(n, 0) == 1. |
|||

msg331337 - (view) | Author: (kellerfuchs) | Date: 2018-12-07 16:52 | |

> I'd personally prefer that floats not be accepted; I think this was a misfeature of `factorial` that we shouldn't compound. 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. |
|||

msg331339 - (view) | Author: (kellerfuchs) | Date: 2018-12-07 17:01 | |

@Steven D'Aprano: > all call this nCr, and nPr for the permutation version. This matches > the notation taught in secondary school maths classes in Australia. > That's common and familiar notation for secondary school students, but > personally I'm not super-keen on it. 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). |
|||

msg331369 - (view) | Author: Steven D'Aprano (steven.daprano) * | Date: 2018-12-07 23:53 | |

Brett, what was the purpose of the title change? > title: The math module should provide a function for computing > binomial coefficients -> Add a function for computing binomial > coefficients to the math module |
|||

msg331402 - (view) | Author: Raymond Hettinger (rhettinger) * | Date: 2018-12-09 01:11 | |

> For what its worth, the colour I prefer for this bikeshed are > "comb" and "perm", which are the names used by the > HP 48GX calculator. Second choice would be to spell the names > out in full, "combinations" and "permutations". +1 These would be my preferences as well :-) |
|||

msg331743 - (view) | Author: Tim Peters (tim.peters) * | Date: 2018-12-13 05:58 | |

My two cents: - Spell it comb(n, k). - TypeError if args aren't ints. - ValueError if not 0 <= k <= n. Not concerned about speed. It's possible to do this with no divisions involving integers bigger than `n` and `k`(*), using O(k) space, but for "practical" arguments I bet that's slower than the dumbest possible loop. (*) 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. |
|||

msg331748 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2018-12-13 08:38 | |

[Tim] > My two cents: > > - Spell it comb(n, k). > - TypeError if args aren't ints. > - ValueError if not 0 <= k <= n. +1 to all of this. |
|||

msg331859 - (view) | Author: Tim Peters (tim.peters) * | Date: 2018-12-14 19:37 | |

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] << shift2 I'm NOT suggesting to adopt this. Just for history in the unlikely case there's worldwide demand for faster `comb` of silly arguments ;-) |
|||

msg332817 - (view) | Author: Yash Aggarwal (FR4NKESTI3N) * | Date: 2018-12-31 13:44 | |

Can I work on C implementation if no-one else is doing it right now? |
|||

msg332826 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2018-12-31 18:26 | |

> 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? |
|||

msg332838 - (view) | Author: Raymond Hettinger (rhettinger) * | Date: 2018-12-31 22:10 | |

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). |
|||

msg332844 - (view) | Author: Yash Aggarwal (FR4NKESTI3N) * | Date: 2019-01-01 08:28 | |

Thanks @mark.dickinson. As @rhettinger suggested, I'll write a basic function that uses division and works in O(k) for now. It's holiday season but hopefully @kellerfuchs will respond by then, and in the meantime I'll write more tests other than pascal's identity and corner cases. |
|||

msg332933 - (view) | Author: Yash Aggarwal (FR4NKESTI3N) * | Date: 2019-01-03 15:38 | |

I have written the function in the latest patch to work only for positive n. Although the definition of combination or nChoosek makes no sense for negative n, negative binomial distribution exists and so binomial coefficient is defined for negative value of n. So my question is should the function be expanded to calculate for negative n or is the function expected to work only in combination sense? |
|||

msg332957 - (view) | Author: Steven D'Aprano (steven.daprano) * | Date: 2019-01-04 02:26 | |

> should the function be expanded to calculate for negative > n or is the function expected to work only in combination sense? If this were my design, I would offer both but in separate functions: def comb(n, k): if n < 0: raise ValueError return bincoeff(n, k) def bincoeff(n, k): if n < 0: return (-1)**k * bincoeff(n+k+1, k) else: # implementation here... I believe we agreed earlier that supporting non-integers was not necessary. Are you also providing a perm(n, k) function? |
|||

msg332960 - (view) | Author: Steven D'Aprano (steven.daprano) * | Date: 2019-01-04 04:10 | |

> return (-1)**k * bincoeff(n+k+1, k) Oops, that's meant to be n+k-1. |
|||

msg332987 - (view) | Author: Yash Aggarwal (FR4NKESTI3N) * | Date: 2019-01-04 17:39 | |

@steven.daprano > Are you also providing a perm(n, k) function? I didn't know it is also being implemented. Should I start on that too? My implementation is based on these requirements: > - Spell it comb(n, k). > - TypeError if args aren't ints. > - ValueError if not 0 <= k <= n. Should the bincoeff function be same with exception of allowing negative n? |
|||

msg332988 - (view) | Author: Tim Peters (tim.peters) * | Date: 2019-01-04 17:56 | |

Please resist pointless feature creep. The original report was about comb(n, k) for integer n and k with 0 <= k <= n and that's all. Everyone who commented appeared to agree they'd find that useful. But nobody has said they'd find generalizing beyond those constraints USEFUL, or that they'd find perm(n, k) USEFUL. They just pointed out that such things are possible. Every bit of new API implies eternal maintenance, porting, testing, doc, and conceptual costs. So please stick to what's (at least nearly) universally agreed to be useful. Complications can be added later if there turns out to be real demand for them. |
|||

msg332990 - (view) | Author: Yash Aggarwal (FR4NKESTI3N) * | Date: 2019-01-04 18:43 | |

@tim.peters Got it. |

History | |||
---|---|---|---|

Date | User | Action | Args |

2019-01-04 18:43:29 | FR4NKESTI3N | set | messages: + msg332990 |

2019-01-04 17:56:09 | tim.peters | set | messages: + msg332988 |

2019-01-04 17:39:22 | FR4NKESTI3N | set | messages: + msg332987 |

2019-01-04 04:10:34 | steven.daprano | set | messages: + msg332960 |

2019-01-04 02:26:15 | steven.daprano | set | messages: + msg332957 |

2019-01-03 15:38:05 | FR4NKESTI3N | set | messages: + msg332933 |

2019-01-02 19:31:23 | FR4NKESTI3N | set | pull_requests: + pull_request10812 |

2019-01-01 08:28:25 | FR4NKESTI3N | set | messages: + msg332844 |

2018-12-31 22:10:52 | rhettinger | set | messages: + msg332838 |

2018-12-31 18:26:09 | mark.dickinson | set | messages: + msg332826 |

2018-12-31 13:44:16 | FR4NKESTI3N | set | nosy:
+ FR4NKESTI3N messages: + msg332817 |

2018-12-14 19:37:40 | tim.peters | set | messages: + msg331859 |

2018-12-13 08:38:30 | mark.dickinson | set | messages: + msg331748 |

2018-12-13 05:58:18 | tim.peters | set | messages: + msg331743 |

2018-12-09 01:11:28 | rhettinger | set | messages: + msg331402 |

2018-12-07 23:53:02 | steven.daprano | set | messages: + msg331369 |

2018-12-07 23:40:31 | brett.cannon | set | title: The math module should provide a function for computing binomial coefficients -> Add a function for computing binomial coefficients to the math module |

2018-12-07 17:01:35 | kellerfuchs | set | messages: + msg331339 |

2018-12-07 16:52:49 | kellerfuchs | set | messages: + msg331337 |

2018-12-07 16:49:39 | kellerfuchs | set | messages: + msg331336 |

2018-12-07 16:47:27 | kellerfuchs | set | messages: + msg331335 |

2018-12-07 14:56:39 | steven.daprano | set | messages: + msg331325 |

2018-12-07 14:52:18 | steven.daprano | set | messages: + msg331323 |

2018-12-07 14:44:09 | steven.daprano | set | messages: + msg331318 |

2018-12-07 13:37:36 | mark.dickinson | set | messages: + msg331312 |

2018-12-07 13:10:08 | mark.dickinson | set | messages: + msg331309 |

2018-12-07 13:08:26 | mark.dickinson | set | messages: + msg331308 |

2018-12-07 11:32:22 | serhiy.storchaka | set | messages: + msg331296 |

2018-12-07 11:30:06 | serhiy.storchaka | set | nosy:
+ serhiy.storchaka messages: + msg331295 |

2018-12-07 11:25:42 | kellerfuchs | set | messages: + msg331293 |

2018-12-07 11:20:51 | kellerfuchs | set | keywords:
+ patch stage: patch review pull_requests: + pull_request10253 |

2018-12-07 09:03:05 | mark.dickinson | set | messages: + msg331281 |

2018-12-07 08:42:05 | mark.dickinson | set | messages: + msg331280 |

2018-12-07 00:04:44 | rhettinger | set | nosy:
+ rhettinger, mark.dickinson, tim.peters messages: + msg331257 |

2018-12-06 23:22:43 | kellerfuchs | set | messages: + msg331256 |

2018-12-06 22:31:53 | steven.daprano | set | nosy:
+ steven.daprano messages: + msg331255 |

2018-12-06 21:18:02 | kellerfuchs | create |