Opened 7 years ago

Last modified 4 years ago

#17123 needs_info enhancement

Extending binomial(n,k) to negative integers n, k.

Reported by: pluschny Owned by:
Priority: minor Milestone: sage-6.4
Component: combinatorics Keywords: binomial
Cc: rws Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: u/chapoton/17123 (Commits, GitHub, GitLab) Commit: 0c994379b150158c84ee583369f2ee4428a35f44
Dependencies: Stopgaps:

Status badges

Description

A simple and coherent extension of the binomial function to negative integers n, k was outlined by M. J. Kronenburg in The Binomial Coefficient for Negative Arguments, http://arxiv.org/abs/1105.3689 (Thanks to John Palmieri for the reference.)

This extension amounts to define

def BINOMIAL(n, k):
    if n in ZZ and k in ZZ: 
        if n >= 0 and k >= 0:
            return binomial(n, k)
        if k >= 0:
            return (-1)^k*binomial(-n+k-1, k)
        if k <= n:
            return (-1)^(n-k)*binomial(-k-1, n-k)
        return 0 
    else:
        return binomial(n, k)

Here 'BINOMIAL' is the targeted version, 'binomial' the implemented version. The targeted behaviour is identical to the behaviour of the Maple and Mathematica function for negative integers n, k.

Change History (65)

comment:1 Changed 7 years ago by vbraun

You'll have our eternal gratitude if you can post a git branch with the change ;-)

comment:2 Changed 7 years ago by SimonKing

Certainly there should not be implemented a modified and renamed binomial function in addition to the existing monomial function. The existing should just be extended to the case of negative input.

comment:3 follow-up: Changed 7 years ago by pluschny

Volker, what is the reason for your sarcasm?

Simon, certainly. But only in this form I could test it against the implemented version without running into a recursion.

comment:4 Changed 7 years ago by vbraun

I'm not sarcastic, just trying to nudge you into the right direction ;-)

comment:5 in reply to: ↑ 3 ; follow-up: Changed 7 years ago by SimonKing

Replying to pluschny:

Volker, what is the reason for your sarcasm?

It could be that he thinks that fixing the problem is a matter of not more than 30 minutes. But if I understand correctly you have no experience with Sage development, yet. And I think Volker should take this into account.

Simon, certainly. But only in this form I could test it against the implemented version without running into a recursion.

The implemented version just raises an error on negative input. Instead of raising the error, the binomial function should simply modify the given negative arguments according to the formula that you present in the ticket description.

Then, the documentation of the binomial function should be modified to contain a description and a test for the new behaviour. And to be on the safe side, the test suite should be run.

When working on a ticket, one ought to create a git branch for the ticket. The documentation will tell you how. For example, if you have the git trac scripts installed, you can do git trac checkout 17123 on the commandline.

How to edit the binomial function? First way: Find out in what module the function is defined.

sage: binomial.__module__
'sage.functions.other'

Hence, take your favourite editor and edit `SAGE_ROOT/src/sage/functions/other.py

Second way: In a running Sage session, do

sage: edit(binomial, 'vim')

where vim is your favourite editor. It will open the correct file in the correct line (if not, report a bug!).

After saving the change, leave Sage and do make (which will also build the changes in the documentation) or better make test in order to also run the test suite of Sage.

If everything works, commit your changes (git commit -a) and push them to this ticket (with the git trac scripts, this is git trac push).

comment:6 in reply to: ↑ 5 Changed 7 years ago by SimonKing

Replying to SimonKing:

The implemented version just raises an error on negative input.

Oops, I stand corrected. No error.

comment:7 Changed 7 years ago by SimonKing

Volker, it actually seems to me that the change is not totally straight forward.

There is the _eval_ method, and the _evalf_ method. Probably both need to be changed.

And there is the _binomial_sym method. Will this need a change as well? It currently returns zero if the second argument is negative:

sage: binomial._binomial_sym(x,-5)
0

Should the return value be a different symbolic expression? After all, the value (according to the formula in the ticket description) is not zero if x evaluates to a negative number that is greater or equal to -5.

So, here I don't know a good answer.

comment:8 Changed 7 years ago by vbraun

I know that there is a learning curve, but once you figure it out it is a very efficient workflow. So what I'm trying to say is: Its worthwhile to figure out how to use git.

_eval_ and _evalf_ punt to either _binomial_sym (for symbolic expressions) or sage.rings.arith.binomial (for numbers), so both of the latter need to be changed.

comment:9 Changed 7 years ago by jhpalmieri

Note the current situation:

sage: binomial(-1, -1)
0

So presumably the change here would require a deprecation warning.

Apparently Maple and Mathematica use the proposed version. Is that sufficient justification to change our version? Note that currently we have binomial(n-1, k-1) + binomial(n-1,k) == binomial(n,k), and I believe this change would break that identity. (So I'm not sure I would agree with the word "coherent" describing this extension of the binomial function.) If we want to make this change, then maybe we should bifurcate and have a separate Pascal's triangle function which still preserves the identity? See also http://en.wikipedia.org/wiki/Pascal's_triangle#Extensions.

comment:10 Changed 7 years ago by vbraun

TFA talks about gamma functions in the abstract, so I would have thought that it matches the extension that we would all write down (and is in the wp article)... though I haven't read it ;-)

comment:11 Changed 7 years ago by pluschny

Volker, acronyms are nice for those who understand. Do I have to understand 'TFA' or 'wp'?

comment:12 Changed 7 years ago by pluschny

John, this is why I offered to discuss things first on sage-devel and hear more arguments.

"Apparently Maple and Mathematica use the proposed version. Is that sufficient justification to change our version?"

Of course not, although for me this is a strong hint that something is missing in the current Sage implementation.

But this is a reason:

sage: binomial(-1, -1)
0

More generally, binomial(z, z) != 1 is absurd, a bug IMO, and definitely a reason to change things. binomial(z, z) = 1 for all z is the first thing to be preserved. Next the relation of the diagonals to the Fibonacci numbers is another important thing which should be preserved (and the proposed extension does this although apparently Kronenburg misses to mention this in his paper). The additive formula is nice in the region where the formula applies, but the multiplicative formula is the more important one, a viewpoint taken by all major modern authors as far as I understand.

comment:13 Changed 7 years ago by darij

Please do not change the current version. No, Peter, your convention is not standard in any way. Read anything by Knuth and you will see that binomial(n, k) == 0 for negative k (no matter what n is) is standard, and this happens to be the current behavior (I don't know why people are saying that it is currently undefined). Having binomial(z, z) != 1 is collateral damage, but there is no way that could be fixed reasonably.

Feel free to add binomial_symmetrized or whatever else you want to call your function, however!

comment:14 Changed 7 years ago by rws

  • Cc rws added

comment:15 Changed 7 years ago by vbraun

We could also have an optional keyword argument binomial(n, k, extension='Knuth')

comment:16 Changed 7 years ago by rws

The pseudocode in the ticket description does not seem to fit the paper given as reference, quote:

For nonnegative integer n and integer k this reduces to [1]:

(1.2) binomial(n,k) = n!/(k!(n-k)!) if 0<=k<=n, 0 otherwise

... This results in the following binomial coefficient identity,
which with identity (1.2) allows computation of the binomial coefficient for all integer
arguments.

(2.1) Theorem 2.1. For negative integer n and integer k:

binomial(n,k) = (-1)^k * binomial(-n+k-1,k) if k>0,
                (-1)^(n-k) * binomial(-k-1,n-k) if k<=n,
                0 otherwise
...

so in my opinion the definition should rather be

def BINOMIAL(n, k):
    if n in ZZ and k in ZZ: 
        if n >= 0:
            return binomial(n, k)
        if k >= 0:
            return (-1)^k*binomial(-n+k-1, k)
        if k <= n:
            return (-1)^(n-k)*binomial(-k-1, n-k)
        return 0 
    else:
        return binomial(n, k)

This leaves a discrepancy between binomial(*,k<0) and BINOMIAL(*,k<0) in the area n<0, k<n and to get to the point I would ask darij to provide a Knuth reference that supports the behaviour of binomial for these values.

comment:17 Changed 7 years ago by pluschny

Ralf, a quick numerical check did not show any difference:

    for n in (-4..4):
        print "plu", [n], [BINOMIAL   (-n , -k) for k in (-6..6)]
        print "rws", [n], [BINOMIALrws(-n , -k) for k in (-6..6)]

So please give one counterexample (n,k) where the definitions differ.

But I do prefer my form over your form. I write

    if n >= 0 and k >= 0:
        return binomial(n, k)

and this makes it perfectly clear that in this region nothing changes, that all changes apply only to negative integers. I think this lowers the burden of understanding the code. Darij: "... your convention is not standard in any way."

(1) If it were standard it would be pointless to call it an 'extension'. (2) As an extension it is standard for more than a quarter of a century due to the fact that both Maple and Mathematica use this extension ever since the beginning.

(3) So even if you prefer not to use this extension in your papers using Sage will give you different results than many of your colleagues will get with different software.

And indeed this was for me the motivation to write this ticket: I found the confusion annoying which shows up in many places in the code on OEIS because the systems behave differently and that this fact is often overlooked. Nobody expects this to happen with a standard function like 'binomial'.

Therefore I would even prefer a "value error" (or what it is called) for arguments not in the range 0<=k<=n for integer k, n, over the current behaviour of Sage.

Darij: "I don't know why people are saying that it is currently undefined."

This is the difference between a formal "... otherwise 0" and a definition guided by mathematical considerations which aims to enlarge relations to the largest possible region of validity. (Instead accepting things like binomial(z, z) != 1 as a collateral damage.)

comment:18 Changed 7 years ago by nbruin

Maxima's documentation on the matter is silent, so I'd be hesitant to take it as authoritative on the matter, but your proposed code is at odds with maxima's for (n,n) with n negative (maxima gives zero in that case). Deviating from maxima requires careful planning because expressions can quite easily end up in maxima.

sage: [ (n,k) for n in [-10..10] for k in [-10..10] if BINOMIAL(n,k) !=maxima_calculus.binomial(n,k)]
[(-10, -10),
 (-9, -9),
 (-8, -8),
 (-7, -7),
 (-6, -6),
 (-5, -5),
 (-4, -4),
 (-3, -3),
 (-2, -2),
 (-1, -1)]

comment:19 Changed 7 years ago by darij

Peter:

(1) If it is an extension, then call it binomial_extension or whatever; do not usurp binomial.

(2) Maple and Mathematica do not set standards in combinatorics.

The 0 definition does enlarge a relation to its largest possible domain of validity -- the recurrence relation, to the whole integer lattice. Unlike the symmetry, the recurrence actually is used all the time in proofs. The 0 definition also matches with the idea that (-1)k (-n choose k) is the number of k-element multisets of elements of {1, 2, ..., n}. It is most definitely guided by mathematical consideration; check, e.g., the identities (5.22) to (5.26) in Graham-Knuth-Patashnik, specifically (5.22) (Vandermonde's convolution, an extremely important identity). It might not be the most interesting definition, but we should not be aiming for that; I don't think we want 1 + 2 + 3 + ... to yield -1/12 just because the most interesting definition for summing powers of positive integers involves the zeta function. We should be going for the definition which is the most reliable and standard in *mathematical* literature. And here it is the 0 one.

comment:20 Changed 7 years ago by pluschny

Darij: "If it is an extension, then call it binomial_extension or whatever; do not usurp binomial."

Volker's proposal above is: "We could also have an optional keyword argument binomial(n, k, extension='Knuth')". This looks like a good idea to me, perhaps better with the name 'symmetrical' instead of 'Knuth'. Darij: "Maple and Mathematica do not set standards in combinatorics."

But you can also not ignore them especially as Sage's mission is to become a viable alternative to theses systems.

"We should be going for the definition which is the most reliable and standard in *mathematical* literature. And here it is the 0 one."

I agree. But this does not mean that we should block interesting developments either. Volker's proposal shows the way.

On a side note: Don Knuth is known to be happily hacking with Mathematica for a long time http://www-cs-faculty.stanford.edu/~uno/screen.jpeg and reading his 'Convolution Polynomials' with it's many implicit and explicit uses of the binomial function I do not see him complain about the way Mathematica defines it.

comment:21 Changed 7 years ago by vbraun

As a string theorist, I approve of 1+2+3+... = -1/12 ;-)

comment:22 follow-ups: Changed 4 years ago by chapoton

I have just been hurt by that issue, and I would like to change the current behaviour. In the usual hypergeometric context, having binomial(-1,-1)=1 would be the only natural thing to do. The current binomial(-1,-1)=0 forces me to do special-casing everywhere in my code, which I take as a very good sign that this is a bad convention.

comment:23 in reply to: ↑ 22 Changed 4 years ago by pluschny

Replying to chapoton:

I have just been hurt by that issue, and I would like to change the current behaviour.

This is great! I have also written on my OEIS blog about it: Extensions of the Binomial

comment:24 in reply to: ↑ 22 ; follow-up: Changed 4 years ago by jhpalmieri

Replying to chapoton:

I have just been hurt by that issue, and I would like to change the current behaviour. In the usual hypergeometric context, having binomial(-1,-1)=1 would be the only natural thing to do. The current binomial(-1,-1)=0 forces me to do special-casing everywhere in my code, which I take as a very good sign that this is a bad convention.

As pointed out elsewhere on this ticket, making the proposed change breaks the identity binomial(n-1, k-1) + binomial(n-1,k) == binomial(n,k), which suggests to me that the proposed change is a bad convention. That is, there are arguments on both sides. I think adding an optional argument to get the behavior you want might be the best solution.

comment:25 Changed 4 years ago by chapoton

  • Branch set to u/chapoton/17123
  • Commit set to 0c994379b150158c84ee583369f2ee4428a35f44
  • Status changed from new to needs_info

I am not convinced at all by the desire (why ?) to get the usual recursion extend to all k and n. Instead, I am doing concrete and complicated hypergeometric computations right now, where I am put in trouble by the fact that sage currently says that binomial(-1,-1) is 0.

But I am worried by the different behaviour in Maxima, saying binomial(-1,-1)=0. We also have giac, that says binomial(-1,-1)=1.

sage: a=ZZ(-1)
sage: a._giac_().binomial(-1)
1
sage: a._maxima_().binomial(-1)
0
sage: a.__pari__().binomial(-1)
0

New commits:

0c99437trac 17123 change values of binomial(n, k) for negative k and n

comment:26 Changed 4 years ago by chapoton

And sympy does not know the answer:

sage: a=ZZ(-1)
sage: a._sympy_().binomial(-1)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-11-c88f55965f06> in <module>()
----> 1 a._sympy_().binomial(-Integer(1))

AttributeError: 'NegativeOne' object has no attribute 'binomial'

*EDIT*

Sympy doc says : For the sake of convenience for negative ‘k’ this function will return zero no matter what valued is the other argument.

Last edited 4 years ago by chapoton (previous) (diff)

comment:27 in reply to: ↑ 24 Changed 4 years ago by pluschny

Replying to jhpalmieri:

As pointed out elsewhere on this ticket, making the proposed change breaks the identity binomial(n-1, k-1) + binomial(n-1,k) == binomial(n,k), which suggests to me that the proposed change is a bad convention.

This is not true. With the solution shown in my blog linked above the identity remains valid. See here.

And it would be far more sensible instead of offering different versions because of some special cases to seek a consistent mathematical solution.

comment:28 follow-up: Changed 4 years ago by chapoton

Summary:

One single natural idea : use the expression of Binomial in term of Gamma functions

==> this gives unique and well-defined values for all integers (n,k)

==> for k >= 0, these values are the same as the ones coming from the polynomiality in n

==> the addition rule works everywhere but for (-1,-1)+(-1,0) != (0,0)

==> this works very well in the context of hypergeometric series, where the natural way to extend binomials is to pass to Gamma products.

==> and this satisfies the symmetry property (n,k) = (n, n-k)

So it seems to me that we have very good reasons to change. I am also personnally totally convinced that this is the only correct convention and that setting blandly to 0 was just done by some kind of lazyness.

Would it be useful to ask the opinion of one of the great masters of hypergeometric series ?

comment:29 Changed 4 years ago by darij

I am strongly against changing the results on the current domain of definition, whatever this is (it's been a while...). I don't remember whether it currently computes binomial(n, k) to 0 for k < 0, but if it does, we cannot just deprecate it. At the very least we need to create a function that replicates this behavior.

There are various sums in combinatorics which sort-of rely on the "zero for k < 0" behavior, because modifying their bounds to get rid of the k < 0 case is impractical (this leads to stupid bounds like floors and ceilings). I suspect that in many of these cases, the n on top of the binomial coefficient is nonnegative (and you are not making changes to this case, right?), but I wouldn't expect this to hold for ALL of these sums.

comment:30 Changed 4 years ago by chapoton

Hello, Darij. I hear your preoccupation, but I do not quite understand it. Do you have one concrete case to show me ?

I have one supporting the change : consider

sum(binomial(n, i) * binomial(n - i - l - 1, i - 1) *
               x ** i * y ** l for i in range(1, floor(n/2) + 1)
               for l in range(n - 2 * i + 1))

This gives for me the result that I expect if and only if binomial(-1,-1)=1.

The proposal would change the value when both n and k are negative, from 0 (now) to some integers obtained as limits of Gamma products.

Of course, we can easily propose either an "old_binomial" function or an "oldstyle" keyword. What would you prefer ?

comment:31 Changed 4 years ago by darij

What would the expected result be?

Here is an example:

binom = binomial

def lr(n, k):
    return sum(binom(n+u-1, u) * binom(n, k-2*u) for u in range(k+1)) - binom(n+k-1, k)

This should yield 0 for all n and k. (See http://www.cip.ifi.lmu.de/~grinberg/QEDMO4P13.pdf for a proof for n and k nonnegative. It then follows for all n because both sides are polynomials in n. Finally, the k < 0 case is simply a 0 = 0 statement.)

If your definition breaks it for negative k, okay. But if it breaks it for positive k and negative n, then your definition is bad. Changing the bound of the sum to something like floor(k/2) is exactly the kind of annoyance that a good definition avoids!

comment:32 Changed 4 years ago by chapoton

Well, this identity of yours stays true for all n and k>=0 if you use the floor bound.

I would like to argue that not using a floor bound is correct if n > 0, but is not the good way to go if n < 0.

Enough for today..

comment:33 Changed 4 years ago by rws

From all this it seems to me the best way to resolve this issue is by providing an additional function. Deprecation can then be discussed in a new ticket. The name C comes to mind. Please give your naming ideas too.

comment:34 in reply to: ↑ 28 ; follow-up: Changed 4 years ago by mantepse

One single natural idea : use the expression of Binomial in term of Gamma functions

==> this gives unique and well-defined values for all integers (n,k)

I must be misunderstanding: the function f(x,y) = gamma(x + 1) / (gamma(y + 1)*gamma(x - y + 1)) is singular when x is a negative integer. See "The Binomial Coefficient Function, David Fowler, The American Mathematical Monthly, Vol. 103, No. 1 (Jan., 1996), pp. 1-17" for an analysis.

Still, I think it would be practical to define the binomial coefficient in this way, throwing an error when the limit does not exist. Then it becomes very easy for the user to extend it whenever necessary, because she can catch the exception. This solution also has the advantage of making the conventions explicit in all methods.

I think the solution of defining the binomial coefficient as binomial(x0, y0)=limit(limit(f(x,y), y=y0), x=x0) (in this order) is inferior, because it yields surprising results, in particular binomial(-1,-1)=0.

comment:35 in reply to: ↑ 34 Changed 4 years ago by pluschny

Replying to mantepse:

One single natural idea : use the expression of Binomial in term of Gamma functions

==> this gives unique and well-defined values for all integers (n,k)

I must be misunderstanding: the function f(x,y) = gamma(x + 1) / (gamma(y + 1)*gamma(x - y + 1)) is singular when x is a negative integer. I think the solution of defining the binomial coefficient as binomial(x0, y0)=limit(limit(f(x,y), y=y0), x=x0) (in this order) is inferior, because it yields surprising results, in particular binomial(-1,-1)=0.

Yes, you are misunderstanding. Try this:

def limit_binomial(n, k):
    return limit(gamma(n + x) / (gamma(k + x)*gamma(n - k + x)), x = 1) 
for n in (-5..5): print [limit_binomial(n, k) for k in (-5..5)]

limit_binomial(-1,-1) = 1

All been said in my blog post linked above.

Last edited 4 years ago by pluschny (previous) (diff)

comment:36 follow-up: Changed 4 years ago by mantepse

OK, I thought that Frederic was referring to the other expression of the binomial in terms of the gamma function. I was further misled by another blog post saying that Mathematica uses the iterated limit definition.

So, one question is whether sage should use

limit(gamma(x0 + h) / (gamma(y0 + h)*gamma(x0-y0+h)), h=1)

or

limit(gamma(x + 1) / (gamma(y + 1)*gamma(x-y+1)), (x,y)=(x0,y0))

as definition, and you are saying that the first is better, because it is works always, right?

comment:37 in reply to: ↑ 36 Changed 4 years ago by pluschny

Replying to mantepse:

you are saying that the first is better, because it is works always, right?

Yes. And I'm sure that Mathematica also uses it contrary to your information because otherwise they would have similar problems as the ones you ran into. But they do not have problems.

Note that this argument (embedding into C in conformance to hypergeometric considerations) is only one argument.

Also from the combinatorial point of view it makes sense to define:

    def Binomial(n, k):
        if n in ZZ and k in ZZ: 
            if n >= 0:
                return binomial(n, k) 
            if k >= 0:
                return (-1)^k*binomial(-n+k-1, k) 
            if k <= n:
                return (-1)^(n-k)*binomial(-k-1, n-k) 
            return 0 
        return binomial(n, k)

You can use this definition with 'binomial' taking Sage's current binomial.

The reflection property of the binomial numbers is similar and related to the reflection property of the Stirling numbers {n,k} = [-k,-n] and of the Lah numbers L|n,k| = L| -k,-n|. This has been observed by Riordan, Knuth and others.

comment:38 follow-up: Changed 4 years ago by mantepse

I think that I dislike the fact that your definition (by the way: is it yours? otherwise, do you have a reference?) is not continuous.

Can we assume that everybody agrees that for x not negative integral, gamma(x+1)/(gamma(y+1)*gamma(x-y+1)) is a good definition of binomial(x,y)?

Do (at least the more important) packages agree on these values, and disagreement is only for negative integral x?

If so, wouldn't it still be best to define binomial with an optional argument, extension, which defaults to None. Then we could even emulate the behaviour of other packages, whenever useful.

I think it is safest to have the optional argument being None, sage could throw a useful error in this case. Then we get an additional check of the usage of binomial in the sage library.

comment:39 in reply to: ↑ 38 ; follow-ups: Changed 4 years ago by pluschny

Replying to mantepse:

"..your definition (by the way: is it yours? otherwise, do you have a reference?).."

I do not claim any originality whatsoever with this definition. For references see what I have written on my blog linked above, and the references given there.

"Can we assume that everybody agrees that for x not negative integral,

gamma(x+1)/(gamma(y+1)*gamma(x-y+1)) is a good definition of binomial(x,y)?"

Pardon me? No, what is g(x,y) for x = 1, y = -1 with this 'definition'?

An extension which does not give the reflection formula for the binomial

binomial(-k, -n) = (-1)^(n-k) binomial(n-1, k-1)

is completely worthless, potentially introducing new errors.

"I think it is safest to have the optional argument being None.. ""

Three years ago Volker Braun suggested above:

"We could also have an optional keyword argument binomial(n, k, extension='Knuth')."

comment:40 in reply to: ↑ 39 ; follow-up: Changed 4 years ago by darij

Replying to pluschny:

An extension which does not give the reflection formula for the binomial

binomial(-k, -n) = (-1)^(n-k) binomial(n-1, k-1)

is completely worthless, potentially introducing new errors.

What's so great about this reflection formula that it has to hold for all integers n and k?

Here's another reason for (n choose k) to be 0 when n < 0 and k < 0: The coefficient in this case counts certain multisubsets of size k. There are none of them when k < 0.

I also don't buy the hypergeometric-functions argument. Graham, Knuth and Patashnik, in Chapter 5 of CM, define binomial coefficients to be 0 when k < 0, and this chapter does a good deal of hypergeometric manipulations. Sure, feel free to introduce another convention, but don't change the defaults!

comment:41 in reply to: ↑ 39 Changed 4 years ago by mantepse

Replying to pluschny:

"Can we assume that everybody agrees that for x not negative integral,

gamma(x+1)/(gamma(y+1)*gamma(x-y+1)) is a good definition of binomial(x,y)?"

Pardon me? No, what is g(x,y) for x = 1, y = -1 with this 'definition'?

Sorry, I was sloppy. I should have added extended by continuity.

comment:42 in reply to: ↑ 40 Changed 4 years ago by pluschny

Replying to darij:

Replying to pluschny:

An extension which does not give the reflection formula for the binomial

binomial(-k, -n) = (-1)^(n-k) binomial(n-1, k-1)

is completely worthless, potentially introducing new errors.

What's so great about this reflection formula that it has to hold for all integers n and k?

What a question!

As a short answer I cite GKP, CM on what they write in the Stirling case:

"In fact, a surprisingly pretty pattern emerges: The two kinds of Stirling numbers are related by an extremely simple law:

[n, k] = {-k, -n}, integers k, n.

We have a "duality", something like the relations between min and max, between floor(x) and ceiling(x), between falling factorial and rising factorial, between gcd and lcm."

If such a relation exists for the binomial numbers, no CAS should presume to conceal it.

comment:43 in reply to: ↑ 39 ; follow-up: Changed 4 years ago by mantepse

Replying to pluschny:

Replying to mantepse:

"..your definition (by the way: is it yours? otherwise, do you have a reference?).."

I do not claim any originality whatsoever with this definition. For references see what I have written on my blog linked above, and the references given there.

OK, I checked and found the definition in Section 3.2 of https://arxiv.org/pdf/math/9502214.pdf. It would be nice to know who came up with this first.

comment:44 in reply to: ↑ 39 ; follow-up: Changed 4 years ago by mantepse

Replying to pluschny:

Three years ago Volker Braun suggested above:

"We could also have an optional keyword argument binomial(n, k, extension='Knuth')."

Yes, and I am saying that sage should have binomial(n, k, extension=None), which should give the continuous version and raise an error for negative integral n, and then binomial(n, k, extension='who_or_whatever') for the various extensions.

(Assuming that all the extensions agree except for negative integral n.)

Is there any downside to this approach?

comment:45 in reply to: ↑ 44 ; follow-up: Changed 4 years ago by pluschny

Replying to mantepse:

Replying to pluschny:

Three years ago Volker Braun suggested above: "We could also have an optional keyword argument binomial(n, k, extension='Knuth')."

Yes, and I am saying that sage should have binomial(n, k, extension=None), which should give the continuous version and raise an error for negative integral n, and then binomial(n, k, extension='who_or_whatever') for the various extensions. Is there any downside to this approach?

This would effectively set the present state forever, which I consider unsatisfactory.

Chapoton reopened this issue yesterday with the words: "I have just been hurt by that issue."

And this issue will continue to hurt other users. I have been hurt by that issue several times while developing software for the OEIS. In the hypergeometric context you 'feel' it, in many other cases it will bite unnoticed.

This is a bug and should be handled like a bug. I therefore propose to follow Johan S. R. Nielsen suggestion on sage-devel:

I vote for having #17123 as the default *in the long run*. If so, what to do in the short run to mitigate user problems. We could

1) Return the value of #17123 but print a deprecation warning on input

with conflicting behaviour. The warning tells the user to explicitly set the optional argument if he wants to disable the warning. Warning is removed in ~1 year.

2) Return the current value, but print a deprecation warning on input

with conflicting behaviour. After one year, change the behaviour with no more deprecation warnings.

comment:46 in reply to: ↑ 45 Changed 4 years ago by mantepse

Yes, and I am saying that sage should have binomial(n, k, extension=None), which should give the continuous version and raise an error for negative integral n, and then binomial(n, k, extension='who_or_whatever') for the various extensions. Is there any downside to this approach?

This would effectively set the present state forever,

This is not true:

currently binomial(-1,-1) gives 0, while I propose that it raises an exception (possibly after a deprecation period), forcing you to use an explicit extension.

You and Frederic might use binomial(-1, -1, extension='supertrooper') whereas Darij might use binomial(-1,-1, extension='themanwhosoldtheworld'), but in any code you write, the extension chosen would be explicit!

Chapoton reopened this issue yesterday with the words: "I have just been hurt by that issue."

Yes, because sage currently silently returns 0, which is not what he expected.

And this issue will continue to hurt other users.

No, please read me proposal again.

Again, what is the downside to my approach?

Martin

comment:47 follow-up: Changed 4 years ago by darij

These look like good options, as long as:

1) the new default agrees with the old default on all values where both are defined (i.e., don't raise an error);

2) the new default includes (n choose k) = n (n-1) ... (n-k+1) / k! for k >= 0 (so the only disagreements are about what happens when k < 0) (e.g., this means that the Roman convention cannot be the default);

3) the case branching and string passing coming from the various conventions does not significantly slow down the function, OR each of the different conventions can be accessed through its own function without passing a string. (I'm not saying these functions need to be in the global namespace; it's enough if code can access them.)

BTW, the paper "A generalization of the binomial coefficients" by Loeb ( http://www.sciencedirect.com/science/article/pii/0012365X92901386 and https://arxiv.org/abs/math/9502218 ) seems to give a good overview of the different conventions (but the zero convention from Concrete Mathematics is not among them -- so there are four now: "zero", "Roman", "classical" and "Gamma"). I am still fairly convinced that the "zero" convention is the one that makes combinatorics the least painful, and that symmetry formulae are less important than the recursion.

comment:48 Changed 4 years ago by mantepse

I think it is calling for trouble not to raise an error for negative integral first argument. As much as I dislike python, I think "explicit is better than implicit" is a sensible rule.

If some functionality relies on binomial(-1,-1)=1, it will be much easier to maintain it if this stated explicitely, for example:

def fun():
    binomial = lambda x,y: binomial(x, y, extension='Gamma')
    return binomial(-1,-1)

comment:49 in reply to: ↑ 47 ; follow-up: Changed 4 years ago by pluschny

Replying to darij: I am still fairly convinced that the "zero" convention is the one that makes combinatorics the least painful ...

pluschny wrote: "I have been hurt by that issue several times while developing software for the OEIS. In the hypergeometric context you 'feel' it, in many other cases it will bite unnoticed."

I will give an example for the latter. Consider:

PartitionCoefficient = lambda p: mul(binomial(p[j], p[j+1]) for j in range(len(p)-1))
for n in (0..6):
    for k in (0..n):
        P = Partitions(n, max_part=k, inner=[k])
        print sum(PartitionCoefficient(p) for p in P), binomial(n-1,k-1) 

What this code shows is for me inconsistency.

darij, you coined the famous sentence in this thread: "Having binomial(z, z) != 1 is collateral damage..." of your favourite definition.

This includes the damage of inconsistency?

comment:50 in reply to: ↑ 49 Changed 4 years ago by darij

Replying to pluschny:

pluschny wrote: "I have been hurt by that issue several times while developing software for the OEIS. In the hypergeometric context you 'feel' it, in many other cases it will bite unnoticed."

I will give an example for the latter. Consider:

PartitionCoefficient = lambda p: mul(binomial(p[j], p[j+1]) for j in range(len(p)-1))
for n in (0..6):
    for k in (0..n):
        P = Partitions(n, max_part=k, inner=[k])
        print sum(PartitionCoefficient(p) for p in P), binomial(n-1,k-1) 

You should be comparing with binomial(n-1, n-k). These are really partitions of n-k you are summing over (the first part is just for convenience, as you clamp it to k).

Sorry, but it makes no sense on earth, in heaven and in other places to diverge from the (n choose k) = n(n-1)...(n-k+1)/k! standard for negative n as long as k >= 0. The binomial coefficients (n choose k) for a fixed k are a polynomial in n when n is nonnegative; why on earth would you want to break that?

Last edited 4 years ago by darij (previous) (diff)

comment:51 Changed 4 years ago by mantepse

Darij, are you sure that Peter proposal changes the default for k>=0?

comment:52 Changed 4 years ago by darij

Oops, maybe not!

Sorry, this thread is a mess (and I'm part of it). Can anyone update me on whether the k >= 0 case is settled?

comment:53 Changed 4 years ago by mantepse

As far as I can see, Peter is proposing the "Gamma" convention from https://arxiv.org/pdf/math/9502218.pdf. I think that this agrees with what is roughly the "classical" convention in the paper whenever the latter is defined, that is limit(gamma(x+1)/(gamma(y+1)*gamma(x-y+1)), (x,y)=(n,k)) whenever the limit exists.

So they agree in particular for k>=0.

However, I still think it would be better to force users to make an explicit choice for negative integral k.

Version 0, edited 4 years ago by mantepse (next)

comment:54 Changed 4 years ago by mantepse

I'm afraid that I also added to the confusion, please let me correct my mistake.

The truth is, limit(gamma(x+1)/(gamma(y+1)*gamma(x-y+1)), (x,y)=(n,k)) does not exist for negative integral n, even if k is positive integral. Which means, that by default binomial(-2, 41) would be undefined instead of -42.

comment:55 in reply to: ↑ 43 Changed 4 years ago by pluschny

Replying to mantepse:

Replying to pluschny:

Replying to mantepse:

"..your definition (by the way: is it yours? otherwise, do you have a reference?).."

I do not claim any originality whatsoever with this definition. For references see what I have written on my blog linked above, and the references given there.

OK, I checked and found the definition in Section 3.2 of https://arxiv.org/pdf/math/9502214.pdf. It would be nice to know who came up with this first.

Meanwhile I have found a reference for the reflexion formula: Martin Aigner, Kombinatorik I, Springer 1975, exercise on page 149. This book became part of his Combinatorial Theory, Grundlehren 234.

Since I learned from this book, it seemed always to me to be self-evident and I do not understand to this day why one should be reluctant to use it (or implement an extension which assures it).

comment:56 follow-up: Changed 4 years ago by mantepse

let me please try to summarize:

  1. all agree that it would be good to have the "gamma" convention from https://arxiv.org/pdf/math/9502218.pdf accessible
  1. also the "zero" convention (I think that's equivalent to the iterated limit expression) should be accessible, and perhaps also other conventions
  1. some formulas involving binomial coefficients depend on the chosen convention, and different CAS make different choices
  1. we also need binomial for non-integral arguments
  1. a continuous definition is desirable, for example for analytic stuff, but limit(gamma(x+1)/(gamma(y+1)*gamma(x-y+1)), (x,y)=(n,k)) only exists when n is not a negative integer, which is very restrictive

possible solutions:

  1. make some convention (very likely "gamma" - limit(gamma(n+h)/(gamma(k+h)*gamma(n-k+h)), h=1)) the default and make other conventions accessible
    • advantages: most users won't notice
    • disadvantages: convention-mismatch will go unnoticed
  1. binomial(n, k) raises an error for negative integral n (even when k is a non-negative integer), and in this situation the convention must be chosen explicitly
    • advantages: convention-mismatch will go unnoticed
    • disadvantages: slightly harder to use and routines calling binomial will probably need adjustment in some places, although these should be easy to find

comment:57 in reply to: ↑ 56 Changed 4 years ago by pluschny

Replying to mantepse:

let me please try to summarize:

I think there should be exactly two functions:

  • one defined on N X N where N = {0,1,2,..},
  • the other defined on C where C are the complex numbers.

This mirrors the case factorial / Gamma.

The first defined on N X N is the classical Pascal function which is zero for all (n,k) which are not 0<=k<=n.

I leave open how the second is precisely defined, but it should reduce on the ZZ-grid to the values described by the Binomial function given above. This means in particular that on the ZZ-grid the reflection formula is valid, all values finite and Binomial(x, x) = 1 for all x.

For the naming: the first one could be called binomial(n,k) and the second Binomial(x,y).

This means in particular that the existing Sage binomial function needs not to be touched. No disadvantages, no deprecations.

comment:58 follow-up: Changed 4 years ago by mantepse

Could you please comment on the advantages and disadvantages I listed.

Moreover, what is the disadvantage of binomial(n,k,extension="Gamma") over Binomial(n,k)?

Finally, how would you accommodate Darij's needs with your proposal?

comment:59 in reply to: ↑ 58 ; follow-up: Changed 4 years ago by pluschny

Replying to mantepse:

Finally, how would you accommodate Darij's needs with your proposal?

Your question amazes me. Darij wrote: "Please do not change the current version." Done. "Feel free to add binomial_symmetrized or whatever else you want to call your function, however!" Fine. So I think this will be OK with him. But he can take position himself, no?

Could you please comment on the advantages and disadvantages I listed.

I concur with what you say.

Moreover, what is the disadvantage of binomial(n,k,extension="Gamma") over Binomial(n,k)?

No big difference for me. I think it is the analogy with factorial/Gamma which I like. Of course we could also drop the Gamma function and introduce factorial(n,extension="Gamma"). Matter of taste.

comment:60 in reply to: ↑ 59 ; follow-up: Changed 4 years ago by mantepse

Replying to pluschny:

Replying to mantepse:

Finally, how would you accommodate Darij's needs with your proposal?

Your question amazes me. Darij wrote: "Please do not change the current version." Done. "Feel free to add binomial_symmetrized or whatever else you want to call your function, however!" Fine. So I think this will be OK with him. But he can take position himself, no?

The current binomial is not defined on N x N, but on a larger domain. Above you proposed at the same time to leave it untouched and to define it on N x N, which confused me.

comment:61 in reply to: ↑ 60 Changed 4 years ago by pluschny

Replying to mantepse:

The current binomial is not defined on N x N, but on a larger domain. Above you proposed at the same time to leave it untouched and to define it on N x N, which confused me.

I see, that was not well expressed. In fact, I am now a little tired of this subject and have nothing else to say. May the issue get a good implementer. Bye.

comment:62 follow-up: Changed 4 years ago by pluschny

I have just read an unpublished manuscript, which is the most detailed and comprehensive analysis of the subject so far. It convinced me that the approach defended by me is not adequate. I withdraw my request and ask to close the subject.

comment:63 Changed 4 years ago by darij

If any of you is interested in a fast review (at least from me), I suggest to implement whatever definitions you like as separate functions:

def binomial_roman...

def binomial_gamma...

etc. (without changing the existing binomial function), document their definitions in full and link them to each other via .. SEEALSO:. This way I'll be able to just review the maths without plunging into a debate I really don't want to take part in (nor do I have the time to).

Last edited 4 years ago by darij (previous) (diff)

comment:64 in reply to: ↑ 62 ; follow-up: Changed 4 years ago by mantepse

Replying to pluschny:

I have just read an unpublished manuscript, which is the most detailed and comprehensive analysis of the subject so far. It convinced me that the approach defended by me is not adequate. I withdraw my request and ask to close the subject.

Could you share this, or maybe just the address of the author?

I am actually digging into the code and the details of the conventions, to see how it goes.

comment:65 in reply to: ↑ 64 Changed 4 years ago by pluschny

Replying to mantepse:

Replying to pluschny:

I have just read an unpublished manuscript,

Could you share this, or maybe just the address of the author?

Several people urged the author to publish the paper. But it is, of course, at the discretion of the author. Nevertheless, it may one day appear in the arXiv and you will recognize it by the excellent plots.

In his opinion, no path leads past the arguments of GKP in CM, with the effect that all integer values on the left half-plane have to be 0. In particular GKP writes:

"The symmetry identity fails for all other negative integers n, too. But unfortunately it's all too easy to forget this restriction, since the expression in the upper index is sometimes negative only for obscure (but legal) values."

Meanwhile I read yet another paper which is much more critical about the presentation of GKP:

David Fowler, The Binomial Coefficient Function, The American Mathematical Monthly, Vol. 103, No. 1 (Jan., 1996), pp. 1-17

Fowler also makes the following interesting suggestion:

"One could generalise the standard binomial coefficients to include an extra argument, the slope at which the directional limit is to be taken, and thus extend such standard identities to negative arguments and perhaps find new ones."

Note: See TracTickets for help on using tickets.