#shifts the index by N
dshift:=proc(E,N,n)
local newE;
newE:=subs(n=n+N,E);
end:

# The procedure varder calculates the kth-order
# variational derivative of a functional E with respect to u.
dvard:=proc(E,k,N,n,u)
local L,m;
L:=diff(add(binomial(m,k)*dshift(E,-m,n),m=k..N),u);
end:

#expands discrete operators
dexpand:=proc(poly,D)
local s,i;
s:=expand(poly);
s:=[op(s)];
s:=[seq([degree(s[i],D),coeffs(s[i])],i=1..nops(s))];
end:

#computes the flux for vecdhom
vecdflux:=proc(E,N,n,u)
local s,D,f,J,j,m,x;
for j from 0 to N-1 do
s:=dexpand((x-1)^j,x);
f[j]:=u*dvard(E,j+1,N,n,u);
D[j]:=add(op(2,s[m])*dshift(f[j],op(1,s[m]),n),m=1..nops(s));
od;
return(op(expand([seq(D[i],i=0..N-1)])));
end:

#modified version of the discrete homotopy operator, a list of terms is returned
vecdhom:=proc(E,N,n,u,a0)
local sh,result,newans,ans,h,newM,Min,MAX,m,a,lambda,r,M,s,F,i;
M:=nops(u);#print(E,N,n,u);
#s is the integrand, to maximize the number of terms, we take this as a vector and a sum
s:=[seq(vecdflux(E,N,n,u[r]),r=1..M)];
s:=subs({seq(seq(op(0,u[i])[n+m]=lambda*op(0,u[i])[n+m],m=-N..N),i=1..nops(u))},s);
sh:=normal(convert(s,`+`));
ans:=expand({seq(int(expand(s[i]/lambda),lambda),i=1..nops(s)),int(expand(sh/lambda),lambda)});
#remove expressions that were not integrated
for s in ans do
 if has(s,int) then
  ans:=ans minus {s}
 fi;
od;
#get as many terms as possible
newans:=[NULL];
for i from 1 to nops(ans) do
 if whattype(ans[i])=`+` then
  newans:=[op(newans),op(ans[i])];
 else
  newans:=[op(newans),ans[i]];
 fi;
od;
#normalize the solution and remove constants
result:={NULL};
for i from 1 to nops(newans) do
 if has(newans[i],n) then
  result:={op(result),SPLIT(newans[i],n)[1]};
 fi;
od;
[op(subs(lambda=1,result))];
end:

#picks out undefined functions in an expression
whatfun:=proc(E,n)
local i,s,sid;
s:={NULL};
sid:=[op(indets(E,`indexed`))];
for i from 1 to nops(sid) do
 if not(type(op(0,sid[i]),`procedure`))
  then s:={op(s),op(0,sid[i])[n]};
 fi;
od;
[op(s)];
end:

#picks out undefined functions and their shifts in an expression
whatfun2:=proc(E,n)
local i,s,sid;
s:={NULL};
sid:=[op(indets(E,`indexed`))];
for i from 1 to nops(sid) do
 if not(type(op(0,sid[i]),`procedure`))
  then s:={op(s),sid[i]};
 fi;
od;
[op(s)];
end:

#maxshift determines the maximum shift in T
maxshift:=proc(T,n)
local ulist;
ulist:=whatfun2(T,n);
ulist:={seq(op(ulist[i])-n,i=1..nops(ulist))};
max(op(ulist));
end:

#minshift determines the minimum shift in T
minshift:=proc(T,n)
local ulist;
ulist:=whatfun2(T,n);
ulist:={seq(op(ulist[i])-n,i=1..nops(ulist))};
min(op(ulist));
end:

#orders a set from large to small
orderset:=proc(s)
    local ss,leftover,k;
    ss:=NULL;
    leftover:=s;
    for k to nops(s) do
        ss:=ss,max(op(leftover));
        leftover:=leftover minus {max(op(leftover))}
    od;
    [ss]
end:

#orders a set according to number of products in a term
torderset:=proc(E)
local w,fin,n,s;
s:={seq(nops(op(1,E[n]))+(n/(nops(E)+1)),n=1..nops(E))};
s:=orderset(s);
w:={seq(nops(op(1,E[n]))+(n/(nops(E)+1))=E[n],n=1..nops(E))};
fin:=subs(w,s);
end:

#orders a list of terms according to the sign of the lowest shift, negative first
negorder:=proc(L,u,n)
local s,neglist,poslist;
neglist:=[NULL];
poslist:=[NULL];
for s in L do
 if minshift(s,n,u)<0 then
  neglist:=[op(neglist),s];
 else
  poslist:=[op(poslist),s];
 fi;
od;
op(neglist),op(poslist);
end:

#orders a set according to the absolute value of the minimum shift, from higest to lowest,
shiftorderlist:=proc(L,u,n)
local finlist,j,diffset,dlist,i,MAX,difflist,ord;
difflist:=[seq(max(abs(minshift(L[i],n,u)),abs(maxshift(L[i],n,u))),i=1..nops(L))];
diffset:=orderset({op(difflist)});
finlist:=[NULL];
for i in diffset do
 dlist:=[NULL];
 for j from 1 to nops(L) do
  if i=difflist[j] then
   dlist:=[op(dlist),L[j]];
  fi;
 od;
finlist:=[op(finlist),dlist];
od;
[seq(op(finlist[i]),i=1..nops(finlist))];
end:

#orders a set according to the difference between the highest and lowest shifts, from higest to lowest.
#if terms are of the same order, then it orders a according to the absolute value of the minimum shift, from higest to lowest,
#with negative shifts first
difforderlist:=proc(L,u,n)
local finlist,j,diffset,dlist,i,MAX,difflist,ord;
difflist:=[seq(maxshift(L[i],n,u)-minshift(L[i],n,u),i=1..nops(L))];
diffset:=orderset({op(difflist)});
finlist:=[NULL];
for i in diffset do
 dlist:=[NULL];
 for j from 1 to nops(L) do
  if i=difflist[j] then
   dlist:=[op(dlist),L[j]];
  fi;
 od;
finlist:=[op(finlist),dlist];
od;
finlist:=[seq(shiftorderlist(finlist[i],u,n),i=1..nops(finlist))];
[seq(negorder(finlist[i],u,n),i=1..nops(finlist))];
end:

#converts a list of functional expressions to strings, orders them according to lexorder, and returns the funtional expressions
funlexorder:=proc(LIST)
local Newlist,i,subrules,s;
subrules:=[seq(convert(LIST[i],`string`)=LIST[i],i=1..nops(LIST))];
Newlist:=[seq(op(1,subrules[i]),i=1..nops(LIST))];
s:=sort(Newlist,`lexorder`);
subs(subrules,s);
end:

#RRED row-reduces A while preserving the order of the rows
#returns the first maxsstep-1 linearly independent rows

RRED:=proc(A,sstep,B,maxsstep)
local step,i,j,newA,flag1,rowdims,newB;

rowdims:=linalg[rowdim](A);

#if rowdims=maxsstep-1 then
 #return(linalg[gausselim](A));
#fi;

if sstep=maxsstep then
 return(B);
fi;

step:=sstep;
newA:=A;
flag1:=0;
for i from 1 to rowdims do
  if newA[i,step]<>0 and flag1=0 then
   flag1:=i;
  fi;
od;
if flag1<>0 then
 newA:=linalg[pivot](newA,flag1,step);
 if step=1 then
  newB:=linalg[matrix]([linalg[row](newA,flag1)]);
 else
  newB:=linalg[stackmatrix](B,linalg[row](newA,flag1));
 fi;
 if rowdims=1 then
  return(newB);
 else
  newA:=linalg[delrows](newA,flag1..flag1);
  RRED(newA,step+1,newB,maxsstep);
 fi;
else
 RRED(newA,step+1,newB,maxsstep);
fi;
end:

#splits the product T into an x-dependent part and a constant part
SPLIT:=proc(T,x)
local xpart,coeffpart,s,opT,xset,coeffset;
xset:={NULL};
coeffset:={NULL};
if whattype(T)=`*` then
 opT:=[op(T)];
else
 opT:=[T];
fi;
for s in opT do
 if has(s,x) then
  xset:={op(xset),s};
 else
  coeffset:={op(coeffset),s};
 fi;
od;
xpart:=convert(xset,`*`);
coeffpart:=convert(coeffset,`*`);
return([xpart,coeffpart]);
end:

#substitutes ssubrules into EE
SUBS:=proc(ssubrules,EE,x)
local E,Elist,subrules;
subrules:=torderset(ssubrules);
#E:=expand(EE);
E:=EE;
if not(whattype(E)=`+`) then
 Elist:=[E];
else
 Elist:=convert(E,`list`);
fi;
Elist:=[seq(SPLIT(Elist[i],x),i=1..nops(Elist))];
Elist:=subs(subrules,Elist);
E:=add(convert(Elist[i],`*`),i=1..nops(Elist));
end:

#sums the expression EE with respect to n, where EE only contains positive shifts
SUM:=proc(EE,n)
local i,ssublist,E,param,Elist,u,N,functionset,nset,k,sumnset,Fnlist,Fn1list,sublist,subrules,Esub1,Fn1listsub,Fnlistsub,R,syst,m,A,Ar,eqnset,slnset,sln,nogood,sumnogood,fresult;

E:=expand(eval(EE));

#convert E to a list
if not(whattype(E)=`+`) then
 Elist:=[E];
else
 Elist:=convert(E,`list`);
fi;

#determine the unknown functions
u:=whatfun(EE,n);
if u=[NULL] then
 return(dshift(sum(EE,n),1,n))
fi;

#maxshift in E
N:=maxshift(E,n,u);

#remove terms which have no functional dependence to be summed directly
functionset:=[seq(seq(dshift(u[r],j,n),j=0..N),r=1..nops(u))];
nset:={NULL};
for k in Elist do
if not(has(k,functionset)) then
 nset:={op(nset),k};
fi;
od;
nset:=expand(convert(nset,`+`));
sumnset:=sum(nset,n);
E:=expand(eval(E-nset));
if not(whattype(E)=`+`) then
 Elist:=[E];
else
 Elist:=convert(E,`list`);
fi;

#Fnlist is a list of terms from the homotopy operator applied to E
Fnlist:=expand(vecdhom(E,N,n,u));#print(Fnlist);
if Fnlist=[NULL] then
 return(combine(sum(E,n)+sum(nset,n)));
fi;
Fn1list:=expand([seq(dshift(Fnlist[i],1,n),i=1..nops(Fnlist))]);

#form substitution rules
ssublist:=expand({op(Elist),op(Fn1list),op(Fnlist)});
sublist:={NULL};
for i in ssublist do
 if not(whattype(i)=`+`) then
  sublist:={op(sublist),i};
 else
  sublist:={op(sublist),op(i)};
 fi;
od;
#remove constant coefficients
sublist:=[op({seq(SPLIT(sublist[i],n)[1],i=1..nops(sublist))})];
#order sublist according to the difference between the highest and lowest shifts, from higest to lowest.
#if terms are of the same order, then it orders according to the absolute value of the minimum shift, from higest to lowest,
#with negative shifts first
sublist:=difforderlist(sublist,u,n);

#perform the substitutions
subrules:=(seq(sublist[i]=b[i],i=1..nops(sublist)));
#order subrules according to number of products in each term
subrules:=torderset({subrules});
Esub1:=SUBS(subrules,E,n);#print(ESUB,Esub1);
Fn1listsub:=[seq(SUBS(subrules,Fn1list[i],n),i=1..nops(Fn1list))];
Fnlistsub:=[seq(SUBS(subrules,Fnlist[i],n),i=1..nops(Fnlist))];

#R is what remains after summing E
#R is optimized
R:=expand(Esub1-sum(a[iii]*(Fn1listsub[iii]-Fnlistsub[iii]),iii=1..nops(Fn1list)));
R:=collect(R,{seq(b[i],i=1..nops(sublist))});
#form the linear system to be solved
syst:={NULL};
for m from 1 to nops(sublist) do
 if has(coeff(R,b[m]),[seq(a[j],j=1..nops(Fn1list))]) then
  syst:=[op(syst),coeff(R,b[m])=0];
 fi;
od;

#A is the overdetermined system for the coefficients, a[n]
#Ar is the first nops(a) linearly independent rows of A
A:=linalg[genmatrix](syst,[seq(a[j],j=1..nops(Fn1list))],flag);#print(AAAA,A);
Ar:=RRED(A,1,A,nops(Fn1list)+1);#print(ArAr,Ar);

slnset:=linalg[backsub](Ar,`false`,`v`);
sln:=seq(a[m]=slnset[m],m=1..nops(Fn1list));#print(solution,sln);
sln:=eval(sln,v=0);
assign(sln);

#all terms which were not summed using homotopy
nogood:=expand(E-add(a[i]*(Fn1list[i]-Fnlist[i]),i=1..nops(Fn1list)));#print(nogood);
sumnogood:=sum(nogood,n);

fresult:=add(a[j]*Fn1list[j],j=1..nops(Fn1list))+combine(sumnogood)+dshift(sumnset,1,n);
unassign(a);

#check if correct
if SUMcheck(EE,fresult,n)<>0 then
 print("No summation done.");
 return(sum(E,n)+dshift(sumnset,1,n))
else
 return(fresult);
fi;
end:

#ouputs the expression G modulo total differences
DifferenceMod:=proc(G,n)
local i,sumpart,Glist;
Glist:=[op(SUM(G,n))];
sumpart:=0;
for i from 1 to nops(Glist) do
 if has(Glist[i],`sum`) then sumpart:=op(1,Glist[i]);
 fi;
od;
sumpart;
end:

#checks if SUM is correct
SUMcheck:=proc(E,SG,n)
local i,sumpart,diffpart,Glist;
if not(whattype(SG)=`+`) then
 Glist:=[SG];
else
 Glist:=[op(SG)];
fi;
sumpart:=0;
diffpart:=SG;
for i from 1 to nops(Glist) do
 if has(Glist[i],`sum`) then
  sumpart:=op(1,Glist[i]);
 fi;
od;
diffpart:=expand(diffpart-sum(sumpart,n));
normal(expand(E-diffpart+dshift(diffpart,-1,n)-sumpart));
end: