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

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

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.

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.

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

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

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

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

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.

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

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

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.

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

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

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

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

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

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

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.

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

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

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

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

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.

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?

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

>         return (-1)**k * bincoeff(n+k+1, k)

Oops, that's meant to be n+k-1.

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

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.

@tim.peters
Got it.

Hi everyone,

Sorry for the lack of reply, I severely underestimated how exhausting the holiday season would be.
Catching up with the comments right now.

> 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

+1 from me on that.

@Yash: Thanks a bunch for starting on the implementation.  I will have a look shortly  :)

So, I rebased Yash's and my branches, and merged them together.  The result is still in PR#11018.

This involved a few changes, which seem to reflect the consensus here:
- raise ValueError if k>n ;
- rename the function to math.combinations.

@Raymond Hettinger
> Let's name this comb() instead of binomial() please (as requested by me, Mark, and Tim).

(Replying here to keep the discussion in a single place.)

As far as I can tell from the discussions here, Steven and you stated a preference for the shortened names, and that's it.
There was also no reply to my comment about `comb` being confusing (due to the collision with an English word).

Since there was, however, pretty clear agreement on calling it after combinations (shortened or not) rather than binomial(), I went with this.

Steven: I'm assuming Brett rearranged the title to put emphasis on the new function and to place it earlier in the title. Especially important if you're reading e-mails with the old subject on an e-mail client with limited subject preview lengths, you end up seeing something like:

"The math module should provide a function for computing..."

rather than the more useful:

"Add a function for computing binomial coefficients to t..."

Sorry for the late reply, I missed Tim's comment when it first came 
through.

> 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 [...] that they'd find perm(n, k) USEFUL.

I'm not going to argue for binomial coefficients with negative n, but I 
find it hard to imagine anyone needing combinations without also needing 
permutations, and I didn't think it was necessary to explicitly say so.

But since you insist, I'll say so: I would find it useful to have a 
function to compute the number of permutations of n taking k at a time.

My perspective may be biased from my experience with secondary school 
maths, where they are taught together, but providing one without the 
other strikes me as weird as providing tan without sin and cos.

There are other functions from combinatorics which I personally use, 
like derangements, but I know when I'm pushing my luck :-)

> This involved a few changes, which seem to reflect the consensus here:
> - raise ValueError if k>n ;
> - rename the function to math.combinations.

I see at least four people (myself, Raymond, Mark and Tim) giving comb 
as first choice, and I didn't see anyone give combinations as their 
first choice.

I don't object to you taking it upon yourself to go with the longer name 
(which is my second choice), but I do object to you claiming concensus 
for the change without evidence of such.

> There was also no reply to my comment about `comb` being confusing 
> (due to the collision with an English word).

To be honest, I didn't think that comment needed a reply.

Collisions between words in different knowledge domains are not unusual. 
I don't think people think that math.tan has anything to do with 
changing the colour of their skin, or math.sin is about being wicked. 
Due to their length, permutation, combination and factorial are 
frequently abbreviated to perm, comb, fact and anyone looking for those 
functions should recognise the abbreviations.

But given the precedent set by itertools and math.factorial, perhaps you 
are right and we ought to stick to the longer name.

> But given the precedent set by itertools and math.factorial,
> perhaps you are right and we ought to stick to the longer name

I disagree.  Let's stick with comb() which is what SciPy uses.  It is tedious to write out combinations in a formula what includes other steps. Also, I really want to differentiate it from the existing combinations() in the itertools module (it is reasonably forseeable that the two will used together.

Agreed, comb sounds much better than combination. And using the name binomial would make it sound like something that would puke out whole binomial series rather than a single coefficient(maybe leave it for that in case is it decided to be useful in the future).

PR 11414 implements simple algorithm that performs slower than using a factorial definition for k>n/3. 
@kellerfuchs I'd prefer if we could work on this since it's conflict free and already reflects the behavior everyone agreed upon. 

Would it be better to create a separate issue for math.perm to discuss its behavior?


If the behavior of comb is satisfactory, can we start with optimizations?

> As far as I can tell from the discussions here, Steven
> and you stated a preference for the shortened names, and
> that's it.

Plus Mark, plus me - all backed "comb()" specifically.

> I find it hard to imagine anyone needing combinations without
> also needing permutations, and I didn't think it was necessary
< to explicitly say so.

Of course it is.  Merely saying something is possible is no reason at all to do it.  The original report didn't say anything about counting partial permutations, and so it's "feature creep" on the face of it to tack that on.

I personally have scant use for `perm()`, but have written my own `comb()` many times.  Far more often than I've written a `factorial()` and a `product()` combined, but I've written each of the latter more than twice, and a `perm()` not even once.  Especially if `prod()` (the topic of a different report) is added, the case for adding a `perm()` gets weaker (rising and falling factorials are special cases of what `prod()` does, and `perm()` is just an instance of falling factorial).

Which doesn't mean `perm()` must not be added ;-)  You're now the first to say it _would_ be useful, which is a start.  Can we get a second?

In any case, I don't want to _see_ this report get bogged down by feature creep:  `comb()` is what it was about from the start, and everyone so far has agreed `comb()` would be useful.

@Steven
> > This involved a few changes, which seem to reflect the consensus here:
> > - raise ValueError if k>n ;
> > - rename the function to math.combinations.
> [...]
> > As far as I can tell from the discussions here, Steven and you stated a preference for the shortened names, and that's it.
> > There was also no reply to my comment about `comb` being confusing (due to the collision with an English word).
> >
> > Since there was, however, pretty clear agreement on calling it after combinations (shortened or not) rather than binomial(), I went with this.
>
> I see at least four people (myself, Raymond, Mark and Tim) giving comb 
as first choice, and I didn't see anyone give combinations as their 
first choice.
>
> I don't object to you taking it upon yourself to go with the longer name 
> (which is my second choice), but I do object to you claiming concensus 

I wasn't claiming consensus on the short-vs.-long name issue, but on the binomial-vs-combinations one.
I thought that would have been clear considering the context quoted above (which was missing from your reply)

Given that people clarified they prefer comb(), and that people conspicuously didn't comment on it being entirely-opaque to people who do not elready know what it is, I guess there is indeed consensus.

> Given that people clarified they prefer comb(), and that people conspicuously didn't comment on it being entirely-opaque to people who do not elready know what it is, I guess there is indeed consensus.

commit afb3d36e82b8d690a410fa9dca8029a8fad42984
Author: The Fox in the Shell <KellerFuchs@hashbang.sh>
Date:   Fri Feb 1 15:42:11 2019 +0100

    Rename math.combinations to math.comb

Can I get a consensus on weather math.perm() is needed?

> Can I get a consensus on weather math.perm() is needed?

It's not, for the purposes of this issue. I think `math.perm` should be the subject of a separate issue and discussion, and a separate PR. That way it doesn't block completion of this issue.

@mark.dickinson
Ok, then I will work on comb for now then.

@kellerfuchs @FR4NKESTI3N

We seem to have two open PRs for this feature, both with some review comments. Which is the one true PR? Can we close the other one?

@mark.dickinson both pr's are more or less same. Keller was offline for some time so I made the new issue. They were merged later. Only difference is in unittests. I'd say it's up to you to decide which one to keep.

New changeset 4a686504eb2bbf69adf78077458508a7ba131667 by Raymond Hettinger (Yash Aggarwal) in branch 'master':
bpo-35431: Implemented math.comb (GH-11414)
https://github.com/python/cpython/commit/4a686504eb2bbf69adf78077458508a7ba131667

Leaving this open for a while in case there is more work that needs to be done or any further comments to be resolved.

Thanks, Raymond. I'm planning to do a post-merge review and test this weekend.

Thanks @rhettinger for cleaning up the code and closing the pr. I didn't get what you meant by long, and sorry for not being much active as well. I am stuck with a pretty time consuming internship.

PR 13725 is a technical fix/optimization/cleanup. Later we can apply algorithmic optimization.

math.comb is leaking, I opened https://bugs.python.org/issue37125 to track it.

New changeset 2b843ac0ae745026ce39514573c5d075137bef65 by Serhiy Storchaka in branch 'master':
bpo-35431: Refactor math.comb() implementation. (GH-13725)
https://github.com/python/cpython/commit/2b843ac0ae745026ce39514573c5d075137bef65

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 `__index__`), in just the same way that the regular math module functions (sin, log, atan, sqrt, ...) accept arbitrary objects defining `__float__`. We already do this for gcd, factorial.

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.

Just a heads up from issue37125: The leak is was fixed by PR13725.

> What are your thoughts?

Sigh. I don't object to extending to `k < 0` and `k > n`, but once we've made that extension it's impossible to undo if we decide we'd rather have had the error checking. I'd really like to see some convincing use-cases. Quotes from Concrete Mathematics (fine book though it is) don't amount to use-cases.

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

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 `comb()` has gotten confused.  It certainly would have been "an error" in any programming use for `comb()` I've ever had.

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

> the sum of the first n positive integers is binomial(n+1, 2), and the formula should still work if n is zero.

I agree that's a convincing use-case.

Can you think of any practical uses for allowing negative `k`? I can't.

> @David Ratcliffe

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 `comb()` _define_ the result via the factorial expression.)

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 `comb(n+1, 2)`.

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 `comb()` at least a dozen times, I've never found that constraint limiting, and enforcing it _has_ caught errors in my code.

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
https://www.wolframalpha.com/input/?i=10+choose+12

I'm going to repeat part of an earlier comment :-)

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

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 `choose` with those arguments would be far better served if the function complained.

I understand that pure mathematics is not the primary use case. But I would
expect that math.comb() would be used to implement mathematical formulas.
As a simple example, the number of edges in a complete graph on n vertices
is binomial(n, 2), and this should be valid if n < 2. Polynomial
interpolation is another problem where this can occur. However, I do agree
that binomial coefficients with negative arguments are relatively
unimportant.

On Sun, Jun 2, 2019 at 8:59 PM Raymond Hettinger <report@bugs.python.org>
wrote:

>
> Raymond Hettinger <raymond.hettinger@gmail.com> added the comment:
>
> 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
> https://www.wolframalpha.com/input/?i=10+choose+12
>
> ----------
>
> _______________________________________
> Python tracker <report@bugs.python.org>
> <https://bugs.python.org/issue35431>
> _______________________________________
>

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 
binomial coefficients with negative arguments. They are used in 
fractional calculus

https://nrich.maths.org/1365

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 
bounds values (e.g. "choose 5 items from 4"), and a seperate binomial() 
function that accepts any positive or negative integer values. Even I 
think that fractional and complex arguments are probably a bit too 
exotic for the std lib -- that's Mathematica territory.

And yes, Mathematica does accept fractional and complex arguments to 
the choose() function.

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 `comb()` is best suited for Python's general audience, and "number of ways to pick k things out of n without regard to order" makes clear sense for any pair of non-negative k and n.  But anything beyond that does not.

> If it were me, I'd drop the k <= n requirement

+1 from me.

Fine by me. The same change should be applied to `math.perm`.

New changeset 963eb0f4738456455b9bef7eb531b46805415208 by Raymond Hettinger in branch 'master':
bpo-35431: Drop the k <= n requirement (GH-13798)
https://github.com/python/cpython/commit/963eb0f4738456455b9bef7eb531b46805415208

New changeset 1b8a46d59734a77cd1f5ffcf3bdfcaafd58a87e7 by Serhiy Storchaka in branch 'master':
bpo-35431: Test math.comb() and math.perm() for OverflowError only on CPython. (GH-14146)
https://github.com/python/cpython/commit/1b8a46d59734a77cd1f5ffcf3bdfcaafd58a87e7

New changeset 914d6b79735e5eabaf4e4d77e3f2ad4eae0beb9a by Serhiy Storchaka in branch '3.8':
[3.8] bpo-35431: Test math.comb() and math.perm() for OverflowError only on CPython. (GH-14146) (#14226)
https://github.com/python/cpython/commit/914d6b79735e5eabaf4e4d77e3f2ad4eae0beb9a