# sum(y,k,a,b,c) is sum of y index k from a to b step c;
# sum is reverse of backward difference;


sum(a_+b_,x_) := sum(a,x)+sum(b,x);
#sum(a_-b_,x_) := sum(a,x)-sum(b,x);
#sum(a_*b_*c_,x_) := If(isfree(a,x), a*sum(b*c,x));
sum(f(x_)-f(x_+(-1)),x_):=f(x);
sum(f(1+x_)-f(x_),x_):=f(x);

sum(binomial(n_,k_), k_):= 2^k;
sum(sin(x_)^2,x_):= x/2+csc(1)*sin(1-2x)/4;
sum(cos(x_)^2,x_):= x/2-csc(1)*sin(1-2x)/4;
sum(sin(a_*x_)^2,x_):= x/2+csc(a)*sin(a-2a*x)/4;
sum(cos(a_*x_)^2,x_):= x/2-csc(a)*sin(a-2a*x)/4;
sum(sin(x_),x_):= -csc(1/2)*cos(x-1/2)/2;
sum(cos(x_),x_):=cot(1/2)*sin(x)/2-cos(x)/2;
sum(sinh(x_),x_):=csch(1/2)*cosh(x-1/2)/2;
sum(cosh(x_),x_):=coth(1/2)*sinh(x)/2-cosh(x)/2;
sum(sin(a_*x_),x_):= -csc(a/2)*cos(a*x-a/2)/2;
sum(cos(a_*x_),x_):=cot(a/2)*sin(a*x)/2-cos(a*x)/2;
sum(sinh(a_*x_),x_):=csch(a/2)*cosh(a/2-a*x)/2;
sum(cosh(a_*x_),x_):=coth(a/2)*sinh(a*x)/2-cosh(a*x)/2;
sum(log(a_*x_),x_):=If(isfree(a,x), log(a^(x-1)*Gamma(x)));
sum(log(x_),x_) :=logGamma(x+1);
sum(exp(x_),x_):= (exp(x+1)-1)/(exp(1)-1);

sum(a_^k_,k_):=if(a==e, exp(k)/(exp(1)-1), If(isfree(a,k), a^(k)/(a-1) ));
sum(a_^(b_+k_),k_):=if(a==e, exp(b+k)/(exp(1)-1), If(isfree(a,k), a^(b+k)/(a-1) ));
sum(a_^(b_*k_),k_):=If(isfree(a,k), a^(b*k-1)/(a^b-1) );

sum(k_^a_, k_) :=If(isfree(a,k), harmonic(-a,k));
sum((c_+k_)^a_, k_) :=If(isfree(a,c,k), harmonic(-a,c+k));
sum((z_+k_)^a_, k_) :=If(isfree(a,z,k), harmonic(-a,z+k));
sum((c_+d_*k_)^a_, k_) :=If(isfree(a,c,k), d^a*harmonic(-a,c/d+k));

sum(a_*x_^(b_*k_),k_):=replace(sum(a*x^k,k),x,x^b);
sum(a_*x_^(b_+k_),k_):=sum(a*x^k,k)*x^b;
sum((-1)^k_*a_*x_^k_,k_):=replace(sum(a*x^k,k),x,-x);

sum(b_^k_*k_,k_):= b/(1-b)^2*GammaQ(k+1,b);
sum(b_^(-k_)*k_,k_):= b/(1-b)^2*GammaQ(k+1,b);

sum(x_^(2k_)/(2k_+1),k_):=atanh(x)*GammaQ(k+1,x)/x;
sum(x_^(2k_+1)/(2k_+1),k_):=atanh(x)*GammaQ(k+1,x);

sum(b_^k_*k_^p_,k_):= -L(b,-p,1+k)*b^(k+1);
sum(k_^p_*x_^k_,k_):= -L(x,-p,1+k)*x^(k+1);
sum((a_+k_)^p_*x_^k_,k_):= -L(x,-p,1+a+k)*x^(k+1);
sum((z_+k_)^p_*x_^k_,k_):= -L(x,-p,1+z+k)*x^(k+1);
sum((a_+b_*k_)^p_*x_^k_,k_):= -L(x,-p,1+a/b+k)*x^(k+1)*b^p;
sum((m_+b_*k_)^p_*x_^k_,k_):= -L(x,-p,1+m/b+k)*x^(k+1)*b^p;

sum((-1)^k_*k_^a_, k_) :=If(isfree(a,k), eta(-a,1+k));
sum((-1)^k_*(c_+k_)^a_, k_) :=If(isfree(a,c,k), eta(-a,c+k+1));
sum((-1)^k_*x_^(2k_)/(2k_+1),k_):=atan(x)*GammaQ(k+1,x)/x;
sum((-1)^k_*x_^(2k_+1)/(2k_+1),k_):=atan(x)*GammaQ(k+1,x);

sum(k_/(k_+1)^3,k_):= -1+psi(2,2)/2-psi(2,k+2)-psi(1,2+k)+pi^2/6;
sum(k_/(k_+1)^2,k_):= psi(k+2)+psi(1,2+k)+gamma-pi^2/6;
sum(k_/(k_+1),k_):= psi(k+2)-gamma+1;

sum(1/(k_+k_^2),k_):=k/(1+k);
sum(1/k_/(k_+1),k_):=k/(1+k);

sum(k_/k_!,k_) := exp(1)*GammaQ(k);
sum(1/k_!,k_) := exp(1)*GammaQ(k+1);
sum(1/(b_+k_)!,k_) := exp(1)*GammaQ(k+1+b);
sum(1/(a_*k_)!,k_) := -mittag(a,k+2,1);
sum(1/(a_*k_+b_)!,k_) :=  -mittag(a,b+2+k,1);
#sum(1/k_!,k_) := -E(1,2+k,1);

sum(x_^k_/(2k_)!!,k_):= exp(x/2)*GammaQ(k+1,x/2);
sum(x_^k_/k_!,k_):=  exp(x)*GammaQ(k+1);
sum(x_^k_/(b_*k_)!,k_):= if(isfree(b,k), -mittag(b,k+2,x));
sum(x_^k_/(b_+k_)!,k_):= if(isfree(b,k), -mittag(1,k+2+b,x));
sum(x_^k_/(a_*k_+b_)!, k_):=if(isfree(a,b,k), -mittag(a,b+2+k,x));
sum(x_^k_/Gamma(a_*k_), k_):=if(isfree(a,k), -mittag(a,1+k,x));
sum(x_^k_/Gamma(k_+b_), k_):=if(isfree(b,k), -mittag(1,b+1+k,x));
sum(x_^k_/Gamma(a_*k_+b_), k_):=if(isfree(a,b,k), -mittag(a,b+1+k,x));
sum(x_^k_/k_!/k_,k_):= -Ein(-x,k);

sum(x_^(a_*k_)/k_!,k_):=-mittag(1,k+2,x^a);
sum(x_^(a_*k_)/(b_*k_)!,k_):=-mittag(b,k+2,x^a);

#sum((-1)^k*x_^k_/k_!,k_):= -mittag(1,k+2,-x);
#sum((-1)^k*x_^(a_*k_)/k_!,k_):=-mittag(1,k+2,-x^a);
#sum((-1)^k_/k_!/k_*x_^k_,k_):= -Ein(x,k);
#sum((-1)^k_/k_*x_^k_,k_):= -log(1+x)*GammaQ(k,x);
#sum((-1)^k_/k_!*x_^k_,k_):= exp(-x)*GammaQ(k+1,-x);

sum(k_^k_,k_):= k*GammaQ(k-1);
sum(k_^(-k_),k_):= -Sophomore(-1,1)*GammaQ(k-1);
sum((-1)^(-k_)*k_^(-k_),k_):= -Sophomore(1)*GammaQ(k-1);
sum(x_,x_) := x^2/2+x/2;


sum(y_,x_,b_) := sum(y, x,0,b);


sum(y_,j_,a_,b_) := If(isfree(y,j), y*b-y*a+y, If(isnumber(b) and b<>inf, block(sum:=0, For(j,a,b,1,sum:=sum+y), sum) ));
sum(y_,k_,n_,n_):=replace(y,k,n);
sum(a_*b_,k_,a_,b_) := If(isfree(a,k), a*sum(b,k,a,b));
sum(a_*b_*c_,k_,a_,b_) := If(isfree(a,k), a*sum(b*c,k,a,b));
sum(-y_,k_,a_,n_) := -sum(y,k,a,n);
sum(y_+z_,x_,a_,b_) := sum(y,x,a,b)+sum(z,x,a,b);
sum(y_-z_,x_,a_,b_) := sum(y,x,a,b)-sum(z,x,a,b);

sum(x_^k_,k_,a_,b_) := if(hasnot(x,k),(x^(1+b)-x^a)/(x-1) );
#sum(x_^k_/k_!,k_,a_,b_) := exp(x)*(Gamma(b+1,x)/Gamma(b+1)-Gamma(a,x)/Gamma(a));
#sum(1/k_!,k_,a_,n_) := exp(1)*Gamma(n+1,1)/Gamma(n+1)-exp(1)*Gamma(a,1)/Gamma(a);

sum(k_^s_*x_^k_,k_) := -lerch(x,-s,k)*x^k;
sum(k_^s_*x_^k_,k_,a_,b_) := x^a*lerch(x,-s,a)-x^(b+1)*lerch(x,-s,b+1);
sum((c_+k_)^s_*x_^k_,k_,a_,b_) := x^a*lerch(x,-s,a+c)-x^(b+1)*lerch(x,-s,b+c+1);
#sum((c_+k_)^s_*x_^k_,k_,a_,b_) := lerch(x,-s,a+c,b);

#sum((c_+k_)^n_,k_,a_,b_) := zeta(-n,a+c)-zeta(-n,1+b+c);
sum(k_^n_,k_,a_,b_) := -H(-n,a-1)+H(-n,b);
#sum(k_^n_,k_,a_,b_) := H(-n,a,b);
#sum((k_+x_)^n_,k_,a_,b_) := H(-n,a+x,b+x);
#sum((c_+k_)^n_,k_,a_,b_) := H(-n,a+c,b+c);
sum(k_^5,k_,a_,b_) := 1/6*b^6+1/2*b^5+5/12*b^4-1/12*b^2-1/6*a^6+1/2*a^5-5/12*a^4+1/12*a^2;
sum(k_^4,k_,a_,b_) := 1/5*b^5+1/2*b^4+1/3*b^3-1/30*b-a^5/5+a^4/2-a^3/3+a/30;
sum(k_^3,k_,a_,b_):= b^4/4+b^3/2+b^2/4-a^4/4+a^3/2-a^2/4;
sum(k_^2,k_,a_,b_):= b^3/3+b^2/2+b/6-a^3/3+a^2/2-a/6;
sum(k_,k_,a_,b_):= b^2/2+b/2-a^2/2+a/2;

sum((-1)^k_*k_^n_,k_,a_,b_) := (-1)^a*eta(-n,a,b);
sum((-1)^k_*(k_+x_)^n_,k_,0,b_) := eta(-n,x,b);
sum((-1)^k_*(c_+k_)^n_,k_,0,b_) := eta(-n,c,b);

sum(exp(k_*x_),x_,a_,b_) := (exp((1+b)*k)-exp(a*k))/(exp(k)-1);
sum(exp(x_),x_,a_,b_) := (exp(1+b)-exp(a))/(exp(1)-1);
sum(exp(-x_),x_,a_,b_) := (exp(1-a)-exp(-b))/(exp(1)-1);

sum(abs(x_),x_,a_,b_):= b^2/2+abs(b)/2+a^2/2+abs(a)/2;

#sum((-1)^k_*k_^2,k_,1,b_) := (-1)^b*b*(b+1)/2;
sum((2*k_+1)^2,k_,a_,b_) := −4/3a^3+4/3b^3+4b^2+1/3a+11/3b+1;
sum((2*k_+1)^3,k_,a_,b_) := −2a^4+2b^4+8b^3+a^2+11b^2+6b+1;
sum(log(x_),x_,a_,b_) := log(b!/(a-1)!);


sum(binomial(n_, k_), k_, 0, n_):=2^n;
sum(x_^k_/k_!, k_,0,b_):=exp(x)-mittag(1,2+b,x);
sum(k_^r_,k_,0,b_) := if(not(isnumber(b)),harmonic(-r,b)+0^r );
sum(k_*k_!, k_,0,b_) := (b+1)!-1;
sum(k_/(k_+1)!, k_,0,b_) := 1-1/(b+1)!;

sum(y_,k_,a_,inf):=lim(sum(y,k,a,n),n=oo);
sum(k_,k_,a_,inf):= infinity;

sum(integrates(y_,x_,k_),k_,c_,inf):= exp(x)*int(y*exp(-x),x);
sum(integrates(y_,x_,k_,a_,b_),k_,c_,inf):= exp(x)*int(y*exp(-x),x,a,b);

sum(b_^(-k_)*k_,k_,1,inf):= if(abs(b)<=1, inf, if(isfree(b,k),b/(b-1)^2 ));
sum(b_^k_*k_,k_,1,inf):= if(abs(b)>=1, inf, if(isfree(b,k),b/(1-b)^2 ));
sum(b_^k_/k_*x_^k_,k_,1,inf) := -log(1-b*x);
sum(b_^(-k_)/k_*x_^k_,k_,1,inf) := -log(1-x/b);
sum(1/k_*x_^k_,k_,1,inf) := -log(1-x);
sum((-1)^k_/k_,k_,1,inf):= -log(2);
sum(k_^x_,k_,1,inf) := zeta(-x);
sum(1/k_^k_,k_,1,oo):=1.99545595750014;
#sum(x_^k_*k_^b_,k_,1,inf):= polylog(-b,x);
#sum((-1)^k_*k_^x_,k_,1,inf):= -eta(-x);
#sum((-1)^k_*k_^b_*x_^k,k_,1,inf):= polylog(-b,-x);
sum(x_^k_*(a_+k_)^b_,k_,1,inf):= polylog(a,-b,x);
sum(x_^(a_*k_)/k_,k_,1,oo):= -log(1-x^a);
sum(x_^(a_+k_)/k_,k_,1,oo):= -log(1-x)*x^a;
sum(x_^(a_+b_*k_)/k_,k_,1,oo):= -log(1-x^b)*x^a;
sum(x_^k_/k_,k_,1,oo):= -log(1-x);
sum((-1)^(-k)*k_^(-k_),k_,1,oo):= -Sophomore(1);

sum(binomial(n_, k_), k_, 0, Infinity):=2^n;
sum(binomial(n_, k_)*x_^k_, k_, 0, Infinity):=(1+x)^n;
sum((-1)^k*binomial(n_, k_)*x_^k_, k_, 0, Infinity):=(1-x)^n;
sum((k_+x_)^n_,k_,0,oo):=zeta(-n,x);

sum((-1)^k_/(1+k_),k_,0,inf):= log(2);
sum((-1)^k_/(1+k_)*x^(1+k_),k_,0,inf):= log(1+x);
sum(a_^k_*x_^k_/k_!,k_,0,inf):= exp(a*x);
sum(1/k_!*x_^k_,k_,0,inf):= exp(x);
sum(x_^(a_+k_)/k_!,k_,0,inf):= exp(x)*x^a;
sum(x_^(a_+b_*k_)/k_!,k_,0,inf):= exp(x^b)*x^a;
#sum((-1)^k_/k_!*x_^k_,k_,0,inf):= exp(-x);
sum(k_/k_!*x_^k_,k_,0,inf):= x*exp(x);
sum((-1)^k_*k_/k_!*x_^k_,k_,0,inf):= x*exp(-x);
sum(x_^k_/(a_*k_)!, k_,0,inf):=mittag(a,x);
sum(x_^k_/(a_*k_+b_)!, k_,0,inf):=mittag(a,b-1,x);
sum(x_^k_/gamma(a_*k_+b_), k_,0,inf):=mittag(a,b,x);
sum(x_^k_/gamma(k_+b_), k_,0,inf):=mittag(1,b,x);

sum((-1)^k_*x_^(2k_)/(2k_+1)! ,k_,0,inf):= sin(x)/x;
sum((-1)^k_*x_^(2k_+1)/(2k_+1)! ,k_,0,inf):= sin(x);
sum((-1)^k_*x_^(2k_)/(2k_)!,k_,0,inf):= cos(x);
sum((-1)^k_*x_^(2k_+1)/(2k_)!,k_,0,inf):= cos(x)*x;
sum((-1)^k_*x_^(2k_+(-1))/(2k_)!,k_,0,inf):= cos(x)/x;
sum(k_*x^(2k_+1)/((2k_+1)!!),k_,0,inf):= tan(x);
sum((-1)^k_/(2k_+1)! ,k_,0,inf):= sin(1);
sum((-1)^k_/(2k_)!,k_,0,inf):= cos(1);
sum(k_/((2k_+1)!!),k_,0,inf):= tan(1);

sum((2k_+(-1))!!/(2k_)!!/(2k_+1)*x_^(2k_+1),k_,0,inf):=asin(x);
sum((-1)^k_*x_^(2k_+1)/(2k_+1),k_,0,inf):= atan(x);
sum((-1)^k_/(2k_+1),k_,0,inf):= atan(1);

sum(x_^(2k_)/(2k_+1)!,k_,0,inf):= sinh(x)/x;
sum(x_^(2k_+a_)/(2k_+1)!,k_,0,inf):= sinh(x)*x^(a-1);
sum(x_^(2k_)/(2k_)!,k_,0,inf):= cosh(x);
sum(x_^(2k_)/(2k_)!/k_,k_,1,inf):= cosh(x)-1;
sum(x_^(2k_+a_)/(2k_)!,k_,0,inf):= cosh(x)*x^a;

sum((-1)^k_*exp(2k_*x_),k_,0,inf):= -1/2*tanh(x);

sum(1/(2k_+1)!,k_,0,inf):= sinh(1);
sum(1/(2k_)!,k_,0,inf):= cosh(1);

sum((-1)^k_*(2k_-1)!!/((2k_)!!*(2k_+1))*x_^(2k_+1),k_,0,oo):=asinh(x);
sum(x_^(2k_+1)/(2k_+1),k_,0,inf):=atanh(x);

sum(k_*x_^(k_-1),k_,0,inf):= if(abs(x)>=1, inf, if(isfree(x,k),1/(1-x)^2 ));
sum(x_^k_,k_,0,inf):= if(abs(x)>=1, inf, if(isfree(x,k),1/(1-x) ));
sum(x_^(-k_),k_,0,inf):= if(abs(x)<=1, inf, if(isfree(x,k),x/(x-1) ));
sum(1/k_!,k_,0,inf):= exp(1);
sum(k_/k_!,k_,0,inf):= exp(1);
sum(k_/k_!,k_,1,inf):= exp(1);
sum(k_^2/k_!,k_,0,inf):= 2 exp(1);
sum(k_^3/k_!,k_,0,inf):= 5 exp(1);
sum(k_^2/k_!,k_,1,inf):= 2 exp(1);
sum(k_^3/k_!,k_,1,inf):= 5 exp(1);
sum((-1)^k_/k_!, k_,0,inf):= exp(-1);
#sum((a_+k_)^p_*x_^k_,k_,0,oo):=L(x,-p,a);
#sum((1+k_)^p_*x_^k_,k_,0,oo):=1+polylog(-p,x);

#sum(exp(a_*x_),x_,0,inf) := when(a<0, exp(-a)/(exp(-a)-1), oo);
#sum(exp(-k_*x_),k_,0,inf):=1/(exp(x)-1);
#sum(exp(-x_),x_,0,inf) := exp(1)/(exp(1)-1);
#sum(exp(-x_),x_,0,inf) := exp(1)/(exp(1)-1);
sum(exp(-x_)*x_,x_,0,inf) := exp(1)/(exp(1)-1)^2;

sum(1/k_,k_,a_,inf):=oo;
sum(1/(c_+k_),k_,a_,inf):=oo;
sum(1/(b_*k_),k_,a_,inf):=oo;
sum(1/(a_+b_*k_),k_,c_,inf):=oo;
sum(1/k_/(2k_+1),k_,1,n_):= 2-2log(2)+H(n)-H(1/2+n);
sum(1/k_/(-1+2k_),k_,1,n_):= 2log(2)-H(n)+H(-1/2+n);


#sum(y_,j_,a_,b_,c_) := If(b>=a, If(isfree(y,j), y*(b-a+1)/c, block(sum:=0, For(j,a,b,c,sum:=sum+y), sum)));
sum(y_,j_,a_,b_,c_) := If(isfree(y,j), y*(b-a+1)/c, If(isnumber(b) and b<>inf, block(sum:=0, For(j,a,b,c,sum:=sum+y), sum) ));

sum(y_) := sum(y,k);