Opened 4 months ago

Last modified 2 months ago

#31844 new defect

Incorrect summation with binomial of non-atomic arguments.

Reported by: charpent Owned by:
Priority: major Milestone: sage-9.5
Component: symbolics Keywords: summation, binomial, maxima
Cc: gh-DaveWitteMorris, slelievre Merged in:
Authors: Reviewers:
Report Upstream: Reported upstream. No feedback yet. Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by slelievre)

From Ask Sage question 57195: Incorrect symbolic sum with binomial coefficients.

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:

Change History (14)

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

  • Cc gh-DaveWitteMorris slelievre added
  • 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: 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 :

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

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

maxima_lib: load("simplify_sum");
\"/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

  • Keywords maxima added

Thanks for making the upstream report!

comment:13 Changed 4 months ago by slelievre

  • Description modified (diff)

Replying to gh-DaveWitteMorris:

Contrary to what is said in the description, the "equivalent sum obtained by shifting the summation index" is not correct.

Apologies for an earlier misleading edit to the ticket description. If I understand correctly, the sum earlier described as an "equivalent sum obtained by shifting the summation index" is indeed correctly computed, but not equivalent to the buggy one. Description amended now.

comment:14 Changed 2 months ago by mkoeppe

  • Milestone changed from sage-9.4 to sage-9.5
Note: See TracTickets for help on using tickets.