######################################################################
##EPDE.txt: Save this file as  EPDE.txt                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `EPDE.txt`<Enter>                             #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################


#Created: 
 
print(`Created:  March 2021`):
print(` This is EPDE.txt `):
print(`It is one of the packages that accompany the article `):
print(`The Enumerative Theory of Partial Differential Equations`):
print(`by AJ Bu and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):


print(`---------------------------------------`):
 print(`For a list of the Generalized procedures type ezraG();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the Bi-variate procedures type ezraB();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the Bi-variate procedures tbe New Way type ezraBN();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

#_EnvAllSolutions := true: 

with(combinat):


ezraG:=proc()

if args=NULL then
 print(` The Generalized procedures are: CompsG, GenPolG, KernG, TotHg, TotHgDim `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: ApplyMonoDO, CheckHC, CheckHC1, Comps, Extract1, Exctract1B, FindBase`):   
 print(`  Hnmd1d2s, GenPol ,  GFup, GuessPol, GuessPolE, IsLC, Kama, Surv, SurvT, VecToHS  `):
 print(``):

else
ezra(args):
fi:

end:


ezraB:=proc()

if args=NULL then
 print(` The bivariate procedures are:  ApplyDOB, ApplyMonoDOB, HPn, HPnm, GenPolB, GenPolBg, KernB, KernBg, TotHB, TotHBdim, TotHBnm,`):
print(`TotHBnmDim, TotHB1, TotHB1br , TotHB1brDim , TotHB1g  `):

 print(``):

else
ezra(args):
fi:

end:



ezraBN:=proc()

if args=NULL then
 print(` The bivariate procedures the New way: DimKernSeqTF,  Hnmd1d2, Hnmd1d2K,  TotHBnmN,  , TotHB1gN, Ynm, YnmF  `):

 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: A, ApplyDO  DimKernSeq,  DimKernSeqF, DimsH, DimsHfast, GuessPolU, HndK, HndKLC, HSS, Kern, SeqH, TotH, TotHdim, Vecs, VecsT `):
 print(` `):


elif nops([args])=1 and op(1,[args])=A then
print(`A(n,d,K): the dimension of the vector space of totally harmonic polynomials in n variables of total degree d and`):
print(`degree in x[n]<=K. It uses the recurrence implied by the lemma. Try:`):
print(`n1:=4:[seq(A(n1,d,d),d=0..binomial(n1,2))]; `):

elif nops([args])=1 and op(1,[args])=ApplyDO then
print(`ApplyDO(Ope,x,D1,n,P): inputs a differential operator (with constant coefficients) in the form`):
print(`Sum(c*D1[1]^a1*...*D1[n]^an) ,c a constants and a1, .., an, non-negative integers and applies it the polynomial P in x[1], ..., x[n]`):
print(`Try:`):
print(`ApplyDO(D1[1]+D1[2],x,D1,2,x[1]^3*x[2]^6);`):


elif nops([args])=1 and op(1,[args])=ApplyDOB then
print(`ApplyDOB(Ope,x,y,D1,D2,n,m,P): inputs a differential operator (with constant coefficients) in the form`):
print(`Sum(c*D1[1]^a1*...*D1[n]^an*D2[1]^b1*...*D2[m]^bm) ,c a constants `):
print(`and a1, .., an, b1, ..., bm non-negative integers and applies it the polynomial P in x[1], ..., x[n], y[1], ..., y[m]`):
print(` Try: `):
print(`ApplyDOB(D1[1]*D2[1]+D1[2]*D2[2],x,y,D1,D2,2,2,x[1]^4*y[1]^4+x[2]^3*y[2]^2);`):

elif nops([args])=1 and op(1,[args])=ApplyMonoDO then
print(`ApplyMonoDO(M,x,D1,n,P): inputs a monomial differential operator (with constant coefficients) in the form`):
print(`c*D1[1]^a1*...*D1[n]^an ,c a constant and a1, .., an, non-negative integers and applies it the polynomial P in x[1], ..., x[n]`):
print(`Try: `):
print(`ApplyMonoDO(3*D1[1]^2*D1[2]^3,x,D1,2,x[1]^3*x[2]^6);`):


elif nops([args])=1 and op(1,[args])=ApplyMonoDOB then
print(`ApplyMonoDOB(M,x,y,D1,D2,n,m,P): inputs a monomial differential operator (with constant coefficients) in the form`):
print(`c*D1[1]^a1*...*D1[n]^an*D2[1]^b1*...*D2[m]^bm ,c a constant and a1, .., an; b1, ..., bm `):
print(`non-negative integers and applies it the polynomial P in x[1], ..., x[n]`):
print(`Try: `):
print(`ApplyMonoDOB(3*D1[1]^2*D1[2]^3*D2[1]^2*D2[2],x,y,D1,D2,2,2,x[1]^3*x[2]^6*y[1]^2*y[2]^3);`):



elif nops([args])=1 and op(1,[args])=CheckHC then
print(`CheckHC(n,d,K): verifies the central lemma that H(n,d,K)/H(n,d,K-1) is the zero space if d>=K>=n, and otherwise`):
print(`H(n-1,d-K,d-K). Try:  `):
print(`CheckHC(4); `):

elif nops([args])=1 and op(1,[args])=CheckHC1 then
print(`CheckHC1(n,d,K): computes dim(HndK(x,n,d,K))- dim(HndK(x,n,d,K-1))-dim(HndK(x,n-1,d-K,d-K). Try:`):
print(`CheckHC1(4,1,1); `):

elif nops([args])=1 and op(1,[args])=Comps then
print(`Comps(n,d): the set of compositions of d into n parts. Try:`):
print(`Comps(4,2);`):


elif nops([args])=1 and op(1,[args])=CompsG then
print(`CompsG(n,d,L): the set of compositions of d into n parts, with upper bound given by the list L. Try:`):
print(`CompsG(4,3,[1,1,1,1]);`):


elif nops([args])=1 and op(1,[args])=DimKernSeq then
print(`DimKernSeq(n,d,D1,S): inputs a pos. integer n, a pos. integer,d, a symbol D1, and a set S of linear differential operators with constant coefficients`):
print(`phrased in terms of D1, where D1[i] is d/dx[i] ,`):
print(`output the first d+1 terms, starting at i=0 of the vector space of polynomials in x[1], ..., x[n] annihilated by the member of S. Try:`):
print(`DimKernSeq(3,4,D1,{seq(add(D1[i]^s,i=1..3),s=1..3)});`):


elif nops([args])=1 and op(1,[args])=DimKernSeqF then
print(`DimKernSeqF(n,d,D1,S): Like DimKernSeq(n,d,D1,S), but using Groebner bases and survivirors. Try: `):
print(`DimKernSeqF(3,4,D1,{seq(add(D1[i]^s,i=1..3),s=1..3)});`):

elif nops([args])=1 and op(1,[args])=DimKernSeqTF then
print(`DimKernSeqTF(m,n,d1,d2,D1,D2,S): The list of listsL such that L[d1][d2] is the dimension of the st of polynomials of b-degree [d1,d2]`):
print(`annihilated by S. Try`):
print(`DimKernSeqTF(3,3,1,1,D1,D2,{seq(seq(add(D1[i]^r*D2[i]^s,i=1..3),r=1..3),s=1..3)});`):

elif nops([args])=1 and op(1,[args])=DimsH then
print(`DimsH(n,R,N): inputs a  pos. integer n,  a set of positive integers R, and a positive integer N`):
print(`outputs the sequence of integers of length N+1 whose i-th entry is the dimesnion of the vector space`):
print(`of homogeneous polynomials of degree i in n variables annihilated by the set of differential operators`):
print(`D1^r+...+Dn^r for r in R. For example, try:`):
print(`DimsH(3,{1,2,3},3);`):


elif nops([args])=1 and op(1,[args])=DimsHfast then
print(`DimsHfast(n,R,N): a faster version of DimsH(n,R,N) (q.v.), using Groebner bases. Try:`):
print(`DimsHfast(5,{1,2,3},6);`):

elif nops([args])=1 and op(1,[args])=Extract1 then
print(`Extract1(f,x,n): the set of coefficients of the polynomial f in the variables x[1], ... x[n]. Try:`):
print(`Extract1(a*x[1]^2+b*x[1]*x[2]*c*x[2]^5,x,2);`):

elif nops([args])=1 and op(1,[args])=Extract1B then
print(`Extract1B(f,x,y,n,m): the set of coefficients of the polynomial f in the variables x[1], ... x[n], y[1], ..., y[m]. Try:`):
print(`Extract1B(a*x[1]^2+b*x[1]*x[2]*c*x[2]^5*y[1]*y[2],x,y,2,2);`):

elif nops([args])=1 and op(1,[args])=FindBase then
print(`FindBase(S,x,n): Given a set S of plolynomials in x[1], ..., x[n], outputs a base for the vector space it generates . Try:`):
print(`FindBase({x[1]+x[2],x[2]+x[3],x[1]+2*x[2]+x[3]},x,3);`):

elif nops([args])=1 and op(1,[args])=HndK then
print(`HndK(x,n,d,K): a basis for the vector space of totally harmonic polynomials in the n variables`):
print(`x[1], ..., x[n], of total degree d in x[1], ..., x[n] and degree in x[n]<=K. Try:`):
print(`HndK(x,4,2,1); `):

elif nops([args])=1 and op(1,[args])=HndKLC then
print(`HndKLC(x,n,d,K): a basis for the vector space of LEADING COEFFICIENT in x[n] for totally harmonic polynomials in the n variables`):
print(`x[1], ..., x[n], of total degree d in x[1], ..., x[n] and degree in x[n]<=K. Try:`):
print(`HndKLC(x,4,2,1);`):


elif nops([args])=1 and op(1,[args])=HSS then
print(`HSS(n1,S,K,q): inputs a positive integer n, a subset S, of {1, ..., n1} a positive integer K, and a variable q,`):
print(`outputs the Hilbert series of the vector space of polynomials in (x[1], ..., x[n1]) annihilated by`):
print(`Sum(D1[i]^s,i=1..n1) for s in S. Using K+1 data points for guessing. Try:`):
print(`HSS(3,{1,2,3},12,q);`):

elif nops([args])=1 and op(1,[args])=GenPol then
print(`GenPol(n,d,x,a): inputs a positive integer n, a positive integer d, a symbol x and a symbol a`):
print(`outputs a generic homog. polynomial in the n variable x[1],...,x[n] of degree d with coefficients of the form a[i]`):
print(`followed by the set of coefficients. Try:`):
print(`GenPol(4,2,x,a);`):

elif nops([args])=1 and op(1,[args])=GenPolG then
print(`GenPolG(n,d,x,a,L): inputs a positive integer n, a positive integer d, a symbol x and a symbol a and a list of length L`):
print(`outputs a generic homog. polynomial in the n variable x[1],...,x[n] of degree d with coefficients and`):
print(`degree sequence <=L , of the form a[i]`):
print(`followed by the set of coefficients. Try:`):
print(`GenPolG(4,2,x,a,[1,1,1,1]);`):

elif nops([args])=1 and op(1,[args])=GenPolB then
print(`GenPolB(n,m,d1,d2,x,y,a): inputs  positive integers n and m , and positive integers d1 and d2, and  symbols x and y and a symbol a`):
print(`outputs a generic bi- homog. polynomial in the n variable x[1],...,x[n]  and m variables y[1], ..., y[m]`):
print(`of degree d1 on the x's and degree d2 in the y's  with coefficients of the form a[i]`):
print(`followed by the set of coefficients. Try:`):
print(`GenPolB(3,3,1,1,x,y,a);`):


elif nops([args])=1 and op(1,[args])=GenPolBg then
print(`GenPolBg(n,m,d1,d2,x,y,a,L1,L2): inputs  positive integers n and m , and positive integers d1 and d2, and  symbols x and y and a symbol a`):
print(`outputs a generic bi- homog. polynomial in the n variable x[1],...,x[n]  and m variables y[1], ..., y[m]`):
print(`of degree d1 on the x's and degree d2 in the y's  and the degree sequence of the x variables is <=L1 and the`):
print(`degree sequence of the y-variables is <=L2`):
print(`with coefficients of the form a[i]`):
print(`followed by the set of coefficients. Try:`):
print(`GenPolBg(3,3,1,1,x,y,a,[0,0,2],[0,0,2]);`):

elif nops([args])=1 and op(1,[args])=GFup then
print(`GFup(L,n,q): Given an ultimate polynomial L=[INI,Pol(n)] in terms of the variable n, and a variable name q, outputs the`):
print(`rational function, in q, that is its  generating function. Try:`):
print(`GFup([[1,1],n^2],n,q); `):

elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(L,n): inputs a list of integers of values starting at n=0, tries to fit a polynomial to the data, or returns FAIL. Try`):
print(`GuessPol([0,1,4,9],n);`):

elif nops([args])=1 and op(1,[args])=GuessPolU then
print(`GuessPolU(L,n): inputs a list of integers of values starting at n=0, tries to fit an ultimate polynomial to the data, or returns FAIL. Try`):
print(`GuessPolU([0,1,4,9,16,25],n);`):

elif nops([args])=1 and op(1,[args])=GuessPolE then
print(`GuessEPol(L,n): inputs a list of integers of values starting at n=0, tries to fit a polynomial to the data, `):
print(`that is eventually a polynomial or returns FAIL. Try`):
print(`GuessPolE([1,1,4,9,16,25],n);`):


elif nops([args])=1 and op(1,[args])=Hnmd1d2 then
print(`Hnmd1d2(x,y,n,m,d1,d2): Like TotHB1(x,y,n,m,d1,d2) but the new way. `):
print(`inputs  symbols x and y and  pos. integers n and m , `):
print(`and integers d1 and d2, outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]`):
print(`and degree d2 in y[1], ..., y[m]. Try:`):
print(`Hnmd1d2(x,y,3,3,1,1);`):


elif nops([args])=1 and op(1,[args])=Hnmd1d2K then
print(`Hnmd1d2K(x,y,n,m,d1,d2,K): inputs  symbols x and y and  pos. integers n and m , `):
print(`and integers d1 and d2, and a non-neg. integer K <=d2,`):
print(`outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]`):
print(`and degree d2 in y[1], ..., y[m] and deg_y[m]<=K. Try:`):
print(`Hnmd1d2K(x,y,3,3,1,2,1);`):

elif nops([args])=1 and op(1,[args])=Hnmd1d2s then
print(`Hnmd1d2s(x,y,n,m,d1,d2): Like Hnmd1d2(x,y,n,m,d1,d2) but the smart way. Try:`):
print(`Hnmd1d2s(x,y,3,3,1,1);`):


elif nops([args])=1 and op(1,[args])=HPn then
print(`HPn(n,p,q): the Hilbert polynomial of diagonally-harmonic polynomials with (n,n) varibles of each kind in terms of p and q. Try:`):
print(`HPn(3,p,q);`):

elif nops([args])=1 and op(1,[args])=HPnm then
print(`HPnm(n,m,p,q): the Hilbert polynomial of diagonally-harmonic polynomials with (n,m) varibles of each kind in terms of p and q. Try:`):
print(`HPnm(3,3,p,q);`):

elif nops([args])=1 and op(1,[args])=IsLC then
print(`IsLC(p,S,x,n): Is the polynomial p in x[1], ..., x[n] a linear combination of the members of S? Try: `):
print(`IsLC (x[1]+x[2],{x[1],x[2]},x,2);`):

elif nops([args])=1 and op(1,[args])=Kama then
print(`Kama(L): The sum of the entries of the list of lists L. Try:`):
print(`Kama([[3,4],[5,3,[1,4],2] ]);`):

elif nops([args])=1 and op(1,[args])=Kern then
print(`Kern(n,d,x,D1, S): inputs pos. integers n and d, a symbol x, and a set S of differential operators with constant coefficients`):
print(`where D1[i] denotes differentiation with respect to x[i], outputs a basis for the`):
print(`set of polynomials in x[1], ..., x[n] of degree d, annihilated by the differential operators in S. Try: `):
print(`Kern(3,3,x,D1,{D1[1]+D1[2]+D1[3],D1[1]^2+D1[2]^2+D1[3]^2,D1[1]^3+D1[2]^3+D1[3]^3 });`):
print(`Kern(4,6,x,D1,{D1[1]+D1[2]+D1[3]+D1[4],D1[1]^2+D1[2]^2+D1[3]^2+D1[4]^2,D1[1]^3+D1[2]^3+D1[3]^3+D1[4]^3,D1[1]^4+D1[2]^4+D1[3]^4+D1[4]^4 });`):

elif nops([args])=1 and op(1,[args])=KernB then
print(`KernB(n,m,d1,d2,x,y,D1, D2, S): inputs pos. integers n and m, non-negative integers`):
print(`d1,d2, symbols x and y, and a set S of differential operators with constant coefficients`):
print(`where D1[i] denotes differentiation with respect to x[i] and D2[i] denotes differentiation`):
print(`with repspect to y[i], outputs a basis for the`):
print(`set of polynomials in (x[1], ..., x[n];y[1],..,y[m]) of bi-degree (d1,d2), annihilated by the differential operators in S. Try: `):
print(`KernB(3,3,2,1,x,y,D1,D2,{D1[1]*D2[1]+D1[2]*D2[2]+D1[3]*D2[3], D1[1]^2*D2[1]^2+D1[2]^2*D2[2]^2+D1[3]^2*D2[3]^2});`):

elif nops([args])=1 and op(1,[args])=KernBg then
print(`KernBg(n,m,d1,d2,x,y,D1, D2, S,L1,L2): inputs pos. integers n and m, non-negative integers`):
print(`d1,d2, symbols x and y, and a set S of differential operators with constant coefficients`):
print(`where D1[i] denotes differentiation with respect to x[i] and D2[i] denotes differentiation`):
print(` with repspect to y[i]. It also inputs lists of integers L1,L2, or length n and m  respectively. `):
print(`It outputs a basis for the`):
print(`set of polynomials in (x[1], ..., x[n];y[1],..,y[m]) of bi-degree (d1,d2),`):
print(`and degree sequence in x<=L1  and degree sequence in y<=L2`):
print(` annihilated by the differential operators in S. Try: `):
print(`KernBg(3,3,2,1,x,y,D1,D2,{D1[1]*D2[1]+D1[2]*D2[2]+D1[3]*D2[3], D1[1]^2*D2[1]^2+D1[2]^2*D2[2]^2+D1[3]^2*D2[3]^2},[3,2,2],[3,2,2]);`):

elif nops([args])=1 and op(1,[args])=KernG then
print(`KernG(n,d,x,D1, S,L): inputs pos. integers n and d, a symbol x, and a set S of differential operators with constant coefficients`):
print(`where D1[i] denotes differentiation with respect to x[i]. It also outputs a list L of length n for upper bounds for the degree sequence`):
print(` outputs a basis for the`):
print(`set of polynomials in x[1], ..., x[n] of degree d and degree sequence <=L, annihilated by the differential operators in S. Try`):
print(`KernG(3,2,x,D1,{D1[1]+D1[2]+D1[3],D1[1]^2+D1[2]^2+D1[3]^2},[1,1,1]);`):
print(`KernG(4,3,x,D1,{D1[1]+D1[2]+D1[3]+D1[4],D1[1]^2+D1[2]^2+D1[3]^2+D1[4]^2,D1[1]^3+D1[2]^3+D1[3]^3+D1[4]^3,D1[1]^4+D1[2]^4+D1[3]^4+D1[4]^4 },[6,6,6,1]);`):

elif nops([args])=1 and op(1,[args])=SeqH then
print(`SeqH(n,x,R,N): inputs a pos. integer n,  a set of positive integers R, and a positive integer N`):
print(`outputs the sequence of integers of length N+1 whose i-th entry is the  sequence of bases for the vector spaces`):
print(`of homogeneous polynomials of degree i in n variables annihilated by the set of differential operators`):
print(`D1^r+...+Dn^r for r in R. For example, try:`):
print(`SeqH(3,x,{1,2,3},3);`):


elif nops([args])=1 and op(1,[args])=Surv then
print(`Surv(x,n,S,d): all the vectors of non-neg. integers of length n, that are never >=then the the vectors in`):
print(`Vecs(x,n,S);. Try:`):
print(`Surv(x,2,{x[1]+x[2],x[1]^2+x[2]^2},1);`):

elif nops([args])=1 and op(1,[args])=SurvT then
print(`SurvT(x,y,m,n,S,d1,d2): all the pair of vectors [v1,v2] of non-neg. integers of length m and n respectively, `):
print(`that are never >=then the the vectors in`):
print(`VecsT(x,y,m,n,S,d1,d2);. Try:`):
print(`SurvT(x,y,2,2,{x[1]+x[2],x[1]^2+x[2]^2,x[1]*y[1]+x[2]*y[2],x[1]^2*y[1]+x[2]^2*y[2]},1,1);`):

elif nops([args])=1 and op(1,[args])=TotH then
print(`TotH(x,n): inputs a symbol x and a pos. integer n, outputs the list of length binomial(n,2)+1 whose`):
print(`(i+1)-th component is a basis for the vector space of totally harmonic homog. polynomials in x[1], .., x[n]`):
print(`of degree i. Try:`):
print(`TotH(x,3); `):

elif nops([args])=1 and op(1,[args])=TotHdim then
print(`TotHdim(n): the dimension sequence of totally harmomic functions. Try:`):
print(`TotHdim(3);`):


elif nops([args])=1 and op(1,[args])=TotHg then
print(`TotHg(x,n): inputs a symbol x and a pos. integer n, outputs the list of length binomial(n,2)+1 whose`):
print(`(i+1)-th component is a  list if the bases for the increasing vector spaces of totally harmonic homog. polynomials in x[1], .., x[n]`):
print(`of degree <=j in x[n] . Try:`):
print(`TotHg(x,3);`):

elif nops([args])=1 and op(1,[args])=TotHgDim then
print(`TotHgDim(n): the broken-up dimension sequence of totally harmomic functions. Try:`):
print(`TotHgDim(3);`):

elif nops([args])=1 and op(1,[args])=TotHBnm then
print(`TotHBnm(x,y,n,m): inputs  symbols x and y and  a pos. integers  n and m, such that n>=m. Ooutputs`):
print(`a list of lists L such that L[d1+1][d2+1] is a basis for the vector space of`):
print(`homog. polynomials in x[1], ..., x[n], y[1], ..., y[m] if bidegree (d1, d2). Try`):
print(`TotHBnm(x,y,3,3);`):


elif nops([args])=1 and op(1,[args])=TotHBnmN then
print(`TotHBnmN(x,y,n,m):  Like TotHBnm(x,y,n,m), but the new way. `):
print(`inputs  symbols x and y and  a pos. integers  n and m, such that n>=m. Ooutputs`):
print(`a list of lists L such that L[d1+1][d2+1] is a basis for the vector space of`):
print(`homog. polynomials in x[1], ..., x[n], y[1], ..., y[m] if bidegree (d1, d2). Try`) :
print(`TotHBnmN(x,y,3,3);`):


elif nops([args])=1 and op(1,[args])=TotHBnmDim then
print(`TotHBnmDim(n,m): the list of lists of dimensions of diagonal harmonics of (n,m) bi-variables.`):
print(`If the output is L then L[d1+1][d2+1] is the dimension of the vector space of`):
print(`totally harmonic bi-polynomials of bi-degree (d1,d2). Try:`):
print(`TotHBnmDim(3,3);`):


elif nops([args])=1 and op(1,[args])=Ynm then
print(`Ynm(n,m):  Like TotHBnmDim(n,m) but the New way`):
print(`the list of lists of dimensions of diagonal harmonics of (n,m) bi-variables.`):
print(`If the output is L then L[d1+1][d2+1] is the dimension of the vector space of`):
print(`totally harmonic bi-polynomials of bi-degree (d1,d2). Try:`):
print(`Ynm(3,3);`):


elif nops([args])=1 and op(1,[args])=YnmF then
print(`YnmF(n,m):  A fast version of Ynm(n,m). Try: `):
print(`YnmF(3,3);`):

elif nops([args])=1 and op(1,[args])=TotHB then
print(`TotHB(x,y,n): inputs  symbols x and y and  a pos. integer  n. Ooutputs`):
print(`a list of lists L such that L[d1+1][d2+1] is a basis for the vector space of`):
print(`homog. polynomials in x[1], ..., x[n], y[1], ..., y[n] if bidegree (d1, d2). Try`):
print(`TotHB(x,y,3);`):


elif nops([args])=1 and op(1,[args])=TotHBdim then
print(`TotHBdim(n): the list of lists of dimensions of diagonal harmonics of n bi-variables.`):
print(`If the output is L then L[d1+1][d2+1] is the dimension of the vector space of`):
print(`totally harmonic bi-polynomials of bi-degree (d1,d2). Try:`):
print(`TotHBdim(3);`):

elif nops([args])=1 and op(1,[args])=TotHB1 then
print(`TotHB1(x,y,n,m,d1,d2): inputs  symbols x and y and  pos. integers n and m , `):
print(`and integers d1 and d2, outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]`):
print(`and degree d2 in y[1], ..., y[m]. Try:`):
print(`TotHB1(x,y,3,3,1,1);`):



elif nops([args])=1 and op(1,[args])=TotHB1br then
print(`TotHB1br(x,y,n,m,d1,d2): Like TotHB1(x,y,n,m,d1,d2) but  outputs a list of lists`):
print(`call it L, such that L[i+1][j+1] is a basic for the subset of`):
print(` TotHB1(x,y,n,m,d1,d2) with degree in x[n]<=i and degree in y[m]<=j. Try:`):
print(`TotHB1br(x,y,3,3,1,1);`):


elif nops([args])=1 and op(1,[args])=TotHB1brDim then
print(`TotHB1brDim(n,m,d1,d2): The dimensios of TotHB1br(x,y,m,n,d1,d2). Try:`):
print(`TotHB1brDim(3,3,1,1);`):

elif nops([args])=1 and op(1,[args])=TotHB1g then
print(`TotHB1g(x,y,n,m,d1,d2,L1,L2): inputs  symbols x and y and  pos. integers n and m , `):
print(`and integers d1 and d2. It also inputs lists L1 and L2 of length m and n respectively.`):
print(`It outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]`):
print(`and degree d2 in y[1], ..., y[m] and x-degree sequence <=L1 and y-degree sequence <=L2. Try:`):
print(`TotHB1g(x,y,3,3,2,2,[3,3,1],[3,3,1]);`):


elif nops([args])=1 and op(1,[args])=TotHB1gN then
print(`TotHB1gN(x,y,n,m,d1,d2,L1,L2): Like TotHB1g(x,y,n,m,d1,d2,L1,L2) but the new way. Try:`):
print(`TotHB1gN(x,y,3,3,2,2,[3,3,1],[3,3,1]);`):

elif nops([args])=1 and op(1,[args])=Vecs then
print(`Vecs(x,n,S), inputs a variable name x, a positive integer n, and a set of polynomials S, in x[1], ..., x[n], outputs`):
print(`the set of exponents of the leading monomial of the Grobner basis of the ideal generated by S with pure lex order x[1], ..., x[n].`):
print(`Try: `):
print(`Vecs(x,2,{x[1]+x[2]});`):


elif nops([args])=1 and op(1,[args])=VecsT then
print(`VecsT(x,y,n,S), inputs  variable names x, and y, a positive integer n, and a set of polynomials S, in x[1], ..., x[n], y[1], ..., y[n], outputs`):
print(`the set of exponents of the leading monomial of the Grobner basis of the ideal generated by S with pure lex order x[1],y[1] ..., x[n], y[n]`):
print(`Try: `):
print(`VecsT(x,y,2,{x[1]+x[2],y[1]+y[2],x[1]*y[1]+x[2]*y[2]});`):

elif nops([args])=1 and op(1,[args])=VecToHS then
print(`VecToHS(S,q,n): Given a set of vectors S in N^n outputs the generating function for the lattice points that are NOT >= and of the vectors of S`):
print(`Try: `):
print(`VecToHS({[1,0],[0,2]},q,2);`):


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


Kama:=proc(L):convert(ListTools[Flatten](L),`+`):end:

#######Start GuessPol
#GuessPol1(L,n,d): inputs a list of integers
#and a pos. integer d, outputs a polynomial of degree d f such that f(i)=L[i+1]
GuessPol1:=proc(L,n,d) local a,i,f,var,eq:

if nops(L)<d+3 then
 print(`List too short`):
 RETURN(FAIL):
fi:

var:={ seq(a[i],i=0..d)}:
f:=add(a[i]*n^i,i=0..d):

eq:={seq( subs(n=i, f)=L[i+1], i=0..d+1)}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

f:=subs(var,f):

if nops(L)>d+3 then
if {seq( subs(n=i, f)-L[i+1], i=d+2..nops(L)-1)}<>{0} then
# print( {seq( subs(n=i, f)-L[i+1], i=d+2..nops(L)-1)}):
# print(f, `did not work out`):
 RETURN(FAIL):
fi:
fi:

f:

end:

#GuessPol(L,n): inputs a list of integers of values starting at n=0, tries to fit a polynomial to the data, or returns FAIL. Try
#GuessPol([0,1,4,9],n);
GuessPol:=proc(L,n) local f,d:

for d from 0 to nops(L)-3 do
 f:=GuessPol1([op(1..d+3,L)] ,n,d):
  if f<>FAIL then
   f:=GuessPol1(L,n,d):
     if f<>FAIL then
       RETURN(factor(f)):
     fi:
   fi:
od:
FAIL:
end:


#GuessPolU(L,n): inputs a list of integers of values starting at n=0, tries to fit an ultimate polynomial to the data, or returns FAIL. Try
#GuessPolU([0,1,4,9],n);
GuessPolU:=proc(L,n) local f1,d,L1,i:

for i from 1 to nops(L)-4 do
 L1:=[op(i..nops(L),L)]:
 f1:=GuessPol(L1,n):
 if f1<>FAIL then
  RETURN([[op(1..i-1,L)],factor(subs(n=n-i+1,f1))]):
 fi:
od:
FAIL:
end:


#GuessEPol(L,n): inputs a list of integers of values starting at n=0, tries to fit a polynomial to the data, 
#that is eventually a polynomial or returns FAIL. Try
#GuessPolE([1,1,4,9,16,25],n);
GuessPolE:=proc(L,n) local i,L1,gu:

for i from 1 to nops(L)-3 do
 L1:=[op(i..nops(L),L)]:
 gu:=GuessPol(L1,n):
 if gu<>FAIL then
  RETURN([[op(1..i-1,L)], factor(expand(subs(n=n-i+1,gu))) ]):
 fi:
od:
FAIL:
end:

#######end GuessPol

#Comps(n,d): the set of compositions of d into n parts. Try:
#Comps(4,2);
Comps:=proc(n,d) local gu,i,lu,lu1:
option remember:

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

gu:={}:

for i from 0 to d do
 lu:=Comps(n-1,d-i):
 gu:=gu union {seq([i,op(lu1)],lu1 in lu)}:
od:

gu:

end:


#GenPol(n,d,x,a): inputs a positive integer n, a positive integer d, a symbol x and a symbol a
#outputs a generic homog. polynomial in the n variable x[1],...,x[n] of degree d with coefficients of the form a[i]
#followed by the set of coefficients. Try:
#GenPol(4,2,x,a);
GenPol:=proc(n,d,x,a) local gu,gu1,P,mu,co,i1:
gu:=Comps(n,d):
P:=0:
mu:={}:
co:=0:
 for gu1 in gu do
 co:=co+1:
  P:=P+a[co]*mul(x[i1]^gu1[i1],i1=1..n):
  mu:=mu union {a[co]}:
 od:
P,mu:
end:


#MyDiff(f,x,i):the i-th derivative of f w.r.t the variable x, including if i=0. Try:
#MyDiff(x^3,x,3);
MyDiff:=proc(f,x,i) :
if i=0 then
 f:
else
expand(diff(f,x$i)):
fi:
end:



#ApplyMonoDO(M,x,D1,n,P): inputs a monomial differential operator (with constant coefficients) in the form
#c*D1[1]^a1*...*D1[n]^an ,c a constant and a1, .., an, non-negative integers and applies it the polynomial P in x[1], ..., x[n]
#Try:
#ApplyMonoDO(3*D1[1]^2*D1[2]^3,x,D1,2,x[1]^3*x[2]^6);
ApplyMonoDO:=proc(M,x,D1,n,P) local c,gu,d,i:

d:=[seq(degree(M,D1[i]),i=1..n)]:
c:=normal(M/mul(D1[i]^d[i],i=1..n)):

if not type(c,numeric) then
  RETURN(FAIL):
fi:

gu:=P:

for i from 1 to n do
 gu:=MyDiff(gu,x[i],d[i]):
od:
expand(c*gu):
end:


#ApplyDO(Ope,x,D1,n,P): inputs a differential operator (with constant coefficients) in the form
#Sum(c*D1[1]^a1*...*D1[n]^an) ,c a constants and a1, .., an, non-negative integers and applies it the polynomial P in x[1], ..., x[n]
#Try:
#ApplyDO(D1[1]+D1[2],x,D1,2,x[1]^3*x[2]^6);
ApplyDO:=proc(Ope,x,D1,n,P) local i:

if not type(Ope,`+`) then
ApplyMonoDO(Ope,x,D1,n,P):
else
expand( add(ApplyMonoDO(op(i,Ope),x,D1,n,P)  , i=1..nops(Ope))):
fi:

end:

#Extract11(f,x): the set of coefficients of the polynomial f in the variable x. Try:
#Extract11(a+2*b*x^3,x);
Extract11:=proc(f,x) local i:
 {seq(coeff(f,x,i),i=0..degree(f,x))} minus {0}:
end:


#Extract1(f,x,n): the set of coefficients of the polynomial f in the variables x[1], ... x[n]. Try:
#Extract1(a*x[1]^2+b*x[1]*x[2]*c*x[2]^5,x,2);
Extract1:=proc(f,x,n) local gu,gu1:
if n=1 then
 RETURN(Extract11(f,x[1])):
fi:

gu:=Extract1(f,x,n-1):

{seq(op(Extract11(gu1,x[n])),gu1 in gu)}:
end:


#Kern(n,d,x,D1, S): inputs pos. integers n and d, a symbol x, and a set S of differential operators with constant coefficients
#where D1[i] denotes differentiation with respect to x[i], outputs a basis for the
#set of polynomials in x[1], ..., x[n] of degree d, annihilated by the differential operators in S. Try
#Kern(3,2,x,D1,{D1[1]+D1[2]+D1[3],D1[1]^2+D1[2]^2+D1[3]^2});
Kern:=proc(n,d,x,D1,S) local gu,P,eq,var,a,P1,s,var1,var11,fr,fr1,gu1:
option remember:

gu:=GenPol(n,d,x,a) :

P:=gu[1]:
var:=gu[2]:

eq:={}:

for s in S do
P1:=ApplyDO(s,x,D1,n,P):
 eq:=eq union Extract1(P1,x,n):
od:

var1:=solve(eq,var):

	
fr:={}:

for var11 in var1 do
 if op(1,var11)=op(2,var11) then
   fr:=fr union {op(1,var11)}:
fi:
od:

P:=expand(subs(var1,P)):

gu:={seq(coeff(P,fr1,1),fr1 in fr)}:

for gu1 in gu do
 if {seq(ApplyDO(s,x,D1,n,gu1),s in S)}<>{0} then
  print(`Soemthing is wrong`):
  RETURN(FAIL):
 fi:
od:
gu:
end:

#TotH(x,n): inputs a symbol x and a pos. integer n, outputs the list of length binomial(n,2)+1 whose
#(i+1)-th component is a basis for the vector space of totally harmonic homog. polynomials in x[1], .., x[n]
#of degree i. Try:
#TotH(x,3):
TotH:=proc(x,n) local S,D1,i,r:
option remember:
S:={seq(add(D1[i]^r,i=1..n),r=1..n)}:
[seq(Kern(n,i,x,D1,S),i=0..binomial(n,2))]:
end:




#DimsH(n,R,N): inputs a pos. integer n,  a set of positive integers R, and a positive integer N
#outputs the sequence of integers of length N+1 whose (i+1)-th entry is the dimesnion of the vector space
#of homogeneous polynomials of degree i in n variables annihilated by the set of differential operators
#D1^r+...+Dn^r for r in R. For example, try:
#DimsH(3,{1,2,3},3);
DimsH:=proc(n,R,N) local gu,x:
option remember:
gu:=SeqH(n,x,R,N):
[seq(nops(gu[i]),i=1..nops(gu))]:
end:


#DimsHfast(n,R,N): inputs a pos. integer n,  a set of positive integers R, and a positive integer N
#outputs the sequence of integers of length N+1 whose (i+1)-th entry is the dimesnion of the vector space
#of homogeneous polynomials of degree i in n variables annihilated by the set of differential operators
#D1^r+...+Dn^r for r in R. For example, try:
#DimsHfast(3,{1,2,3},3);
DimsHfast:=proc(n,R,N) local gu,x,S1,D1:
option remember:
S1:={seq(add(D1[i]^r,i=1..n), r in R)}:
DimKernSeqF(n,N,D1,S1):
end:


#SeqH(n,x,R,N): inputs a pos. integer n,  a set of positive integers R, and a positive integer N
#outputs the sequence of integers of length N+1 whose i-th entry is the  sequence of  vector spaces
#of homogeneous polynomials of degree i in n variables annihilated by the set of differential operators
#D1^r+...+Dn^r for r in R. For example, try:
#SeqH(3,x,{1,2,3},3);
SeqH:=proc(n,x,R,N) local S,D1,i,r:
option remember:
S:={seq(add(D1[i]^r,i=1..n),r in R)}:
[seq(Kern(n,i,x,D1,S),i=0..N)]:
end:




###Start Bi-variate



#GenPolB(n,m,d1,d2,x,y,a): inputs  positive integers n and m , and positive integers d1 and d2, and  symbols x and y and a symbol a
#outputs a generic bi- homog. polynomial in the n variable x[1],...,x[n]  and m variables y[1], ..., y[m]
#of degree d1 on the x's and degree d2 in the y's  with coefficients of the form a[i]
#followed by the set of coefficients. Try:
#GenPolB(3,3,1,1,x,y,a);
GenPolB:=proc(n,m,d1,d2,x,y,a) local guX,guY, guX1,guY1, P,mu,co,i1:
guX:=Comps(n,d1):
guY:=Comps(m,d2):
P:=0:
mu:={}:
co:=0:
 for guX1 in guX do
  for guY1 in guY do

 co:=co+1:
  P:=P+a[co]*mul(x[i1]^guX1[i1],i1=1..n)*mul(y[i1]^guY1[i1],i1=1..m):
  mu:=mu union {a[co]}:
 od:
od:
P,mu:
end:


#ApplyMonoDOB(M,x,y,D1,D2,n,m,P): inputs a monomial differential operator (with constant coefficients) in the form
#c*D1[1]^a1*...*D1[n]^an*D2[1]^b1*...*D2[m]^bm ,c a constant and a1, .., an; b1, ..., bm 
#non-negative integers and applies it the polynomial P in x[1], ..., x[n]
#Try:
#ApplyMonoDOB(3*D1[1]^2*D1[2]^3*D2[1]^2*D2[2],x,y,D1,D2,2,2,x[1]^3*x[2]^6*y[1]^2*y[2]^3);
ApplyMonoDOB:=proc(M,x,y,D1,D2,n,m,P) local c,gu,d1,d2, i:

d1:=[seq(degree(M,D1[i]),i=1..n)]:
d2:=[seq(degree(M,D2[i]),i=1..m)]:
c:=normal(M/mul(D1[i]^d1[i],i=1..n)/mul(D2[i]^d2[i],i=1..m)):

if not type(c,numeric) then
  RETURN(FAIL):
fi:

gu:=P:

for i from 1 to n do
 gu:=MyDiff(gu,x[i],d1[i]):
od:


for i from 1 to m do
 gu:=MyDiff(gu,y[i],d2[i]):
od:

expand(c*gu):
end:



#ApplyDOB(Ope,x,y,D1,D2,n,m,P): inputs a differential operator (with constant coefficients) in the form
#Sum(c*D1[1]^a1*...*D1[n]^an*D2[1]^b1*...*D2[m]^bm) ,c a constants 
#and a1, .., an, b1, ..., bm non-negative integers and applies it the polynomial P in x[1], ..., x[n], y[1], ..., y[m]
#Try:
#ApplyDOB(D1[1]*D2[1]+D1[2]*D2[2],x,y,D1,D2,2,2,x[1]^4*y[1]^4+x[2]^3*y[2]^2);
ApplyDOB:=proc(Ope,x,y,D1,D2,n,m,P) local i:

if not type(Ope,`+`) then
ApplyMonoDOB(Ope,x,y,D1,D2,n,m,P):
else
expand( add(ApplyMonoDOB(op(i,Ope),x,y,D1,D2,n,m,P)  , i=1..nops(Ope))):
fi:

end:






#Extract1B(f,x,y,n,m): the set of coefficients of the polynomial f in the variables x[1], ... x[n], y[1], ..., y[m]. Try:
#Extract1B(a*x[1]^2+b*x[1]*x[2]*c*x[2]^5*y[1]*y[2],x,y,2,2);
Extract1B:=proc(f,x,y,n,m) local gu,gu1:

if m=0 then
RETURN(Extract1(f,x,n)):
fi:

gu:=Extract1B(f,x,y,n,m-1):
{seq(op(Extract11(gu1,y[m])),gu1 in gu)}:
end:

#KernB(n,m,d1,d2,x,y,D1, D2, S): inputs pos. integers n and m, non-negative integers
#d1,d2, symbols x and y, and a set S of differential operators with constant coefficients
#where D1[i] denotes differentiation with respect to x[i] and D2[i] denotes differentiation
#with repspect to y[i], outputs a basis for the
#set of polynomials in (x[1], ..., x[n];y[1],..,y[m]) of bi-degree (d1,d2), annihilated by the differential operators in S. Try
#KernB(3,3,2,1,x,y,D1,D2,{D1[1]*D2[1]+D1[2]*D2[2]+D1[3]*D2[3], D1[1]^2*D2[1]^2+D1[2]^2*D2[2]^2+D1[3]^2*D2[3]^2});
KernB:=proc(n,m,d1,d2,x,y,D1,D2,S) local gu,P,eq,var,a,P1,s,var1,var11,fr,fr1,gu1:


gu:=GenPolB(n,m,d1,d2,x,y,a) :



P:=gu[1]:
var:=gu[2]:

eq:={}:


for s in S do
P1:=ApplyDOB(s,x,y,D1,D2,n,m,P):
 eq:=eq union Extract1B(P1,x,y,n,m):
od:


var1:=solve(eq,var):

	
fr:={}:

for var11 in var1 do
 if op(1,var11)=op(2,var11) then
   fr:=fr union {op(1,var11)}:
fi:
od:

P:=expand(subs(var1,P)):

gu:={seq(coeff(P,fr1,1),fr1 in fr)}:

for gu1 in gu do
 if {seq(ApplyDOB(s,x,y,D1,D2,n,m,gu1),s in S)}<>{0} then
  print(`Soemthing is wrong`):
  RETURN(FAIL):
 fi:
od:

gu:

end:



#TotHB1(x,y,n,m,d1,d2): inputs  symbols x and y and  pos. integers n and m , 
#and integers d1 and d2, outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]
#and degree d2 in y[1], ..., y[m]. Try:
#TotHB1(x,y,3,3,1,1);
TotHB1:=proc(x,y,n,m,d1,d2) local S,D1,D2,i,r,s:
option remember:
if n<m then
print(`n>=m`):
 RETURN(FAIL):
fi:

S:={}:


for r from 1 to n do
 S:= S union {add(D1[i]^r,i=1..n)}:
od:


for s from 1 to m do
 S:= S union {add(D2[i]^s,i=1..m)}:
od:


for r from 1 to n do
for s from 1 to m do
 S:= S union {add(D1[i]^r*D2[i]^s,i=1..m)+add(D1[i]^r,i=m+1..n)}:
od:
od:



KernB(n,m,d1,d2,x,y,D1,D2,S):

end:



#TotHBnm(x,y,n,m): inputs  symbols x and y and  a pos. integer  n outputs
#a list of lists L such that L[d1+1][d2+1] is a basis for the vector space of
#homog. polynomials in x[1], ..., x[n], y[1], ..., y[n] if bidegree (d1, d2). Try
#TotHBnm(x,y,3,2);
TotHBnm:=proc(x,y,n,m) local d1,d2,T:
option remember:

for d1 from 0 to binomial(n,2) do
 for d2 from 0 to binomial(m,2) do
  T[d1,d2]:=TotHB1(x,y,n,m,d1,d2):
 od:
od:

[seq([seq(T[d1,d2],d2=0..binomial(m,2))],d1=0..binomial(n,2))]:
end:


#TotHB(x,y,n): inputs  symbols x and y and  a pos. integer  n outputs
#a list of lists L such that L[d1+1][d2+1] is a basis for the vector space of
#homog. polynomials in x[1], ..., x[n], y[1], ..., y[n] if bidegree (d1, d2). Try
#TotHB(x,y,3);
TotHB:=proc(x,y,n) local d1,d2,T:
option remember:

for d1 from 0 to binomial(n,2) do
 for d2 from 0 to binomial(n,2)-d1 do
  T[d1,d2]:=TotHB1(x,y,n,n,d1,d2):
 od:
od:

[seq([seq(T[d1,d2],d2=0..binomial(n,2)-d1)],d1=0..binomial(n,2))]:
end:



#HPnm(n,m,p,q): the Hilbert polynomial of diagonally-harmonic polynomials with n varibles of each kind in terms of p and q. Try:
#HPnm(3,2,p,q);
HPnm:=proc(n,m,p,q) local gu,x,y,i,j:
option remember:
gu:=TotHBnm(x,y,n,m):
add(add(nops(gu[i+1][j+1])*p^i*q^j,j=0..nops(gu[i+1])-1),i=0..nops(gu)-1):
end:


#HPn(n,p,q): the Hilbert polynomial of diagonally-harmonic polynomials with n varibles of each kind in terms of p and q. Try:
#HPn(3,p,q);
HPn:=proc(n,p,q) local gu,x,y,i,j:
option remember:
gu:=TotHB(x,y,n):
add(add(nops(gu[i+1][j+1])*p^i*q^j,j=0..nops(gu[i+1])-1),i=0..nops(gu)-1):
end:

##Start generalized procedures

#CompsG(n,d,L): the set of compositions of d into n parts <=L. Try:
#CompsG(4,3,[2,2,2,2]);
CompsG:=proc(n,d,L) local gu,i,lu,lu1:
option remember:

if not (type(L,list) and nops(L)=n) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


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

gu:={}:

for i from 0 to min(d,L[n]) do
 lu:=CompsG(n-1,d-i,[op(1..n-1,L)]):
 gu:=gu union {seq([op(lu1),i],lu1 in lu)}:
od:

gu:

end:



#GenPolG(n,d,x,a,L): inputs a positive integer n, a positive integer d, a symbol x and a symbol a and a list of length L
#outputs a generic homog. polynomial in the n variable x[1],...,x[n] of degree d with coefficients and
#degree sequence <=L , of the form a[i]
#followed by the set of coefficients. Try:
#GenPolG(4,2,x,a,[1,1,1,1]);
GenPolG:=proc(n,d,x,a,L) local gu,gu1,P,mu,co,i1:
gu:=CompsG(n,d,L):
P:=0:
mu:={}:
co:=0:
 for gu1 in gu do
 co:=co+1:
  P:=P+a[co]*mul(x[i1]^gu1[i1],i1=1..n):
  mu:=mu union {a[co]}:
 od:
P,mu:
end:


#KernG(n,d,x,D1, S,L): inputs pos. integers n and d, a symbol x, and a set S of differential operators with constant coefficients
#where D1[i] denotes differentiation with respect to x[i]. It also outputs a list L of length n for upper bounds for the degree sequence
# outputs a basis for the
#set of polynomials in x[1], ..., x[n] of degree d and degree sequence <=L, annihilated by the differential operators in S. Try
#KernG(3,2,x,D1,{D1[1]+D1[2]+D1[3],D1[1]^2+D1[2]^2+D1[3]^2},[1,1,1]);

KernG:=proc(n,d,x,D1,S,L) local gu,P,eq,var,a,P1,s,var1,var11,fr,fr1,gu1:
option remember:

gu:=GenPolG(n,d,x,a,L) :

P:=gu[1]:
var:=gu[2]:

eq:={}:

for s in S do
P1:=ApplyDO(s,x,D1,n,P):
 eq:=eq union Extract1(P1,x,n):
od:

var1:=solve(eq,var):

	
fr:={}:

for var11 in var1 do
 if op(1,var11)=op(2,var11) then
   fr:=fr union {op(1,var11)}:
fi:
od:

P:=expand(subs(var1,P)):

gu:={seq(coeff(P,fr1,1),fr1 in fr)}:

for gu1 in gu do
 if {seq(ApplyDO(s,x,D1,n,gu1),s in S)}<>{0} then
  print(`Soemthing is wrong`):
  RETURN(FAIL):
 fi:
od:
gu:
end:

##End generalized procedures



#TotHg(x,n): inputs a symbol x and a pos. integer n, outputs the list of length binomial(n,2)+1 whose
#(i+1)-th component is a  list if the bases for the increasing vector spaces of totally harmonic homog. polynomials in x[1], .., x[n]
#of degree <=j in x[n] . Try:
#TotHg(x,3):
TotHg:=proc(x,n) local S,D1,i,r,j:
option remember:
S:={seq(add(D1[i]^r,i=1..n),r=1..n)}:
[seq([seq(KernG(n,i,x,D1,S,[n$(n-1),j]),j=0.. min(i,n))],i=0..binomial(n,2))]:
end:

#TotHdim(n): the dimension sequence of totally harmomic functions. Try:
#TotHdim(3);
TotHdim:=proc(n) local gu,x,i:
gu:=TotH(x,n):
[seq(nops(gu[i]),i=1..nops(gu))]:
end:


#TotHgDim(n): the broken-up dimension sequence of totally harmomic functions. Try:
#TotHgDimG(3);
TotHgDim:=proc(n) local gu,x,i,i1,lu,mu,mu1:
option remember:
gu:=TotHg(x,n):
mu:=[]:

for i from 1 to nops(gu) do
 lu:=gu[i]:
 mu1:=[nops(lu[1]),seq(nops(lu[i1])-nops(lu[i1-1]),i1=2..nops(lu))]:
 mu:=[op(mu),mu1]:
od:

mu:

end:



#TotHBdim(n): the list of lists of dimensions of diagonal harmonics of n bi-variables.
#If the output is L then L[d1+1][d2+1] is the dimension of the vector space of
#totally harmonic bi-polynomials of bi-degree (d1,d2). Try:
#TotHBdim(3);
TotHBdim:=proc(n) local gu,x,y,i,j:
option remember:
gu:=TotHB(x,y,n):
[seq([seq(nops(gu[i+1][j+1]),j=0..nops(gu[i+1])-1)],i=0..nops(gu)-1)]:
end:



#GenPolBg(n,m,d1,d2,x,y,a,L1,L2): inputs  positive integers n and m , and positive integers d1 and d2, and  symbols x and y and a symbol a
#outputs a generic bi- homog. polynomial in the n variable x[1],...,x[n]  and m variables y[1], ..., y[m]
#of degree d1 on the x's and degree d2 in the y's  and the degree sequence of the x variables is <=L1 and the
#degree sequence of the y-variables is <=L2
#with coefficients of the form a[i]
#followed by the set of coefficients. Try:
#GenPolBg(3,3,1,1,x,y,a,[0,0,2],[0,0,2]);
GenPolBg:=proc(n,m,d1,d2,x,y,a,L1,L2) local guX,guY, guX1,guY1, P,mu,co,i1:
guX:=CompsG(n,d1,L1):
guY:=CompsG(m,d2,L2):
P:=0:
mu:={}:
co:=0:
 for guX1 in guX do
  for guY1 in guY do

 co:=co+1:
  P:=P+a[co]*mul(x[i1]^guX1[i1],i1=1..n)*mul(y[i1]^guY1[i1],i1=1..m):
  mu:=mu union {a[co]}:
 od:
od:
P,mu:
end:



#KernBg(n,m,d1,d2,x,y,D1, D2, S,L1,L2): inputs pos. integers n and m, non-negative integers
#d1,d2, symbols x and y, and a set S of differential operators with constant coefficients
#where D1[i] denotes differentiation with respect to x[i] and D2[i] denotes differentiation
#with repspect to y[i]. It also inputs lists of integers L1,L2, or length n and m  respectively.
#It outputs a basis for the
#set of polynomials in (x[1], ..., x[n];y[1],..,y[m]) of bi-degree (d1,d2),
#and degree sequence in x<=L1  and degree sequence in y<=L2
# annihilated by the differential operators in S. Try
#KernBg(3,3,2,1,x,y,D1,D2,{D1[1]*D2[1]+D1[2]*D2[2]+D1[3]*D2[3], D1[1]^2*D2[1]^2+D1[2]^2*D2[2]^2+D1[3]^2*D2[3]^2},[3,2,2],[3,2,2]);

KernBg:=proc(n,m,d1,d2,x,y,D1,D2,S,L1,L2) local gu,P,eq,var,a,P1,s,var1,var11,fr,fr1,gu1:


gu:=GenPolBg(n,m,d1,d2,x,y,a,L1,L2) :



P:=gu[1]:
var:=gu[2]:

eq:={}:


for s in S do
P1:=ApplyDOB(s,x,y,D1,D2,n,m,P):
 eq:=eq union Extract1B(P1,x,y,n,m):
od:


var1:=solve(eq,var):

	
fr:={}:

for var11 in var1 do
 if op(1,var11)=op(2,var11) then
   fr:=fr union {op(1,var11)}:
fi:
od:

P:=expand(subs(var1,P)):

gu:={seq(coeff(P,fr1,1),fr1 in fr)}:

for gu1 in gu do
 if {seq(ApplyDOB(s,x,y,D1,D2,n,m,gu1),s in S)}<>{0} then
  print(`Soemthing is wrong`):
  RETURN(FAIL):
 fi:
od:
gu:
end:



#TotHB1g(x,y,n,m,d1,d2,L1,L2): inputs  symbols x and y and  pos. integers n and m , 
#and integers d1 and d2, and lists L1 and L2 of length n and m
#respectively. Outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]
#and degree d2 in y[1], ..., y[m] and x-degree sequence <=L1 and y-degree sequence <=L2  Try:
#TotHB1g(x,y,3,3,1,1,[3,3,1],[3,3,1]);
TotHB1g:=proc(x,y,n,m,d1,d2,L1,L2) local S,D1,D2,i,r,s:
option remember:
if n<m then
print(`n>=m`):
 RETURN(FAIL):
fi:

S:={}:


for r from 1 to n do
 S:= S union {add(D1[i]^r,i=1..n)}:
od:


for s from 1 to m do
 S:= S union {add(D2[i]^s,i=1..m)}:
od:


for r from 1 to n do
for s from 1 to m do
 S:= S union {add(D1[i]^r*D2[i]^s,i=1..m)+add(D1[i]^r,i=m+1..n)}:
od:
od:



KernBg(n,m,d1,d2,x,y,D1,D2,S,L1,L2):

end:


#TotHB1br(x,y,n,m,d1,d2): Like TotHB1(x,y,n,m,d1,d2) but  outputs a list of lists
#call it L, such that L[i+1][j+1] is a basic for the subset of
# TotHB1(x,y,n,m,d1,d2) with degree in x[n]<=i and degree in y[m]<=j. Try:
#TotHB1br(x,y,3,3,1,1);

TotHB1br:=proc(x,y,n,m,d1,d2) local T,i,j:

for i from 0 to d1 do
 for j from 0 to d2 do
 T[i,j]:=TotHB1g(x,y,n,m,d1,d2,[d1$(n-1),i],[d2$(m-1),j]):
od:
od:

[seq([seq(T[i,j],j=0..d2)],i=0..d1)]:
end:


#TotHB1brDim(n,m,d1,d2): The dimensios of TotHB1br(x,y,m,n,d1,d2). Try:
#TotHB1brDim(3,3,1,1);
TotHB1brDim:=proc(n,m,d1,d2) local x,y,gu,i,j,T:

gu:=TotHB1br(x,y,n,m,d1,d2):

gu:=[seq([seq(nops(gu[i][j]),j=1..nops(gu[i]))],i=1..nops(gu))]:



T[1,1]:=gu[1][1]:

for j from 2 to nops(gu[1]) do
T[1,j]:=gu[1][j]-gu[1,j-1]:
od:


for i from 2 to nops(gu) do
T[i,1]:=gu[i][1]-gu[i-1][1]:
od:

for i from 2 to nops(gu)  do
 for j from 2 to nops(gu[i]) do
  T[i,j]:=gu[i][j]-gu[i][j-1]-gu[i-1][j]+gu[i-1][j-1]:
od:
od:

[seq([seq(T[i,j],j=1.. nops(gu[i]))],i=1..nops(gu))]:

end:





#TotHBnmDim(n,m): the list of lists of dimensions of diagonal harmonics of (n,m) bi-variables.
#If the output is L then L[d1+1][d2+1] is the dimension of the vector space of
#totally harmonic bi-polynomials of bi-degree (d1,d2). Try:
#TotHBnmDim(3,3);
TotHBnmDim:=proc(n,m) local gu,x,y,i,j:
option remember:
gu:=TotHBnm(x,y,n,m):
[seq([seq(nops(gu[i+1][j+1]),j=0..nops(gu[i+1])-1)],i=0..nops(gu)-1)]:
end:


#Hnmd1d2(x,y,n,m,d1,d2): inputs  symbols x and y and  pos. integers n and m , 
#and integers d1 and d2, outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]
#and degree d2 in y[1], ..., y[m]. Try:
#Hnmd1d2(x,y,3,3,1,1);
Hnmd1d2:=proc(x,y,n,m,d1,d2) local S,D1,D2,i,r,s:
option remember:
if n<m then
print(`n>=m`):
 RETURN(FAIL):
fi:

S:={}:


for r from 1 to n do
 S:= S union {add(D1[i]^r,i=1..n)}:
od:



for r from 0 to n do
for s from 1 to m do
 S:= S union {add(D1[i]^r*D2[i]^s,i=1..m)}:
od:
od:



KernB(n,m,d1,d2,x,y,D1,D2,S):

end:


#TotHBnmN(x,y,n,m): inputs  symbols x and y and  a pos. integer  n outputs
#a list of lists L such that L[d1+1][d2+1] is a basis for the vector space of
#homog. polynomials in x[1], ..., x[n], y[1], ..., y[n] if bidegree (d1, d2). Try
#TotHBnmN(x,y,3,2);
TotHBnmN:=proc(x,y,n,m) local d1,d2,T:
option remember:

for d1 from 0 to binomial(n,2) do
 for d2 from 0 to binomial(m,2) do
  T[d1,d2]:=Hnmd1d2(x,y,n,m,d1,d2):
 od:
od:

[seq([seq(T[d1,d2],d2=0..binomial(m,2))],d1=0..binomial(n,2))]:
end:



#TotHBnmDimN(n,m): the list of lists of dimensions of diagonal harmonics of (n,m) bi-variables.
#If the output is L then L[d1+1][d2+1] is the dimension of the vector space of
#totally harmonic bi-polynomials of bi-degree (d1,d2). Try:
#TotHBnmDimN(3,3);
TotHBnmDimN:=proc(n,m) local gu,x,y,i,j:
option remember:
gu:=TotHBnmN(x,y,n,m):
[seq([seq(nops(gu[i+1][j+1]),j=0..nops(gu[i+1])-1)],i=0..nops(gu)-1)]:
end:



#Ynm(n,m): the list of lists of dimensions of diagonal harmonics of (n,m) bi-variables.
#If the output is L then L[d1+1][d2+1] is the dimension of the vector space of
#totally harmonic bi-polynomials of bi-degree (d1,d2). Try:
#Ynm(3,3);
Ynm:=proc(n,m) local gu,x,y,i,j:
option remember:
gu:=TotHBnmN(x,y,n,m):
[seq([seq(nops(gu[i+1][j+1]),j=0..nops(gu[i+1])-1)],i=0..nops(gu)-1)]:
end:



#HndK(x,n,d,K): a basis for the vector space of totally harmonic polynomials in the n variables
#x[1], ..., x[n], of total degree d in x[1], ..., x[n] and degree in x[n]<=K. Try:
#HndK(x,4,2,1);
HndK:=proc(x,n,d,K) local i,r,D1:
option remember:
 KernG(n,d,x,D1,{seq(add(D1[i]^r,i=1..n),r=1..n)},[d$(n-1),K]):
end:


#CheckHC1(n,d,K): computes dim(HndK(x,n,d,K))- dim(HndK(x,n,d,K-1))-dim(HndK(x,n-1,d-K,d-K)
CheckHC1:=proc(n,d,K) local x:

if K=0 then
 nops(HndK(x,n,d,0))-nops(HndK(x,n-1,d,d)),  [nops(HndK(x,n,d,0)),nops(HndK(x,n-1,d,d))]:
elif K>=n and d>=n then
 nops(HndK(x,n,d,K))- nops(HndK(x,n,d,K-1)), [nops(HndK(x,n,d,K)), nops(HndK(x,n,d,K-1))]:
else
 nops(HndK(x,n,d,K))- nops(HndK(x,n,d,K-1))-nops(HndK(x,n-1,d-K,d-K)), [nops(HndK(x,n,d,K)), nops(HndK(x,n,d,K-1)),nops(HndK(x,n-1,d-K,d-K))]:
fi:
end:

CheckHC:=proc(n) local d,K:
evalb({seq(seq(CheckHC1(n,d,K)[1],K=0..d),d=0..binomial(n,2))}={0}):
end:

#A(n,d,K): the dimension of the vector space of totally harmonic polynomials in n variables of total degree d and
#degree in x[n]<=K. Try:
#n1:=4:[seq(A(n1,d,d),d=0..binomial(n1,2))]:
A:=proc(n,d,K)
option remember:

if n=1 then
  if d=0 and K=0 then
    RETURN(1):
  else
   RETURN(0):
 fi:
fi:

if K<0 or K>d then 
0:
elif K=0 then
A(n-1,d-K,d-K):
elif
K>=n then
A(n,d,K-1):
else
A(n,d,K-1)+A(n-1,d-K,d-K):
fi:
end:




#TotHB1gN(x,y,n,m,d1,d2,L1,L2): New way
#inputs  symbols x and y and  pos. integers n and m , 
#and integers d1 and d2, and lists L1 and L2 of length n and m
#respectively. Outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]
#and degree d2 in y[1], ..., y[m] and x-degree sequence <=L1 and y-degree sequence <=L2  Try:
#TotHB1gN(x,y,3,3,1,1,[3,3,1],[3,3,1]);
TotHB1gN:=proc(x,y,n,m,d1,d2,L1,L2) local S,D1,D2,i,r,s:
option remember:
if n<m then
print(`n>=m`):
 RETURN(FAIL):
fi:

S:={}:


for r from 1 to n do
 S:= S union {add(D1[i]^r,i=1..n)}:
od:



for r from 0 to n do
for s from 1 to m do
 S:= S union {add(D1[i]^r*D2[i]^s,i=1..m)}:
od:
od:

KernBg(n,m,d1,d2,x,y,D1,D2,S,L1,L2):

end:





#Hnmd1d2s(x,y,n,m,d1,d2,L1,L2): Smart New way, like Hnmd1d2(x,y,n,m,d1,d2,L1,L2) but smart
Hnmd1d2s:=proc(x,y,n,m,d1,d2) local S,D1,D2,i,j,r,s,gu1,gu2,a,lu,var,var1,fr,fr1,P,P1, eq, gu, var11:
option remember:
if n<m then
print(`n>=m`):
 RETURN(FAIL):
fi:

S:={}:


for r from 1 to n do
 S:= S union {add(D1[i]^r,i=1..n)}:
od:



for r from 0 to n do
for s from 1 to m do
 S:= S union {add(D1[i]^r*D2[i]^s,i=1..m)}:
od:
od:

gu1:=TotH(x,n)[d1+1]:
gu2:=TotH(y,n)[d2+1]:

lu:=add(add(a[i,j]*gu1[i]*gu2[j],i=1..nops(gu1)),j=1..nops(gu2)):
var:={seq(seq(a[i,j],i=1..nops(gu1)),j=1..nops(gu2))}:

eq:={}:
for s in S do
P1:=ApplyDOB(s,x,y,D1,D2,n,m,lu):
 eq:=eq union Extract1B(P1,x,y,n,m):
od:

var1:=solve(eq,var):
lu:=subs(var1,lu):

	
fr:={}:

for var11 in var1 do
 if op(1,var11)=op(2,var11) then
   fr:=fr union {op(1,var11)}:
fi:
od:

P:=expand(subs(var1,lu)):

gu:={seq(coeff(P,fr1,1),fr1 in fr)}:

for gu1 in gu do
 if {seq(ApplyDOB(s,x,y,D1,D2,n,m,gu1),s in S)}<>{0} then
  print(`Soemthing is wrong`):
  RETURN(FAIL):
 fi:
od:

gu:


end:






#IsLC(p,S,x,n): Is the polynomial p in x[1], ..., x[n] a linear combination of the members of S?
#Try (x[1]+x[2],{x[1],x[2]},x,2);
IsLC:=proc(p,S,x,n) local S1,eq,var,a,gu,i:
S1:=convert(S,list):
var:={seq(a[i],i=1..nops(S1))}:
gu:=expand(add(a[i]*S1[i],i=1..nops(S1) )-p):
eq:=Extract1(gu,x,n) minus {0}:
var:=solve(eq,var):
if var=NULL then
 false:
else
 true:
fi:
end:

#KickOut(S,x,n): tries to kick  a member of S that is a linear combination of the other ones, or returns FAIL
KickOut:=proc(S,x,n) local s:
for s in S do
if IsLC(s,S minus {s},x,n) then
 RETURN(S minus {s}):
fi:
od:
S:
end:


#FindBase(S,x,n): Given a set S of plolynomials in x[1], ..., x[n], outputs a base for the vector space it generates . Try:
#FindBase({x[1]+x[2],x[2]+x[3],x[1]+2*x[2]+x[3]},x,3);
FindBase:=proc(S,x,n) local S1,S2,S3:
S1:=S:
S2:=KickOut(S1,x,n):
while S1<>S2 do
 S3:=KickOut(S2,x,n):
 S1:=S2:
 S2:=S3:
od:
S1:
end:




#HndKLC(x,n,d,K): a basis for the vector space of LEADING COEFFICIENT in x[n] for totally harmonic polynomials in the n variables
#x[1], ..., x[n], of total degree d in x[1], ..., x[n] and degree in x[n]<=K. Try:
#HndKLC(x,4,2,1);
HndKLC:=proc(x,n,d,K) local gu,gu1:

gu:=HndK(x,n,d,K):

gu:={seq(coeff(gu1,x[n],K),gu1 in gu)}:
FindBase(gu,x,n-1):
end:




#Hnmd1d2K(x,y,n,m,d1,d2,K): inputs  symbols x and y and  pos. integers n and m , 
#and integers d1 and d2, and a non-neg. integer K <=d2,
#outputs a basis for the set of doubly-harmonic polynomials of degree d1 in x[1], ..., x[n]
#and degree d2 in y[1], ..., y[m] and deg_y[m]<=K. Try:
#Hnmd1d2K(x,y,3,3,1,2,1);
Hnmd1d2K:=proc(x,y,n,m,d1,d2,K) local S,D1,D2,i,r,s:
option remember:
if n<m then
print(`n>=m`):
 RETURN(FAIL):
fi:

S:={}:


for r from 1 to n do
 S:= S union {add(D1[i]^r,i=1..n)}:
od:



for r from 0 to n do
for s from 1 to m do
 S:= S union {add(D1[i]^r*D2[i]^s,i=1..m)}:
od:
od:



KernBg(n,m,d1,d2,x,y,D1,D2,S,[d1$n],[d2$(m-1),K]):

end:






#IsLCB(p,S,x,y,n,m): Is the polynomial p in x[1], ..., x[n] a linear combination of the members of S?
#Try (x[1]+x[2],{x[1],x[2]},x,2);
IsLCB:=proc(p,S,x,y,n,m) local S1,eq,var,a,gu,i:
S1:=convert(S,list):
var:={seq(a[i],i=1..nops(S1))}:
gu:=expand(add(a[i]*S1[i],i=1..nops(S1) )-p):
eq:=Extract1B(gu,x,y,n,m) minus {0}:
var:=solve(eq,var):
if var=NULL then
 false:
else
 true:
fi:
end:


#KickOutB(S,x,y,n,m): tries to kick  a member of S that is a linear combination of the other ones, or returns FAIL
KickOutB:=proc(S,x,y,n,m) local s:
for s in S do
if IsLCB(s,S minus {s},x,y,n,m) then
 RETURN(S minus {s}):
fi:
od:
S:
end:

#FindBaseB(S,x,y,n,m): Given a set S of plolynomials in x[1], ..., x[n], y[1], ..., y[m]
#outputs a base for the vector space it generates . Try:
#FindBaseB({x[1]+x[2],x[2]+x[3],x[1]+2*x[2]+x[3]},x,y,3,0);
FindBaseB:=proc(S,x,y,n,m) local S1,S2,S3:
S1:=S:
S2:=KickOutB(S1,x,y,n,m):
while S1<>S2 do
 S3:=KickOutB(S2,x,y,n,m):
 S1:=S2:
 S2:=S3:
od:
S1:
end:

#Hnmd1d2KLC(x,y,n,m,d1,d2,K): a basis for the vector space of LEADING COEFFICIENT in y[m] for totally bi-harmonic 
#polynomials in the n x variables and m y-variables for the vector space
#of total bi-degree (d1,d2)  and degree in y[m]<=K. Try:
#HndKLC(x,y,4,3,2,2,1);
Hnmd1d2KLC:=proc(x,y,n,m,d1,d2,K) local gu,gu1:

gu:=Hnmd1d2K(x,y,n,m,d1,d2,K):

gu:={seq(coeff(gu1,y[m],K),gu1 in gu)}:
FindBaseB(gu,x,y,n,m-1):
end:


#DimKernSeq(n,d,D1,S), inputs a pos. integer n, a pos. integer,d, a symbol D1, and a set S of linear differential operators with constant coefficients
#output the first d+1 terms, starting at i=0 or the vector space of polynomials in x[1], ..., x[n] annihilated by the member of S. Try:
#DimKernSeq(3,4,D1,{seq(add(D1[i]^s,i=1..3),s=1..3)});
DimKernSeq:=proc(n,d,D1,S)  local x,d1:
[seq(nops(Kern(n,d1,x,D1,S)),d1=0..d)]:
end:
ez:=proc(): print(`Vecs(x,n,S), Comps(n,d), IsMore1(v1,v2) , IsMore(v1,v2), Surv(x,n,S,d)`): end:

#Using Grobner bases

#Vecs(x,n,S), inputs a variable name x, a positive integer n, and a set of polynomials S, in x[1], ..., x[n], outputs
#the set of exponents of the leading monomial of the Grobner basis of the ideal generated by S with pure lex order x[1], ..., x[n].
#Try:
#Vecs(x,2,{x[1]+x[2]});
Vecs:=proc(x,n,S) local gu,i,j:
option remember:
gu:=Groebner[Basis](S,plex(seq(x[i],i=1..n))):
gu:=Groebner[LeadingMonomial](gu,plex(seq(x[i],i=1..n)) ):

[seq(
[seq(degree(gu[j],x[i]),i=1..n)],
   j=1..nops(gu)
    )
]:
end:

#Comps(n,d): The set of vectors of length n with non-negative integers that add-up to d. Try:
#Comps(3,4);
Comps:=proc(n,d) local gu,d1,mu,mu1:
option remember:

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

gu:={}:

for d1 from 0 to d do
 mu:=Comps(n-1,d-d1):
 gu:=gu union {seq([op(mu1),d1], mu1 in mu)}:
od:
gu:
end:

#IsMore1(v1,v2) Is v1>=v2?
IsMore1:=proc(v1,v2) local i:
if nops(v1)<>nops(v2) then
 RETURN(FAIL):
fi:

if min(seq(v1[i]-v2[i],i=1..nops(v1)))<0 then
 RETURN(false):
fi:
true:
end:




# IsMore(v1,V): is v1 more than at least one member of V
IsMore:=proc(v1,V)  local v2:

for v2 in V do
if IsMore1(v1,v2) then
 RETURN(true):
fi:
od:
false:
end:

#Surv(x,n,S,d): all the vectors of non-neg. integers of length n, that are never >=then the the vectors in
#Vecs(x,n,S);. Try:
#Surv(x,2,{x[1]+x[2],x[1]^2+x[2]^2},1);
Surv:=proc(x,n,S,d) local lu,mu,mu1,gu:

lu:=Vecs(x,n,S):

mu:=Comps(n,d):

gu:={}:

for mu1  in mu do
 if not IsMore(mu1,lu) then
  gu:=gu union {mu1}:
 fi:
od:
gu:
end:




#DimKernSeqF(n,d,D1,S): Like  DimKernSeq(n,d,D1,S),  but using Grobner bases and survivors. Try:
#DimKernSeqF(3,4,D1,{seq(add(D1[i]^s,i=1..3),s=1..3)});
DimKernSeqF:=proc(n,d,D1,S)  local   d1:

[seq(nops(Surv(D1,n,S,d1)),d1=0..d)]:
end:

#End using Grobner bases





#VecsTold(x,y,n,S), inputs  variable names x, and y, a positive integer n, and a set of polynomials S, in x[1], ..., x[n], y[1], ..., y[n], outputs
#the set of exponents of the leading monomial of the Grobner basis of the ideal generated by S with pure lex order x[1],y[1] ..., x[n], y[n]
#Try:
#VecsTold(x,y,2,{x[1]+x[2],y[1]+y[2],x[1]*y[1]+x[2]*y[2]});
VecsTold:=proc(x,y,n,S) local gu,i,j:
option remember:
gu:=Groebner[Basis](S,plex(seq(op([x[i],y[i]]),i=1..n))):
gu:=Groebner[LeadingMonomial](gu,plex(seq(op([x[i],y[i]]),i=1..n)) ):

[seq([seq([degree(gu[j],x[i]),degree(gu[j],y[i])],i=1..n)],   j=1..nops(gu)    )]:
end:



#VecsT(x,y,m,n,S), inputs  variable names x, and y, a positive integers m and n, and a set of polynomials S, in x[1], ..., x[m], y[1], ..., y[n], outputs
#the set of exponents of the leading monomial of the Grobner basis of the ideal generated by S with pure lex order x[1], ..., x[m], y[1], ..., y[n],
#Try:
#VecsT(x,y,2,2,{x[1]+x[2],y[1]+y[2],x[1]*y[1]+x[2]*y[2]});
VecsT:=proc(x,y,m,n,S) local gu,i,j:
option remember:
gu:=Groebner[Basis](S,plex(seq(x[i],i=1..m), seq(y[i],i=1..n) )):
gu:=Groebner[LeadingMonomial](gu,plex(seq(x[i],i=1..m), seq(y[i],i=1..n) ) ) :

[seq([[seq(degree(gu[j],x[i]),i=1..m)], [seq(degree(gu[j],y[i]),i=1..n)]],j=1..nops(gu))]:
end:





# IsMoreT(v1,v2,V): is [v1,v2] more than at least one member of the set of pairs V
IsMoreT:=proc(v1,v2,V)  local w:

for w in V do
if IsMore1(v1,w[1]) and IsMore1(v2,w[2]) then
 RETURN(true):
fi:
od:
false:
end:

#SurvT(x,y,m,n,S,d1,d2): all the pair of vectors [v1,v2] of non-neg. integers of length m and n respectively, 
#that are never >=then the the vectors in
#VecsT(x,y,m,n,S,d1,d2);. Try:
#SurvT(x,y,2,2,{x[1]+x[2],x[1]^2+x[2]^2,x[1]*y[1]+x[2]*y[2],x[1]^2*y[1]+x[2]^2*y[2]},1,1);

SurvT:=proc(x,y,m,n,S,d1,d2) local lu,mu1,mu2,mu1a, mu2a, gu:

lu:=VecsT(x,y,m,n,S):

mu1:=Comps(m,d1):
mu2:=Comps(n,d2):

gu:={}:

for mu1a  in mu1 do
 for mu2a  in mu2 do
  if not IsMoreT(mu1a,mu2a,lu) then
   gu:=gu union {[mu1a,mu2a]}:
  fi:
 od:
od:

gu:
end:




#DimKernSeqTF(m,n,d1,d2,D1,D2,S): The list of listsL such that L[d1][d2] is the dimension of the set of polynomials of b-degree [d1,d2]
#annihilated by S. Try
#DimKernSeqTF(3,3,1,1,D1,D2,{seq(seq(add(D1[i]^r*D2[i]^s,i=1..3),r=1..3),s=1..3)});
DimKernSeqTF:=proc(n,m,d1,d2,D1,D2,S)  local   d1a,d2a:
[seq([seq(nops(SurvT(D1,D2,n,m,S,d1a,d2a)),d2a=0..d2)],d1a=0..d1)]:
end:


#YnmF(n,m): Faster version of Ynm(n,m), using Groebner bases
YnmF:=proc(n,m)  local  D1,D2,S,r,s,i:
 S:={seq(add(D1[i]^r,i=1..n),r=1..n), seq(seq(add(D1[i]^r*D2[i]^s,i=1..m),r=0..n),s=1..m)}:
DimKernSeqTF(n,m,binomial(n,2),binomial(m,2),D1,D2,S):

end:


#gadol(S,n): given a set of vectors  S, each of length n, outputs the pointwise max
gadol:=proc(S,n) local s,i:

if S={} then
[0$n]:
else
 [seq(max(seq(s[i],s in S)),i=1..n)]:
fi:
end:


#VecToHS(S,q,n): Given a set of vectors S in N^n outputs the generating function for the lattice points that are NOT >= and of the vectors of S
#Try:
#VecToHS({[1,0],[0,2]},x,2);
VecToHS:=proc(S,q,n) local gu,S1,s1,lu,i,x:
S1:=powerset(convert(S,set)):
gu:=0:

for s1 in S1 do

 lu:=gadol(s1,n):
 lu:=mul(x[i]^lu[i],i=1..n):
 gu:=gu+(-1)^nops(s1)*lu:
od:

gu:=gu/mul(1-x[i],i=1..n):

factor(subs({seq(x[i]=q,i=1..n)},gu)):

end:


#GFup(L,n,q): Given an ultimate polynomial L=[INI,Pol(n)] in terms of the variable n, outputs the
#rational function that is its  genearting function. Try:
#GFup([[1,1],n^2],q):
GFup:=proc(L,n,q) local gu,s,mu,i,P:

gu:=add(L[1][i+1]*q^i,i=0..nops(L[1])-1):

s:=nops(L[1]):


mu:=q^s/(1-q):

P:=L[2]:

gu:=normal(gu+coeff(P,n,0)*mu):

for i from 1 to degree(P,n) do
 mu:=normal(q*diff(mu,q)):
 gu:=normal(gu+coeff(P,n,i)*mu):
od:

gu:

end:


#HSS(n1,S,K,q): inputs a positive integer n, a subset S, of {1, ..., n1} a positive integer K, and a variable q,
#outputs the Hilbert series of the vector space of polynomials in (x[1], ..., x[n1]) annihilated by
#Sum(D1[i]^s,i=1..n1) for s in S. Using K+1 data points for guessing. Try:
#HSS(3,{1,2,3},12,q);
HSS:=proc(n1,S,K,q) local gu,D1,i1,n,mu,s:

gu:=DimKernSeqF(n1,K,D1, {seq(add(D1[i1]^s,i1=1..n1), s in S)}):
mu:=GuessPolU(gu,n):

if mu=FAIL then
 RETURN(gu):
fi:

GFup(mu,n,q):

end: