Opened 4 months ago

Incorrect summation with binomial of non-atomic arguments.

Reported by: Owned by: charpent major sage-9.5 symbolics summation, binomial, maxima gh-DaveWitteMorris, slelievre Reported upstream. No feedback yet.

This sum of binomials (expressed in two ways) is incorrect:

```sage: t, n, k = SR.var('t, n, k')
sage: sum(t^(k - 1)*binomial(n - 1, k - 1), k, 1, n)  # incorrect
0
sage: sum(t^(k)*binomial(n - 1, k), k, 0, n - 1)  # incorrect
0
```

Expected in both cases: `(t + 1)^(n - 1)`.

These somewhat similar (but different) sums are correct:

```sage: sum(t^(k - 1)*binomial(n, k - 1), k, 1, n + 1)
(t + 1)^n
sage: sum(t^(k)*binomial(n, k), k, 0, n)
(t + 1)^n
sage: sum(t^(k - 1)*binomial(n, k - 1), k, 1, n)
(t + 1)^n - t^n
sage: sum(t^(k)*binomial(n, k), k, 0, n - 1)
(t + 1)^n - t^n
```

Details and workaround in Emmanuel Charpentier's answer to Ask Sage question 57195.

Upstream report:

comment:1 Changed 4 months ago by charpent

My hunch is that the problem lies with the management of non-atomic arguments by `binomial`, but I haven't a clue on how to explore this.

Suggestions ?

comment:2 Changed 4 months ago by slelievre

• Description modified (diff)

comment:3 Changed 4 months ago by gh-DaveWitteMorris

Contrary to what is said in the description, the "equivalent sum obtained by shifting the summation index" is not correct. In the equivalent sum, the top of the binomial coefficient needs to stay at `n-1`, not change to `n`:

```sage: sum(t^(k)*binomial(n-1, k), k, 0, n-1)
0
```

I don't know what is going on, but here is some additional data. In this example, only `j=-1` is wrong:

```sage: for j in range(-3, 4):
....:     j, sum(t^k * binomial(n + j, k), k, 0, n + j)
....:
(-3, sum(t^k*binomial(n - 3, k), k, 0, n - 3))
(-2, sum(t^k*binomial(n - 2, k), k, 0, n - 2))
(-1, 0)
(0, (t + 1)^n)
(1, (t + 1)^(n + 1))
(2, (t^2 + 2*t + 1)*(t + 1)^n)
(3, (t^3 + 3*t^2 + 3*t + 1)*(t + 1)^n)
```

But in this example, all of the negative values of `j` are wrong:

```sage: for j in range(-3, 4):
....:     j, sum(t^k * binomial(n + j, k), k, 0, n)
....:
(-3, (t + 1)^n)
(-2, (t + 1)^n)
(-1, (t + 1)^n)
(0, (t + 1)^n)
(1, (t + 1)^(n + 1) - t^(n + 1))
(2, (t^2 + 2*t + 1)*(t + 1)^n - ((n + 2)*t + t^2)*t^n)
(3, (t^3 + 3*t^2 + 3*t + 1)*(t + 1)^n
- 1/2*(2*(n + 3)*t^2 + 2*t^3 + (n^2 + 5*n + 6)*t)*t^n)
```

comment:4 Changed 4 months ago by charpent

Contrary to what is said in the description, the "equivalent sum obtained by shifting the summation index" is not correct. In the equivalent sum, the top of the binomial coefficient needs to stay at n-1, not change to n:

Oh yeah ?

```sage: var("n, t, j, k, n1, k1")
(n, t, j, k, n1, k1)
sage: foo = sum(t^k*binomial(n-1, k-1), k, 1, n, hold=True)
sage: foo.unhold() # doubleplus ungood...
0
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).unhold() # doubleplus ungood again...
0
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1).unhold() # better
(t + 1)^n1*t
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1).unhold().subs(n1==n-1) # better, and expressed in terms of the original variable
(t + 1)^(n - 1)*t
```

A little numerical check won't hurt :

```sage: foo.subs(n==5).unhold().factor() # Binomial has no problem with numerical values ==> reference
(t + 1)^4*t
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).unhold().subs(n==5) # 0.subs(n==5) is still 0...
0
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1).unhold().subs(n1==n-1).subs(n=5) # Correct...
(t + 1)^4*t
```

As for the extension to other offset values :

```sage: sum(t^(k-j)*binomial(-j + n, -j + k), k, j, n) # Sage doesn't compute this
sum(t^(-j + k)*binomial(-j + n, -j + k), k, j, n)
sage: sum(t^(k-j)*binomial(-j + n, -j + k), k, j, n, hold=True).maxima_methods().changevar(k1==k-j, k1,  k).subs(n==n1+j).unhold().subs(n1==n-j) # same transformation as before
(t + 1)^(-j + n)
sage: sum(t^(k-j)*binomial(-j + n, -j + k), k, j, n, algorithm="giac") # giac concurs
(t + 1)^(-j + n)
sage: sum(t^(k-j)*binomial(-j + n, -j + k), k, j, n, algorithm="mathematica") # so does Mathematica
(t + 1)^(-j + n)
```

Numerical checks left to the reader as an exercise...

It turns out that Sage can (awkwardly but) correctly solve the general case, and gets the special cases wrong. Premature optimization ?

Note : the quantities you have computed are not those sought for by the original author. Re-read your code...

Last edited 4 months ago by charpent (previous) (diff)

comment:5 Changed 4 months ago by gh-DaveWitteMorris

I'm sorry if I wasn't clear. Replacing `n` with `6` in the original sum gives:

```sage: sum(t^(k - 1)*binomial(6 - 1, k - 1), k, 1, 6)
t^5 + 5*t^4 + 10*t^3 + 10*t^2 + 5*t + 1
```

But replacing `n` with `6` in the "equivalent" sum gives:

```sage: sum(t^(k)*binomial(6, k), k, 0, 6-1)
6*t^5 + 15*t^4 + 20*t^3 + 15*t^2 + 6*t + 1
```

This is not the same answer, which means that the two sums are mathematically different. On the other hand, replacing `n` with `6` in the sum I suggested gives:

```sage: sum(t^(k)*binomial(6-1, k), k, 0, 6-1)
t^5 + 5*t^4 + 10*t^3 + 10*t^2 + 5*t + 1
```

This is the same answer as the original sum, and (unless I made a mistake) it will give the same answer for all values of `n`. I thought that was what we wanted, so that is why I said that `n-1` should not be changed to `n`. I apologize if I was not clear about what I meant, or if I misunderstood what you were trying to say.

My other examples were not intended to be directly related to the original question. However, since they are slightly different, and present some additional situations, I thought they might (but only might) help locate the bug. I don't think I need to re-read my code -- it acts exactly as I intended -- but I'm sorry if it wasted your time because I didn't explain my intentions properly.

Last edited 4 months ago by gh-DaveWitteMorris (previous) (diff)

comment:6 Changed 4 months ago by charpent

You :

But replacing n with 6 in the "equivalent" sum gives:

```sage: sum(t^(k)*binomial(6, k), k, 0, 6-1)
6*t^5 + 15*t^4 + 20*t^3 + 15*t^2 + 6*t + 1
```

Nope. Step by step :

```sage: foo
sum(t^k*binomial(n - 1, k - 1), k, 1, n)
sage: foo.maxima_methods().changevar(k1==k-1, k1, k)
sum(t^(k1 + 1)*binomial(n - 1, k1), k1, 0, n - 1)
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1)
sum(t^(k1 + 1)*binomial(n1, k1), k1, 0, n1)
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1).unhold()
(t + 1)^n1*t
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1).unhold().subs(n1==n-1)
(t + 1)^(n - 1)*t
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1).unhold().subs(n1==n-1).subs(n==6)
(t + 1)^5*t
sage: foo.maxima_methods().changevar(k1==k-1, k1, k).subs(n==n1+1).unhold().subs(n1==n-1).subs(n==6).expand()
t^6 + 5*t^5 + 10*t^4 + 10*t^3 + 5*t^2 + t
```

You:

This is not the same answer, which means that the two sums are mathematically different. On the other hand, replacing n with 6 in the sum I suggested gives:

```sage: sum(t^(k)*binomial(6-1, k), k, 0, 6-1)
t^5 + 5*t^4 + 10*t^3 + 10*t^2 + 5*t + 1
```

Nope :

```             +-- By your logic, this should be k-1,
|                           +-- and this should be 6 !
V                           V
sage: sum(t^(k)*binomial(6-1, k), k, 0, 6-1)
t^5 + 5*t^4 + 10*t^3 + 10*t^2 + 5*t + 1
```

By the way :

```sage: sum(t^(k-1)*binomial(6-1, k), k, 0, 6).expand()
t^4 + 5*t^3 + 10*t^2 + 10*t + 1/t + 5
```

I don't think I need to re-read my code

I'm devastated to be of a different opinion...

comment:7 Changed 4 months ago by gh-DaveWitteMorris

I apologize for the confusion that my comments have caused, and that I got us sidetracked from finding (and fixing) bugs in the sage code. Let's focus on finding the problem.

comment:8 Changed 4 months ago by gh-DaveWitteMorris

I also apologize for saying "Contrary to what is said in the description ...". I should have chosen my words more carefully when trying to call attention to a closely related sum that still exhibited the bug. I certainly agree that the second sum in the ticket description can be obtained from the original sum by substitutions, and that sage evaluates it correctly (and that this was a very helpful observation!).

comment:9 follow-up: ↓ 10 Changed 4 months ago by gh-DaveWitteMorris

I think this is a maxima bug, after all:

```sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='maxima')
0
```

I confirmed this by testing directly in maxima:

```(%i8)	load("simplify_sum");
(%o8)	"/Applications/Maxima.app/Contents/Resources/opt/share/maxima/5.43.0/share/solve_rec/simplify_sum.mac"
(%i9)	simplify_sum(sum(t^(k-1)*binomial(n-1, k-1), k, 1, n));
(%o9)	0
```

For comparison:

```sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='maple')
(t + 1)^(n - 1)

sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='giac')
(t + 1)^(n - 1)

sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='sympy')
...
AttributeError: Unable to convert SymPy result (=Piecewise((t*(t + 1)**(n - 1), ((Abs(t) <= 1) & (re(n) - 1 > 0)) | ((Abs(t) < 1) & (re(n) - 1 <= -1)) | (Ne(t, -1) & (Abs(t) <= 1) & (re(n) - 1 <= 0) & (re(n) - 1 > -1))), (Sum(t**k*binomial(n - 1, k - 1), (k, 1, n)), True))/t) into Sage
```

I don't have mathematica, so I didn't test `algorithm='mathematica'`, but it gave the correct result in comment:4, so it must be fine.

comment:10 in reply to: ↑ 9 Changed 4 months ago by charpent

• Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.

You may be right :

I think this is a maxima bug, after all:

```sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='maxima')
0
```

I confirmed this by testing directly in maxima:

```(%i8)	load("simplify_sum");
(%o8)	"/Applications/Maxima.app/Contents/Resources/opt/share/maxima/5.43.0/share/solve_rec/simplify_sum.mac"
(%i9)	simplify_sum(sum(t^(k-1)*binomial(n-1, k-1), k, 1, n));
(%o9)	0
```

By the way, addressing the general case :

```sage: var("t, n, k, j");
sage: sum(t^(k - j)*binomial(n - j, k-j), k, j, n)
sum(t^(-j + k)*binomial(-j + n, -j + k), k, j, n)
sage: maxima_calculus.interact()

--> Switching to Maxima_lib <--

\"/usr/local/sage-9/local/share/maxima/5.44.0/share/solve_rec/simplify_sum.mac\"
maxima_lib: simplify_sum(sum(binomial(n - j, k - j)*t^(k - j), k, j, n));
'sum(binomial(n-j,k-j)*t^(k-j),k,j,n)
maxima_lib: subst([nj = n-j], simplify_sum(subst([n = nj+j],changevar(sum(binomial(n - j, k - j)*t^(k - j), k, j, n), kj = k - j, kj, k))));
(t+1)^(n-j)
maxima_lib:
--> Exiting back to Sage <--
```

For comparison:

```sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='maple')
(t + 1)^(n - 1)

sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='giac')
(t + 1)^(n - 1)

sage: symbolic_sum(t^(k-1)*binomial(n - 1, k - 1), k, 1, n, algorithm='sympy')
...
AttributeError: Unable to convert SymPy result (=Piecewise((t*(t + 1)**(n - 1), ((Abs(t) <= 1) & (re(n) - 1 > 0)) | ((Abs(t) < 1) & (re(n) - 1 <= -1)) | (Ne(t, -1) & (Abs(t) <= 1) & (re(n) - 1 <= 0) & (re(n) - 1 > -1))), (Sum(t**k*binomial(n - 1, k - 1), (k, 1, n)), True))/t) into Sage
```

The latter can be (penibly) generalized as :

```sage: from sympy import sympify, Sum as SSum
sage: %time SSum(*map(sympify, (t^(k - j)*binomial(n - j, k-j), (k, j, n)))).doit().simplify().args[0][0]._sage_()
CPU times: user 13.6 s, sys: 0 ns, total: 13.6 s
Wall time: 13.6 s
(t + 1)^(-j + n)
```

I don't have mathematica,

which isn' that hard to steal ;-)

so I didn't test `algorithm='mathematica'`, but it gave the correct result in comment:4, so it must be fine.

Last edited 4 months ago by charpent (previous) (diff)

comment:11 Changed 4 months ago by charpent

• Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.

This is also Maxima bug 3788.

comment:12 Changed 4 months ago by gh-DaveWitteMorris

Thanks for making the upstream report!

comment:13 Changed 4 months ago by slelievre

• Description modified (diff)