kill(all)$

Maxima: Cauchysum explained with a proof ,version 1.0
peter.vlasschaert@gmail.com,21/02/2022

load(lrats)$

finite sum :

pa:sum(x^i/i!,i,0,5);

infinite sum :

p1:sum(x^i/i!,i,0,inf);

rem : i ,dummy

simpsum;

1e) step : build , exp(x)
----------------------------

p1a:p1,simpsum;

load(simplify_sum)$

p1b:simplify_sum(p1a);

p1c:p1a;

definition : exp(x)
=============

p1c1:p1b = p1c;

2e) step : build , exp(y)
----------------------------

p1d:subst([x=y,i=j],p1a);

rem : j ,dummy

p1e:simplify_sum(p1d);

p1e1:p1e = p1d;

3e)step :We make product : double sum both are to infinity

p2:p1c1*p1e1;

p2a:lhs(p2)=subst([reverse(p1c1),reverse(p1e1)],rhs(p2));

p2a1:part(lhs(p2),2);

4e)step : construct infinite sum lhs(p2)

part a:

p3:subst([x=a],p1c1);

part b:

p3a:subst([a=p2a1],p3);

part c:(y+x)^i, Binomium van Newton

ref:https://en.wikipedia.org/wiki/Binomial_theorem (English)
     https://nl.wikipedia.org/wiki/Binomium_van_Newton
ref :https://maxima.sourceforge.io/docs/manual/maxima_50.html

use :Function: binomial (x, y)

p4:(p2a1)^i = sum(binomial (i, l)*x^(i-l)*y^l,l,0,i);

l:dummy

p4a:binomial (i, l);

p4b1:makefact(binomial (i, l));

p4b2:makegamma(binomial (i, l));

rem : choose between makefact,makegamma (advanced : mathematics)
ref :https://en.wikipedia.org/wiki/Gamma_function

p4c:p4a=p4b1;

use :p4c into p4

p4d:subst([p4c],p4);

p4d1:lhs(p4d)=intosum(rhs(p4d));

use : p2 & p3a & p4d1 or p4d 

p4e1:rhs(p2) = subst([p4d1],rhs(p3a));

double sum both to infiniy → sum to infinity & finite sum.
===========================================

p4e2:rhs(p2) = ratsimp(subst([p4d],rhs(p3a)));

set : c[ii]=x^ii/ii! , d[jj]=y^jj/jj!

look to rhs "p4e2"

p5a: c[ii]=x^ii/ii!;
p5b: d[jj]=y^jj/jj!;

lhs:change : ii =i,jj =j,{p5a,p5b}

p5a1:reverse(subst( [ii =i,jj =j],p5a));
p5b1:reverse(subst( [ii =i,jj =j],p5b));

rhs:change : ii =i-l,jj =l,{p5a,p5b}

p5a2:reverse(subst([ii =i-l,jj =l],p5a));
p5b2:reverse(subst([ii =i-l,jj =l],p5b));

final general result
==============

p6:subst([p5a1,p5b1],lhs(p4e2))=lratsubst([p5a2,p5b2],rhs(p4e2));

How to do this Maxima : with cauchysum
===============================

1e) step : default values (=check)
-------------------------------

p7a1:sumexpand;
p7a2:cauchysum;

2e) step :make product of double infinite sum

s1:sum(f1[i],i,0,inf);/*i : dummy*/
s2:sum(f2[j],j,0,inf);

sd:s1*s2;

3e) step : change default values
---------------------------------------

p7a2:sumexpand:true;
p7b2:cauchysum:true;

4e) step : make cauchysum
---------------------------------------

p7a31:''sd;

5e) step : see above proof i1=i,i2=l
-------------------------------------------

p7a32:subst([i1=i,i2=l],p7a31);

6e) step : see above proof f2=c,f1=d
-------------------------------------------

p7a33:subst([f2=c,f1=d],p7a32);

7e) step : check against rhs(p6)
-------------------------------------------

is(p7a33=rhs(p6));