rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

simplify_sum on summation of binomials incorrect #2872

Open rtoy opened 2 months ago

rtoy commented 2 months ago

Imported from SourceForge on 2024-07-07 00:11:49 Created by zmth on 2021-09-11 01:41:18 Original: https://sourceforge.net/p/maxima/bugs/3852


(assume_pos:true,
file_output_append:true, 
ratprint:false,
showtime:true,
load(simplify_sum),
simpsum:true,  
load("lrats"),
letrat:true,
ratfac:true,
algebraic:true );   

(factorial_expand:true,
declare([n,a],integer),
assume(n>2,a>0),
t:simplify_sum(sum((-1)^k*binomial(k+a,a)*binomial(k+n,2*(k+a))*4^k,k,0,n)),
disp(["t",t]),
fl:true,
nn:2+3,
kill(a),
for ns:3 thru nn
do for as thru ns+1 
do(sm:sum((-1)^k*binomial(k+as,as)*binomial(k+ns,2*(k+as))*4^k,k,0,ns),
sm1:sum((-1)^k*binomial(k+as,as)*binomial(k+ns,2*(k+as))*4^k,k,0,ns-2*as),
if 2*as-ns>-1
then tg:0 
else(tg:2*as*(2*ns-2*as+1)*ns!/(2*as-1)/(ns-2*as)/(2*as)!^2,
if -2*as<2*as-ns
then tm:prod(2*as- ns-i,i,0,4*as-ns-1) 
else tm:1/prod(- 2*as-i,i,0,ns-4*as-1),
tg:tg*tm),
if sm#sm1 
then(fl:false,
disp(["ns",ns,"as",as,"sm",sm,"sm1",sm1])),
if sm#tg 
then
(disp(["maxima wrong ans,ns",ns,"as",as,"sm",sm,"tg",tg])) ),
fl);

Maxima gives wrong answer for this sum. It asks the question Is n-a+2 pos., neg. or 0 and i ans. positive then what is labled t is the maxima ans. out. You can verify that it is wrong for if n=4,a=2 t is 0 but the sum is actually equal to 1.

Note in the trial runs i used as,ns for a,n.

rtoy commented 2 months ago

Imported from SourceForge on 2024-07-07 00:11:50 Created by robert_dodier on 2021-09-17 05:32:29 Original: https://sourceforge.net/p/maxima/bugs/3852/#4e91


I'm not sure I understand the validation test shown above, but anyway, when I try just

load (simplify_sum);
sum((-1)^k*binomial(k+a,a)*binomial(k+n,2*(k+a))*4^k,k,0,n);
simplify_sum (%);

I get (after answering "positive" to asksign questions about a and n - a)

-((2*((-n)+2*a-1)!*n+(1-2*a)*((-n)+2*a-1)!)*n!)
 /((4*(-2*a)!*a^2-2*(-2*a)!*a)*(2*a-1)!^2)

which is problematic for integer a since then it involves factorials of negative integers. Even for noninteger a, the result seems incorrect; evaluating the summation for a = 5/3, n = 8 yields something different from evaluating the simplify_sum result for those values (116.0677 versus -371.4579, respectively).

verbose_level: 2 before calling simplify_sum shows the result is coming from the Zeilberger method, so the Zeilberger function thinks it knows how to solve it, and it's returning an apparently-invalid result.

rtoy commented 2 months ago

Imported from SourceForge on 2024-07-07 00:11:53 Created by robert_dodier on 2022-01-30 21:54:00 Original: https://sourceforge.net/p/maxima/bugs/3852/#f0f0


Looking at this again. I wasn't able to make any further progress, but anyway here is a simplified form of the original problem statement which might be easier to study for additional analysis.

(load(simplify_sum),
 ratfac:true,
 factorial_expand:true,
 declare([n,a],integer),
 assume(n>2,a>0));

t:simplify_sum(sum((-1)^k*binomial(k+a,a)*binomial(k+n,2*(k+a))*4^k,k,0,n));

for ns:3 thru 5
    do for as thru ns+1 
           do (sm:sum((-1)^k*binomial(k+as,as)*binomial(k+ns,2*(k+as))*4^k,k,0,ns),
               errcatch (subst (['a = as, 'n = ns], t)),
               tg: if %% = [] then false else %%[1],
               if sm=tg 
                   then disp (["maxima apparently correct ans,ns",ns,"as",as,"sm",sm,"tg",tg])
                   else disp (["maxima wrong ans,ns",ns,"as",as,"sm",sm,"tg",tg]));
rtoy commented 2 months ago

Imported from SourceForge on 2024-07-07 00:11:57 Created by robert_dodier on 2022-01-31 08:30:58 Original: https://sourceforge.net/p/maxima/bugs/3852/#99aa


I think I've identified the problem. The summand contains a term binomial(k + n, 2*(k + a)), and somewhere along the way, simplify_sum calls makefact, which turns that into (k + n)!/(2*(k + a))!/(n - k - 2*a)!. For k > n - 2*a, that last term is going to be the factorial of a negative integer, which is undefined. simplify_sum doesn't notice that, and eventually produces an invalid result.

Replacing the upper limit of summation by n - 2*a to ensure the makefact result is valid, and then renaming n - 2*a to m, and applying simplify_sum to that, seems to yield a valid result. (The sum is the same because binomial(n, k) = 0 for k > n, it's just the makefact output which is invalid.)

I think it's the responsibility of simplify_sum to ensure the output from makefact is valid. If so, then simplify_sum has to inspect any binomial and either figure out how to adjust the limits of summation, or refuse to apply makefact.