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

read `ShapiroData.txt`:

Digits:=30: 
print(`Created:  April-May 2016`):
print(` This is HaroldSilentShapiro.txt `):
print(`It is a package that accompanies 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, in the following url:`):
print(``):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/HaroldSilentShapiro.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 the paper`):
 print(` is available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/hss.html  .`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, 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 Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


with(combinat):

ezraS:=proc()

if args=NULL then
 print(` The Story procedures are:  GDv, HaroldSMv,HaroldSMvv,  MamarD, MamarDseq, MamarDseqG, MamarH, MamarHseq, MedioPv1, MedioPv, `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  AsyRat, BigAsy, BigAsySlow, BigMates,  BigMatrSlow,  Canon,Canon1,  CheckHarold, CheckNuMedio, `):
 print(`  CheckCP1, Coe, `):
  print(` CollectA, CollectB, CPF, CtoR, DHMnSeq, EvalExp, `):
 print(`  findrec GD, Gr, Gr1,  GuessA, GuessB, GuessRec,  Harold, HaroldF, HaroldM, HaroldP, HaroldSM,`):
 print(` K, Mates,  Monos, MtoSeq, Nake1, Pn, RtoS,  Sid, SidE, SidP`):
 print(`Shorashim, ShorashimH `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  Anr, AnrE,   BigGF, BigGFpf, BigMatr, CheckCP, CheckEV, Medio, MedioP, RS, SeqRS, Spc `):
 print(` `):


elif nops([args])=1 and op(1,[args])=AsyRat then
print(`AsyRat(R,t): inputs a rational function R of t, and outputs the asymptotic formula (with exponential error)`):
print(`it it exists, otherwise the approximate formula in the form [a, [C1,C2]] `):
print(`where C1 is the min and C2 is the max of the coefficients of t^i divided by a^i from 100 to 1000`):
print(`If successful, it returns `):
print(`[a,C], where the asymptotics is C*a^n *(1+ O(beta)^n)), for some beta strictly less than a`):


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

print(`Anr(n,r,t,J1): the beginning of the expression for the continued-fraction expression for rational function A_n^{(r)}(t)`):
print(`In includes all the largest J1 positive eignevalues and the smallest J1 negative eigenvalues`):
print(`Try: `):
print(`Anr(n,r,t,1); `):

elif nops([args])=1 and op(1,[args])=AnrE then
print(`AnrE(n,r,k,J1): Inputs SYMBOLS n,r, k, and a non-negatite integer J1, outputs `):
print(`the  expression, in n,r,k,  for the first J1+1 and last J1+1 terms  exact formula for`):
print(` the "Big Guys" approximation to  `):
print(`Int((P[k](exp(2*Pi*I*t))*P[k](exp(-2*Pi*I*t)))^(n-r)*(P[k](-exp(2*Pi*I*t))*P[k](-exp(-2*Pi*I*t)))^r,t=0..1)):`):
print(`Try:`):
print(`AnrE(n,r,k,1);`):

elif nops([args])=1 and op(1,[args])=AnrE1 then
print(`AnrE1(n,r1,k,J1): Inputs SYMBOLS n,k, a non-neg. integer r1, and a non-negatite integer J1, outputs `):
print(`the  expression, in n and k, for the first J1+1 and last J1+1 terms  exact formula for`):
print(` the "Big Guys" approximation to  `):
print(`Int((P[k](exp(2*Pi*I*t))*P[k](exp(-2*Pi*I*t)))^(n-r1)*(P[k](-exp(2*Pi*I*t))*P[k](-exp(-2*Pi*I*t)))^r1,t=0..1)):`):
print(`Try:`):
print(`AnrE1(n,r,k,1);`):

elif nops([args])=1 and op(1,[args])=BigAsySlow then
print(`BigAsySlow(n): Inputs a positive integer n, and outputs the list of constants C, such that`):
print(`CT_z (P[k](z)*P[k](1/z))^(n-i)*(P[k](-z)*P[k](-1/z))^i) is asymptotic to C[i+1]*(2^m)^k`):
print(`as k goes to infinity for i from 0 to n, under the "little guys don't count" assumption `):
print(` Try: `):
print(`BigAsySlow(4); `):

elif nops([args])=1 and op(1,[args])=BigAsy then
print(`BigAsy(n): A fast version of BigAsySlow(n) (q.v.), not using generating functions. `):
print(` Try: `):
print(`BigAsy(4); `):

elif nops([args])=1 and op(1,[args])=BigGF then
print(`BigGF(n,t): The asymptotic generating functions for the monomials a^(n-i)*A^(n-i)*b^i*B^i for i from 0 to n. Under the "little guys do not count"`):
print(`assumptiom. Try:`):
print(`BigGF(4,t);`):

elif nops([args])=1 and op(1,[args])=BigGFpf then
print(`BigGFpf(n): inputs a positive integer n, and outputs BigGF(n,t) in compact, partial-fraction format given by CPF(R,t) (q.v).`):
print(`Try:`):
print(`BigGFpf(6);`):

elif nops([args])=1 and op(1,[args])=BigMates then
print(`BigMates(n,a,b,A,B,z,K) : inputs a positive integer n and symbols a,b,A,B,z,k`):
print(`returns the set of monomials in {a,b,A,B,z} that have the same rate of growth as a^n*A^n`):
print(`in the scheme that is used to compute its generating function`):
print(`is to investigate a route to proving Saffari's conjecture. Try:`):
print(`BigMates(2,a,b,A,B,z,1000);`):

elif nops([args])=1 and op(1,[args])=BigMatr then
print(`BigMatr(n): The n+1 by n+1 transition matrix for the important monomials (a*A)^(n-i)*(b*B)^i, 0<=i<=n. Try:`):
print(`BigMatr(4);`):

elif nops([args])=1 and op(1,[args])=BigMatrSlow then
print(`BigMatrSlow(n):  slow (original) version of BigMatr(n) (q.v.) using the explicit formula. Try: `):
print(`BigMatrSlow(4);`):


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])=Canon1 then
print(`Canon1(w,a,b,A,B,z): the Canonical form of the word w. Try:`):
print(`Canon1(b*B,a,b,A,B,z);`):


elif nops([args])=1 and op(1,[args])=CheckEV then
print(`CheckEV(N): checks empirically that 2^n is indeed an eigenvalue of the matrix M^(n) and`):
print(`[seq(1/(binomial(n,r),r=0..n)] is an eigenvector for all n<=N. Try`);
print(` CheckEV(20); `):


elif nops([args])=1 and op(1,[args])=CheckCP1 then
print(`CheckCP1(n): checks empirically that the characteristic polynomial of BigMatr(n) is indeed`):
print(`mul((x-2^(n-4*j)*binomial(4*j,2*j),j=0..trunc(n/4)))*mul(x+2^(n-4*j-2)*binomial(4*j+2,2*j+1),j=0..trunc((n-2)/4))*x^(trunc(n+3)/2)`):
print(`Try: `):
print(`CheckCP1(5);`):

elif nops([args])=1 and op(1,[args])=CheckCP then
print(`CheckCP(N): checks empirically, for n from 2 to N, that the characteristic polynomial of BigMatr(n) is indeed`):
print(`mul((x-2^(n-4*j)*binomial(4*j,2*j),j=0..trunc(n/4)))*mul(x+2^(n-4*j-2)*binomial(4*j+2,2*j+1),j=0..trunc((n-2)/4))*x^(trunc(n+3)/2)`):
print(`Try: `):
print(`CheckCP(10);`):

elif nops([args])=1 and op(1,[args])=CheckHarold then
print(`CheckHarold(w,a,b,A,B,z,K,m): checks procedure Harold for the first m+1 terms. Try:`):
print(`CheckHarold(a*A,a,b,A,B,z,10,6);`):
print(`CheckHarold(a^2*A^2,a,b,A,B,z,10,6);`):
print(`CheckHarold(a^2*B^2*z,a,b,A,B,z,10,6);`):


elif nops([args])=1 and op(1,[args])=CheckNuMedio then
print(`CheckNuMedio(N,K): checks that the number of mediocre terms for n equals the sequence Ak(n)`):
print(`  mentioned in the website http://mrob.com/pub/math/MCS-10.html, MCS42828,`):
print(`defined by Ak(0)=-1, Ak(1)=0, Ak(k)=Ak(k-2)+k-1, that in fact equals`):
print(`n^2/4-5/8-3*(-1)^n/8`):
print(`Try:`):
print(`CheckNuMedio(8,1000);`):

elif nops([args])=1 and op(1,[args])=Coe then
print(`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:`):
print(`Coe(15*a^2*b^3*A*B*z,a,b,A,B,z);`):



elif nops([args])=1 and op(1,[args])=CollectA then
print(`CollectA(n1,j1,k): inputs   numeric n1 and j1 such that j1<=trunc(n1/4) outputs `):
print(`the polynomial in k for the coefficients of 1/(1-2^(n-4*j)*binomial(4*j,2*j)*t) in R_{n1,k} times `):
print(`binomial(n1,k)/(2^n1)*(n1+4*j1+1)!/n1!)`):
print(` Try: `):
print(` CollectA(24,0,k); `):

elif nops([args])=1 and op(1,[args])=CollectB then
print(`CollectB(n1,j1,k): inputs   numeric n1 and j1 such that j1<=trunc((n1-2)/4) and a variable k, outputs`):
print(`the polynomial in k for the coefficients of 1/(1+2^(n-4*j-2)*binomial(4*j+2,2*j+1)*t) in R_{n1,k} times`):
print(`binomial(n1,k)/(2^n1)*(n1+4*j1+1)!/n1!)`):
print(`Try: `):
print(` CollectB(24,1,k): `):

elif nops([args])=1 and op(1,[args])=CtoR then
print(`CtoR(S,t), the rational function, in t, whose coefficients`):
print(`are the members of the C-finite sequence S. For example, try:`):
print(`CtoR([[1,1],[1,1]],t);`):


elif nops([args])=1 and op(1,[args])=CPF then
print(`CPF(R,t): inputs a rational function whose denominator factorizes into a product of linear factors`):
print(`outputs its partial fraction decomosition in as a list [[C1,a1],[C2,a2], ..., [Ck,ak]] such that`):
print(`a1>a2>...>ak and R=C1/(1-a1*t)+ ...+ Ck/(1-ak*t). Try:`):
print(`CPF((t^2+1)/((1-2*t)*(1-3*t)*(1+2*t)),t);`):

elif nops([args])=1 and op(1,[args])=DHMnSeq then
print(`DHMnSeq(N,q): the sequence 2q-th moment of the Rudin-Shapiro polynomial P_n(z) from n=0 to n=N`):
print(`Should give the same output as  as SeqRS(q,q,0,0,10000,N);`):
print(`Try: `):
print(` DHMnSeq(20,3); `):


elif nops([args])=1 and op(1,[args])=EvalExp then
print(`EvalExp(n,w,a,b,A,B,z): given a non-negative integer n, a polynomial expression w, in a,b,A,B, and z and a positive integer n, evaluates what happens`):
print(`when a is replaced by Pn(z), A by Pn(1/z), b by Pn(-z), andB by Pn(-1/z), and then the CONSTANT TERM is taken`):
print(`Try: `):
print(`EvalExp(4,(a*A)^2,a,b,A,B,z);`):


elif nops([args])=1 and op(1,[args])=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(` the sequence f of degree DEGREE and order ORDER `):
print(` For example, try: findrec([seq(i,i=1..10)],0,2,n,N); `):



elif nops([args])=1 and op(1,[args])=GuessA then
print(`GuessA(n,k,j1):inputs  symbols n and k and a non-negative integer j1`):
print(`outputs a polynomial in k and n `):
print(`of degree 4*j1 in k,  that is the the numerator of 1/(1-binomial(4*j1,2*j1)*2^(n-4*j1)*t) in`):
print(`the partial fraction expansion of R_{n,k}(t) described in the paper `):
print(`Try:  `):
print(` GuessA(n,k,1); `):

elif nops([args])=1 and op(1,[args])=GuessB then
print(`GuessB(n,k,j1):inputs  symbols n and k and a non-negative integer j, and a positive integer N`):
print(`outputs a polynomial in k and n`):
print(`of degree 4*j in k let's call it`):
print(` B_j(n,k), such that the numerator of 1/(1+binomial(4*j+2,2*j+1)*2^(n-4*j-2)*t) is`):
print(`A_j(n,k)`):
print(`Try: `):
print(`GuessB(n,k,1);`):

elif nops([args])=1 and op(1,[args])=GuessRec then
print(`GuessRec(L): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant coefficients `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):


elif nops([args])=1 and op(1,[args])=GD then
print(`GD(F,a,b,A,B,z): The Going Down recurrence for an expression F.`):
print(`Let P_n(z) be the Rudin-Shapiro polynomials called Pn(n,z) in this program.`):
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 defining recurrence`):
print(` P_n(z)=P_{n-1}(z^2)+ z*P_n(-z^2), 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(` GD((a*A)^2,a,b,A,B,z); `):

elif nops([args])=1 and op(1,[args])=GDv then
print(`GDv(F,a,b,A,B,z): Verbose form of GD(F,a,b,A,B,z) (q.v.). Try: `):
print(` GDv((a*A)^2,a,b,A,B,z); `):

elif nops([args])=1 and op(1,[args])=Gr then
print(`Gr(L,n) guesses a rational function in n that fits the list L`):
print(`Try: `):
print(`Gr([seq((i+1)/(i+2),i=1..10)],n); `):

elif nops([args])=1 and op(1,[args])=Gr1 then
print(`Gr1(L,n,d1,d2) guesses a rational function of numerator degree d1 and denominator degree d2, for the list L`):
print(`Try: `):
print(`Gr1([seq((i+1)/(i+2),i=1..10)],n,1,1); `):

elif nops([args])=1 and op(1,[args])=Harold then
print(`Harold(w,a,b,A,B,z,K,t): inputs a word w in symbols a,b,A,B,z, a positive integer K, and  variable t outputs`):
print(`the ordiidnary generating function,  in the variable t, whose coefficient of t^n is  for the constant term of w where `):
print(` a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).`):
print(`It uses a recursive scheme with at most K member. If none os found it returns FAIL.`):
print(`Try: `):
print(` Harold(a*A,a,b,A,B,z,10,t); `):
print(`Harold(a^2*A^2,a,b,A,B,z,10,t);`):
print(`Harold(a^2*B^2*z,a,b,A,B,z,10,t);`):


elif nops([args])=1 and op(1,[args])=HaroldF then
print(`HaroldF(w,a,b,A,B,z,K,t): A (hopefully) faster version of Harold(w,a,b,A,B,z,K,t): (q.v.)`):
print(`it does it via HaroldSM(w,a,b,A,B,z,K) (q.v.), generating sufficiently many terms and then`):
print(`(rigorously!) guessing the linear recurrence equation with constant coefficients satisfied by`):
print(`the sequence, and then finding the rational function expression for the generating function`):
print(`Try: `):
print(`HaroldF(a*A,a,b,A,B,z,10,t);`):
print(`HaroldF(a^2*A^2,a,b,A,B,z,10,t);`):
print(`HaroldF(a^2*B^2*z,a,b,A,B,z,10,t);`):

elif nops([args])=1 and op(1,[args])=HaroldM then
print(`HaroldM(w,a,b,A,B,z,K): inputs a word w in symbols a,b,A,B,z, a positive integer K, and  outputs`):
print(`the  transition  matrix, given as a list of lists, that enables to compute`):
print(`ordinary generating function,  or many terms, or the eigenvalues`):
print(` a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).`):
print(`It uses a recursive scheme with at most K member. If none os found it returns FAIL.`):
print(`It also outputs the initial vector, and the "ledend", i.e. the list whose i-th item is the corresponding word `):
print(`Try: `):
print(`HaroldM(a*A,a,b,A,B,z,10);`):
print(`HaroldM(a^2*A^2,a,b,A,B,z,10);`):
print(`HaroldM(a^2*B^2*z,a,b,A,B,z,10);`):

elif nops([args])=1 and op(1,[args])=HaroldP then
print(`HaroldP(w,a,b,A,B,z,K,t): inputs a word w in symbols a,b,A,B,z, a positive integer K, and  variable t outputs`):
print(`the ordinary generating function,  in the variable t, whose coefficient of t^n is  for the constant term of w where `):
print(` a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).`):
print(` It also returns, as  extra gifts, the gernerating functions for other words `):
print(`It uses a recursive scheme with at most K member. If none os found it returns FAIL.`):
print(`Try: `):
print(` HaroldP(a*A,a,b,A,B,z,10,t); `):
print(`HaroldP(a^2*A^2,a,b,A,B,z,10,t);`):
print(`HaroldP(a^2*B^2*z,a,b,A,B,z,10,t);`):



elif nops([args])=1 and op(1,[args])=HaroldSM then
print(`HaroldSM(w,a,b,A,B,z,K): inputs a word w in symbols a,b,A,B,z, a positive integer K, 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(`It uses a recursive scheme with at most K member. If none os found it returns FAIL.`):
print(`It also outputs the initial vector, and the "ledend", i.e. the list whose i-th item is the corresponding word `):
print(`Try:`):
print(`HaroldSM(a*A,a,b,A,B,z,10);`):
print(`HaroldSM(a^2*A^2,a,b,A,B,z,10);`):
print(`HaroldSM(a^2*B^2*z,a,b,A,B,z,10);`):

elif nops([args])=1 and op(1,[args])=HaroldSMv then
print(`HaroldSMv(w,a,b,A,B,z,K,t): Verbose version of HaroldSM(w,a,b,A,B,z,K) (q.v.)`):
print(`Try: `):
print(`HaroldSMv(a^2*A^2,a,b,A,B,z,100,t); `):


elif nops([args])=1 and op(1,[args])=HaroldSMvv then
print(`HaroldSMvv(w,a,b,A,B,z,K,t): A VERY Verbose version of HaroldSM(w,a,b,A,B,z,K) (q.v.)`):
print(`it includes explanation of the going down formula . `):
print(`Try: `):
print(`HaroldSMv(a^2*A^2,a,b,A,B,z,100,t); `):


elif nops([args])=1 and op(1,[args])=K then
print(`K(n,k,r) Inputs positive integers 0<=k<=n, and 0<=r<=n: and outputs (add((-1)^a*binomial(n-k,r-a)*binomial(k,a),a=0..r))^2:`):
print(`Try:`):
print(`K(6,3,2);`):


elif nops([args])=1 and op(1,[args])=MamarD then
print(`MamarD(N,K,t,EE): inputs a positive integer N <=10, and a positive integer K and a variable name t`):
print(`and EE that is either 0 or 1 (if EE=0 there is no error estimate. if EE=1 then there is)`):
print(`outputs an article with explicit expressions as rational functions for the sequences`):
print(`of the first N even moments. It also gives the first K+1 terms of the actual sequence,`):
print(`and asymptotics. Try:`):
print(`MamarD(4,30,t,1):`):
print(`MamarD(10,100,t,0):`):

elif nops([args])=1 and op(1,[args])=MamarDseq then
print(`MamarDseq(N,L,K): inputs a positive integer N <=10, and a positive integer K (take it to be 10000), `):
print(`(it is for internal purposes, for the upper bound for the size of the scheme (see main (human) article)).`):
print(`Outputs an article that lists the first L+1 terms (from k=0 to k=L)of `):
print(` the even Moments, from 2 to 2N, of the Rudin-Shapiro polynomials. With sufficiently large L, `):
print(` this data can be used to find the generating function, if desired, but for large N it may take a long time. `):
print(` Try: `):
print(` MamarDseq(4,100,1000): `):

elif nops([args])=1 and op(1,[args])=MamarDseqG then
print(`MamarDseqG(N,L,K): inputs a positive integer N <=10, and a positive integer K (take it to be 10000), `):
print(`(it is for internal purposes, for the upper bound for the size of the scheme (see main (human) article)).`):
print(`Outputs an article that lists the first L+1 terms (from k=0 to k=L)of`):
print(`the sequences CT_z (P[k](z)*P[k](1/z))^(n-i)*(P[k](-z)*P[k](-1/z))^i) is asymptotic to C[i+1]*(2^m)^k`):
print(` n from 0 to N and i from 0 to n, where P[k](z) is the k-th Rudin-Shapiro polynomial. With sufficiently large L, `):
print(`this data can be used to find the generating function, if desired, but for large N it may take a long time. `):
print(` Try: `):
print(` MamarDseqG(4,100,1000): `):

elif nops([args])=1 and op(1,[args])=MamarH then
print(`MamarH(N,K,t): inputs a positive integer N, and a positive integer K and a variable name t`):
print(`outputs an article with explicit expressions as rational functions, in t, for the sequences`):
print(`of M[m,n](k)=Int(P_k(exp(2*Pi*I*theta)^m*P_k(exp(-2*Pi*I*theta)^n,theta=0..1) , k=0..infinity`):
print(`for 1<=m<n<=N, and testing for these values a conjecture of Hugh Montgomery made in Oberwolfach`):
print(`20 October – 26 October 2013, published on-line, and available from the url`):
print(`https://www.mfo.de/document/1343/OWR_2013_51.pdf (see pp. 67-68)`):
print(`It also gives the first K+1 terms of the actual sequence,`):
print(`and confirms Hugh Montgomery's conjecture for these m and n. Try`):
print(`MamarH(4,100,t):`):


elif nops([args])=1 and op(1,[args])=MamarHseq then
print(`MamarHseq(N,K,L): inputs a positive integer N, and a positive integer K (take it to be 10000), `):
print(`and a positive integer L`):
print(`outputs an article with the first L+1 term of the sequence`):
print(`M[m,n](k):=Int(P_k(exp(2*Pi*I*theta)^m*P_k(exp(-2*Pi*I*theta)^n,theta=0..1) , k=0..infinity`):
print(`for 1<=m<n<=N, and testing that the absolute  values satisfy the upper bound  conjectured by Hugh Montgomery made in Oberwolfach`):
print(`20 October – 26 October 2013, published on-line, and available from the url `):
print(`https://www.mfo.de/document/1343/OWR_2013_51.pdf (see pp. 67-68). Try:`):
print(`MamarHseq(4,10000,40);`):

elif nops([args])=1 and op(1,[args])=Mates then
print(`Mates(w,a,b,A,B,z): given a word in a,b,A,B,z where a=P_n(z), b=P_n(-z), A=P_n(1/z),B=P_n(-1/z) finds`):
print(`all equivalent monomials that trivially give the same Constant Term, under the action of the group z->{z,-z,1/z,-1/z}`):
print(`Try: `):
print(`Mates(a*A,a,b,A,B,z); `):

elif nops([args])=1 and op(1,[args])=Medio then
print(`Medio(n,a,b,A,B,z,K): inputs a positive integer n, symbols a,b,A,B,z, and a large positive integer K`):
print(`Let the monomials a^i*A^i*b^(n-i)*B^(n-i), i<=n/2 be called important.`):
print(`It outputs the set of unimportant monomials that have important childen. Try:`):
print(`Medio(3,a,b,A,B,z,1000);`):

elif nops([args])=1 and op(1,[args])=MedioP then
print(`MedioP(n,a,b,A,B,z,K): inputs a positive integer n, symbols a,b,A,B,z, and a large positive integer K`):
print(`Let the monomials a^i*A^i*b^(n-i)*B^(n-i), i<=n/2 be called important.`):
print(`It outputs the set of unimportant monomials that have important childen. Followed by the same set but with their`):
print(` expression in terms of important monomials for each.`):
print(`Try: `):
print(`MedioP(3,a,b,A,B,z,1000);`):

elif nops([args])=1 and op(1,[args])=MedioPv1 then
print(`MedioPv1(n,K):  Inputs a positive integer n, and a large positive integer K (for the upper bound for the scheme)`):
print(` tests the Important Monomial Hypothesis for the`):
print(`Saffari conjecture for Int(|P_n(exp(2*Pi*t))|^(2*n),t=0..1)`):
print(` that all the Non-Important Children of important monomials `):
print(`Try:` ):
print(`MedioPv1(3,1000);`):

elif nops([args])=1 and op(1,[args])=MedioPv then
print(`MedioPv(N,K):  Inputs a positive integer N, and a large positive integer K (for the upper bound for the scheme)`):
print(` tests the Important Monomial Hypothesis for the`):
print(`Saffari conjecture for Int(|P_n(exp(2*Pi*t))|^(2*n),t=0..1)`):
print(` that all the Non-Important Children of important monomials `):
print(`for all n between 1 and N`):
print(`Try:` ):
print(`MedioPv(10,1000);`):

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])=MtoSeq then
print(`MtoSeq(V1,M,N): inputs an initial vector V1 a sparse matrix M, and a positive integer N`):
print(`outputs the first N+1 terms, starting at n=0 of the sequence M^n(v1)[1]`):
print(`Try:`):
print(`MtoSeq([1],[{[2,1]}],10);`):



elif nops([args])=1 and op(1,[args])=Nake1 then
print(`Nake1(f,t): inputs a reciprocal of an affine linear function and outputs the pair [C,alpha] such that it equals`):
print(`C/(1-alpha*t). Try:`):
print(`Nake1(4/(11+6*t),t);`):

elif nops([args])=1 and op(1,[args])=Pn then
print(`Pn(n,z): the n-th Shapiro polynomial. Try: `):
print(` Pn(4,z); `):


elif nops([args])=1 and op(1,[args])=RS then
print(`RS(m,n,t,K): Inputs non-negative integers m and n, and a variable t, (and a large pos. integer K),outputs`):
print(`the (ordinary) generating function, in t, of the sequence `):
print(`Int(P_k(e(theta))^m*(P_k(e(-theta))^n,theta=0..1)`):
print(`(i.e. the rational function whose coefficient of t^k is) is the above.`):
print(`K is the upper size of the scheme. Try:`):
print(`RS(2,2,t,1000);`):
print(`RS(2,5,t,1000);`):

elif nops([args])=1 and op(1,[args])=RtoS then
print(`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:`):
print(`RtoS(1/(1-t-t^3),t,30);`):


elif nops([args])=1 and op(1,[args])=SeqRS then
print(`SeqRS(m,n,m1,n1,K,N): Inputs non-negative integers m and n, and a large pos. integer K,`):
print(`and a positive integer N, outputs`):
print(`the first N terms, starting at k=0`):
print(` Int(P_k(e(theta))^m*(P_k(e(-theta))^n*P_k(-e(theta))^m1*(P_k(-e(-theta))^n1 ,theta=0..1) `):
print(` SeqRS(2,2,0,0,1000,100); `):
print(` SeqRS(5,5,0,0,1000,100); `):


elif nops([args])=1 and op(1,[args])=Sid then
print(`Sid(N,r,s): the first N+1 terms of the coeff. of z^s in (Pn(n,z)*Pn(n,1/z))^r.`):
print(`Try: `):
print( ` Sid(6,2,0); `):

elif nops([args])=1 and op(1,[args])=SidE then
print(`SidE(N,w,a,b,A,B,z): the first N+1 terms of the constant term in  EvalExp(n,w,a,b,A,B,z) (q.v.)`):
print(`when a is replaced by Pn(z), A by Pn(1/z), b by Pn(-z), and B by Pn(-1/z), and then the coefficient of z^s is taken.`):
print(`Try: `):
print(`SidE(6,(a*A)^2,a,b,A,B,z);`):

elif nops([args])=1 and op(1,[args])=SidP then
print(`SidP(N,r,s,z): the first N+1 terms of the central part, from z^(-s) to z^s  in (Pn(n,z)*Pn(n,1/z))^r.`):
print(` Try: `):
print(` SidP(6,2,2,z); `):

elif nops([args])=1 and op(1,[args])=Shorashim then
print(`Shorashim(n): all the roots that show up in BigGFpf(n). Try: `):
print(`Shorashim(8); `):

elif nops([args])=1 and op(1,[args])=ShorashimH then
print(`ShorashimH(n): the conjectured roots that show up in BigGFpf(n).`):
print(`Should be the same as Shorashim(n); Try:`):
print(`ShorashimH(8); `):

elif nops([args])=1 and op(1,[args])=Spc then
print(`Spc(k,t): If k is from 1 to 10, the pre-computed generating function for the (2k)-th moment of the Rudin-Shapiro polynomial`):
print(`Same output as RS(k,k,t,10000) but, much faster, since precomputed `):
print(`You need to also download, into the same directory, as this Maple package, the data file: ShapiroData.txt `):
print(`available from`):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/ShapiroData.txt  .`):
print(``):
print(`Try:`):
print(` Spc(5,t); `):

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

ka:=a,b,A,B,z:

##Start from Findrec and GuessRat
#Gr1(L,n,d1,d2) guesses a rational function for the list L
Gr1:=proc(L,n,d1,d2) local a,b,eq,var,R,i:

if nops(L)<=d1+d2+2 then
 print(`List must have at least `, d1+d2+2, `terms `):
 RETURN(FAIL):
fi:

R:=(n^d1+add(a[i]*n^i,i=0..d1-1))/add(b[i]*n^i,i=0..d2):
var:={seq(a[i],i=0..d1-1),seq(b[i],i=0..d2)}:
eq:={seq(numer(subs(n=i,R)-L[i]),i=1..d1+d2+2)}:




var:=solve(eq,var):

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


R:=subs(var,R):

if {seq(subs(n=i,R)-L[i],i=1..nops(L))}<>{0} then
  RETURN(FAIL):
fi:

R:

end:


##Start from Findrec and GuessRat
#Gr(L,n) guesses a rational function for the list L of sums of degrees<=d
Gr:=proc(L,n) local C,gu,d1:

 for C from 0 to nops(L)-4 do
  for d1 from 0 to C do
     gu:=Gr1(L,n,d1,C-d1):
   if gu<>FAIL then
    RETURN(factor(gu)):
   fi:
 od:
od:
FAIL:
end:
   

#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:

if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:

###end from Findrec
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:


#Pn(n,z): the n-th Shapiro polynomial. Try: Pn(4,z);
Pn:=proc(n,z) local gu:
option remember:

if n=0 then
 RETURN(1):
else
gu:=Pn(n-1,z):
RETURN(expand( subs(z=z^2,gu) +z*subs(z=-z^2,gu) )):
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


#Sid(N,r,s): the first N+1 terms of the coefficient of z^s in (Pn(n,z)*Pn(n,1/z))^r.
#Try:
#Sid(6,2,0);
Sid:=proc(N,r,s) local n,z, gu,lu:

lu:=[]:

for n from 0 to N do
gu:=Pn(n,z):
gu:=expand((gu*subs(z=1/z,gu))^r):
lu:=[op(lu),coeff(gu,z,s)]:
od:
lu:

end:

#SidP(N,r,s,z): the first N+1 terms of the central part, from z^(-s) to z^s  in (Pn(n,z)*Pn(n,1/z))^r.
#Try:
#SidP(6,2,2,z);
SidP:=proc(N,r,s,z) local n, gu,lu,i:

lu:=[]:

for n from 0 to N do
gu:=Pn(n,z):
gu:=expand((gu*subs(z=1/z,gu))^r):
lu:=[op(lu),add(coeff(gu,z,i)*z^i,i=-s..s)]:
od:
lu:

end:




#EvalExp(n,w,a,b,A,B,z): given a non-negative integer n, a polynomial expression w, in a,b,A,B, and z and a positive integer n, evaluates what happens
#when a is replaced by Pn(z), A by Pn(1/z), b by Pn(-z), andB by Pn(-1/z), and then the CONSTANT TERM is taken
#Try:
#EvalExp(4,(a*A)^2,a,b,A,B,z);
EvalExp:=proc(n,w,a,b,A,B,z) local gu,mu:
gu:=Pn(n,z):
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): the first N+1 terms of the coeff. of the constant term in  EvalExp(n,w,z,a,b,A,B,s) (q.v.)
#Try:
#SidE(6,(a*A)^2/z^2,a,b,A,B,z);
SidE:=proc(N,w,a,b,A,B,z) local n:

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

end:


#Mates(w,a,b,A,B,z): given a word in a,b,A,B,z where a=P_n(z), b=P_n(-z), A=P_n(1/z),B=P_n(-1/z) finds
#all equivalent monomials that trivially give the same Constant Term, undel the group z->{z,-z,1/z,-1/z}
#Try:
#Mates(a*A,a,b,A,B,z);
Mates:=proc(w,a,b,A,B,z)
{w,subs({z=-z,a=b,b=a,A=B,B=A},w),subs({z=1/z,a=A,A=a,b=B,B=b},w),subs({z=-1/z,a=B,B=a,b=A,A=b},w)}:
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:





#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:


#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:




#GD(F,a,b,A,B,z): The Going Down recurrence for an expression F.
#Let P_n(z) be the Rudin-Shapiro polynomials called Pn(n,z) in this program.
#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 defining recurrence
# P_n(z)=P_{n-1}(z^2)+ z*P_n(-z^2), 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:
#GD((a*A)^2,a,b,A,B,z);
GD:=proc(F,a,b,A,B,z) local gu,i:

gu:=expand(subs({a=a+z*b,b=a-z*b, A=A+B/z, B=A-B/z},F));

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

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

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:

#CheckHarold(w,a,b,A,B,z,K,m): checks procedure Harold for the first m+1 terms. Try:
#CheckHarold(a*A,a,b,A,B,z,10,6);
#CheckHarold(a^2*A^2,a,b,A,B,z,10,6);
#CheckHarold(a^2*B^2*z,a,b,A,B,z,10,6);

CheckHarold:=proc(w,a,b,A,B,z,K,m) local t,gu,mu1,mu2,i:
gu:=Harold(w,a,b,A,B,z,K,t):

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

mu1:=SidE(m,w,a,b,A,B,z):

mu2:=[seq(coeff(taylor(gu,t=0,m+1),t,i),i=0..m)]:

if mu1=mu2 then
  RETURN(true);
else
 RETURN(false):
fi:

end:


#Harold(w,a,b,A,B,z,K,t): inputs a word w in symbols a,b,A,B,z, a positive integer K, and  variable t outputs
#the ordinary generating function,  in the variable t, whose coefficient of t^n is  for the constant term of w where 
# a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).
#It uses a recursive scheme with at most K member. If none os found it returns FAIL.
#Try
#Harold(a*A,a,b,A,B,z,10,t);
#Harold(a^2*A^2,a,b,A,B,z,10,t);
#Harold(a^2*B^2*z,a,b,A,B,z,10,t);
Harold:=proc(w,a,b,A,B,z,K,t) local eq,eq1,var,S,f,StillToDo,adam,yeled,c,i,AlreadyDone,var1,lu:
option remember:
var:={}:
eq:={}:

S:={w}:

StillToDo:={w}:

AlreadyDone:={}:

while  nops(S)<=K and StillToDo<>{} do
 adam:=StillToDo[1]:
 AlreadyDone:=AlreadyDone union {adam}:
  var:=var union {f[adam]}:
  

 StillToDo:= StillToDo minus {adam}:
 yeled:=Monos(GD(adam,a,b,A,B,z),[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
    c:=1:
   else
    c:=0:
   fi:

 eq1:=f[adam]-c-t*add(yeled[i][1]*f[yeled[i][2]],i=1..nops(yeled)):
 eq:=eq union {eq1}:
od:



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



var1:=solve(eq,var):



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

lu:=normal(subs(var1,eq)):
	  
if lu<>{0} then
  RETURN(FAIL):
fi:

factor(subs(var1,f[w])):

end:





 

#HaroldP(w,a,b,A,B,z,K,t): inputs a word w in symbols a,b,A,B,z, a positive integer K, and  variable t outputs
#the ordinary generating function,  in the variable t, whose coefficient of t^n is  for the constant term of w where 
# a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).
#It uses a recursive scheme with at most K member. If none os found it returns FAIL.
#It also returns the expressions for other words needed for the proof
#Try
#HaroldP(a*A,a,b,A,B,z,10,t);
#HaroldP(a^2*A^2,a,b,A,B,z,10,t);
#HaroldP(a^2*B^2*z,a,b,A,B,z,10,t);
HaroldP:=proc(w,a,b,A,B,z,K,t) local eq,eq1,var,S,f,StillToDo,adam,yeled,c,i,AlreadyDone,var1,lu,ka,w1:
option remember:
var:={}:
eq:={}:

S:={w}:

StillToDo:={w}:

AlreadyDone:={}:

while  nops(S)<=K and StillToDo<>{} do
 adam:=StillToDo[1]:
 AlreadyDone:=AlreadyDone union {adam}:
  var:=var union {f[adam]}:
  

 StillToDo:= StillToDo minus {adam}:
 yeled:=Monos(GD(adam,a,b,A,B,z),[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
    c:=1:
   else
    c:=0:
   fi:

 eq1:=f[adam]-c-t*add(yeled[i][1]*f[yeled[i][2]],i=1..nops(yeled)):
 eq:=eq union {eq1}:
od:



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



var1:=solve(eq,var):



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

lu:=normal(subs(var1,eq)):
	  
if lu<>{0} then
  RETURN(FAIL):
fi:

ka:=factor(subs(var1,f[w])):

S:=S minus {w}:

lu:={}:

for w1 in S do
lu:=lu union {[w1,factor(subs(var1,f[w1]))]}:
od:

ka,lu:
end:


#HaroldM(w,a,b,A,B,z,K): inputs a word w in symbols a,b,A,B,z, a positive integer K, and  outputs
#the  transition  matrix, given as a list of lists, that enables to compute
#ordinary generating function,
# a <-P_n(z),  b <-P_n(-z),  A <-P_n(1/z),  B <-P_n(-1/z).
#It uses a recursive scheme with at most K member. If none os 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
#HaroldM(a*A,a,b,A,B,z,10);
#HaroldM(a^2*A^2,a,b,A,B,z,10);
#HaroldM(a^2*B^2*z,a,b,A,B,z,10);
HaroldM:=proc(w,a,b,A,B,z,K) local S,StillToDo,adam,yeled,i,AlreadyDone,T1,V1,C1,co,i1,j1,T2,shap:
option remember:
shap:={}:
S:={w}:

StillToDo:={w}:

co:=0:
AlreadyDone:={}:

while  nops(S)<=K and StillToDo<>{} do
 adam:=StillToDo[1]:
 co:=co+1:
 C1[co]:=adam:
  
 AlreadyDone:=AlreadyDone union {adam}:
 StillToDo:= StillToDo minus {adam}:
 yeled:=Monos(GD(adam,a,b,A,B,z),[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:

  for i1 from 1 to nops(yeled) do
    T1[adam,yeled[i1][2]]:=yeled[i1][1]:
    shap:=shap union {[adam,yeled[i1][2]]}:
  od:
od:

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

for i1 from 1 to nops(S) do
 for j1 from 1 to nops(S) do
   if member([C1[i1],C1[j1]],shap) then
    T2[i1,j1]:=T1[C1[i1],C1[j1]]:
   else
     T2[i1,j1]:=0:
  fi:
 od:
od:

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

end:


#HaroldSM(w,a,b,A,B,z,K): inputs a word w in symbols a,b,A,B,z, a positive integer K, 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).
#It uses a recursive scheme with at most K member. If none os 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
#HaroldSM(a*A,a,b,A,B,z,10);
#HaroldSM(a^2*A^2,a,b,A,B,z,10);
#HaroldSM(a^2*B^2*z,a,b,A,B,z,10);
HaroldSM:=proc(w,a,b,A,B,z,K) 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(GD(adam,a,b,A,B,z),[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:


#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:



#HaroldF(w,a,b,A,B,z,K,t): A (hopefully) faster version of Harold(w,a,b,A,B,z,K,t): (q.v.)
#it does it via HaroldSM(w,a,b,A,B,z,K) (q.v.), generating sufficiently many terms and then
#(rigorously!) guessing the linear recurrence equation with constant coefficients satisfied by
#the sequence, and then finding the rational function expression for the generating function
#Try
#HaroldF(a*A,a,b,A,B,z,10,t);
#HaroldF(a^2*A^2,a,b,A,B,z,10,t);
#HaroldF(a^2*B^2*z,a,b,A,B,z,10,t);
HaroldF:=proc(w,a,b,A,B,z,K,t) local lu,mu,ku:
lu:=HaroldSM(w,a,b,A,B,z,K):

 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:

CtoR(ku,t):
end:




#RS(m,n,t,K): Inputs non-negative integers m and n, and a variable t, (and a large pos. integer K),outputs
#the (ordinary) generating function, in t, of the sequence 
#Int(P_k(e(theta))^m*(P_k(e(-theta))^n,theta=0..1)
#(i.e. the rational function whose coefficient of t^k is) is the above.
#K is the upper size of the scheme. Try:
#RS(2,2,t,1000);
#RS(2,5,t,1000);

RS:=proc(m,n,t,K) local a,b,A,B,z:

HaroldF(a^m*A^n,a,b,A,B,z,K,t):

end:


#AsyRat(R,t): inputs a rational function R of t, and outputs the asymptotic formula (with exponential error)
#it it exists, otherwise the approximate formula in the form [a, [C1,C2]] 
#where C1 is the min and C2 is the max of the coefficients of t^i divided by a^i from 100 to 1000
#If successful, it returns 
#[C,a], where the asymptotics is C*a^n *(1+ O(beta)^n)), for some beta strictly less than a

#AsyRat(R,t): inputs a rational function R of t, and outputs the asymptotic formula (with exponential error)
#it it exists, otherwise the approximate formula in the form [FAIL, a^n]. 
#If successful, it returns 
#[C*alpha^n ,beta], where the asymptotics is C*alpha^n *(1+ O(beta)^n))
AsyRat:=proc(R,t) local vu,i,lu,a,c:
lu:=denom(R):

vu:=[solve(lu,t)]:

a:=1/min(seq(evalf(abs(vu[i])),i=1..nops(vu))):

lu:=[seq(coeff(taylor(R,t=0,1001),t,i),i=990..1000)]:


if abs(lu[11]/lu[10]-a)<1/10^12 and abs(lu[10]/lu[9]-a)<1/10^12   and abs(lu[9]/lu[8]-a)<1/10^12  
then
c:=identify(evalf(lu[11]/a^1000,10)):
if evalf(abs(a-round(a)))<1/10^10 then
   a:=round(a):
 else
  a:=identify(a):
fi:

 RETURN([a,[c]]):
else

lu:=evalf([seq(abs(coeff(taylor(R,t=0,1001),t,i))/a^i,i=100..1000)]):
 RETURN([a,[min(op(lu)), max(op(lu))] ]):
fi:

end:
 

   



#SeqRS(m,n,m1,n1,K,N): Inputs non-negative integers m,n, m1,n1, and a large pos. integer K,
#and a positive integer N, outputs
#the first N terms, starting at k=0
#Int(P_k(e(theta))^m*(P_k(e(-theta))^n*P_k(-e(theta))^m1*(P_k(-e(-theta))^n1,theta=0..1)
#SeqRS(2,2,0,0,1000,100);
#SeqRS(5,5,0,0,1000,100);
SeqRS:=proc(m,n,m1,n1,K,N) local a,b,A,B,z,lu:

lu:=HaroldSM(a^m*A^n*b^m1*B^n1,a,b,A,B,z,K):
MtoSeq(lu[1],lu[2],N):

end:


#MamarD(N,K,t,EE): inputs a positive integer N <=10, and a positive integer K and a variable name t
#and EE that is either 0 or 1 (if EE=0 there is no error estimate. if EE=1 then there is):
#outputs an article with explicit expressions as rational functions for the sequences
#of the first N even moments. It also gives the first K terms of the actual sequence,
#and asymptotics. Try:
#MamarD(4,100,t):
MamarD:=proc(N,K,t,EE)  local t0,P,gu,lu,i,A,k,m,f,gu1,c,mu,z:
t0:=time():
if N>10 then
 print(`The first argument must be <=10`):
 RETURN(FAIL):
fi:
print(`Explicit Generating Functions, Asymptotics, and More for the first`, N, `even moments of the Rudin-Shapiro polynomials`):
print(``):
print(`By Shalosh B. Ekhad (ShaloshBEkhad@gmail.com )`):
print(``):
print(` Let `, P[k](z), ` be the Rudin-Shapiro polynomial as described for example, in:`):
print(``):
print(` https://en.wikipedia.org/wiki/Shapiro_polynomials .` ):
print(``):
print(`The best way to define them is via the recurrence `):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(``):
print(`and the initial condition`):
print(` P[0](z)=1 `):
print(``):

print(`Let m and k be a positive integers, and define the even moment`):
print(``):
print(A[k,m]=Int(abs(P[k](exp(2*Pi*I*t)))^(2*m),t=0..1)):
print(`It was first proved by Christoph Doche and Laurent Habsieger: `):
print(` J. of Fourier Analysis and Applications  Vol. 10 Issue 5, 497-505, (2004)`):
print(``):
print(`available from:  http://web.science.mq.edu.au/~doche/049.pdf `):
print(``):
print(`that for each specific positive integer m, the generating function`):
print(``):
print(f[m](t)=Sum(A[k,m]*t^k,k=0..infinity)):
print(``):
print(`is always a rational function of t. This fact seems also to follow from our approach.`):
print(``):
print(`Since the original Pari-program referred to by these`):
print(`authors is no longer available, and neither is the output, and since these polynomials, and associated sequences are so important`):
print(`we decided to recompute them, and also confirm, for m from 1 to`, N, `the famous conjecture  of Saffari that  for any specific m`):
print( A[k,m], ` is asymptotic to`,  (2^m/(m+1))*(2^m)^k, `as k goes to infinity. `):
print(``):
print(`For the sake of our beloved OEIS, we also supply the first`, K+1, `terms of the actual sequences`):
print(``):
print(`---------------------------------------------------------------------`):
print(``):

for m from 1 to N do
gu:=Spc(m,t):
lu:=taylor(gu,t=0,K+1):
lu:=[seq(coeff(lu,t,i),i=0..K)]:
print(``):
print(`-------------------------------------------------`):
print(``):
print(`Theorem Number`, m  ):
print(``):
lprint(f[m](t)=gu):
print(``):

gu1:=denom(gu):
gu1:=normal(gu1/(1-2^m*t)):
if degree(denom(gu1),t)=0 then
print(`Note that the smallest root of the denominator is`, 1/2^m):
c:=subs(t=1/2^m,normal((1-2^m*t)*gu)):
print(`and the residue there is`, c, `that implies that the asymptotics is indeed`, c*(2^m)^k,`as conjectured by Saffari.`):
fi:

if EE=1 then
gu:=normal((1-2^m*t)*gu):
mu:= [solve(denom(gu),t)]:
mu:=min(seq(evalf(abs(mu[i])),i=1..nops(mu))):
print(`The largest modulo of the second largest root is`, mu):
 if mu>1/2^m then
  print(`that implies that the error term is exponentially small, of the order of`, (1/2^m/mu)^k ):
 fi:
fi:

print(`The first`, K+1, `terms are `):
print(``):
print(lu):


od:

print(``):
print(`---------------------------------------------------------------`):
print(`This concludes this article that took`, time()-t0, `seconds to generate (but using pre-computed values for the generating functions).`):
print(` I hope that you enjoyed it as much as I did generating it.  `):
end:



#MamarH(N,K,t): inputs a positive integer N, and a positive integer K and a variable name t
#outputs an article with explicit expressions as rational functions, in t, for the sequences
#of M[m,n](k)=Int(P_k(exp(2*Pi*I*theta)^m*P_k(exp(-2*Pi*I*theta)^n,theta=0..1) , k=0..infinity`):
#for 1<=m<n<=N, and testing for these values a conjecture of Hugh Montgomery made in Oberwolfach
#20 October – 26 October 2013, published on-line, and available from the url
#https://www.mfo.de/document/1343/OWR_2013_51.pdf (see pp. 67-68)
#It also gives the first K terms of the actual sequence,
#and confirms Hugh Montgomery's conjecture for these m and n. Try
#MamarH(4,100,t):
MamarH:=proc(N,K,t)  local t0,P,lu,i,k,m,n,f,gu,mu,co,z:
t0:=time():
co:=0:
print(`Explicit Generating Functions, Asymptotics, and More for the "Mixed" Moments of Rudin-Shapiro polynomials up to`, N):
print(``):
print(`By Shalosh B. Ekhad (ShaloshBEkhad@gmail.com )`):
print(``):
print(`Let `, P[k](z), ` be the Rudin-Shapiro polynomial as described for example, in: https://en.wikipedia.org/wiki/Shapiro_polynomials .` ):
print(`The best way to define them is via the recurrence `):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(`Let m, n and k be a positive integers, and define, with Hugh Montgomery (in Oberwolfach)`):
print(`https://www.mfo.de/document/1343/OWR_2013_51.pdf (see pp. 67-68)`):
print(``):
print(M[m,n](k)=Int(P[k](exp(2*Pi*I*t))^m*P[k](exp(-2*Pi*I*t))^n,t=0..1)):
print(`It turns out (from our approach) that, for every specific positive integers m and n, the generating function of the sequence M[m,n](k)`):
print(``):
print(f[m,n](t)=Sum(M[m,n](k)*t^k,k=0..infinity)):
print(``):
print(`is always a rational function of t.`):
print(`Here we will derive these generating functions, fully rigorously, for all 1<m<n <=`, N ):
print(` and also confirm Hugh Montgomery's conjecture that the growth`):
print(`is less than`, (2^((m+n)/2))^k, `in fact, except for (m,n)=(1,2) it is turns out fo be way smaller (at least for m<n<=`, N):
print( `considered in this article`):
print(``):
print(`For the sake of our beloved OEIS, we also supply the first`, K+1, `terms of the actual sequences`):
print(``):
print(`---------------------------------------------------------------------`):
print(``):

for m from 1 to N do
for n from m+1 to N do
co:=co+1:
gu:=RS(m,n,t,10000):
lu:=taylor(gu,t=0,K+1):
lu:=[seq(coeff(lu,t,i),i=0..K)]:
print(``):
print(`-------------------------------------------------`):
print(``):
print(`Theorem Number`, co  ):
print(``):
lprint(f[m,n](t)=gu):
print(``):

mu:=evalf(AsyRat(gu,t)):

if nops(mu[2])=1 then

 print(`The smallest root of the denominator is`, 1/mu[1] , `and all other roots have strictly larger absolute value, hence`):
 print(`the asymptotic expression is`, mu[2][1]*mu[1]^k ):

  if mu[1]<evalf(2^((m+n)/2)) then
   print(`This is indeed smaller than 2^((m+n)/2)=`, 2^((m+n)/2), `as conjectured by Montgomery `):
  fi:

else
 print(`The smallest real root of the denominator is`, 1/mu[1] , ` but there are other roots with the same absolute value, hence`):
 print(`the sequence is oscillating, but its absolute value is bounded by a constant times`, mu[1]^k):
 print(`The constant in front,for the first 1000 terms, does not exceed`, mu[2][2] ):
  if mu[1]<evalf(2^((m+n)/2)) then
   print(`The growth constant `	, mu[1], `is indeed smaller than 2^((m+n)/2)=`, evalf(2^((m+n)/2)), `as conjectured by Montgomery. `):
  fi:

fi:

  
print(`The first`, K+1, `terms are `):
print(``):
print(lu):


od:
od:

print(``):
print(`---------------------------------------------------------------`):
print(`This concludes this article that took`, time()-t0, `seconds to generate.`):
print(` I hope that you enjoyed it as much as I did generating it.  `):
end:










#BigMates(n,a,b,A,B,z,K) : inputs a positive integer n and symbols a,b,A,B,z,k
#returns the set of monomials in {a,b,A,B,z} that have the same rate of growth as a^n*A^n
#in the scheme that is used to compute its generating function
#is to investigate a possible route to proving Saffari's conjecture. Try:
#BigMates(2,a,b,A,B,z,1000);
BigMates:=proc(n,a,b,A,B,z,K)  local t,lu,i,kv:

lu:=HaroldP(a^n*A^n,a,b,A,B,z,K,t)[2]: 

kv:={}:

for i from 1 to nops(lu) do
 if subs(t=1/2^n,normal((1-2^n*t)*lu[i][2]))<>0 then
   kv:=kv union {lu[i][1]}:
 fi:
od:

kv:

end:




#MamarDseq(N,L,K): inputs a positive integer N <=10, and a positive integer K (take it to be 10000), 
#(it is for internal purposes, for the upper bound for the size of the scheme (see main (human) article)).
#Outputs an article that lists the first L+1 terms (from k=0 to k=L)of
#the even Moments, from 2 to 2N, of the Rudin-Shapiro polynomials. With sufficiently large L,
#this data can be used to find the generating function, if desired, but for large N it may take a long time.
#Try:
#MamarDseq(4,100,1000):
MamarDseq:=proc(N,L,K)  local t0,P,lu,i,A,k,m,f,z,t:
t0:=time():
print(`The First,`,L+1 , `terms of the Sequences for the even moments of the Rudin-Shapiro polynomials from the Second to the`, 2*N, `-th moment. `):
print(``):
print(`By Shalosh B. Ekhad (ShaloshBEkhad@gmail.com )`):
print(``):
print(` Let `, P[k](z), ` be the Rudin-Shapiro polynomial as described for example, in:`):
print(``):
print(` https://en.wikipedia.org/wiki/Shapiro_polynomials .` ):
print(``):
print(`The best way to define them is via the recurrence `):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(``):
print(`and the initial condition`):
print(` P[0](z)=1 `):
print(``):

print(`Let m and k be a positive integers, and define the even moment`):
print(``):
print(A[k,m]=Int(abs(P[k](exp(2*Pi*I*t)))^(2*m),t=0..1)):
print(`It was first proved by Christoph Doche and Laurent Habsieger: `):
print(` J. of Fourier Analysis and Applications  Vol. 10 Issue 5, 497-505, (2004)`):
print(``):
print(`available from:  http://web.science.mq.edu.au/~doche/049.pdf `):
print(``):
print(`that for each specific positive integer m, the sequence A[k,m] satisfies a linear recurrence equation with`):
print(`constant coefficients, or equivalently, that the generating function`):
print(``):
print(f[m](t)=Sum(A[k,m]*t^k,k=0..infinity)):
print(``):
print(`is always a rational function of t. This fact seems also to follow from our approach.`):
print(``):
print(`Since the original Pari-program referred to by these`):
print(`authors is no longer available, and neither is the output, and since these polynomials, and associated sequences are so important`):
print(`we decided to recompute the first`, L+1, `terms of these sequences for even moments from the second through the`, 2*N, `-th . `):
print(`We will also test Saffari's conjecture that `):
print( A[k,m], ` is asymptotic to`,  (2^m/(m+1))*(2^m)^k, `as k goes to infinity. `):
print(``):
print(`---------------------------------------------------------------------`):
print(``):

for m from 1 to N do
print(`Let `):
print(``):
print(A[k,m]=Int(abs(P[k](exp(2*Pi*I*t)))^(2*m),t=0..1)):
print(``):
lu:=SeqRS(m,m,0,0,K,L):

print(`The first`, L+1, `terms are `):
print(``):
lprint(lu):
print(``):
print(`The last 10 terms of this sequence, divided by (2^m)^k*2^m/(m+1), i.e.`, (2^m)^k*2^m/(m+1) , `are `):
print(``):
print(evalf([seq(lu[i+1]/2^(m*(i+1))*(m+1),i=nops(lu)-11..nops(lu)-1)])):
od:

print(``):
print(`---------------------------------------------------------------`):
print(`This concludes this article that took`, time()-t0, `seconds to generate.`):
print(` I hope that you enjoyed it as much as I did generating it.  `):

end:




#MamarHseq(N,K,L): inputs a positive integer N, and a positive integer K (take it to be 10000), 
#and a positive integer L
#outputs an article with the first L+1 term of the sequence
#M[m,n](k):=Int(P_k(exp(2*Pi*I*theta)^m*P_k(exp(-2*Pi*I*theta)^n,theta=0..1) , k=0..infinity`):
#for 1<=m<n<=N, and testing that the absolute  values satisfy the upper bound  conjectured by Hugh Montgomery made in Oberwolfach
#20 October – 26 October 2013, published on-line, and available from the url
#https://www.mfo.de/document/1343/OWR_2013_51.pdf (see pp. 67-68). Try:
#MamarHseq(4,10000,40);
MamarHseq:=proc(N,K,L)  local t0,P,i,k,m,n,f,gu,co,z,lu1,t:
t0:=time():
co:=0:
print(`The first`, L+1, `terms for the "Mixed" Moments of Rudin-Shapiro polynomials up to`, N):
print(``):
print(`By Shalosh B. Ekhad (ShaloshBEkhad@gmail.com )`):
print(``):
print(`Let `, P[k](z), ` be the Rudin-Shapiro polynomial as described for example, in: https://en.wikipedia.org/wiki/Shapiro_polynomials .` ):
print(`The best way to define them is via the recurrence `):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(``):
print(`and the initial condition`):
print(` P[0](z)=1 `):
print(``):

print(`Let m, n and k be a positive integers, and define, with Hugh Montgomery (in Oberwolfach)`):
print(`https://www.mfo.de/document/1343/OWR_2013_51.pdf (see pp. 67-68)`):
print(``):
print(M[m,n](k)=Int(P[k](exp(2*Pi*I*t))^m*P[k](exp(-2*Pi*I*t))^n,t=0..1)):
print(`It turns out (from our approach) that, for every specific positive integers m and n, the generating function of the sequence M[m,n](k)`):
print(``):
print(f[m,n](t)=Sum(M[m,n](k)*t^k,k=0..infinity)):
print(``):
print(`is always a rational function of t.`):
print(`Equivalently, that the sequence satisfies a  linear recurrence equation with constant coefficients `):
print(``):
print(`Here we state the first`, L+1, `terms of these sequences for 1<=m<n <=`, N ):
print(` We also confirm Hugh Montgomery's conjecture that the growth`):
print(`is less than`, (2^((m+n)/2))^k, `in fact, except for (m,n)=(1,2) it is turns out fo be way smaller (at least for m<n<=`, N):
print( `considered in this article`):
print(``):
print(``):
print(`---------------------------------------------------------------------`):
print(``):

for m from 1 to N do
for n from m+1 to N do
co:=co+1:
gu:=SeqRS(m,n,0,0,K,L):
print(`Fact Number`, co  ):
print(``):
lprint(`The first`, L+1, `terms of the sequence`, M[m,n](k), `are ` ):
print(``):
lprint(gu):
print(``):
if nops(gu)>20 then
lu1:=evalf([seq(  (abs(gu[i+1]))^(1/i),  i=nops(gu)-21..nops(gu)-1)]):


print(`The minimum and maximum of the k-th root of the absolute values for the last 20 terms are`):
print(min(op(lu1)),max(op(lu1))):
print(``):
print(`compare this with Hugh Montgomery's conecjtured upper bound of 2^((m+n)/2), i.e.`, evalf(2^((m+n)/2))):
print(``):
fi:

od:
od:

print(``):
print(`---------------------------------------------------------------`):
print(`This concludes this article that took`, time()-t0, `seconds to generate.`):
print(` I hope that you enjoyed it as much as I did generating it.  `):
end:


#Medio(n,a,b,A,B,z,K): inputs a positive integer n, symbols a,b,A,B,z, and a large positive integer K
#Let the monomials a^i*A^i*b^(n-i)*B^(n-i), i<=n/2 be called important.
#It outputs the set of unimportant monomials that have important childen. Try:
#Medio(3,a,b,A,B,z,1000);
Medio:=proc(n,a,b,A,B,z,K) local gu,gadol,M,mu,Gadol,Katan,i,yeladim,y,Beno:
gu:=HaroldSM(a^n*A^n,a,b,A,B,z,K):
M:=gu[2]:
mu:=gu[3]:

gadol:={seq((a*A)^(n-i)*(b*B)^i,i=0..trunc(n/2))}:

Gadol:={}:
Katan:={}:

for i from 1 to nops(mu) do
 if member(mu[i],gadol) then
  Gadol:=Gadol union {i}:
 else
  Katan:=Katan union {i}:
 fi:
od:


Beno:={}:

for i from 1 to nops(M) do

 if member(i,Katan) then
  yeladim:={seq(y[2], y in M[i])}:

 if yeladim intersect Gadol<>{} then
   Beno:=Beno union {mu[i]}:
   
 fi:
 fi:
od:
  
Beno:

end:




#BigMatrSlow(n): The n+1 by n+1 transition matrix for the important monomials (a*A)^(n-i)*(b*B)^i, 0<=i<=n. Try:
#BigMatrSlow(4);
BigMatrSlow:=proc(n) local i,L,L1,gu,a,b,A,B,z,i1:

L:=[]:

for i from 0 to n do
 gu:=expand(subs({a=a+z*b,b=a-z*b, A=A+B/z, B=A-B/z},(a*A)^(n-i)*(b*B)^i));
 L1:=[seq(coeff(coeff(coeff(coeff(coeff(gu,z,0),a,n-i1),A,n-i1),B,i1),b,i1),i1=0..n)]:
 L:=[op(L),L1]:
od:

L:

end:

#BigGF(n,t): The asymptotic generating functions for the monomials a^(n-i)*A^(n-i)*b^i*B^i for i from 0 to n. Under the "little guys do not count"
#assumptiom. Try:
#BigGF(4,t);
BigGF:=proc(n,t) local x,i,j,L,var,eq:
option remember:
L:=BigMatr(n):

var:={seq(x[i],i=0..n)}:
eq:={seq(x[i]-1-t*add(L[i+1][j+1]*x[j],j=0..n),i=0..n)}:
var:=solve(eq,var):

factor(subs(var,[seq(x[i],i=0..n)])):
end:

#BigAsySlow(n): Inputs a positive integer n, and outputs the list of constants C, such that
#CT_z (P[k](z)*P[k](1/z))^(n-i)*(P[k](-z)*P[k](-1/z))^i) is asymptotic to C[i+1]*(2^m)^k
#as k goes to infinity for i from 0 to n, under the "little guys don't count" assumption
#Try:
#BigAsySlow(4);
BigAsySlow:=proc(n) local gu,t,i:

gu:=BigGF(n,t):

[seq(subs(t=1/2^n,normal((1-2^n*t)*gu[i])),i=1..nops(gu))]:
end:


#BigAsy(n): A fast version of BigAsy(n), just using the transtion matrix, BigMatr(n),  and
#the fact that the word (a*A)^n is asymptotic to (2^n)^k* 2^n/(n+1)
#Try:
#BigAsy(4);
BigAsy:=proc(n) local x,i,j,L,eq,var:

L:=BigMatr(n):
var:={seq(x[i],i=0..n)}:
eq:={seq(x[i]*2^n-add(L[i+1][j+1]*x[j],j=0..n),i=0..n),x[0]=2^n/(n+1)}:
var:=solve(eq,var):

subs(var,[seq(x[i],i=0..n)]):

end:




#MamarDseqG(N,L,K): inputs a positive integer N <=10, and a positive integer K (take it to be 10000), 
#(it is for internal purposes, for the upper bound for the size of the scheme (see main (human) article)).
#Outputs an article that lists the first L+1 terms (from k=0 to k=L)of
#the sequences CT_z (P[k](z)*P[k](1/z))^(n-i)*(P[k](-z)*P[k](-1/z))^i) is asymptotic to C[i+1]*(2^m)^k
#n from 0 to N and i from 0 to n, where P[k](z) is the k-th Rudin-Shapiro polynomial. With sufficiently large L,
#this data can be used to find the generating function, if desired, but for large N it may take a long time.
#Try:
#MamarDseqG(4,100,1000):
MamarDseqG:=proc(N,L,K)  local t0,P,lu,i,A,k,m,f,z,t,i1:
t0:=time():
print(``):
print(`The First,`,L+1 , `terms of the Sequences for`):
print(Int((P[k](exp(2*Pi*I*t))*P[k](exp(-2*Pi*I*t)))^(m-i)*(P[k](-exp(2*Pi*I*t))*P[k](-exp(-2*Pi*I*t)))^i,t=0..1)):
print(`where P[k](z) are the Rudin-Shapiro polynomials. For 0m<=`, N, `and i<=m/2 `):
print(``):
print(`By Shalosh B. Ekhad (ShaloshBEkhad@gmail.com )`):
print(``):
print(` Let `, P[k](z), ` be the Rudin-Shapiro polynomial as described for example, in:`):
print(``):
print(` https://en.wikipedia.org/wiki/Shapiro_polynomials .` ):
print(``):
print(`The best way to define them is via the recurrence `):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(``):
print(`and the initial condition`):
print(` P[0](z)=1 `):
print(``):

print(`Let k be a positive integer, and let 0<=i<=m, be positive integers  define the squence, in k,`):
print(``):
print(A[m,i](k)=Int((P[k](exp(2*Pi*I*t))*P[k](exp(-2*Pi*I*t)))^(m-i)*(P[k](-exp(2*Pi*I*t))*P[k](-exp(-2*Pi*I*t)))^i,t=0..1)):
print(``):
print(`For every specific pos. integer m, and a pos. integer i between 0 and m, the sequence A[m,i](k) satisfies a linear recurrence equation with`):
print(`constant coefficients, or equivalently,  the generating function`):
print(``):
print(f[m,i](t)=Sum(A[m,i](k)*t^k,k=0..infinity)):
print(``):
print(`is always a rational function of t.`):
print(``):
print(`In this article we will give the first`, L+1, `terms of these sequences for m<=`, N, `and i<=m/2 `):
print(`Note that by symmetry, A[m,i]=A[m,m-i], so it suffices to consider i<=m/2 `):
print(``):
print(`We will also test the generalization of Saffari's conjecture that `):
print( A[m,i](k), ` is asymptotic to`,  2^m/((m+1)*binomial(m,i))*(2^m)^k, `as k goes to infinity. `):
print(``):
print(`---------------------------------------------------------------------`):
print(``):

for m from 1 to N do
 for i from 0 to trunc(m/2) do
print(`Let `):
print(``):
print(A[m,i](k)=Int((P[k](exp(2*Pi*I*t))*P[k](exp(-2*Pi*I*t)))^(m-i)*(P[k](-exp(2*Pi*I*t))*P[k](-exp(-2*Pi*I*t)))^i,t=0..1)):
print(``):
lu:=SeqRS(m-i,m-i,i,i,K,L):

print(`The first`, L+1, `terms are `):
print(``):
lprint(lu):
print(``):
print(`The last 10 terms of this sequence, divided by (2^m)^k*2^m/((m+1)*binomial(m,i)), i.e.`, (2^m)^k*2^m/(m+1).binomial(m,i) , `are `):
print(``):
print(evalf([seq(lu[i1+1]/2^(m*(i1+1))*(m+1)*binomial(m,i),i1=nops(lu)-11..nops(lu)-1)])):
od:
od:

print(``):
print(`---------------------------------------------------------------`):
print(`This concludes this article that took`, time()-t0, `seconds to generate.`):
print(` I hope that you enjoyed it as much as I did generating it.  `):

end:


#K(n,k,r) Inputs positive integers 0<=k<=n, and 0<=r<=n: and outputs (add((-1)^a*binomial(n-k,r-a)*binomial(k,a),a=0..r))^2:
#Try:
#K(6,3,2);
K:=proc(n,k,r) local a:
(add((-1)^a*binomial(n-k,r-a)*binomial(k,a),a=0..r))^2:
end:


#BigMatr(n): Fast version, ising the explicit formula for the entriesof BigMatr(n) (q.v.). Try:
#BigMatr(4);
BigMatr:=proc(n) local k,r:
[seq([seq(K(n,k,r),r=0..n)],k=0..n)]:
end:

#Just for checking
Ak:=proc(k) option remember: if k=0 then -1: elif k=1 then  0: else Ak(k-2)+k-1: fi: end:  Bk:=proc(k): k^2/4-5/8-3*(-1)^k/8:end:
#End Just for checking

#CheckNuMedio(N,K): checks that the number of mediocre terms for n equals the sequence Ak(n)
#  mentioned in the website http://mrob.com/pub/math/MCS-10.html, MCS42828,
#defined by Ak(0)=-1, Ak(1)=0, Ak(k)=Ak(k-2)+k-1, that in fact equals
#n^2/4-5/8-3*(-1)^n/8
CheckNuMedio:=proc(N,K) local n,a,b,A,B,z:
evalb({seq(nops(Medio(n,a,b,A,B,z,K))-(n^2/4-5/8-3*(-1)^n/8),n=1..N)}={0});
end:


#Nake1(f,t): inputs a reciprocal of an affine linear function and outputs the pair [C,alpha] such that it equals
#C/(1-alpha*t). Try:
#Nake1(4/(11+6*t),t);
Nake1:=proc(f,t) local gu,a,b:
gu:=normal(1/f):
if degree(gu,t)<>1 then
  RETURN(FAIL):
fi:

a:=coeff(gu,t,0):
b:=coeff(gu,t,1):

if a=0 then
  RETURN(FAIL):
fi:

if normal((1/a)/(1+(b/a)*t)-f)<>0 then
  RETURN(FAIL):
fi:

[1/a,-b/a]:
end:


#CPF(R,t): inputs a rational function whose denominator factorizes into a product of linear factors
#outputs its partial fraction decomosition in as a list [[C1,a1],[C2,a2], ..., [Ck,ak]] such that
#a1>a2>..>ak and R=C1/(1-a1*t)+ ...+ Ck/(1-ak*t). Try:
#CPF((t^2+1)/((1-2*t)*(1-3*t)*(1+2*t)),t);
CPF:=proc(R,t) local mu,mu1,gu,lu,i,fu,T1:

gu:=[]:

mu:=convert(R,parfrac):

for i from 1 to nops(mu) do
mu1:=op(i,mu):

 lu:=Nake1(mu1,t):

 if lu=FAIL then

  RETURN(FAIL):

  else
  gu:=[op(gu),lu]:
 
 fi:
od:

fu:={seq(gu[i][2],i=1..nops(gu))}:

if nops(fu)<nops(gu) then
 RETURN(FAIL):
fi:

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

fu:=sort([seq(gu[i][2],i=1..nops(gu))]):

gu:=[seq([T1[fu[i]],fu[i]],i=1..nops(gu))]:
gu:=[seq(gu[nops(gu)-i+1],i=1..nops(gu))]:

if normal(R-add(gu[i][1]/(1-t*gu[i][2]),i=1..nops(gu)))<>0 then
 FAIL:
else
 gu:
fi:

end:

#BigGFpf(n): inputs a positive integer n, and outputs BigGF(n,t) in compact, partial-fraction format given by CPF(R,t) (q.v).
#Try:
#BigGFpf(6);
BigGFpf:=proc(n) local t,gu,i:
option remember:
gu:=BigGF(n,t):

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

end:

#Shorashim(n): all the roots that show up in BigGFpf(n)
Shorashim:=proc(n) local lu,i:
lu:=BigGFpf(n)[1]:
[seq(lu[i][2],i=1..nops(lu))]:
end:

#ShorashimH(n):  the conjectured roots that show up in BigGFpf(n)
ShorashimH:=proc(n) local gu,j,i:

gu:=[seq(2^(n-4*j)*binomial(4*j,2*j),j=0..trunc(n/4)),seq(-2^(n-4*j-2)*binomial(4*j+2,2*j+1),j=0..trunc((n-2)/4))]:
gu:=sort(gu):

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

end:








#Anr(n,r,t,J1): the conjectured expression for R_n^{(r)}(t) up to the first J1 largest postive roots and the J1 smallest negative roots
Anr:=proc(n,r,t,J1) local j1:


add( GuessA(n,r,j1) /(1-2^(n-4*j1)*binomial(4*j1,2*j1)*t),j1=0..J1)+
add(GuessB(n,r,j1) /(1+2^(n-4*j1-2)*binomial(4*j1+2,2*j1+1)*t),j1=0..J1):

end:

#AnrE(n,r,k,J1): Inputs SYMBOLS n,r, k, and an integer J1
#outputs the  expression for the beginning and end for the exact formula for
#the "Big Guys" approximation to 
#Int((P[k](exp(2*Pi*I*t))*P[k](exp(-2*Pi*I*t)))^(n-r)*(P[k](-exp(2*Pi*I*t))*P[k](-exp(-2*Pi*I*t)))^r,t=0..1)):
#Try:
#AnrE(n,r,k,1);
AnrE:=proc(n,r,k,J1) local j1:
add( GuessA(n,r,j1) *(2^(n-4*j1)*binomial(4*j1,2*j1))^k,j1=0..J1)+
add(GuessB(n,r,j1) *(-2^(n-4*j1-2)*binomial(4*j1+2,2*j1+1))^k,j1=0..J1):
end:

#AnrE1(n,r,k,J1): Inputs SYMBOLS n  k, and a positive integer r1, and a positive  integer J1
#the conjectured expression for the beginning and end for the exact formula for
#the "Big Guys" approximation to 
#Int((P[k](exp(2*Pi*I*t))*P[k](exp(-2*Pi*I*t)))^(n-r)*(P[k](-exp(2*Pi*I*t))*P[k](-exp(-2*Pi*I*t)))^r,t=0..1)):
#Try:
#AnrE(n,r,k,1);
AnrE1:=proc(n,r1,k,J1) local j1:
add( simplify(eval(subs(r=r1,GuessA(n,r,j1)))) *(2^(n-4*j1)*binomial(4*j1,2*j1))^k,j1=0..J1)+
add(simplify(eval(subs(r=r1,GuessB(n,r,j1)))) *(-2^(n-4*j1-2)*binomial(4*j1+2,2*j1+1))^k,j1=0..J1):
end:


#CollectA(n1,j1,k): inputs   numeric n1 and j1 such that j1<=trunc(n1/4) and a variable k, outputs
#the polynomial in k for the coefficients of 1/(1-2^(n-4*j)*binomial(4*j,2*j)*t) in R_{n1,k} times
#binomial(n1,k)/(2^n1)*(n1+4*j1+1)!/n1!)
#Try:
#CollectA(24,1,k):
CollectA:=proc(n1,j1,k) local gu,k1,lu,efes,a1:

if not (j1>=0 and j1<=trunc(n1/4)) then
 print(`j1 must be between 0 and trunc(n1/4)`):
 RETURN(FAIL):
fi:

efes:=2^(n1-4*j1)*binomial(4*j1,2*j1):

lu:=BigGFpf(n1):

gu:=[]:

for k1 from 0 to n1 do
 if not member(efes,{seq(lu[k1+1][a1][2],a1=1..nops(lu[k1+1]))}) then
   gu:=[op(gu),0]:
 else

  for a1 from 1 to nops(lu[k1+1]) while lu[k1+1][a1][2]<>efes do od:
   gu:=[op(gu),lu[k1+1][a1][1]]:
 fi:
od:

gu:=[seq(gu[a1+1]*binomial(n1,a1)*(n1+4*j1+1)!/n1!/(2^n1),a1=0..n1)]:

expand(subs(k=k+1,Gr(gu,k))):

end:




#CheckEV(N): checks empirically that 2^n is indeed an eigenvalue of the matrix M^(n) and
#[seq(1/(binomial(n,r),r=0..n)] is an eigenvector for all n<=N. Try
#CheckEV(20);
CheckEV:=proc(N) local r,a,n,k:
evalb({
seq(seq(2^(n)/binomial(n,k)-add( add((-1)^a*binomial(n-k,r-a)*binomial(k,a),a=0..r)^2/binomial(n,r),r=0..n),
k=0..n),n=1..N)}={0}):
end:




#GuessA(n,k,j1):inputs  symbols n and k and a non-negative integer j, and a positive integer N
#outputs a polynomial in k and n
#of degree 4*j in k let's call it
# A_j(n,k), such that the numerator of 1/(1-binomial(4*j,2*j)*2^(n-4*j)*t) is
#A_j(n,k)
#Try: 
#GuessA(n,k,1);
GuessA:=proc(n,k,j1) local lu,n1,mu,lu1,r,a1,N:
N:=6+8*j1:
lu:=[seq(CollectA(n1,j1,k),n1=4*j1+3..N)]:

mu:=0:

for r from 0 to 4*j1 do
 lu1:=[seq(coeff(lu[a1],k,r),a1=1..nops(lu))]:
 lu1:=Gr(lu1,n):
  if lu1=FAIL then
     RETURN(FAIL):
  else
    lu1:=factor(subs(n=n-4*j1-2,lu1)):
    mu:=mu+lu1*k^r:
  fi:
od:
2^n*k!*(n-k)!/(n+4*j1+1)!* add(factor(coeff(mu,k,a1))*k^a1,a1=0..degree(mu,k)):
end:




#CollectB(n1,j1,k): inputs   numeric n1 and j1 such that j1<=trunc((n1-2)/4) and a variable k, outputs
#the polynomial in k for the coefficients of 1/(1+2^(n-4*j-2)*binomial(4*j+2,2*j+1)*t) in R_{n1,k} times
#binomial(n1,k)/(2^n1)*(n1+4*j1+1)!/n1!)
#Try:
#CollectB(24,1,k):
CollectB:=proc(n1,j1,k) local gu,k1,lu,efes,a1:

if not (j1>=0 and j1<=trunc((n1-2)/4)) then
 print(`j1 must be between 0 and trunc(n1/4)`):
 RETURN(FAIL):
fi:

efes:=-2^(n1-4*j1-2)*binomial(4*j1+2,2*j1+1):

lu:=BigGFpf(n1):

gu:=[]:

for k1 from 0 to n1 do
 if not member(efes,{seq(lu[k1+1][a1][2],a1=1..nops(lu[k1+1]))}) then
   gu:=[op(gu),0]:
 else

  for a1 from 1 to nops(lu[k1+1]) while lu[k1+1][a1][2]<>efes do od:
   gu:=[op(gu),lu[k1+1][a1][1]]:
 fi:
od:

gu:=[seq(gu[a1+1]*binomial(n1,a1)*(n1+4*j1+1)!/n1!/(2^n1),a1=0..n1)]:

expand(subs(k=k+1,Gr(gu,k))):

end:

#GuessB(n,k,j1):inputs  symbols n and k and a non-negative integer j, and a positive integer N
#outputs a polynomial in k and n
#of degree 4*j in k let's call it
# B_j(n,k), such that the numerator of 1/(1+binomial(4*j+2,2*j+1)*2^(n-4*j-2)*t) is
#A_j(n,k)
#Try: 
#GuessB(n,k,1);
GuessB:=proc(n,k,j1) local lu,n1,mu,lu1,r,a1,N:
N:=6+8*j1+7:
lu:=[seq(CollectB(n1,j1,k),n1=4*j1+5..N)]:


mu:=0:

for r from 0 to 4*j1+2 do
 lu1:=[seq(coeff(lu[a1],k,r),a1=1..nops(lu))]:
 lu1:=Gr(lu1,n):
  if lu1=FAIL then
     RETURN(FAIL):
  else
    lu1:=factor(subs(n=n-4*j1-4,lu1)):
    mu:=mu+lu1*k^r:
  fi:
od:
2^n*k!*(n-k)!/(n+4*j1+1)!* add(factor(coeff(mu,k,a1))*k^a1,a1=0..degree(mu,k)):
end:



#MamarD(N,K,t,EE): inputs a positive integer N <=10, and a positive integer K and a variable name t
#and EE that is either 0 or 1 (if EE=0 there is no error estimate. if EE=1 then there is):
#outputs an article with explicit expressions as rational functions for the sequences
#of the first N even moments. It also gives the first K terms of the actual sequence,
#and asymptotics. Try:
#MamarD(4,100,t):
MamarD:=proc(N,K,t,EE)  local t0,P,gu,lu,i,A,k,m,f,gu1,c,mu,z:
t0:=time():
if N>10 then
 print(`The first argument must be <=10`):
 RETURN(FAIL):
fi:
print(`Explicit Generating Functions, Asymptotics, and More for the first`, N, `even moments of the Rudin-Shapiro polynomials`):
print(``):
print(`By Shalosh B. Ekhad (ShaloshBEkhad@gmail.com )`):
print(``):
print(` Let `, P[k](z), ` be the Rudin-Shapiro polynomial as described for example, in:`):
print(``):
print(` https://en.wikipedia.org/wiki/Shapiro_polynomials .` ):
print(``):
print(`The best way to define them is via the recurrence `):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(``):
print(`and the initial condition`):
print(` P[0](z)=1 `):
print(``):

print(`Let m and k be a positive integers, and define the even moment`):
print(``):
print(A[k,m]=Int(abs(P[k](exp(2*Pi*I*t)))^(2*m),t=0..1)):
print(`It was first proved by Christoph Doche and Laurent Habsieger: `):
print(` J. of Fourier Analysis and Applications  Vol. 10 Issue 5, 497-505, (2004)`):
print(``):
print(`available from:  http://web.science.mq.edu.au/~doche/049.pdf `):
print(``):
print(`that for each specific positive integer m, the generating function`):
print(``):
print(f[m](t)=Sum(A[k,m]*t^k,k=0..infinity)):
print(``):
print(`is always a rational function of t. This fact seems also to follow from our approach.`):
print(``):
print(`Since the original Pari-program referred to by these`):
print(`authors is no longer available, and neither is the output, and since these polynomials, and associated sequences are so important`):
print(`we decided to recompute them, and also confirm, for m from 1 to`, N, `the famous conjecture  of Saffari that  for any specific m`):
print( A[k,m], ` is asymptotic to`,  (2^m/(m+1))*(2^m)^k, `as k goes to infinity. `):
print(``):
print(`For the sake of our beloved OEIS, we also supply the first`, K+1, `terms of the actual sequences`):
print(``):
print(`---------------------------------------------------------------------`):
print(``):

for m from 1 to N do
gu:=Spc(m,t):
lu:=taylor(gu,t=0,K+1):
lu:=[seq(coeff(lu,t,i),i=0..K)]:
print(``):
print(`-------------------------------------------------`):
print(``):
print(`Theorem Number`, m  ):
print(``):
lprint(f[m](t)=gu):
print(``):

gu1:=denom(gu):
gu1:=normal(gu1/(1-2^m*t)):
if degree(denom(gu1),t)=0 then
print(`Note that the smallest root of the denominator is`, 1/2^m):
c:=subs(t=1/2^m,normal((1-2^m*t)*gu)):
print(`and the residue there is`, c, `that implies that the asymptotics is indeed`, c*(2^m)^k,`as conjectured by Saffari.`):
fi:

if EE=1 then
gu:=normal((1-2^m*t)*gu):
mu:= [solve(denom(gu),t)]:
mu:=min(seq(evalf(abs(mu[i])),i=1..nops(mu))):
print(`The largest modulo of the second largest root is`, mu):
 if mu>1/2^m then
  print(`that implies that the error term is exponentially small, of the order of`, (1/2^m/mu)^k ):
 fi:
fi:

print(`The first`, K+1, `terms are `):
print(``):
print(lu):


od:

print(``):
print(`---------------------------------------------------------------`):
print(`This concludes this article that took`, time()-t0, `seconds to generate (but using pre-computed values for the generating functions).`):
print(` I hope that you enjoyed it as much as I did generating it.  `):
end:





#CheckCP1(n): checks empirically that the characteristic polynomial of BigMatr(n) is indeed
#mul((x-2^(n-4*j)*binomial(4*j,2*j),j=0..trunc(n/4)))*mul(x+2^(n-4*j-2)*binomial(4*j+2,2*j+1),j=0..trunc((n-2)/4))*x^(trunc(n+3)/2):
#Try:
#CheckCP1(5);
CheckCP1:=proc(n) local gu1,gu2,j:

gu1:=factor(linalg[charpoly](BigMatr(n),x)):

gu2:=mul((x-2^(n-4*j)*binomial(4*j,2*j),j=0..trunc(n/4)))*mul(x+2^(n-4*j-2)*binomial(4*j+2,2*j+1),j=0..trunc((n-2)/4))*x^(trunc((n+1)/2)):

evalb(normal(gu1/gu2)=1):

end:


#CheckCP(N): checks empirically that the characteristic polynomial of BigMatr(n) is indeed true from n=2 to n=N
#mul((x-2^(n-4*j)*binomial(4*j,2*j),j=0..trunc(n/4)))*mul(x+2^(n-4*j-2)*binomial(4*j+2,2*j+1),j=0..trunc((n-2)/4))*x^(trunc(n+3)/2):
#Try:
#CheckCP(5);
CheckCP:=proc(N) local n: evalb({seq(CheckCP1(n),n=2..N)}={true}):end:



#MedioP(n,a,b,A,B,z,K): inputs a positive integer n, symbols a,b,A,B,z, and a large positive integer K
#Let the monomials a^i*A^i*b^(n-i)*B^(n-i), i<=n/2 be called important.
#It outputs the set of unimportant monomials that have important childen. Followed bu their expression in terms of important monomials
#Try:
#MedioP(3,a,b,A,B,z,1000);
MedioP:=proc(n,a,b,A,B,z,K) local gu,gadol,M,mu,Gadol,Katan,i,yeladim,y,Beno, Geno,lu,lu1, w,co:
gu:=HaroldSM(a^n*A^n,a,b,A,B,z,K):
M:=gu[2]:
mu:=gu[3]:

gadol:={seq((a*A)^(n-i)*(b*B)^i,i=0..trunc(n/2))}:

Gadol:={}:
Katan:={}:

for i from 1 to nops(mu) do
 if member(mu[i],gadol) then
  Gadol:=Gadol union {i}:
 else
  Katan:=Katan union {i}:
 fi:
od:


Beno:={}:

for i from 1 to nops(M) do

 if member(i,Katan) then
  yeladim:={seq(y[2], y in M[i])}:

 if yeladim intersect Gadol<>{} then
   Beno:=Beno union {mu[i]}:
   
 fi:
 fi:
od:
  
Beno:

Geno:={}:

for w in Beno do

 lu:=Monos(GD(w,a,b,A,B,z),[a,b,A,B,z]):
  co:=0:
  mu:=0:

  for lu1 in lu do
    if member(lu1[2],gadol) then
     mu:=mu+lu1[1]*lu1[2]:
    co:=co+lu1[1]/binomial(n,degree(lu1[2],a)):
    fi:
  od:
Geno:=Geno union {[w,mu,co]}:
      
od:

Beno, Geno:

end:



#MedioPv1(n,K):  Inputs a positive integer n, and a large positive integer K (for the upper bound for the scheme)
# tests the Important Monomial Hypothesis for the
#Saffari conjecture for Int(|P_n(exp(2*Pi*t))|^(2*n),t=0..1)
# that all the Non-Important Children of important monomials
#Try:
#MedioPv1(3,1000);
MedioPv1:=proc(n,K) local gu,a,b,A,B,z,r,k,gu1,t:

gu:=MedioP(n,a,b,A,B,z,K):

print(`We will assume the hypothesis that the monomial `, (a*A)^(n-r)*(b*B)^r, `is proportional to`, 1/binomial(n,r) ):
print(``):
print(`The scheme that yields the generating functions for `):
print( Int( (P[k](exp(2*Pi*t))*P[k](exp(-2*Pi*t)) )^(n-r)* (P[k](-exp(2*Pi*t))*P[k](-exp(-2*Pi*t)) )^(r) ,t=0..1)):
print(`has the following mediocre monomials, who have at least one important monomial as a child`):
print(``):
print(gu[1]):
print(``):
gu:=gu[2]:

for gu1 in gu do
 print(`The important part in the children of`, gu1[1], `are `):
 print(``):
 print(gu1[2]):
 print(``):
 print(`and under the above hypothesis the sum total is`, gu1[3]):
od:

print(``):

end:


#MedioPv(N,K):  Inputs a positive integer N, and a large positive integer K (for the upper bound for the scheme)
# tests the Important Monomial Hypothesis for the
#Saffari conjecture for Int(|P_n(exp(2*Pi*t))|^(2*n),t=0..1)
# that all the Non-Important Children of important monomials
#for all n between 1 and N
#Try:
#MedioPv(10,1000);
MedioPv:=proc(N,K) local n:

for n from 1 to N do
 MedioPv1(n,K):
od:
end:



#Begin Doche-Habsieger approach


DHSn:=proc(n,z,a,a1,b,b1,q) local y,lu:option remember:
if n=0 then RETURN(expand(((a+b)*(a1+b1))^q)): fi:
lu:=DHSn(n-1,z,a,a1,b,b1,q): lu:=subs({z=y,a=a+b, a1=a1+b1,b=(a-b)*y,b1=(a1-b1)/y},lu):
lu:=expand((lu+subs(y=-y,lu))/2):  simplify(subs(y=sqrt(z),lu),symbolic):end:

#DHMn(n,q): the 2q-th moment of the Rudin-Shapiro polynomial P_n(z), using Christophe Doche and Laurent Habsieger's approach
DHMn:=proc(n,q) local lu,z,a,a1,b,b1: lu:=coeff(DHSn(n,z,a,a1,b,b1,q),z,0):subs({a=1,b=0,a1=1,b1=0},lu):end:

#DHMnSeq(N,q): the sequence 2q-th moment of the Rudin-Shapiro polynomial P_n(z) from n=0 to n=N
#Should give the same output as  as SeqRS(q,q,0,0,10000,N);
DHMnSeq:=proc(N,q) local n:[seq(DHMn(n,q),n=0..N)]: end:

#End Doche-Habsieger approach




#HaroldSMv(w,a,b,A,B,z,K): Verbose version of HaroldSM(w,a,b,A,B,z,K) (q.v.)
#Try
#HaroldSMv(a^2*A^2,a,b,A,B,z,10);
HaroldSMv:=proc(w,a,b,A,B,z,K,t)	local S,StillToDo,adam,yeled,i,AlreadyDone,T1,V1,C1,C1i,co,i1,j1,T2,ka,E,k,c1,c2,c3,c4,c5,kha,F,eq,var,t0,i2:
				 
t0:=time():
print(`A Linear Recurrence Scheme for the Constant Term of `,subs({a=P[k](z),A=P[k](1/z),b=P[k](-z), B=P[k](-1/z)},w)):
print(`Where P[k](z) is the Rudin-Shapiro polynomial `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(` Let `, P[k](z), `be the Rudin-Shapiro polynomial, that may be defined by the recurrence`):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(``):
print(`and the initial condition`):
print(` P[0](z)=1 `):
print(``):
print(`Definition: For any monomial w, in a,b,A,B,z,`,  a^c1*A^c2*b^c3*B^c4*z^c5, `let `):
print(E[a^c1*A^c2*b^c3*B^c4*z^c5](k), `be the coefficient of z^0  in the polynomial `):
print(``):
print(P[k](z)^c1*P[k](1/z)^c2*P[k](-z)^c3*P[k](-1/z)^c4*z^c5):
print(``):
print(``):
print(`We are interested in getting a scheme for computing the sequence`,  E[w](k)):

S:={w}:

StillToDo:={w}:

co:=0:
AlreadyDone:={}:

while  nops(S)<=K and StillToDo<>{} do
 adam:=StillToDo[1]:
 print(`Let's express`, E[adam](k), `in terms of the values of other, related, sequences at` , k-1):	
 co:=co+1:
 C1[co]:=adam:
 C1i[adam]:=co:
  
 AlreadyDone:=AlreadyDone union {adam}:
 StillToDo:= StillToDo minus {adam}:

 yeled:=Monos(GD(adam,a,b,A,B,z),[a,b,A,B,z]):
 print(`Applying the defining recurrence, and using trivial equivalences`):
 print(E[adam](k)=add(yeled[i1][1]*E[yeled[i1][2]](k-1),i1=1..nops(yeled))):


 kha:={seq(yeled[i][2],i=1..nops(yeled))}:
 
 if kha intersect AlreadyDone<>{} then
  print(`Note that, in the left side, the following  monomials`, op(kha intersect AlreadyDone), `are already treated, `):
 	     print(``):

 if kha minus AlreadyDone<>{} then
 print(`but we have to handle the new arrivals`, op(kha minus AlreadyDone), ` . `):
 fi:

 else
   if kha minus AlreadyDone<>{} then
     print(`We have to handle the new arrivals`, op(kha minus AlreadyDone)):
   fi:
 
fi:
 	     print(``):

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

if StillToDo<>{} then
 print(`Note that we still have to do handle the monomials in the set`, StillToDo ):	
else
 print(`Nothing left to do.`):
fi:


 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
print(`The scheme has more than`, K, `members `):	
 RETURN(FAIL):
fi:


print(`Summing up  we found the following scheme for the sequences `):
print(``):
print(seq(E[adam](k), adam in S)):
print(``):
print(`as follows `):
for adam in S do
  yeled:=T1[adam]:
  print(E[adam](k)=add(yeled[i1][1]*E[yeled[i1][2]](k-1),i1=1..nops(yeled))):
od:


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:

print(`In order to simplify notation, let `):
print(``):
print(seq(E[i1](k)=E[C1[i1]](k),i1=1..nops(S)))	:
print(``):
print(`Our scheme becomes `):

for i1 from 1 to nops(S) do
  yeled:=T2[i1]:
   print(E[i1](k)=add(yeled[i2][1]*E[yeled[i2][2]](k-1),i2=1..nops(yeled))):
  print(``):
od:

print(`Subject to the obvious initial conditions`):
print(``):
print(seq(E[i1](0)=V1[i1],i1=1..nops(S)))	:
print(``):

print(`Now let's try to find explicit expressions for the (ordinary) generating functions, in the variable `):
print(``):
print(F[i](t)=Sum(E[i](k)*t^k,k=0..infinity)):
print(``):
print(`For i from 1 to`, nops(S)):
print(`The above recurrences for the sequences`, E[i](k), `translate to the following system of linear equations for`):
print(F[i](t), `for i from 1 to`, nops(S)):
	       
for i1 from 1 to nops(S) do
 yeled:=T2[i1]:
print(F[i1](t)=V1[i1]+ t*(add(yeled[i2][1]*F[yeled[i2][2]](t),i2=1..nops(yeled)))):
od:


var:={seq(F[i1],i1=1..nops(S))}:

eq:={}:
for i1 from 1 to nops(S) do
 yeled:=T2[i1]:
eq:=eq union {F[i1]=V1[i1]+ t*(add(yeled[i2][1]*F[yeled[i2][2]],i2=1..nops(yeled)))}:
od:


var:=solve(eq,var):

print(`Solving this system gives explicit expressions for each of`, F[i](t), `and in particular, we find that our object of desire` ):
print(F[1](t)=factor(normal(subs(var,F[1])))):

print(``):
print(`This concludes the article that took`, time()-t0, `seconds to generate. `):
print(``):

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

end:


#GDv(F,a,b,A,B,z): Verbose version of GD(F,a,b,A,B,z) (q.v.). Try:
#GDv((a*A)^2,a,b,A,B,z);

GDv:=proc(F,a,b,A,B,z) local gu,i,k:

gu:=expand(subs({a=a+z*b,b=a-z*b, A=A+B/z, B=A-B/z},F));

print(`The monomial`, F, `is shortand for `):
print(subs({a=P[k](z),A=P[k](1/z),b=P[k](-z), B=P[k](-1/z)},F)):
print(``):
print(`Using the defining recurrence for the Shapiro polynomials, we can replace`):
print(a=a+z*b,b=a-z*b, A=A+B/z, B=A-B/z):
print(`where now a,b,A,B, stand for `, P[k-1](z^2),  P[k-1](1/z^2), P[k-1](-z^2),  P[k-1](-1/z^2), `respectively `):
print(``):
print(`we get`):
print(gu):
print(``):
print(`Discarding all odd powers of z, and replacing z^2 by z we get that it is`):
gu:=add(coeff(gu,z,2*i)*z^i,i=-trunc(-ldegree1(gu,z)/2)-1..trunc(degree1(gu,z)/2)+1):
print(``):
print(gu):
print(``):
print(`Replacing each term by  its canonical form (that is trivially the same, due to the symmetry of the dihedral group we have `):

gu:=Canon(gu,a,b,A,B,z):
print(gu):
gu:

end:

#HaroldSMvv(w,a,b,A,B,z,K): Super Verbose version of HaroldSM(w,a,b,A,B,z,K) (q.v.), including explanation of the Going Down process
#Try
#HaroldSMvv(a^2*A^2,a,b,A,B,z,10);
HaroldSMvv:=proc(w,a,b,A,B,z,K,t) local S,StillToDo,adam,yeled,i,AlreadyDone,T1,V1,C1,C1i,co,i1,j1,T2,ka,E,k,c1,c2,c3,c4,c5,kha,F,eq,var,t0,i2:
t0:=time():
print(`A Linear Recurrence Scheme for the Constant Term of `,subs({a=P[k](z),A=P[k](1/z),b=P[k](-z), B=P[k](-1/z)},w)):
print(`Where P[k](z) is the Rudin-Shapiro polynomial `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(` Let `, P[k](z), `be the Rudin-Shapiro polynomial, that may be defined by the recurrence`):
print(P[k](z)=P[k-1](z^2)+z*P[k-1](-z^2)):
print(``):
print(`and the initial condition`):
print(` P[0](z)=1 `):
print(``):
print(`Definition: For any monomial w, in a,b,A,B,z,`,  a^c1*A^c2*b^c3*B^c4*z^c5, `let `):
print(E[a^c1*A^c2*b^c3*B^c4*z^c5](k), `be the coefficient of z^0  in the polynomial `):
print(``):
print(P[k](z)^c1*P[k](1/z)^c2*P[k](-z)^c3*P[k](-1/z)^c4*z^c5):
print(``):
print(``):
print(`We are interested in getting a scheme for computing the sequence`,  E[w](k)):

S:={w}:

StillToDo:={w}:

co:=0:
AlreadyDone:={}:

while  nops(S)<=K and StillToDo<>{} do
 adam:=StillToDo[1]:
 print(`Let's express`, E[adam](k), `in terms of the values of other, related, sequences at` , k-1):	
 co:=co+1:
 C1[co]:=adam:
 C1i[adam]:=co:
  
 AlreadyDone:=AlreadyDone union {adam}:
 StillToDo:= StillToDo minus {adam}:

 yeled:=Monos(GD(adam,a,b,A,B,z),[a,b,A,B,z]):
  
 print(`Applying the defining recurrence, and using trivial equivalences, we get the following`):
  GDv(adam,a,b,A,B,z):
 print(``):
 print(`that implies`):
 print(``):
 print(E[adam](k)=add(yeled[i1][1]*E[yeled[i1][2]](k-1),i1=1..nops(yeled))):


 kha:={seq(yeled[i][2],i=1..nops(yeled))}:
 
 if kha intersect AlreadyDone<>{} then
  print(`Note that, in the left side, the following  monomials`, op(kha intersect AlreadyDone), `are already treated, `):
 	     print(``):

 if kha minus AlreadyDone<>{} then
 print(`but we have to handle the new arrivals`, op(kha minus AlreadyDone), ` . `):
 fi:

 else
   if kha minus AlreadyDone<>{} then
     print(`We have to handle the new arrivals`, op(kha minus AlreadyDone)):
   fi:
 
fi:
 	     print(``):

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

if StillToDo<>{} then
 print(`Note that we still have to do handle the monomials in the set`, StillToDo ):	
else
 print(`Nothing left to do.`):
fi:


 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
print(`The scheme has more than`, K, `members `):	
 RETURN(FAIL):
fi:


print(`Summing up  we found the following scheme for the sequences `):
print(``):
print(seq(E[adam](k), adam in S)):
print(``):
print(`as follows `):
for adam in S do
  yeled:=T1[adam]:
  print(E[adam](k)=add(yeled[i1][1]*E[yeled[i1][2]](k-1),i1=1..nops(yeled))):
od:


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:

print(`In order to simplify notation, let `):
print(``):
print(seq(E[i2](k)=E[C1[i2]](k),i2=1..nops(S)))	:
print(``):
print(`Our scheme becomes `):

for i1 from 1 to nops(S) do
  yeled:=T2[i1]:
   print(E[i1](k)=add(yeled[i2][1]*E[yeled[i2][2]](k-1),i2=1..nops(yeled))):
  print(``):
od:

print(`Subject to the obvious initial conditions`):
print(``):
print(seq(E[i1](0)=V1[i1],i1=1..nops(S)))	:
print(``):

print(`Now let's try to find explicit expressions for the (ordinary) generating functions, in the variable `):
print(``):
print(F[i](t)=Sum(E[i](k)*t^k,k=0..infinity)):
print(``):
print(`For i from 1 to`, nops(S)):
print(`The above recurrences for the sequences`, E[i](k), `translate to the following system of linear equations for`):
print(F[i](t), `for i from 1 to`, nops(S)):
	       
for i1 from 1 to nops(S) do
 yeled:=T2[i1]:
print(F[i1](t)=V1[i1]+ t*(add(yeled[i2][1]*F[yeled[i2][2]](t),i2=1..nops(yeled)))):
od:


var:={seq(F[i1],i1=1..nops(S))}:

eq:={}:
for i1 from 1 to nops(S) do
 yeled:=T2[i1]:
eq:=eq union {F[i1]=V1[i1]+ t*(add(yeled[i2][1]*F[yeled[i2][2]],i2=1..nops(yeled)))}:
od:


var:=solve(eq,var):

print(`Solving this system gives explicit expressions for each of`, F[i](t), `and in particular, we find that our object of desire` ):
print(F[1](t)=factor(normal(subs(var,F[1])))):

print(``):
print(`This concludes the article that took`, time()-t0, `seconds to generate. `):
print(``):


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

end: