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: |
Description (last modified by )
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
comment:2 Changed 4 months ago by
- Cc gh-DaveWitteMorris slelievre added
- Description modified (diff)
comment:3 Changed 4 months ago by
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
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...
comment:5 Changed 4 months ago by
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.
comment:6 Changed 4 months ago by
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
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
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
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
- 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') 0I 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.
comment:11 Changed 4 months ago by
- Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.
This is also Maxima bug 3788.
comment:13 Changed 4 months ago by
- 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
- Milestone changed from sage-9.4 to sage-9.5
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 ?