######################################################################
##LEON: Save this file as  LEON                                      #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read LEON<Enter>                                   #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Aug./Sept. 2010
with(combinat): 
print(`Created:  Aug./Sept. 2010`):
print(`updated: July 2011.`):
print(`First Posted: Sept. 15, 2010 `):
print(` This is LEON, a Maple package `):
print(`to compute fundamental solutions `):
print(`of PARTIAL DIFFERENCE OPERATORS with constant coefficients`):
print(`and polynomial solutions of partial differential and partial `):
print(`difference (recurrence) equations as well as computing`):
print(`Ehrenpreis' multiplicity varities in the case of a`):
print(`maximially over-determined systems of constant-coefficient PDEs`):
print(`These are implementations, and `):
print(` discrete analogs of some of the seminal results of`):
print(`Professor Rabbi Leon Ehrenpreis, zekher tsadik lebracha. `):
print(`It accompanies Doron Zeilberger's talk: `):
print(`"Leon Ehrenpreis(1930-2010), a truly FUNDAMENTAL mathemtician"`):
print(` delivered at the Rutgers University Experimental Mathematics `):
print(`seiminar on Sept. 16, 2010 (5:00-5:48pm) and viewable from`):
print(`http://www.math.rutgers.edu/~zeilberg/pj.html `):
print(`where one can also find sample input and output files. `):
print(`It also accompanies the short article`):
print(`The Discrete Analog of the Ehrenpreis-Malgrange Theorem`):
print(`submitted to the Ehrenpreis memorial volume edited by`):
print(`Marvin Knopp and other people.`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the main procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(`For a list of the supporting procedures type ezra1(); `):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` ApplyOper, ApplyOperDis, ChekFS, Chop, MV1, Solomon, SV `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: CanBasPol `):
 print(` DAP, FS, FS1, FS2D, , HilSer, HilSerH `):
 print(` MV, PolSols, PolSolsDis `):
  print(` PolSolsH, PolSolsDisH `):
 print(` `):



elif nops([args])=1 and op(1,[args])=ApplyOper then
print(`ApplyOper(Q,D,x,P): applies the differential operator`):
print(`Q(D) to the polynomial P(x) where x is a list of variables`):
print(`and D is the corresponding list of differentiions for example`):
print(`try ApplyOper(D1+D2,[D1,D2],[x1,x2],x1^2+x2^2);`):

elif nops([args])=1 and op(1,[args])=ApplyOperDis then
print(`ApplyOperDis(Q,E,m,P): applies the difference operator`):
print(`Q(E) to the polynomial P(m) where m is a list of discrete variables`):
print(`and E is the corresponding list of shift operators for example`):
print(`try ApplyOperDis(E1+E2,[E1,E2],[m1,m2],m1^2+m2^2);`):

elif nops([args])=1 and op(1,[args])=CanBasPol then
print(`CanBasPol(ope,M,N,m,n,K): Outputs a list of length K+1`):
print(`whose (i+1)-th entry (i=0..K)`):
print(`gives a base for `):
print(`the vector space of polynomial`):
print(`solutions of the the homogeneous partial-difference`):
print(` (alias recurrence)`):
print(`equation with constant coefficient`):
print(`ope(M,N) f(m,n)=0  `):
print(`of degree <=i for i=0..K`):
print(`polynomias. At present either the order in m `):
print(` (i.e. degree of M in ope)`):
print(`or the order in n (i.e. degree of N in ope) have to be 1.`):
print(`For example, try:`):
print(`CanBasPol(1+I*M-I*N-M*N,M,N,m,n,5);`):

elif nops([args])=1 and op(1,[args])=CheckFS then
print(`CheckFS(P,s,x,R): Given a proposed fundamental solution,`):
print(`s, for a Laurent polynomial P, checks whether indeed`):
print(`P*s=1 up to R .`):
print(` For example, try: `):
print(`CheckFS(1-(x+1/x+y+1/y)/4,FS1(1-(x+1/x+y+1/y)/4,[x,y],6),[x,y],6);`):

elif nops([args])=1 and op(1,[args])=Chop then
print(`Chop(f,x,R): given a Laurent polynomial f in the`):
print(`list of variables x and a pos. integer R`):
print(`only retains the terms whose exponents are`):
print(`between -R and R.`):
print(`For example, try`):
print(` Chop((x+1/x+y+1/y)^5,[x,y],2); `):

elif nops([args])=1 and op(1,[args])=DAP then
print(`DAP(m,n,K): The first K elements of the basis for Discrete Analytic `):
print(`polynomias. For example, try:`):
print(`DAP(m,n,5);`):

elif nops([args])=1 and op(1,[args])=FS then
print(`FS(f,x,R): given a Laurent polynomial f in the`):
print(`list of variables x, and a positive integer R,`):
print(`finds the truncation (up to order R) of a SYMMETRIZED`):
print(` reciprocal in the ring of Laurent formal power series. `):
print(`For example, try:`):
print(`FS(x+y,[x,y],R);`):

elif nops([args])=1 and op(1,[args])=FS1 then
print(`FS1(f,x,R): given a Laurent polynomial f in the`):
print(`list of variables x, finds a reciprocal in the`):
print(`ring of formal Laurent series with the ordering induced`):
print(`by x (i.e. the powers of x[1] are >= some finite integer)`):
print(`For example, try:`):
print(`FS1(x+y,[x,y],R);`):

elif nops([args])=1 and op(1,[args])=FS2D then
print(`FS2D(f,x,y,R): 2D display of the Fundamental `):
print(`solution of a linear partial difference equation in the plane`):
print(`whose symbol is f. For example, for the discrete-diagonal`):
print(`(Duffin-style) Laplacian, type:`):
print(`FS2D(4-x-1/x-y-1/y,x,y,10);`):

elif nops([args])=1 and op(1,[args])=HilSer then
print(`HilSer(P,D,K): The first K terms of the Hilbert Series of`):
print(`the space of polynomials annihilated by the operators in`):
print(`the list P, for example, try:`):
print(`HilSer([D1+D2,D1^2+D2^2],[D1,D2],6);`):

elif nops([args])=1 and op(1,[args])=HilSerH then
print(`HilSerH(P,D,K): The first K terms of the Hilbert Series of`):
print(`the space of polynomials annihilated by the operators in`):
print(`the list of HOMOG. diff. operators with polynomial `):
print(`coefficients, P, for example, try:`):
print(`HilSerH([D1+D2,D1^2+D2^2],[D1,D2],6);`):

elif nops([args])=1 and op(1,[args])=MV then
print(`MV(P,x,D): Given a set of polynomials P in the list of variables x`):
print(`finds a basis for the polynomials annihilated by P(D_x).`):
print(`of degree<=K`):
print(`For example try:`):
print(`MV([D1^2-D2,D2^2],[x1,x2],[D1,D2],4);`):


elif nops([args])=1 and op(1,[args])=MV1 then
print(`MV1(P,x,D,K): Given a set of polynomials P in the list of variables x`):
print(`finds a basis for the polynomials annihilated by P(D_x).`):
print(`of degree<=K.`):
print(`For example try:`):
print(`MV1([D1^2-D2,D2^2],[x1,x2],[D1,D2],2);`):

elif nops([args])=1 and op(1,[args])=PolSols then
print(`PolSols(P,x,D,K): Given a set of diff. operators with const. coeffs.`):
print(` P in the list of`):
print(`differentiantions  D in the list of variables x`):
print(`finds a basis for the polynomials annihilated by all p(D_x) in P`):
print(`of degree<=K.`):
print(`For example try:`):
print(`PolSols([D1^2+D2^2,D1+D2],[x1,x2],[D1,D2],2);`):


elif nops([args])=1 and op(1,[args])=PolSolsDisH then
print(`PolSolsDisH(P,m,E,K): Given a set of polynomials P in the list of `):
print(`shift-operators E in the list of discrete variables m`):
print(`finds a basis for the HOMOGENEOUS polynomials annihilated by all P[i](E)`):
print(`of (homog.) degree=K.`):
print(`For example try:`):
print(`PolSolsDisH([E1^2-2*E1+1],[m1],[E1],2);`):

elif nops([args])=1 and op(1,[args])=PolSolsDis then
print(`PolSolsDis(P,m,E,K): Given a set of polynomials P in the list of `):
print(`shift-operators E in the list of discrete variables m`):
print(`finds a basis for the polynomials, in m, annihilated by all P[i](E)`):
print(`of degree<=K.`):
print(`For example try:`):
print(`PolSolsDis([1+I*M-I*N-M*N],[m,n],[M,N],2);`):

elif nops([args])=1 and op(1,[args])=PolSolsH then
print(`PolSolsH(P,x,D,K): Given a HOMOGENEOUS set of differential `):
print(` operators with polynomial coefficients, P in the list of`):
print(`differentiantions  D in the list of variables x`):
print(`finds a basis for the HOMOGENEOUS `):
print(` polynomials annihilated by all p(D_x) in P`):
print(`of degree=K.`):
print(`For example try:`):
print(`PolSolsH([D1^2+D2^2,D1+D2],[x1,x2],[D1,D2],2);`):

elif nops([args])=1 and op(1,[args])=Solomon then
print(`Solomon(n,K,x): A basis for the polynomials in x[1], ..., x[n]`):
print(`of degree K, that are totally harmonic. For example, try:`):
print(`Solomon(4,6):`):

elif nops([args])=1 and op(1,[args])=SV then
print(`SV(n): all the sign-vectors of length n. For example, try:`):
print(`Sv(4);`):


else
print(`There is no ezra for`,args):
fi:
 
end:





#Chop(f,x,R): given a Laurent polynomial f in the
#list of variables x and a pos. integer R
#only retains the terms whose exponents are
#between -R and R
#For example, try
#Chop((x+1/x+y+1/y)^5,[x,y],2);
Chop:=proc(f,x,R) local f1,xn,X,i:
f1:=expand(f):

if not type(x,list) or x=[] then
RETURN(FAIL):
fi:

xn:=x[nops(x)]:
if nops(x)=1 then
RETURN(add(coeff(f1,xn,i)*xn^i,i=-R..R)):
fi:

X:=[op(1..nops(x)-1,x)]:

RETURN(expand(add(Chop(coeff(f1,xn,i),X,R)*xn^i,i=-R..R))):
end:



#FS1(f,x,R): given a Laurent polynomial f in the
#list of variables x, finds a reciprocal in the
#ring of formal Laurent series with the ordering induced
#by x (i.e. the powers of x[1] are >= some finite integer)
#For example, try:
#FS1(x+y,[x,y],R);
FS1:=proc(f,x,R) local n,d,g1,a0,Y,gu,j,i,Z,f1,xn:

if not type(x,list) or x=[] then
RETURN(FAIL):
fi:



n:=nops(x):
xn:=x[n]:

if n=1 then
d:=ldegree(f,xn):
f1:=normal(f/xn^d):
g1:=taylor(1/f1,xn=0,max(R+d+1,1)):
g1:=add(coeff(g1,xn,i)*xn^i,i=0..R+d):
g1:=expand(xn^(-d)*g1):
g1:=Chop(g1,[xn],R):
RETURN(g1):
fi:


d:=ldegree(f,xn):
a0:=coeff(f,xn,d):
f1:=normal(f/xn^d):
g1:=FS1(a0,[op(1..n-1,x)],R) :
Y:=add(coeff(f1,xn,i)*g1*xn^i,i=1..degree(f1,xn)):

gu:=1:
Z:=1:

for j from 1 to R+d do
Z:=expand(Z*Y):
Z:=expand(add(Chop(coeff(Z,xn,i),[op(1..n-1,x)],R)*xn^i,i=0..R+d)):
gu:=gu+(-1)^j*Z:
od:
gu:=Chop(xn^(-d)*g1*gu,x,R):
end:


#CheckFS(P,s,x,R): Given a proposed fundamental solution,
#s, for a Laurent polynomial P, checks whether indeed
#P*s=1 up to R
#For example, try:
#CheckFS(1-(x+1/x+y+1/y)/4,FS1(1-(x+1/x+y+1/y)/4,[x,y],6),[x,y],6);
CheckFS:=proc(P,s,x,R) local d,i,gu,R1:
d:=max(seq(degree(P,x[i]),i=1..nops(x)),
 seq(-ldegree(P,x[i]),i=1..nops(x))):
R1:=max(R-d,0):
gu:=expand(P*s-1):
gu:=Chop(gu,x,R1):
gu:

end:


#SV(n): all the sign-vectors of length n. For example, try:
#Sv(4);
SV:=proc(n) local gu,gu1:
option remember:

if n<0 then
 RETURN({}):
fi:

if n=0 then
 RETURN({[]}):
fi:

gu:=SV(n-1):

{seq([op(gu1),1],gu1 in gu),seq([op(gu1),-1],gu1 in gu)}:

end:


#hofkhi(tau): the inverse of the permutation tau
hofkhi:=proc(tau)
local i,sig1,sig:
 
for i from 1 to nops(tau) do
sig1[tau[i]]:=i:
od:
 
sig:=[]:
 
for i from 1 to nops(tau) do
sig:=[op(sig),sig1[i]]:
od:
sig:
end:


#FSs(f,x,R): given a Laurent polynomial f in the
#set of variables x, finds a symmetrized
#(with respect to the group permutations)
#reciprocal in the
#ring of formal Laurent series with the ordering induced
#by x (i.e. the powers of x[1] are >= some finite integer)
#For example, try:
#FSs(x+y,[x,y],R);
FSs:=proc(f,x,R) local gu,mu1,n,f1,gu1,pi,pi1,i:
n:=nops(x):
mu1:=permute(n):

gu:=0:

for pi in mu1 do
  f1:=subs({seq(x[i]=x[pi[i]],i=1..n)},f):
  gu1:=FS1(f1,x,R):
  pi1:=hofkhi(pi):
  gu1:=subs({seq(x[i]=x[pi1[i]],i=1..n)},gu1):
  gu:=gu+gu1:
od:


expand(gu/n!):
end:


#FS(f,x,R): given a Laurent polynomial f in the
#set of variables x, finds a symmetrized
#(with respect to the group of signed-permutations)
#reciprocal in the
#ring of formal Laurent series with the ordering induced
#by x (i.e. the powers of x[1] are >= some finite integer)
#For example, try:
#FS(x+y,[x,y],R);
FS:=proc(f,x,R) local gu,mu2,s,n,f1,gu1,i:
n:=nops(x):
mu2:=SV(n):

gu:=0:

for s in mu2 do
  f1:=subs({seq(x[i]=x[i]^s[i],i=1..n)},f):
  gu1:=FSs(f1,x,R):
  gu1:=subs({seq(x[i]=x[i]^s[i],i=1..n)},gu1):
  gu:=gu+gu1:
od:


expand(gu/2^n):
end:





 
 
#FS2D(f,x,y,R): 2D display of the Fundamental solution of a linear
#partial difference equation whose symbol is f in the plane
#For example, try:
#FS2D(4-x-1/x-y-1/y,x,y,10);
FS2D:=proc(f,x,y,R) local gu,i,j:
gu:=FS(f,[x,y],R):
[seq([seq(coeff(coeff(gu,x,i),y,j),j=-R..R)],i=-R..R)]:
end:



############GenPOl

#IV1(d,n): all the vectors of non-negative integres of length d whose
#sum is n
IV1:=proc(d,n)
local gu,i,gu1,i1:
if d=0 then
 if n=0 then
   RETURN({[]}):
 else
    RETURN({}):
fi:
fi:
 
gu:={}:

for i from 0 to n do
 gu1:=IV1(d-1,n-i):
 gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}:
od:
gu:
end:

#IV(d,n): all the vectors of non-negative integres of length d whose
#sum is <=n
IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end:


#GenPol(kList,a,deg): The generic polynomial in kList of
#degree deg,  using the indexed variable a, followed by the set
#of coeffs.
GenPol:=proc(kList,a,deg)
local gu,i,i1:
gu:=IV(nops(kList),deg):
convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`),
i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}:
end:

#GenPolH(kList,a,deg): The generic homog. polynomial in kList of
#degree deg,  using the indexed variable a, followed by the set
#of coeffs.
GenPolH:=proc(kList,a,deg)
local gu,i,i1:
gu:=IV1(nops(kList),deg):
convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`),
i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}:
end:

#Extract1(POL,kList): extracts all the coeffs. of a POL in kList
Extract1:=proc(POL,kList) local k1,kList1,i:
k1:=kList[nops(kList)]:
kList1:=[op(1..nops(kList)-1,kList)]:
if nops(kList)=1 then
 RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))} minus {0}):

else
 RETURN({seq(
op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}):
fi:

end:

#ApplyOper(Q,D,x,P): applies the differential operator
#Q(D) to the polynomial P(x) where x is a list of variables
#and D is the correponding list of differentiiosn for example
#try ApplyOper(D1+D2,[D1,D2],[x1,x2],x1^2+x2^2);
ApplyOper:=proc(Q,D,x,P) local i,Qn:
Qn:=expand(Q):

if not type(Qn,`+`) then
RETURN(ApplyOperM(Qn,D,x,P)):
fi:

expand(add(ApplyOperM(op(i,Qn),D,x,P),i=1..nops(Qn))):
end:


#ApplyOperM(Q,D,x,P): applies the monomial differential operator
#Q(D) to the polynomial P(x) where x is a list of variables
#and D is the correponding list of differentiiosn for example
#try ApplyOpeMr(3*D1^2*D2,[D1,D2],[x1,x2],x1^2+x2^2);
ApplyOperM:=proc(Q,D,x,P) local gu,d,Q1,i:
d:=[seq(degree(Q,D[i]),i=1..nops(D))];
Q1:=normal(Q/mul(D[i]^d[i],i=1..nops(d))):
if {seq(diff(Q1,D[i]),i=1..nops(D))}<>{0} then
    RETURN(FAIL):
fi:

gu:=P:

for i from 1 to nops(D) do
if d[i]>0 then
  gu:=diff(gu,x[i]$d[i]):
fi:
od:
Q1*gu:
end:

####End GenPol


#MV1(P,x,D,K): Given a set of polynomials P in the list of variables x
#finds a basis for the polynomials annihilated by P(D_x).
#of degree<=K.
#For example try:
#MV1([D1^2-D2,D2^2],[x1,x2],[D1,D2],2);
MV1:=proc(P,x,D,K) local gu,Q,var,i,kha,var1,a,kha1:
Q:=GenPol(x,a,K):

var:=convert(Q[2],list):
Q:=Q[1]:
gu:={}:

for i from 1 to nops(P) do

gu:= gu union Extract1(ApplyOper(P[i],D,x,Q),x):

od:

var1:=PtorSheli(gu,var):
Q:=subs(var1[1],Q):

kha:=var1[2]:
{seq(coeff(Q,kha1,1),kha1 in kha)}:

end:







  




#MV(P,x,D,MaxK): Given a set of polynomials P in the list of variables x
#finds a basis for the polynomials annihilated by P(D_x).
#of degree<=K
#For example try:
#MV([D1^2-D2,D2^2],[x1,x2],[D1,D2],6);
MV:=proc(P,x,D,MaxK) local gu,P1,i,memad,K,mu:

if expand(subs({seq(D[i]=0,i=1..nops(D))},P))<>[0$nops(P)] then
 print(` origin is not a solution`):
 RETURN(FAIL):
fi:


mu:=Groebner[Solve](P):

if nops(mu)<>1 then
print(`Origin is not the ONLY solution`):
RETURN(FAIL):
fi:

P1:=Groebner[Basis](P,plex(op(D))):

if not  Groebner[IsZeroDimensional](P1) then
 print(`Not Zero-Dimensional`):
 RETURN(FAIL):
fi:



P1:=Groebner[NormalSet](P1,plex(op(D))):

memad:=nops(P1[1]):


for K from 1 to MaxK do 
gu:=MV1(P,x,D,K):

if nops(gu)=memad then
 RETURN(gu):
fi:
od:
print(`MaxK is not big enought, here is an incomplete list`):
gu:
end:





#PolSols(P,x,D,K): Given a set of polynomials P in the list of variables D
#differentiantions  D in the list of variables x
#finds a basis for the polynomials annihilated by all P(D_x).
#of degree<=K.
#For example try:
#PolSols([D1^2+D2^2,D1+D2],[x1,x2],[D1,D2],2);
PolSols:=proc(P,x,D,K) local gu,Q,var,i,kha,var1,a,kha1:
Q:=GenPol(x,a,K):

var:=convert(Q[2],list):
Q:=Q[1]:
gu:={}:

for i from 1 to nops(P) do

gu:= gu union Extract1(ApplyOper(P[i],D,x,Q),x):

od:



var1:=PtorSheli(gu,var):

Q:=subs(var1[1],Q):

for i from 1 to nops(P) do
if ApplyOper(P[i],D,x,Q)<>0 then
print(`Something is wrong`):
 RETURN(FAIL):
fi:
od:

kha:=var1[2]:

{seq(coeff(Q,kha1,1),kha1 in kha)}:

end:


#HilSer(P,D,K): The first K terms of the Hilbert Series of
#the space of polynomials annihilated by the operators in
#the list P, for example, try:
#HilSer([D1+D2,D1^2+D2^2],[D1,D2],6);
HilSer:=proc(P,D,K) local K1,gu,gu1,x,i,gu2,ku:

gu1:=nops(PolSols(P,[seq(x[i],i=1..nops(D))],D,0)):
gu:=[gu1]:

for K1 from 1 to K do
ku:=PolSols(P,[seq(x[i],i=1..nops(D))],D,K1):
if ku=FAIL then
 RETURN(FAIL):
fi:
gu2:=nops(ku):
gu:=[op(gu),gu2-gu1]:
gu1:=gu2:
od:
gu:
end:


ez:=proc(): print(` PtorSheli(Eq,Vars) `): end:

#PtorSheli(Eq,Vars): My solve
PtorSheli:=proc(Eq,Vars) local i,Eq1,xn,n,m,lu,gu,Vars1:

n:=nops(Vars):
Eq1:=Eq minus {0}:
m:=nops(Eq1):
if Eq1={} then
 RETURN({},Vars):
fi:

if Vars=[] and Eq1<>{} then
 RETURN(FAIL):
fi:

xn:=Vars[n]:

Vars1:=[op(1..n-1,Vars)]:

for i from 1 to m while coeff(Eq1[i],xn,1)=0 do od:

if i=m+1 then
 lu:=PtorSheli(Eq1,Vars1):
  if lu=FAIL then
   RETURN(FAIL):
  else
   RETURN(lu[1],[op(lu[2]),xn]):
  fi:
fi:

gu:=-coeff(Eq1[i],xn,0)/coeff(Eq1[i],xn,1):
Eq1:=expand(subs(xn=gu,Eq1)) minus {0}:
 lu:=PtorSheli(Eq1,Vars1):
 if lu=FAIL then
   RETURN(FAIL):
 fi:
 lu[1] union {xn=subs(lu[1],gu)},lu[2]:
end:

#PolSolsH(P,x,D,K): Given a set of polynomials P in the list of 
#differentiantions  D in the list of variables x
#finds a basis for the HOMOGENEOUS polynomials annihilated by all P(D_x).
#of degree=K.
#For example try:
#PolSolsH([D1^2+D2^2,D1+D2],[x1,x2],[D1,D2],2);
PolSolsH:=proc(P,x,D,K) local gu,Q,var,i,kha,var1,a,kha1,t,P0:

P0:=expand(subs({seq(D[i]=t,i=1..nops(D))},P)):
if degree(P0,t)<>ldegree(P0,t) then
 print(`Not homog. diff. opers.`):
 RETURN(FAIL):
fi:

Q:=GenPolH(x,a,K):


var:=convert(Q[2],list):
Q:=Q[1]:
gu:={}:

for i from 1 to nops(P) do

gu:= gu union Extract1(ApplyOper(P[i],D,x,Q),x):

od:

gu:=expand(gu):

var1:=PtorSheli(gu,var):

Q:=expand(subs(var1[1],Q)):

for i from 1 to nops(P) do
if ApplyOper(P[i],D,x,Q)<>0 then
print(`Something is wrong`):
 RETURN(FAIL):
fi:
od:


kha:=var1[2]:


{seq(coeff(Q,kha1,1),kha1 in kha)}:

end:


#HilSerH(P,D,K): The first K terms of the Hilbert Series of
#the space of polynomials annihilated by the HOMOG. operators in
#the list P, for example, try:
#HilSerH([D1+D2,D1^2+D2^2],[D1,D2],6);
HilSerH:=proc(P,D,K) local K1,gu,gu1,x,i,ku,P0,t:

P0:=expand(subs({seq(D[i]=t,i=1..nops(D))},P)):
if degree(P0,t)<>ldegree(P0,t) then
 print(`Not homog. diff. opers.`):
 RETURN(FAIL):
fi:


gu1:=nops(PolSolsH(P,[seq(x[i],i=1..nops(D))],D,0)):
gu:=[gu1]:

for K1 from 1 to K do
ku:=PolSolsH(P,[seq(x[i],i=1..nops(D))],D,K1):
if ku=FAIL then
 RETURN(FAIL):
fi:
gu1:=nops(ku):
gu:=[op(gu),gu1]:
od:
gu:
end:







#Solomon(n,K,x): A basis for the polynomials in x[1], ..., x[n]
#of degree K, that are totally harmonic. For example, try:
#Solomon(4,6):
Solomon:=proc(n,K,x) local D,i,r:
PolSolsH([seq(add(D[i]^r,i=1..n),r=1..n)],[seq(x[i],i=1..n)],
[seq(D[i],i=1..n)],K):
end:





#ApplyOperDis(Q,E,m,P): applies the difference operator
#Q(E) to the polynomial P(m) where m is a list of discrete variables
#and E is the correponding list of shift operators. For example
#try ApplyOperDis(E1+E2,[E1,E2],[m1,m2],m1+m2);
ApplyOperDis:=proc(Q,E,m,P) local i,Qn:

Qn:=expand(Q):
if not type(Qn,`+`) then
RETURN(ApplyOperDisM(Qn,E,m,P)):
fi:

expand(add(ApplyOperDisM(op(i,Qn),E,m,P),i=1..nops(Qn))):
end:


#ApplyOperDisM(Q,E,m,P): applies the monomial difference operator
#Q(E) to the polynomial P(m) where m is a list of discrete variables
#and Eis the correponding list of shift operators for example
#try ApplyOperDisM(3*E1^2*E2,[E1,E2],[m1,m2],m1^2+m2^2);
ApplyOperDisM:=proc(Q,E,m,P) local gu,d,Q1,i:
d:=[seq(degree(Q,E[i]),i=1..nops(E))];
Q1:=normal(Q/mul(E[i]^d[i],i=1..nops(d))):
if {seq(diff(Q1,E[i]),i=1..nops(E))}<>{0} then
    RETURN(FAIL):
fi:

gu:=P:

for i from 1 to nops(E) do
  gu:=expand(subs(m[i]=m[i]+d[i],gu)):
od:
Q1*gu:
end:








#PolSolsDisH(P,m,E,K): Given a set of polynomials P in the list of 
#shift-operators E in the list of discrete variables m
#finds a basis for the HOMOGENEOUS polynomials annihilated by all P[i](E)
#of (homog.) degree=K.
#For example try:
#PolSolsDisH([E1^2-2*E1+1],[m1],[E1],2);
PolSolsDisH:=proc(P,m,E,K) local gu,Q,var,i,kha,var1,a,kha1:

Q:=GenPolH(m,a,K):


var:=convert(Q[2],list):
Q:=Q[1]:
gu:={}:

for i from 1 to nops(P) do

gu:= gu union Extract1(ApplyOperDis(P[i],E,m,Q),m):

od:

gu:=expand(gu):

var1:=PtorSheli(gu,var):

Q:=expand(subs(var1[1],Q)):

for i from 1 to nops(P) do
if ApplyOperDis(P[i],E,m,Q)<>0 then
print(`Something is wrong`):
 RETURN(FAIL):
fi:
od:


kha:=var1[2]:


{seq(coeff(Q,kha1,1),kha1 in kha)}:

end:




#PolSolsDis(P,m,E,K): Given a set of polynomials P in the list of 
#shift-operators E in the list of discrete variables m
#finds a basis for the polynomials annihilated by all P[i](E)
#of (homog.) degree<=K.
#For example try:
#PolSolsDis([E1^2-2*E1+1],[m1],[E1],2);
PolSolsDis:=proc(P,m,E,K) local gu,Q,var,i,kha,var1,a,kha1:

Q:=GenPol(m,a,K):


var:=convert(Q[2],list):
Q:=Q[1]:
gu:={}:

for i from 1 to nops(P) do

gu:= gu union Extract1(ApplyOperDis(P[i],E,m,Q),m):

od:

gu:=expand(gu):

var1:=PtorSheli(gu,var):

Q:=expand(subs(var1[1],Q)):

for i from 1 to nops(P) do
if ApplyOperDis(P[i],E,m,Q)<>0 then
print(`Something is wrong`):
 RETURN(FAIL):
fi:
od:


kha:=var1[2]:


{seq(coeff(Q,kha1,1),kha1 in kha)}:

end:


#DAP(m,n,K): The first K elements of the basis for Discrete Analytic 
#polynomias. For example, try:
#DAP(m,n,5);
DAP:=proc(m,n,K) local z,gu,i:
gu:=z^m*((1+I*z)/(z+I))^n:
gu:=subs(z=z+1,gu):
gu:=taylor(gu,z=0,K+1):
[seq(expand(coeff(gu,z,i)),i=1..K)]:
end:




#CanBasPol(ope,M,N,m,n,K): Outputs a list of length K+1
#whose (i+1)-th entry (i=0..K)
#gives a base for 
#the vector space of polynomial
#solutions of the the homogeneous partial-difference (alias recurrence)
#equation with constant coefficient
#ope(M,N) f(m,n)=0  
#of degree <=i for i=0..K
#polynomias. At present either the order in m (i.e. degree of M in ope)
#or the order in n (i.e. degree of N in ope) have to be 1
#For example, try:
#CanBasPol(1+I*M-I*N-M*N,M,N,m,n,5);
CanBasPol:=proc(ope,M,N,m,n,K) local z,gu,i,ope1,lu,w:
ope1:=numer(normal(ope)):

if degree(ope1,M)>1 and degree(ope1,N)>1 then
 print(`Not yet implemented `):
 RETURN(FAIL):
fi:

if expand(subs({M=1,N=1},ope1))<>0 then
 RETURN({}):
fi:

if degree(ope1,M)=1 then
  lu:=solve(ope1,M):
  lu:=subs(N=w,lu):
  gu:=lu^m*w^n:
gu:=subs(w=w+1,gu):
gu:=taylor(gu,w=0,K+1):
RETURN([seq(expand(coeff(gu,w,i)),i=0..K)]):
fi:

if degree(ope1,N)=1 then
  lu:=solve(ope1,N):
  lu:=subs(M=z,lu):
  gu:=z^m*lu^n:
gu:=subs(z=z+1,gu):
gu:=taylor(gu,z=0,K+1):
RETURN([seq(expand(coeff(gu,z,i)),i=0..K)]):
fi:

end: