#########################################################################
##ShapiroGeneral.txt: Save this file as   ShapiroGeneral.txt            #
## To use it, stay in the                                               #
##same directory, get into Maple (by typing: maple <Enter> )            #
##and then type:  read `ShapiroGeneral.txt`<Enter>                      #
##Then follow the instructions given there                              #
##                                                                      #
##Written by Doron Zeilberger, Rutgers University ,                     #
#zeilberg at math dot rutgers dot edu                                   #
#########################################################################
 
#Created: May 2016



print(`Created:  May 2016`):
print(` This is ShapiroGeneral.txt `):
print(`It is one of the packages that accompany the article `):
print(``):
print(` Integrals Involving Rudin-Shapiro Polynomials and Sketch of a Proof of Saffari's Conjecture`):
print(``):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(`and also available from Zeilberger's website`):
print(``):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/ShapiroGeneral.txt  .`):
print(``):
print(`Dedicated TO  Krishnaswami "Krishna" Alladi , on his 60th birthday, without whom this project would not have to be.`):
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(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezraOri(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 Original Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):



ezraS:=proc()

if args=NULL then
 print(` The story procedures are: GFwS `):
else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Canon, CheckGFw, Eval1,  GDg, GFwD, Monos,  Sch, SidE, SolSch  `):
else
ezra(args):
fi:

end:

ezra:=proc() :
if args=NULL then
 print(`The main procedures are:  GFw, PnG `):
 print(` `):



elif nops([args])=1 and op(1,[args])=Canon then
print(`Canon(F,a,b,A,B,z): the Canonical form of an expression F in a,b,A,B,z Try:`):
print(`Canon(a*A+b*B,a,b,A,B,z);`):


elif nops([args])=1 and op(1,[args])=CheckGFw then
print(`CheckGFw(w,a,b,A,B,z,K,t,C,N,r): `):
print(` checks the N first terms of GFw(w,a,b,A,B,z,K,t,C,N,r) against`):
print(` direct calculations. `):
print(` Try: `):
print(`CheckGFw(a^2*B^2,a,b,A,B,z,100,t,[1,z,z,1],6,2);`):
print(`CheckGFw(a^2*B^2,a,b,A,B,z,100,t,[1,z,z,1],6,3);`):

elif nops([args])=1 and op(1,[args])=Eval1 then
print(`Eval1(n,w,a,b,A,B,z,C,r): given a non-negative integer n, a polynomial expression w, `):
print(`in a,b,A,B, and z and a recurrence `):
print(` templace C, of degree r, and positive integer n, evaluates what happens `):
print(` when a is replaced by f(z):=PnG(n,C,z,r), A by f(1/z), b by f(-z) and B by f(-1/z) `):
print(` and then the CONSTANT TERM is taken `):
print(` Try: `):
print(` Eval1(4,(a*A)^2,a,b,A,B,z,[1,z,0,0],2); `):

elif nops([args])=1 and op(1,[args])=GDg then
print(`GDg(F,a,b,A,B,z,C,r): The Going Down recurrence for an expression F in generalized`):
print(`Shapiro polynomials given by capsule C of degree r `):
print(`Let PnG(n,z,C) be the generalized Shapiro polynomials with capsule `):
print(`(where a,b,c,d are Laruent polynomials of degree 1 and lowedgree -1)`):
print(`i.e. definining the recurrence`):
print(`P_n=C[1](z)*P_{n-1}(z^r)+C[2](z)*P_{n-1}(-z^r)+ C[3](z)*P_{n-1}(1/z^r)+ C[4](z)*P_{n-1}(-1/z^r).`):
print(`This procedure inputs a polynomial in a,b,A,B,z and outputs another polynomial such that`):
print(`if we replace in F`):
print(` a <- P_n(z),  b <- P_n(-z),  A <- P_n(1/z),  B <- P_n(-1/z), and use the abobe defining recurrence`):
print(`the ConstanTerm of F is the same as the constant term of`):
print(`the output when we do the above analogous replacements `):
print(` a <-P_(n-1)(z),  b <-P_(n-1)(-z),  A <-P_(n-1)(1/z),  B <-P_(n-1)(-1/z). Try: `):
print(`GDg((a*A)^2,a,b,A,B,z,[1,z,0,0],2);`):


elif nops([args])=1 and op(1,[args])=GFw then
print(`GFw(w,a,b,A,B,z,K,t,C,r): The generating function of the monomial w in`):
print(`in a,b,A,B,z with limit of scheme K, variable t, and capsule C of degree r`):
print(` Try: `):
print(`GFw(a*A,a,b,A,B,z,10,t,[1,z,0,0],2); `):
print(` GFw(a^2*A^2,a,b,A,B,z,10,t,[1,z,0,0],2); `):
print(`GFw(a^2*B^2*z,a,b,A,B,z,100,t,[1,z,z+1/z,z],2); `):
print(`GFw(a^2*B^2*z,a,b,A,B,z,100,t,[1+z^2,1+1/z,z+1/z,z],3); `):

elif nops([args])=1 and op(1,[args])=GFwD then
print(`GFwD(w,a,b,A,B,z,K,t,C,r): Like GFwe(w,a,b,A,B,z,K,t,C,r) (q.v.) but w/o guessing, but rather by solving the scheme `):
print(`A little slower, but perhaps "safer" `):
print(` Try: `):
print(`GFwD(a*A,a,b,A,B,z,10,t,[1,z,0,0],2); `):
print(` GFwD(a^2*A^2,a,b,A,B,z,10,t,[1,z,0,0],2); `):
print(`GFwD(a^2*B^2*z,a,b,A,B,z,100,t,[1,z,z+1/z,z],2); `):
print(`GFwD(a^2*B^2*z,a,b,A,B,z,100,t,[1+z^2,1+1/z,z+1/z,z],3); `):

elif nops([args])=1 and op(1,[args])=GFwS then
print(`GFwS(w,a,b,A,B,z,K,t,C,r): Verbose version of GFw(w,a,b,A,B,z,K,t,C,r) (q.v.)`):
print(` Try: `):
print(` GFwS(a*A,a,b,A,B,z,10,t,[1,z,0,0],2); `):
print(` GFwS(a^2*A^2,a,b,A,B,z,10,t,[1,z,0,0],2); `):
print(` GFwS(a^2*B^2*z,a,b,A,B,z,100,t,[1/z,z,z,1],2); `):
print(` GFwS(a^2*B^2,a,b,A,B,z,100,t,[1+z+z^2,1+1/z+z,z^2+1/z^2,3+z+1/z^2],3); `):


elif nops([args])=1 and op(1,[args])=Monos then
print(`Given a polynomial F in the list of variables var, finds the set {[coe,mono]}  such that the`):
print(`polynomial is a the sum of coe*mono for [coe,mono] in S. Try:`):
print(`Monos(3*x*y*z,[x,y,z]);`):
print(`Monos((x+y)*(x+3*z),[x,y,z]);`):

elif nops([args])=1 and op(1,[args])=PnG then

print(`PnG(n,C,z,r): Inputs a non-negative integer n, a "capsule" C that is a list of four Laurent polynomials`):
print(`of degree <r and low-degree larger than -r, C=[a,b,c,z] , and r is the degree of the Capsule`):
print(`and outputs the n-th Laurent polynomial in the sequence`):
print(`of Laurent polynomials defined by the recurrence`):
print(`P_n=a(z)*P_{n-1}(z^r)+b(z)*P_{n-1}(-z^r)+ c(z)*P_{n-1}(1/z^r)+ d(z)*P_{n-1}(-1/z^r).`):
print(`For example to get the fourth original Rudin-Shapiro polynomial, type:`):
print(`PnG(4,[1,z,0,0],z,2):`):
print(``):


elif nops([args])=1 and op(1,[args])=Sch then
print(` Sch(w,a,b,A,B,z,K,C,r): inputs a word w in symbols a,b,A,B,z, a positive integer K,  `):
print(`and a capsule C, of degree r, of length 4 whose entries are of the form e1/z+e2+e3*z (for e1,e2,e3 numbers)`):
print(`and  outputs `):
print(`the  transition  matrix, in SPARSE format, as a list, where the i-th entry {[c1,j1], ...}   `):
print(`means that  M[i,j1]=c1, and M[i,j]=0 if not mentioned.`):
print(`that enables to compute`):
print(`ordinary generating function,  or computing many terms,`):
print(` a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).`):
print(`where P_n(z) are the generalized Shapiro polynomials with capsule C, PnG(n,z,C,r);`):
print(`It uses a recursive scheme with at most K member. If none is found it returns FAIL.`):
print(`It also outputs the initial vector, and the "ledend", i.e. the list whose `):
print(`i-th item is the corresponding word `):
print(`Try `):
print(`Sch(a*A,a,b,A,B,z,10,[1,z,0,0],2);`):
print(`Sch(a^2*A^2,a,b,A,B,z,10,[1,z,0,0],2); `):
print(`Sch(a^2*B^2*z,a,b,A,B,z,10,[1,z,-1,-1/z],2);`):
print(`Sch(a^2*B^2*z,a,b,A,B,z,10,[1,z,-1,-1/z],3);`):

elif nops([args])=1 and op(1,[args])=SidE then
print(`SidE(N,w,a,b,A,B,z,C,r): the first N+1 terms of the coeff. of the constant term in  `):
print(`Eval1(n,w,z,a,b,A,B,z,C,r) (q.v.) `):
print(` Try: `):
print(` SidE(6,(a*A)^2/z^2,a,b,A,B,z,[1,z,0,0],2); `):

elif nops([args])=1 and op(1,[args])=SolSch then
print(`SolSch(M,t): solves the scheme M with respect of t, only finds the value of the first component`):
print(`Try:`):
print(`SolSch([[1],[{[2,1]}],[a]],t);`):

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


##From Cfinite

#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:

#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant coefficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d then
 print(`The list must be of size >=`, 2*d+1 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:


eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..2*d)}:

var:=solve(eq,var):


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

if {seq(L[n]-add(subs(var,a[i])*L[n-i],i=1..d),n=2*d+1..nops(L))}<>{0} then
  RETURN(FAIL):
fi:


[[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]:


end:



#RtoS(f,t,N): the first N+1 coefficients, stating at t^0 of the Maclaurin expansion of the function f of t. For example, try:
#RtoS(1/(1-t-t^3),t,30);
RtoS:=proc(f,t,N) local gu,i:
gu:=taylor(f,t=0,N+1):
[seq(coeff(gu,t,i),i=0..N)]:

end:

#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant coefficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc((nops(L)-1)/2) do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:

#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree1(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(factor(f)):
fi:
end:

##end From Cfinite


#MtoSeq(V1,M,N): inputs an initial vector V1 a sparse matrix M, and a positive integer N
#outputs the first N+1 terms, starting at n=0 of the sequence M^n(v1)[1]
#Try:
#MtoSeq([1],[{[2,1]}],10);
MtoSeq:=proc(V1,M,N) local gu,mu,n,i,j:

mu:=V1:

gu:=[mu[1]]:

for n from 1 to N do
mu:=[seq(add(M[i][j][1]*mu[M[i][j][2]],j=1..nops(M[i])),i=1..nops(M))]:
gu:=[op(gu),mu[1]]:
od:

gu:

end:



#Given a polynomial F in the list of variables var, finds the set {[coe,mono]}  such that the
#polynomial is a the sum of coe*mono for [coe,mono] in S. Try:
#Monos(3*x*y*z,[x,y,z]);
#Monos((x+y)*(x+3*z),[x,y,z]);
Monos:=proc(F,var) local F1,mono,c,gu,F11,i1,i:
F1:=expand(F):

if not type(F1,`+`) then
 mono:=mul(var[i1]^degree1(F1,var[i1]),i1=1..nops(var)):
 c:=normal(F1/mono):
  if not type(c,numeric) then
   RETURN(FAIL):
   else
    RETURN({[c,mono]}):
  fi:

fi:

gu:={}:

for i from 1 to nops(F1) do
 F11:=op(i,F1):

 mono:=mul(var[i1]^degree1(F11,var[i1]),i1=1..nops(var)):
 c:=normal(F11/mono):
  if not type(c,numeric) then
   RETURN(FAIL):
   else
    gu:=gu union {[c,mono]}:
  fi:
od:

gu:

end:




#PnG(n,C,z,r): Inputs a non-negative integer n, a "capsule" C that is a list of four Laurent polynomials
#of degree <r and low-degree larger than -r C=[a,b,c,z] and outputs the n-th Laurent polynomial in the sequence
#of Laurent polynomials defined by the recurrence
#P_n=a(z)*P_{n-1}(z^r)+b(z)*P_{n-1}(-z^r)+ c(z)*P_{n-1}(1/z^r)+ d(z)*P_{n-1}(-1/z^r).
#For example to get the fourth original Rudin-Shapiro polynomial, type:
#PnG(4,[1,z,0,0],z,2):
PnG:=proc(n,C,z,r) local gu:
option remember:

if n=0 then
 RETURN(1):
else
gu:=PnG(n-1,C,z,r):
RETURN(expand(C[1]* subs(z=z^r,gu) +C[2]*subs(z=-z^r,gu) +C[3]*subs(z=1/z^r,gu) +C[4]*subs(z=-1/z^r,gu) )):
fi:

end:




#Coe(w,a,b,A,B,z): given a monomial w, in a,b,A,B,z, finds it coefficients or 
#returns FAIL if not a monomial. Try:
#Coe(15*a^2*b^3*A*B*z,a,b,A,B,z);
Coe:=proc(w,a,b,A,B,z) local d1,d2,d3,d4,d5,coe:

d1:=degree1(w,a): d2:=degree1(w,b): d3:=degree1(w,A): d4:=degree1(w,B): d5:=degree1(w,z):

coe:=normal(w/a^d1/b^d2/A^d3/B^d4/z^d5):

if not type(coe,integer) then
 FAIL:
else
 coe:
fi:

end:




degree1:=proc(F,z)
if F=0 then
 RETURN(0):
else
 RETURN(degree(F,z)):
fi:
end:

ldegree1:=proc(F,z)
if F=0 then
 RETURN(0):
else
 RETURN(ldegree(F,z)):
fi:
end:


#Canon1(w,a,b,A,B,z): the Canonical form of the word w. Try:
#Canon1(b*B,a,b,A,B,z);
Canon1:=proc(w,a,b,A,B,z) local w1:

w1:=w:

if degree1(w1,z)<0 then
 w1:=subs({z=1/z,a=A,A=a,b=B,B=b},w1):
fi:

if degree1(w1,b)>degree1(w1,a) or (degree1(w1,b)=degree1(w1,a) and degree1(w1,B)>degree1(w1,A)) then
 w1:=subs({z=-z,a=b,b=a,A=B,B=A},w1):
fi:

w1:
end:

#Canon(F,a,b,A,B,z): the Canonical form of the expression F in a,b,A,B,z. Try:
#Canon(a*A+b*B,a,b,A,B,z);
Canon:=proc(F,a,b,A,B,z) local F1,gu,w,c,i:


F1:=expand(F):

if F1=0 then
  RETURN(0):
fi:

if not type(F1,`+`) then
 c:=Coe(F1,a,b,A,B,z) :
  if c=FAIL then
    RETURN(FAIL):
  fi:
 F1:=Canon1(normal(F1/c),a,b,A,B,z):
  RETURN(c*F1):

fi:

gu:=0:

for i from 1 to nops(F1) do
 w:=op(i,F1):


 c:=Coe(w,a,b,A,B,z) :


  if c=FAIL then
    RETURN(FAIL):
  fi:
 w:=Canon1(normal(w/c),a,b,A,B,z):

  gu:=gu+c*w:
od:

gu:
end:

#GDg(F,a,b,A,B,z,C,r): The Going Down recurrence for an expression F in generalized
#Shapiro polynomials given by capsule C with degree-r recurrence
#Let PnG(n,z,C,r) be the generalized Shapiro polynomials with capsule C of degree r
#(where a,b,c,d are Laruent polynomials of degree 1 and lowedgree -1)
#i.e. definining the recurrence
#P_n=C[1](z)*P_{n-1}(z^r)+C[2](z)*P_{n-1}(-z^r)+ C[3](z)*P_{n-1}(1/z^r)+ C[4](z)*P_{n-1}(-1/z^r).
#This procedure inputs a polynomial in a,b,A,B,z and outputs another polynomial such that
#if we replace in F
# a <- P_n(z),  b <- P_n(-z),  A <- P_n(1/z),  B <- P_n(-1/z), and use the abobe defining recurrence
#the ConstanTerm of F is the same as the constant term of
#the output when we do the above analogous replacements 
#a <-P_(n-1)(z),  b <-P_(n-1)(-z),  A <-P_(n-1)(1/z),  B <-P_(n-1)(-1/z). Try:
#GDg((a*A)^2,a,b,A,B,z,[1,z,0,0],2);
GDg:=proc(F,a,b,A,B,z,C,r) local gu,i:

if r mod 2=0 then
gu:=expand(subs({
a=C[1]*a+C[2]*b+C[3]*A+C[4]*B,
b=subs(z=-z,C[1])*a+subs(z=-z,C[2])*b+subs(z=-z,C[3])*A+subs(z=-z,C[4])*B,
A=subs(z=1/z,C[1])*A+subs(z=1/z,C[2])*B+subs(z=1/z,C[3])*a+subs(z=1/z,C[4])*b,
B=subs(z=-1/z,C[1])*A+subs(z=-1/z,C[2])*B+subs(z=-1/z,C[3])*a+subs(z=-1/z,C[4])*b},F));

gu:=add(coeff(gu,z,r*i)*z^i,i=-trunc(-ldegree1(gu,z)/r)-1..trunc(degree1(gu,z)/r)+1):


else

gu:=expand(subs({
a=C[1]*a+C[2]*b+C[3]*A+C[4]*B,
b=subs(z=-z,C[1])*b+subs(z=-z,C[2])*a+subs(z=-z,C[3])*B+subs(z=-z,C[4])*A,
A=subs(z=1/z,C[1])*A+subs(z=1/z,C[2])*B+subs(z=1/z,C[3])*a+subs(z=1/z,C[4])*b,
B=subs(z=-1/z,C[1])*B+subs(z=-1/z,C[2])*A+subs(z=-1/z,C[3])*b+subs(z=-1/z,C[4])*a},F));

gu:=add(coeff(gu,z,r*i)*z^i,i=-trunc(-ldegree1(gu,z)/r)-1..trunc(degree1(gu,z)/r)+1):


fi:
Canon(gu,a,b,A,B,z):

end:


#Sch(w,a,b,A,B,z,K,C,r): inputs a word w in symbols a,b,A,B,z, a positive integer K, 
#and a capsule C, of degree r, of length 4 whose entries are of the form e1/z+e2+e3*z (for e1,e2,e3 numbers)
#and  outputs
#the  transition  matrix, in SPARSE format, as a list, where the i-th entry {[c1,j1], ...}  
#means that  M[i,j1]=c1, and M[i,j]=0 if not mentioned.
#that enables to compute
#ordinary generating function,  or computing many terms,
# a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).
#where P_n(z) are the generalized Shapiro polynomials with capsule C, PnG(n,z,C);
#It uses a recursive scheme with at most K member. If none is found it returns FAIL.
#It also outputs the initial vector, and the "ledend", i.e. the list whose i-th item is the corresponding word 
#Try
#Sch(a*A,a,b,A,B,z,10,[1,z,0,0],2);
#Sch(a^2*A^2,a,b,A,B,z,10,[1,z,0,0],2);
#Sch(a^2*B^2*z,a,b,A,B,z,10,[1,z,-1,-1/z],2);
Sch:=proc(w,a,b,A,B,z,K,C,r) local S,StillToDo,adam,yeled,i,AlreadyDone,T1,V1,C1,C1i,co,i1,j1,T2,ka:
option remember:

S:={w}:

StillToDo:={w}:

co:=0:
AlreadyDone:={}:

while  nops(S)<=K and StillToDo<>{} do
 adam:=StillToDo[1]:
 co:=co+1:
 C1[co]:=adam:
 C1i[adam]:=co:
  
 AlreadyDone:=AlreadyDone union {adam}:
 StillToDo:= StillToDo minus {adam}:
 yeled:=Monos(GDg(adam,a,b,A,B,z,C,r),[a,b,A,B,z]):

 StillToDo:=StillToDo union ({seq(yeled[i][2],i=1..nops(yeled))} minus AlreadyDone):

 S:=S union {seq(yeled[i][2],i=1..nops(yeled))}:

   if degree1(adam,z)=0 then
    V1[co]:=1:
   else
    V1[co]:=0:
   fi:

 T1[adam]:=yeled:
od:

if StillToDo<>{} then
 RETURN(FAIL):
fi:

for i1 from 1 to nops(S) do
  ka:=T1[C1[i1]]:
  ka:={seq([ka[j1][1],C1i[ka[j1][2]]],j1=1..nops(ka))}:
 T2[i1]:=ka:
od:

[[seq(V1[i1],i1=1..nops(S))], [seq(T2[i1],i1=1..nops(S))],[seq(C1[i1],i1=1..nops(S))]]:

end:



#GFw(w,a,b,A,B,z,K,t,C,r): The generating function of the monomial w in
#in a,b,A,B,z with limit of scheme K, variable t, and capsule C of degree r
#Try:
#GFw(a*A,a,b,A,B,z,10,t,[1,z,0,0],2);
#GFw(a^2*A^2,a,b,A,B,z,10,t,[1,z,0,0],2);
#GFw(a^2*B^2*z,a,b,A,B,z,40,t,[1,z,z,1],2);
GFw:=proc(w,a,b,A,B,z,K,t,C,r) local lu,mu,ku:
option remember:
lu:=Sch(w,a,b,A,B,z,K,C,r):

 if lu=FAIL then
  print(` Make `, K, ` larger `):
  RETURN(FAIL):
 fi:

mu:=MtoSeq(lu[1],lu[2],2*nops(lu[3])+2 ):

ku:=GuessRec(mu):

if ku=FAIL then
 RETURN(FAIL):
fi:

factor(CtoR(ku,t)):
end:


#Eval1(n,w,a,b,A,B,z,C,r): given a non-negative integer n, a polynomial expression w, 
#in a,b,A,B, and z and a recurrence 
#templace C, of degree r, and positive integer n, evaluates what happens
#when a is replaced by f(z):=PnG(n,C,z,r), A by f(1/z), b by f(-z) and B by f(-1/z)
#and then the CONSTANT TERM is taken
#Try:
#Eval1(4,(a*A)^2,a,b,A,B,z,[1,z,0,0],2);
Eval1:=proc(n,w,a,b,A,B,z,C,r) local gu,mu:
gu:=PnG(n,C,z,r):
mu:=subs({a=gu,b=subs(z=-z,gu),A=subs(z=1/z,gu),B=subs(z=-1/z,gu)},w):
coeff(mu,z,0):
end:


#SidE(N,w,a,b,A,B,z,C,r): the first N+1 terms of the coeff. of the constant term in  
#Eval1G(n,w,z,a,b,A,B,z,C,r) (q.v.)
#Try:
#SidE(6,(a*A)^2/z^2,a,b,A,B,z,[1,z,0,0],2);
SidE:=proc(N,w,a,b,A,B,z,C,r) local n:

[seq(Eval1(n,w,a,b,A,B,z,C,r),n=0..N)]:

end:



#CheckGFw(w,a,b,A,B,z,K,t,C,N,r): 
# checks the N first terms of GFw(w,a,b,A,B,z,K,t,C,N,r) against
#direct calulations
#in a,b,A,B,z with limit of scheme K, variable t, and capsule C
#Try:
#CheckGFw(a^2*B^2*z,a,b,A,B,z,40,t,[1,z,z,1],6,2);
CheckGFw:=proc(w,a,b,A,B,z,K,t,C,N,r) local lu1,lu2,i,gu:
gu:=GFw(w,a,b,A,B,z,K,t,C,r):

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

gu:=taylor(gu,t=0,N+1):


lu1:=[seq(coeff(gu,t,i),i=0..N)]:

lu2:=SidE(N,w,a,b,A,B,z,C,r):

lu1,lu2,evalb(lu1=lu2):

end:




#GFwS(w,a,b,A,B,z,K,t,C,r): Verbose version of GFw(w,a,b,A,B,z,K,t,C,r) (q.v.)
#Try:
#GFwS(a*A,a,b,A,B,z,10,t,[1,z,0,0],2);
#GFwS(a^2*A^2,a,b,A,B,z,10,t,[1,z,0,0],2);
#GFwS(a^2*B^2*z,a,b,A,B,z,40,t,[1,z,z,1],2);
GFwS:=proc(w,a,b,A,B,z,K,t,C,r) local lu, gu,P,k,f,F,t0:
t0:=time():
lu:=GFw(w,a,b,A,B,z,K,t,C,r):
if lu=FAIL then
 RETURN(FAIL):
fi:

print(``):
gu:=subs({a=P[k](z),A=P[k](1/z),b=P(k)(-z),B=P(k)(-1/z)},w):
print(`A rational generating function for the Constant Term of `, gu):
print(``):
print(`where the sequence of Laurent polynomials`, P[k](z) , `is defined by the recurrence`):
print(``):
print(P[k](z)=C[1]*P[k-1](z^r)+C[2]*P[k-1](-z^r)+C[3]*P[k-1](1/z^r)+ C[4]*P[k-1](-1/z^r), P[0](z)=1):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, P[k](z) , `be the sequence of Laurent polynomials defined by the`):
print(` recurrence `):
print(``):
print(P[k](z)=C[1]*P[k-1](z^r)+C[2]*P[k-1](-z^r)+C[3]*P[k-1](1/z^r)+ C[4]*P[k-1](-1/z^r), P[0](z)=1):
print(``):
print(`subject to the initial condition`, P[0](z)=1 ):
print(``):
print(`For any non-negative integer, k, Let  f(k) be the constant term of `):
print(``):
print(gu):
print(``):
print(`Let F(t) be the (ordinary) generating function of the sequence f(k), in other words`):
print(``):
print(F(t)=Sum(f(k)*t^k,k=0..infinity)):
print(``):
print(`Then F(t) equals the rational function`):
print(``):
print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):
print(``):
print(`-----------------------`):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds. `):

end:






#SolSch(M,t): solves the scheme M with respect of t, only finds the value of the first component
SolSch:=proc(M,t) local eq,var,i,x,lu,lu1:
lu:=M[2]:
var:={seq(x[i],i=1..nops(lu))}:
eq:={seq(x[i]=M[1][i]+t*add(lu1[1]*x[lu1[2]],lu1 in lu[i]),i=1..nops(lu))}:
var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN(normal(subs(var,x[1]))):
fi:

end:


#GFwD(w,a,b,A,B,z,K,t,C,r): Like GFw(w,a,b,A,B,z,K,t,C,r): , but w/o guessing
#in a,b,A,B,z with limit of scheme K, variable t, and capsule C of degree r
#Try:
#GFwD(a*A,a,b,A,B,z,10,t,[1,z,0,0],2);
#GFwD(a^2*A^2,a,b,A,B,z,10,t,[1,z,0,0],2);
#GFwD(a^2*B^2*z,a,b,A,B,z,40,t,[1,z,z,1],2);
GFwD:=proc(w,a,b,A,B,z,K,t,C,r) local lu:
option remember:
lu:=Sch(w,a,b,A,B,z,K,C,r):
if lu=FAIL then
  RETURN(FAIL):
fi:

factor(SolSch(lu,t)):
end: