######################################################################
##Cfinite.txt  Save this file as  Cfinite.txt                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Cfinite.txt<Enter>                            #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: July 2011
#Previous versions Feb. 8, 2013, Sept. 30, 2015
#This version: Nov. 1, 2022
Digits:=100:
 
print(`Created:  July 2011`):
print(`Previous version of  Feb. 2013`):
print(`This version: Sept. 30, 2015.` ):
print(``):
print(` This is Cfinite.txt `):
print(`It accompanyies the article `):
print(`The Algebra of C-finite Sequences`):
print(`by Doron Zeilberger`):
print(`dedicatd to Mourad Ismail and Dennis Stanton.`):
print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):


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


print(`-----------------------`):
 print(`For a list of the  procedures generating Diphantine equations type`):
print(` ezraDio();`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`-------------------------------`):

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

with(combinat):
with(numtheory):


ezraDio:=proc():
print(`The procdures to generate and solve Diophantine equations are:`):
print(` PellG, GoodQ2, SolPG, Theo2, findPells`):
end:

ezraFamous:=proc():
print(`The procedures coding famous sequences are`):
print(`fn,Fn, Ln, Pn, Ux,Trn, Tx`):
end:

ezraStory:=proc()

if args=NULL then
 print(` The Story procedures are: BTv,FactorizeI1v, FindCHv, GuessNLRv, `):
 print(` KefelV, mSectV, SeferGW `):
else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Agel, CartProd, CartProd2 `):
 print(` CHseq, Coe, CtoR, Curt1,`):
 print(` etop, GenHomPol, `):
 print(` GenPol,  FindCH1, FindCHg, GuessInNLR1, GuessNLR1, GuessNLRg1, GuessPR1`):
 print(` GuessRec1, GuessRec1mini,IISF, ISF, Krovim, ptoe `):
 print(`  ProdIndicator, ProdIndicatorG, ProdProfile, `):
 print(` ProdProfileE, RtoC, TDB`):

else
ezra(args):
fi:

end:


ezra:=proc()

if args=NULL then
 print(`The main procedures are: BT,  CH, DimerData, Factorize, FactorizeI1`):
 print(` FindCH,GuessNLR, GuessNLRg, GuessPR, GuessRec, IsProd, IsProdG `):
 print(` Kast, Kefel, KefelR, KefelS, KefelSk, KefelSm, Khibur, KhiburS`):
 print(` mSect, SeqFromRec `):
 print(` `):



elif nops([args])=1 and op(1,[args])=Agel then
print(`Agel(p,z): inputs a polynomial p in floating-point, rounds it. Try:`):
print(`Agel(1.0000001+3.000000001*z,z);`):

elif nops([args])=1 and op(1,[args])=BT then
print(`BT(C): The Binomial Transform of the C-finite sequence C,`):
print(`in other words, the sequence add(binomial(n,i)*C[i],i=0..n):`):
print(`For example, try:`):
print(`BT([[0,1],[1,1]]);`):

elif nops([args])=1 and op(1,[args])=BTv then
print(`BTv(C): Verbose form of BT(C)`):
print(`For example, try:`):
print(`BTv([[0,1],[1,1]]);`):

elif nops([args])=1 and op(1,[args])=CartProd then
print(`CartProd(L): The Cartesian product of the lists`):
print(`in the list of lists L`):
print(`for example, try:`):
print(`CartProd([[2,5],[4,9],[-2,-7]]);`):

elif nops([args])=1 and op(1,[args])=CartProd2 then
print(`CartProd2(L1,L2): The Cartesian product of lists L1 and L2,`):
print(`for example, try:`):
print(`CartProd2([2,5],[4,9]);`):

elif nops([args])=1 and op(1,[args])=CH then
print(`CH(L,n,j): inputs  a list of pairs [C,e], where C is a C-finite`):
print(`sequence, and e is an affine-linear expression in the symbols`):
print(`n and j, of the form a*n+b*j+c, where a and a+b are non-negative`):
print(`integers. Outputs the C-finite description of`):
print(`Sum(L[1][1](L[1][2])*L[2][1](L[2][2]),j=0..n-1)`):
print(`Implements the sequences discussed by Curtis Greene and`):
print(`Herb Wilf ("Closed form summation of C-finite sequences",`):
print(`Trans. AMS.)`):
print(`For example, try:`):
print(`CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j);`):

elif nops([args])=1 and op(1,[args])=CHseq then
print(`CHseq(L,n,j,K): inputs  a list of pairs [C,e], where C is a C-finite`):
print(`sequence, and e is an affine-linear expression in the symbols`):
print(`n and j, of the form a*n+b*j+c, where a and a+b are non-negative`):
print(`integers. Outputs the first K terms of`):
print(`Sum(L[1][1](L[1][2])*L[2][1](L[2][2]),j=0..n-1)`):
print(`using up to K values. For example, try:`):
print(`CHseq([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j,20);`):

elif nops([args])=1 and op(1,[args])=Coe then
print(`Coe(f,xList) : the set of coefficients of the polynomial f in the list of variables xList. Try:`):
print(`Coe((a*x+y*y)^5,[x,y]);`):

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])=Curt1 then
print(`Curt1(j,n,M): outputs `):
print(`all {a1n+b1j} with 0<=a1<=M and 0<=a1+b1<=M`):
print(`Type: Curt1(j,n,4);`):


elif nops([args])=1 and op(1,[args])=DimerData then
print(`DimerData(z): the list of length 10 whose i-th item is the C-finite representation (only the recurrence part)`):
print(` of the sequence of weight-enumerators of tilings of an i by m rectangle, with weight z^(NumberOfVerticalTiles). `):
print(` Try: `):
print(` DimerData(z); `):

elif nops([args])=1 and op(1,[args])=etop then
print(`etop(e): If e is a list of numbers (or expressions)`):
print(`representing the elementary symmetric functions`):
print(`outputs the power-symmetric functions`):
print(`For example, try`):
print(`etop([1,3,6]);`):

elif nops([args])=1 and op(1,[args])=Factorize then
print(`Factorize(C,L): Given a C-finite recurrence `):
print(`(without the initial conditions)`):
print(`tries to find a factorization into C-finite sequences of order`):
print(`L[i],i=1..nops(L).`):
print(`It returns all the factoriziations, if found, followed by `):
print(`the free parameters and FAIL otherwise`):
print(` For example, try:`):
print(`Factorize([-3,106,-72,-576],[2,2]);`):
print(`Factorize([2,7,2,-1],[2,2]);`):

elif nops([args])=1 and op(1,[args])=FactorizeI1 then
print(`FactorizeI1(C,k): Given a C finite integer sequence,C`):
print(`tries to factorize it into a C finite sequence of`):
print(`order k and order nops(C[1])-k. For example, try:`):
print(`FactorizeI1([[2,10,36,145],[2,7,2,-1]],2);`):
print(`FactorizeI1([[1,8,95,1183],[15,-32,15,-1]],2);`):
print(`FactorizeI1(RtoC(TDB(t)[3],t),2);`):

elif nops([args])=1 and op(1,[args])=FactorizeI1v then
print(`FactorizeI1v(C,k): verbose version of FactorizeI1(C,k)`):
print(` For example, try:`):
print(`FactorizeI1v([[2,10,36,145],[2,7,2,-1]],2);`):
print(`FactorizeI1v([[1,8,95,1183],[15,-32,15,-1]],2);`):
print(`FactorizeI1v(RtoC(TDB(t)[3],t),2);`):

elif nops([args])=1 and op(1,[args])=Fn then
print(`Fn(): The C-finite description of the (usual) Fibonacci sequence`):

elif nops([args])=1 and op(1,[args])=fn then
print(`fn(): The C-finite description of the fibonacci sequence `):
print(`fibonacci(n+1)`):

elif nops([args])=1 and op(1,[args])=GenHomPol then
print(`GenHomPol(Resh,a,deg): Given a list of variables Resh, a symbol a,`):
print(`and a non-negative integer deg, outputs a generic  HOMOGENEOUS `):
print(` polynomial`):
print(`of total degree deg,in the variables Resh, indexed by the letter a of`):
print(`followed by the set of coeffs. For example, try:`):
print(`GenHomPol([x,y],a,1);`):

elif nops([args])=1 and op(1,[args])=GenPol then
print(`GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a,`):
print(`and a non-negative integer deg, outputs a generic polynomial`):
print(`of total degree deg,in the variables Resh, indexed by the letter a of`):
print(`followed by the set of coeffs. For example, try:`):
print(`GenPol([x,y],a,1);`):


elif nops([args])=1 and op(1,[args])=FindCH then
print(`FindCH(C1,r,d, x,L,n,j): Inputs a C-finite sequence C1 and pos. integers`):
print(`r and d, and a letter x, and a list L of affine-linear expressions`):
print(`in n,j, and`):
print(`and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that `):
print(`P(C1[n],C1[n-1],C1[n-r])=C2[n]`):
print(`where C2(n) is the sum(C1[L[1]]*C1[L[2]*...C1[L[nops(L)],j=0..n-1):`):
print(`If none found it returns FAIL. For example to reproduce`):
print(`Eq. (12) of the Greene-Wilf article.`):
print(`If there is none of degree <=d then it returns FAIL`):
print(` type:`):
print(`FindCH(Fn(),1,3,x,[j,j,2*n-j],n,j);`):

elif nops([args])=1 and op(1,[args])=FindCHv then
print(`FindCHv(C1,r,d, x,L,n,j):  verbose version of`):
print(`FindCH(C1,r,d, x,L,n,j)`):
print(` type:`):
print(`FindCHv(Fn(),1,3,[j,j,2*n-j],n,j);`):

elif nops([args])=1 and op(1,[args])=FindCHg then
print(`FindCHg(C1,r,d,x,C2): Inputs a C-finite sequence C1 and`):
print(` a pos. integer`):
print(`d and outputs a polynomial P(x[1],..., x[r]) of minimal`):
print(` degree  such that `):
print(`P(C1[n],C1[n-1],C1[n-r])=C2[n]`):
print(`for all n. If none found it returns FAIL. For example to reproduce`):
print(`Eq. (12) of the Greene-Wilf article,  type:`):
print(`FindCHg([[0,1],[1,1]],1,3,x,CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j));`):


elif nops([args])=1 and op(1,[args])=GoodQ2 then
print(`GoodQ2(K): finds all quadruples [a0,a1,b0,b1] such that a0*b1-a1*b0=1 and 1<=a0,a1,b0<=K. Try:`):
print(`GoodQ2(K)`):

elif nops([args])=1 and op(1,[args])=GuessCH1 then
print(`GuessCH1(C1,r,d,x,C2): Inputs a C-finite sequence C1 and`):
print(` a pos. integer`):
print(`d and outputs a polynomial P(x[1],..., x[r]) of degree d such that `):
print(`P(C1[n],C1[n-1],C1[n-r])=C2[n]`):
print(`for all n. If none found it returns FAIL. For example to reproduce`):
print(`Eq. (12) of the Greene-Wilf article,  type:`):
print(`GuessCH1([[0,1],[1,1]],1,3,x,CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j));`):

elif nops([args])=1 and op(1,[args])=GuessInNLR1 then
print(`GuessInNLR1(C1,r,d,x): Inputs a C-finite sequence C and `):
print(` a pos. integer d and outputs an `):
print(` interesting polynomial`):
print(` P(x[0],x[1],..., x[r]) of degree d such that  `):
print(` of the form x[0]*P1(x[1], ..., x[r])+P2(x[1], ...,x[r])+c `):
print(` where P1,P2 are homog. of degree d, d-1 resp. such that `):
print(` P(C1[n],C1[n-1],C1[n-r])=0 `):
print(`Try: `):
print(` GuessInNLR1([[0,1],[3,-1]],2,2,x); `):

elif nops([args])=1 and op(1,[args])=GuessNLRv then
print(`GuessNLRv(C1,r,d): Verbose form of GuessNLR(C1,r,d,x)`):
print(`Inputs a C-finite sequence C and `):
print(`a pos. integer d and outputs a polynomial relation`):
print(` of todal degree<= d such that P(C1[n],C1[n-1],C1[n-r])=0`):
print(`for all n. If none found it returns FAIL. For example to solve`):
print(`the Hugh Montgomery proposed (but not selected) Putnam problem,`):
print(`(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,`):
print(`where F(n) is the Fibonacci sequence type`):
print(`GuessNLRv([[0,1],[1,1]],1,10);`):


elif nops([args])=1 and op(1,[args])=GuessNLRg then
print(`GuessNLRg(INI, REC,d,x): Inputs a a list of lists of length r (say) each list of lentgh k`):
print(`and a list REC of length r, say`):
print(`d and outputs a polynomial P(x[1],..., x[r]) of degree d such that `):
print(`P(C1[n],C2[n],Cr[n])=0`):
print(`for all n. `):
print(`where C1[n]=SeqFromRec([INI[1],REC]) etc.`):
print(`If none found it returns FAIL. For example to solve`):
print(`the Hugh Montgomery proposed (but not selected) Putnam problem,`):
print(`(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,`):
print(`where F(n) is the Fibonacci sequence type`):
print(`GuessNLRg([[0,1],[1,1]],[1,1],4,x);`):

elif nops([args])=1 and op(1,[args])=GuessNLR then
print(`GuessNLR(C1,r,d,x): Inputs a C-finite sequence C and `):
print(`a pos. integer d and outputs a polynomial P(x[0],..., x[r])`):
print(` of todal degree<= d such that P(C1[n],C1[n-1],C1[n-r])=0`):
print(`for all n. If none found it returns FAIL. For example to solve`):
print(`the Hugh Montgomery proposed (but not selected) Putnam problem,`):
print(`(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,`):
print(`where F(n) is the Fibonacci sequence type`):
print(`GuessNLR([[0,1],[1,1]],1,10,x);`):

elif nops([args])=1 and op(1,[args])=GuessNLR1 then
print(`GuessNLR1(C1,r,d,x): Inputs a C-finite sequence C and `):
print(`a pos. integer d and outputs a polynomial P(x[0],..., x[r])`):
print(` of degree d such that P(C1[n],C1[n-1],C1[n-r])=0`):
print(`for all n. If none found it returns FAIL. For example to solve`):
print(`the Hugh Montgomery proposed (but not selected) Putnam problem,`):
print(`(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,`):
print(`where F(n) is the Fibonacci sequence type`):
print(`GuessNLR1([[0,1],[1,1]],1,4,x);`):


elif nops([args])=1 and op(1,[args])=GuessNLRg then
print(`GuessNLRg(INI,REC,d,x): Inputs a C-finite recurrence REC and a list of length nops(REC), INI and a positive integer d`):
print(`d and outputs a polynomial P(x[1],..., x[r]) of degree <=d such that `):
print(`if x[1](n)=SeqFromRec([INI[1],REC]), ...,  x[r](n)=SeqFromRec([INI[r],REC]),  then`):
print(`P(x[1](n),.., x[r](n))=0`):
print(`for all n. If none found it returns FAIL. For example try:`):
print(`GuessNLRg([[1,k],[0,b]],[2*k,-1],2,x);`):
print(`GuessNLRg([[1,1,2], [0,1,1],[0,0,1]],[1,1,1],3,x);`):

elif nops([args])=1 and op(1,[args])=GuessPR then
print(`GuessPR(C1,C2,d,x,y): `):
print(`Inputs two C-finite sequences and a pos. integer`):
print(`d and outputs a polynomial P(x,y), of lowest possible degree`):
print(` <=d such that P(C1[n],C2[n])=0`):
print(`for all n. If none found it returns FAIL. For example to solve`):
print(`the Hugh Montgomery proposed (but not selected) Putnam problem,`):
print(`(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,`):
print(`where F(n) is the Fibonacci sequence type`):
print(`GuessPR([[0,1],[1,1]],[[1,1],[1,1]],10,x,y);`):

elif nops([args])=1 and op(1,[args])=GuessPR1 then
print(`GuessPR1(C1,C2,d,x,y): `):
print(`Inputs two C-finite sequences and a pos. integer`):
print(`d and outputs a polynomial P(x,y)`):
print(` of degree d such that P(C1[n],C2[n])=0`):
print(`for all n. If none found it returns FAIL. For example to solve`):
print(`the Hugh Montgomery proposed (but not selected) Putnam problem,`):
print(`(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,`):
print(`where F(n) is the Fibonacci sequence type`):
print(`GuessPR1([[0,1],[1,1]],[[1,1],[1,1]],4,x,y);`):

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 cofficients`):
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])=GuessRec1 then
print(`GuessRec1(L,d): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant cofficients of order d`):
print(`satisfying it. It returns the initial d values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec1([1,1,1,1,1,1],1);`):

elif nops([args])=1 and op(1,[args])=GuessRec1mini then
print(`GuessRec1mini(L,d): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant cofficients of order d`):
print(`satisfying it. It returns the initial d values and the operator`):
print(`It is like GuessRec1(L,d) but allows shorter lists`):
print(` as a list. For example try:`):
print(`GuessRec1mini([1,1,1,1,1,1],1);`):


elif nops([args])=1 and op(1,[args])=IISF then
print(`IISF(L,k): given a sequence of integers finds all its`):
print(`Increasing integer`):
print(`factors of length k. For example try:`):
print(`IISF([1,4,9,16,32],3);`):

elif nops([args])=1 and op(1,[args])=ISF then
print(`ISF(L,k): given a sequence of integers finds all its`):
print(`factors of length k. For example try:`):
print(`ISF([1,4,9,16,32],3);`):

elif nops([args])=1 and op(1,[args])=IsProd then
print(`IsProd(f,t,m,n,eps): Given a rational function f of t`):
print(`with NUMERICAL coefficients`):
print(`decides, empirically,`):
print(` whether the C-finite sequence of its coeeficients`):
print(`is a product of a C-finite sequences of order m and n.`):
print(`For example, try:`):
print(`IsProd(1/((1-t)*(1-2*t)*(1-3*t)*(1-6*t)),t,2,2,1/1000);`):


elif nops([args])=1 and op(1,[args])=IsProdG then
print(`IsProdG(f,t,L,eps): Given a rational function f of t`):
print(`with NUMERICAL coefficients`):
print(`and a list of positive integers, L`):
print(`decides, empirically,`):
print(`decides whether the C-finite sequence of its coeeficients`):
print(`is a product of a C-finite sequences of orders given by the`):
print(`members of L`):
print(`For example, try:`):
print(`IsProdG(1/((1-t)*(1-2*t)*(1-3*t)*(1-6*t)),t,[2,2],1/1000);`):
print(`IsProdG(CtoR([[2,10,36,145],[2,7,2,-1]],t),t,[2,2],1/100000);`):
print(`IsProdG(TDB(t)[1],t,[2,2],1/100000);`):
print(`IsProdG(TDB(t)[2],t,[2,2,2],1/100000);`):
print(`IsProdG(TDB(t)[3],t,[2,2,2,2],1/100000);`):
print(`IsProdG(TDB(t)[4],t,[2,2,2,2,2],1/100000);`):


elif nops([args])=1 and op(1,[args])=KastF then
print(`KastF(n,z): The C-finite representation (using numeric-floating point)`):
print(`of the Kastelian solution to the Dimer problem for n by m rectangles`):
print(`where z->1, z'->z. Try:`):
print(`KastF(6,z); `):


elif nops([args])=1 and op(1,[args])=KastF1 then
print(`KastF1(n,z0):  Like KastF1(n,z) (q.v) but now z0 is a number, not a symbol. Try: `):
print(`KastF1(6,3); `):

elif nops([args])=1 and op(1,[args])=Kefel then
print(`Kefel(S1,S2): inputs two C-finite sequences S1,S2, and  outputs`):
print(`their product, For example, try:`):
print(`Kefel([[1,2],[3,4]],[[2,5],[-1,6]]);`):

elif nops([args])=1 and op(1,[args])=KefelV then
print(`KefelV(C1,C2): Verbose version of Kefel(S1,S2);`):
print(`For example try:`):
print(`KefelV([[1,2],[3,4]],[[2,5],[-1,6]]);`):

elif nops([args])=1 and op(1,[args])=KefelR then
print(`KefelR(R1,R2,t): inputs rational functions of t,`):
print(`R1 and R2,`):
print(`and outputs the rational function `):
print(`whose coeffs. is their Hadamard product, for example try:`):
print(`KefelR(1/(1-t-t^2),1/(1-t-t^3),t);`):


elif nops([args])=1 and op(1,[args])=KefelS then
print(`KefelS(a,b,d1,d2): The recurrence operator`):
print(`satisfied by the product of any solution of`):
print(`the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]`):
print(`and any solution of the recurrence `):
print(` g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1] `):
print(`For example, try:`):
print(`KefelS(a,b,2,2);`):

elif nops([args])=1 and op(1,[args])=KefelSk then
print(`KefelSk(a,b): Inputs lists of numbers (or symbols), a, and b, and outputs the recurrence operator`):
print(` with constant coefficients satisfied by the product of any solution of`):
print(`the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]`):
print(`and any solution of the recurrence `):
print(` g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1] `):
print(`For example, try:`):
print(`KefelSk([a[1],a[2]],[b[1],b[2]]);`):

elif nops([args])=1 and op(1,[args])=KefelSm then
print(`KefelSm(a,L): The recurrence operator`):
print(`satisfied by the product of nops(L) recurrences`):
print(`the recurrence f[i][n]=a[i,1]*f[n-1]+...+a[i,L[i]]*f[n-L[i]]`):
print(`It also returns the set of generic coefficients.`):
print(`For example, try:`):
print(`KefelSm(a,[2,2,2]);`):

elif nops([args])=1 and op(1,[args])=Khibur then
print(`Khibur(S1,S2): inputs two C-finite sequences S1,S2, and  outputs`):
print(`their sum, For example, try:`):
print(`Khibur([[1,2],[3,4]],[[2,5],[-1,6]]);`):


elif nops([args])=1 and op(1,[args])=KhiburS then
print(`KhiburS(a,b,d1,d2): The recurrence operator`):
print(`satisfied by the sum of any solution of`):
print(`the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]`):
print(`and any solution of the recurrence `):
print(` g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1]`):
print(`For example, try:`):
print(`KhiburS(a,b,2,2);`):

elif nops([args])=1 and op(1,[args])=Krovim then
print(`Krovim(L,i,eps): Given a list L and an index between`):
print(`1 and nops(L), finds the set of indices (including i)`):
print(`such that abs(L[j]-L[i])<eps. For example, try:`):
print(`Krovim([1.0001,1.01,1.00011],1,1/1000);`):

elif nops([args])=1 and op(1,[args])=Ln then
print(`Ln(): The C-finite description of the Lucas number sequence`):

elif nops([args])=1 and op(1,[args])=mSect then
print(`mSect(C1,m,i): Given a C-finite sequence C1(n), finds`):
print(`the C-finite description of C(m*n+i) where m is a pos. integer`):
print(`and 0<=i<m. For example, try:`):
print(`mSect([[0,1],[1,1]],2,0);`):

elif nops([args])=1 and op(1,[args])=mSectV then
print(`mSectV(C1,m,i): Verbose form of mSect(C1,m,i)`):
print(`Given a C-finite sequence C1(n), finds`):
print(`the C-finite description of C(m*n+i) where m is a pos. integer`):
print(`and 0<=i<m. For example, try:`):
print(`mSectV([[0,1],[1,1]],2,0);`):

elif nops([args])=1 and op(1,[args])=PellG then
print(`PellG(a0,a1,b0,b1,x,y): The Pell-like Diophantine equation with parameters a0,a1,b0,b1. Try:`):
print(`PellG(5,6,6,7,x,y);`):
print(`PellG(a0,a1,b0,b1,x,y);`):


elif nops([args])=1 and op(1,[args])=Pn then
print(`Pn(): The C-finite description of the Pell Numbers `):

elif nops([args])=1 and op(1,[args])=ptoe then
print(`ptoe(p): If p is a list of numbers (or expressions)`):
print(`representing the power-sum symmetric functions`):
print(`outputs the elementary symmetic functions`):
print(`For example, try`):
print(`ptoe([1,3,6]);`):

elif nops([args])=1 and op(1,[args])=ProdIndicator then
print(`ProdIndicator(m,n): The profile of multiplicity of`):
print(`the (m*n)^2 ratios of a list of m*n numbers to`):
print(`be the Cartesian product of a set of m numbers with`):
print(`a set of n numbers. For example, try:`):
print(`ProdIndicator(2,2);`):

elif nops([args])=1 and op(1,[args])=ProdIndicatorG then
print(`ProdIndicatorG(L): The profile of multiplicity of`):
print(`the convert(L,`*`)^2 ratios of a list of `):
print(`be the Cartesian product of lists of size given by L`):
print(`For example, try:`):
print(`ProdIndicatorG([2,2]);`):

elif nops([args])=1 and op(1,[args])=ProdProfile then
print(`ProdProfile(L,eps): Given a list of complex numbers finds its`):
print(`(approximate) profile to within eps. For example, try:`):
print(`ProdProfile([1,2,3,6],1/1000);`):

elif nops([args])=1 and op(1,[args])=ProdProfileE then
print(`ProdProfileE(L): Given a list of numbers (or expressions) finds its`):
print(`(exact) profile to within eps. For example, try:`):
print(`ProdProfileE([1,2,3,6]);`):



elif nops([args])=1 and op(1,[args])=RtoC then
print(`RtoC(R,t): Given a rational function R(t) finds the`):
print(`C-finite description of its sequence of coefficients`):
print(`For example, try:`):
print(`RtoC(1/(1-t-t^2),t);`):


elif nops([args])=1 and op(1,[args])=SeferGW then
print(`SeferGW(m,d,M): outputs a book of Fibonacci identities with`):
print(`right hand sides of degree at most d in F(n),F(n+1)`):
print(`for Greene-Wilf type sums where the summand is`):
print(`of the form F(j)*F(a1n+b1j)*... with m terms, for example, try:`):
print(`SeferGW(2,5,2);`):


elif nops([args])=1 and op(1,[args])=SeqFromRec then
print(`SeqFromRec(S,N0): Given a C-finite sequence S, written as`):
print(`S=[INI,ope], where INI is the list of initial conditions`):
print(`ope is the  recurrence operator (linear with constant coefficients)`):
print(`codes a list of the same size as INI`):
print(`(the recurrence is: f(n)=ope[1]*f(n-1)+ope[2]*f(n-2)+ ...) , `):
print(` finds the first N0 terms of the sequence satisfying`):
print(`the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).`):
print(`For example, for the first 20 Fibonacci numbers,`):
print(`try: SeqFromRec([[1,1],[1,1]],20);`):


elif nops([args])=1 and op(1,[args])=SolPG then
print(`SolPG(a0,a1,b0,b1,t): The pair of rational functions [f1(t),f2(t)] such that its respective Maclaurin coefficients solve the generalized Pell equation. Try:`):
print(`SolPG(3,5,11,17,t);`):

elif nops([args])=1 and op(1,[args])=TDB then
print(`TDB(t): the list of rational functions in t`):
print(`of length 4, consisting of the generating functions`):
print(`for the number of perfect matchings in a `):
print(`4 by n, 6 by n, 8 by n and 10 by n`):
print(`Do TDB(t);`):

elif nops([args])=1 and op(1,[args])=Theo2 then
print(`Theo2(a0,a1,b0,b1,x,y): outputs a theorem about the generalized Pell equation coming from C-finite sequenes of order 2. Try:`):
print(`Theo2(5,6,6,7,x,y);`):
print(`Theo2(a0,a1,b0,b1,x,y):`):

elif nops([args])=1 and op(1,[args])=findPells then
print(`findPells(K, x, y): find Pell equations of the form x^2 - d y^2 = m coming from the C-finite ansatz.`):
print(`Looks at all parameters in GoodQ3(K)`):
print(`Output is a set of lists [ics, rec, P].`):
print(`Try:`):
print(`findPells(8, x, y)`):


elif nops([args])=1 and op(1,[args])=Trn then
print(`Trn(): The C-finite description of the Triboanacci numbers Numbers`):

elif nops([args])=1 and op(1,[args])=Tx then
print(`Tx(x): The C-finite description of the Chebychev polynomials of`):
print(`the first kind. Try: Tx(x);`):

elif nops([args])=1 and op(1,[args])=Ux then
print(`Ux(x): The C-finite description of the Chebychev polynomials of`):
print(`the second kind. Try: Ux(x);`):

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


TDB:=proc(t):
[-(1+t)/(-5*t^2-t+1-t^3+t^4)*(t-1), -(t^6-8*t^4+2*t^3+8*t^2-1)/(1+t)/(t-1)/(t^6+t^5-19*t^4+11*t^3+19*t^2+t-1), -(-1-1033*t^8+110*t^9+26*t^3+43*t^2-360*t^4-26*t^11+t
^14+1033*t^6-43*t^12-110*t^5+360*t^10)/(t^16-t^15-76*t^14-69*t^13+921*t^12+584*t^11-4019*t^10-829*t^9+7012*t^8-829*t^7-4019*t^6+584*t^5+921*t^4-69*t^3-76*t^2-t+1),
-(-1+t^30-2064705*t^8+156848*t^9+10324398*t^13+214*t^3+197*t^2-54699758*t^16-15137114*t^15-9741*t^4-2972710*t^11+54699758*t^14+202037*t^6+56736*t^7-32425754*t^12-\
7262*t^5+11058754*t^10-197*t^28+214*t^27+9741*t^26-202037*t^24-7262*t^25+10324398*t^17-2972710*t^19+32425754*t^18-11058754*t^20+2064705*t^22+156848*t^21+56736*t^23)
/(1-t+t^32+411*t^29-285*t^30+t^31+6149853*t^8+471319*t^9-58545372*t^13-411*t^3-285*t^2+429447820*t^16+123321948*t^15+18027*t^4+10402780*t^11-335484428*t^14-472275*t
^6-271027*t^7+157353820*t^12+20689*t^5-42303393*t^10+18027*t^28-20689*t^27-472275*t^26+6149853*t^24+271027*t^25-123321948*t^17+58545372*t^19-335484428*t^18+
157353820*t^20-42303393*t^22-10402780*t^21-471319*t^23)]:
end:

#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 cofficients 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+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 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..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#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)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:


#Khibur(S1,S2): inputs two C-finite sequences S1,S2, and  outputs
#their sum, For example, try:
#Khibur([[1,2],[3,4]],[[2,5],[-1,6]]);
Khibur:=proc(S1,S2) local d1,d2,L1,L2,L,i,K:

if nops(S1[1])<>nops(S1[2]) or nops(S2[1])<>nops(S2[2]) then
print(`Bad input`):
RETURN(FAIL):
fi:

d1:=nops(S1[1]): d2:=nops(S2[1]):
K:=2*(d1+d2)+4:

L1:=SeqFromRec(S1,K):
L2:=SeqFromRec(S2,K):

L:=[seq(L1[i]+L2[i],i=1..K)]:

GuessRec(L):

end:



#Kefel(S1,S2): inputs two C-finite sequences S1,S2, and  outputs
#their product, For example, try:
#Kefel([[1,2],[3,4]],[[2,5],[-1,6]]);
Kefel:=proc(S1,S2) local d1,d2,L1,L2,L,i,K:

if nops(S1[1])<>nops(S1[2]) or nops(S2[1])<>nops(S2[2]) then
print(`Bad input`):
RETURN(FAIL):
fi:

d1:=nops(S1[1]): d2:=nops(S2[1]):
K:=2*(d1*d2)+4:

L1:=SeqFromRec(S1,K):
L2:=SeqFromRec(S2,K):

L:=[seq(L1[i]*L2[i],i=1..K)]:

GuessRec(L):

end:



#KhiburSslow(a,b,d1,d2): The recurrence operator
#satisfied by the sum of any solution of
#the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]
#and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1]
#For example, try:
#KhiburSslow(a,b,2,2);
KhiburSslow:=proc(a,b,d1,d2) local S1,S2,i:

S1:=[[seq(ithprime(5+i),i=1..d1)],[seq(a[i],i=1..d1)]]:

S2:=[[seq(ithprime(9+d1+i),i=1..d2)],[seq(b[i],i=1..d2)]]:

Khibur(S1,S2)[2]:

end:



#KefelSslow(a,b,d1,d2): The recurrence operator
#satisfied by the product of any solution of
#the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]
#and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1]
#For example, try:
#KefelSslow(a,b,2,2);
KefelSslow:=proc(a,b,d1,d2) local S1,S2,i:

S1:=[[seq(ithprime(5+i),i=1..d1)],[seq(a[i],i=1..d1)]]:

S2:=[[seq(ithprime(9+d1+i),i=1..d2)],[seq(b[i],i=1..d2)]]:

Kefel(S1,S2)[2]:

end:


#KhiburS(a,b,d1,d2): The recurrence operator
#satisfied by the sum of any solution of
#the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]
#and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1]
#For example, try:
#KhiburS(a,b,2,2);
KhiburS:=proc(a,b,d1,d2) local S1,S2,i,x,y,L,C,var,eq,L1,L2,gu:

S1:=[[seq(x[-d1+i],i=0..d1-1)],[seq(a[i],i=1..d1)]]:

S2:=[[seq(y[-d2+i],i=0..d2-1)],[seq(b[i],i=1..d2)]]:

L1:=SeqFromRec(S1,d1+d2+1):
L2:=SeqFromRec(S2,d1+d2+1):

L:=expand([seq(L1[i]+L2[i],i=1..d1+d2+1)]):

gu:=expand(L[d1+d2+1]-add(C[i]*L[d1+d2+1-i],i=1..d1+d2)):

var:={seq(C[i],i=1..d1+d2)}:
eq:={seq(coeff(gu,x[-i],1),i=1..d1),seq(coeff(gu,y[-i],1),i=1..d2)}:
var:=solve(eq,var):

[seq(subs(var,C[i]),i=1..d1+d2)]:


end:




#KefelS(a,b,d1,d2): The recurrence operator
#satisfied by the product of any solution of
#the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]
#and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1]
#For example, try:
#KefelS(a,b,2,2);
KefelS:=proc(a,b,d1,d2) local S1,S2,i,x,y,L,C,var,eq,L1,L2,gu,j:

S1:=[[seq(x[-d1+i],i=0..d1-1)],[seq(a[i],i=1..d1)]]:

S2:=[[seq(y[-d2+i],i=0..d2-1)],[seq(b[i],i=1..d2)]]:

L1:=SeqFromRec(S1,d1*d2+1):
L2:=SeqFromRec(S2,d1*d2+1):

L:=expand([seq(L1[i]*L2[i],i=1..d1*d2+1)]):

gu:=expand(L[d1*d2+1]-add(C[i]*L[d1*d2+1-i],i=1..d1*d2)):

var:={seq(C[i],i=1..d1*d2)}:
eq:={seq(seq(coeff(coeff(gu,x[-i],1),y[-j],1),j=1..d2),i=1..d1)}:
var:=solve(eq,var):

[seq(subs(var,C[i]),i=1..d1*d2)]:


end:


#ProdIndicatorN(m,n): The profile of multiplicity of
#the (m*n)^2 ratios of a list of m*n numbers to
#be the Cartesian product of a set of m numbers with
#a set of n numbers. For example, try:
#ProdIndicatorN(2,2);
ProdIndicatorN:=proc(m,n) local a,b,i,j,gu,x,lu,mu:
option remember:
gu:=[seq(seq(a[i]*b[j],j=1..n),i=1..m)]:
ProdProfileE(gu):
end:

#ProdIndicator(m,n): The profile of multiplicity of
#the (m*n)^2 ratios of a list of m*n numbers to
#be the Cartesian product of a set of m numbers with
#a set of n numbers. For example, try:
#ProdIndicator(2,2);
ProdIndicator:=proc(m,n) local a,b,i,j,gu,x,lu,mu:
option remember:
gu:=[seq(seq(a[i]*b[j],j=1..n),i=1..m)]:
lu:=[seq(seq(gu[i]/gu[j],j=1..m*n),i=1..m*n)]:
lu:=add(x[lu[i]],i=1..nops(lu)):

mu:=[]:

for i from 1 to nops(lu) do
 if type(op(1,op(i,lu)),integer) then
   mu:=[op(mu),op(1,op(i,lu))]:
 else:
   mu:=[op(mu),1]:
 fi:
od:

sort(mu):

end:


#ProdProfile(L,eps): Given a list of complex numbers finds its
#(approximate profile). For example, try:
#ProdProfile([1,2,3,6],10^(-4));
ProdProfile:=proc(L,eps) local i,j,lu,mu,S,gu:

gu:=[seq(seq(L[i]/L[j],j=1..nops(L)),i=1..nops(L))]:
mu:=[]:
S:={seq(i,i=1..nops(gu))}:

while S<>{} do
 lu:=Krovim(gu,S[1],eps):
 mu:=[op(mu),nops(lu)]:
 S:=S minus lu:
od:
sort(mu):
end:

#IsProd(f,t,m,n,eps): Given a rational function f of t
#decides whether the C-finite sequence of its coeeficients
#is a product of a C-finite sequences of order m and n.
#For example, try:
#IsProd(1/((1-t)*(1-2*t)*(1-3*t)*(1-6*t)),t,2,2);
IsProd:=proc(f,t,m,n,eps) local lu,gu:
lu:=denom(f):

if degree(lu,t)<>m*n then
 print(`The degree of the denominator of`, f, `should be`, m*n):
 RETURN(FAIL):
fi:

gu:=evalf([solve(lu,t)]):


if ProdProfile(gu,eps)=ProdIndicator(m,n) then
 RETURN(true):
else

 RETURN(false,gu):
fi:

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:=degree(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(f):
fi:
end:


#RtoC(R,t): Given a rational function R(t) finds the
#C-finite description of its sequence of coefficients
#For example, try:
#RtoC(1/(1-t-t^2),t);
RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a:
R1:=normal(R):
D1:=denom(R1):
d:=degree(D1,t):
if degree(numer(R1),t)>=d then
print(`The degree of the top must be less than the degree of the bottom`):
 RETURN(FAIL):
fi:

a:=coeff(D1,t,0):
S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]:
f1:=taylor(R,t=0,d+1):
S1:=[seq(coeff(f1,t,i),i=0..d-1)]:
[S1,S2]:
end:



#KefelR(R1,R2,t): inputs rational functions of t,
#R1 and R2,
#and outputs the rational function 
#whose coeffs. is their Hadamard product, for example try:
#KefelR(1/(1-t-t^2),1/(1-t-t^3),t);
KefelR:=proc(R1,R2,t) local S1,S2,S:
S1:=RtoC(R1,t):
S2:=RtoC(R2,t):
S:=Kefel(S1,S2):
normal(CtoR(S,t)):
end:





#Krovim(L,i,eps): Given a list L and an index between
#1 and nops(L), finds the set of indices (including i)
#such that abs(L[j]-L[i])<eps. For example, try:
#Krovim([1.0001,1.01,1.00011],1,1/1000);
Krovim:=proc(L,i,eps) local gu,j,a:
gu:={}:
a:=L[i]:

for j from 1 to nops(L) do
 if abs(L[j]-a)<eps then
  gu:=gu union {j}:
 fi:
od:
gu:
end:

#KrovimE(L,i): Given a list L and an index between
#1 and nops(L), finds the set of indices (including i)
#such that L[i]=L[j] For example, try:
#KrovimE([1,2,1,3],1);
KrovimE:=proc(L,i) local gu,j,a:
gu:={}:
a:=L[i]:

for j from 1 to nops(L) do
 if L[j]=a then
  gu:=gu union {j}:
 fi:
od:
gu:
end:


#ProdProfileE(L): Given a list of numbers (or symbols) finds its
#exact profile. For example, try:
#ProdProfileE([1,2,3,6]);
ProdProfileE:=proc(L) local i,j,lu,mu,S,gu:

gu:=[seq(seq(normal(L[i]/L[j]),j=1..nops(L)),i=1..nops(L))]:
mu:=[]:
S:={seq(i,i=1..nops(gu))}:

while S<>{} do
 lu:=KrovimE(gu,S[1]):
 mu:=[op(mu),nops(lu)]:
 S:=S minus lu:
od:
sort(mu):
end:


#etop(e): If e is a list of numbers (or expressions)
#representing the elementary symmetric functions
#outputs the power-symmetric functions
#For example, try
#etop([1,3,6]);
etop:=proc(e) local k,p,i:

for k from 1 to nops(e) do
 p[k]:=expand(add((-1)^(i-1)*e[i]*p[k-i],i=1..k-1)+(-1)^(k-1)*k*e[k]):
od:

[seq(p[k],k=1..nops(e))]:

end:


#ptoe(p): If p is a list of numbers (or expressions)
#representing the power-sum symmetric functions
#outputs the elementary symmetric functions
#For example, try
#ptoe([1,3,6]);
ptoe:=proc(p) local k,e,i:

for k from 1 to nops(p) do
 e[k]:=expand(add((-1)^(i-1)*p[i]*e[k-i]/k,i=1..k-1)+(-1)^(k-1)*p[k]/k):
od:

[seq(e[k],k=1..nops(p))]:

end:




#IsProdN(f,t,m,n): Given a rational function f of t
#decides whether the C-finite sequence of its coeeficients
#is a product of a C-finite sequences of order m and n.
#For example, try:
#IsProd(1/((1-t)*(1-2*t)*(1-3*t)*(1-6*t)),t,2,2);
IsProdN:=proc(f,t,m,n) local lu,k,e,p,i:
lu:=denom(f):
k:=degree(lu,t):
if k<>m*n then
print(`The denominator of `, f, `should have been of degree`, m*n ):
print(`but it is in fact of degree`, k):
RETURN(FAIL):
fi:

lu:=expand(lu/coeff(lu,t,k)):
e:=[seq((-1)^i*coeff(lu,t,k-i),i=1..k)]:

p:=etop(e):
RETURN(p):
ProdProfileE(p), ProdIndicator(m,n):
end:



#CartProd2(L1,L2): The Cartesian product of lists L1 and L2,
#for example, try:
#CartProd2([2,5],[4,9]);
CartProd2:=proc(L1,L2) local i,j:
[seq(seq(L1[i]*L2[j],i=1..nops(L1)),j=1..nops(L2))]:
end:


#CartProd(L): The Cartesian product of the lists in the list
#of lists L
#for example, try:
#CartProd([[2,5],[4,9],[4,7]]);
CartProd:=proc(L) local L1,L2,i:
if not type(L,list) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {seq(type(L[i],list),i=1..nops(L))}<>{true} then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(L)=1 then
RETURN(L[1]):
fi:


L1:=CartProd([op(1..nops(L)-1,L)]):
L2:=L[nops(L)]:
CartProd2(L1,L2):

end:



#ProdIndicatorG(L): The profile of multiplicity of
#the convert(L,`*`)^2 ratios of a list of 
#be the Cartesian product of lists of size given by L
#For example, try:
#ProdIndicatorG([2,2]);
ProdIndicatorG:=proc(L) local a,i,j,gu,x,lu,mu,M:
option remember:
M:=[seq([seq(a[i][j],j=1..L[i])],i=1..nops(L))]:
gu:=CartProd(M):

lu:=[seq(seq(gu[i]/gu[j],j=1..nops(gu)),i=1..nops(gu))]:
lu:=add(x[lu[i]],i=1..nops(lu)):

mu:=[]:

for i from 1 to nops(lu) do
 if type(op(1,op(i,lu)),integer) then
   mu:=[op(mu),op(1,op(i,lu))]:
 else:
   mu:=[op(mu),1]:
 fi:
od:

sort(mu):

end:



#IsProdG(f,t,L,eps): Given a rational function f of t
#decides whether the C-finite sequence of its coeeficients
#is a product of a C-finite sequences of the orders given
#by the list of positive integers L
#For example, try:
#IsProdG(1/((1-t)*(1-2*t)*(1-3*t)*(1-6*t)),t,[2,2],10^(-10));
IsProdG:=proc(f,t,L,eps) local lu,gu,ku1,ku2:
lu:=denom(f):

if degree(lu,t)<>convert(L,`*`) then
 print(`The degree of the denominator of`, f, `should be`, convert(L,`*`)):
 RETURN(FAIL):
fi:

gu:=evalf([solve(lu,t)]):


ku1:=ProdProfile(gu,eps):
ku2:=ProdIndicatorG(L):

if ku1=ku2 then
 RETURN(true):
else

 RETURN(false,[ku1,ku2]):
fi:

end:




#KefelSm(a,L): The recurrence operator
#satisfied by the product of nops(L) recurrences
#the recurrence f[i][n]=a[i,1]*f[n-1]+...+a[i,L[i]]*f[n-L[i]]
#It also returns the set of generic coefficients
#For example, try:
#KefelSm(a,[2,2,2]);
KefelSm:=proc(a,L) local gu,L1,mu,A,B,i,i1,j:

if not (type(a,symbol) and type(L,list) ) then
print(`Bad input`):
 RETURN(FAIL):
fi:

if not {seq(type(L[i],integer),i=1..nops(L))}={true} then
print(`Bad input`):
 RETURN(FAIL):
fi:

if not {seq(evalb(L[i]>=2),i=1..nops(L))}={true} then
print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(L)=1 then
 RETURN([seq(a[1,j],j=1..L[1])]):
fi:

L1:=[op(1..nops(L)-1,L)]:
gu:=expand(KefelSm(a,L1)):


mu:=KefelS(A,B,convert(L1,`*`),L[nops(L)]):
mu:=subs({seq(A[i1]=gu[i1],i1=1..nops(gu)),
seq(B[i1]=a[nops(L),i1],i1=1..L[nops(L)])},mu):
mu:=expand(mu):
mu,{seq(seq(a[i,j],j=1..L[i]),i=1..nops(L))}:

end:

#Factorize(C,L): Given a C-finite recurrence (without the initial conditions)
#tries to find a factorization into C-finite sequences of order
#L[i],i=1..nops(L). 
#It returns the factoriziation, if found, followed by the free parameters
#and FAIL otherwise
#For example, try:
#Factorize([-3,106,-72,-576],[2,2]);
Factorize:=proc(C,L) local gu,a,eq,var,i,var1,j,lu,k,lu1,mu:
if nops(C)<>convert(L,`*`) then
 print(`The number of elements of `, C):
 print(`should be`, convert(L,`*`)):
  ERROR(`Bad input`):
fi:

gu:=KefelSm(a,L):
var:=gu[2]:
gu:=gu[1]:
eq:={seq(gu[i]=C[i],i=1..nops(gu))}:
var1:=[solve(eq,var)]:

if var1=[] then
print(`There is no  solution`):
RETURN(FAIL):
fi:

for k from 1 to nops(var1) do
var1:=var1[k]:
lu:={}:

for i from 1 to nops(var1) do
 if op(1,op(i,var1))=op(2,op(i,var1)) then
   lu:=lu union {op(1,op(i,var1))}:
 fi:
od:



if nops(lu)>0 then
 mu:=[[seq([seq(subs(var1,a[i,j]),j=1..L[i])],i=1..nops(L))]]:
 mu:=subs({seq(lu1=1,lu1 in lu)},mu):
 RETURN(mu):
fi:
od:

FAIL:
end:




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

if nops(L)<2*d then
 print(`The list must be of size >=`, 2*d ):
 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..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#ISF(L,k): given a sequence of integers finds all its
#factors of length k. For example try:
#ISF([1,4,9,16,32],3);
ISF:=proc(L,k) local gu,gu1,lu,lu1:
if k<0 then
 RETURN({}):
fi:

if k>nops(L) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if k=0 then
 RETURN({[]}):
fi:
gu:=ISF(L,k-1):
lu:=divisors(L[k]):

{seq(seq([op(gu1),lu1],lu1 in lu),gu1 in gu)}:


end:



#IISF(L,k): given a sequence of integers finds all its
#increasing sequences
#factors of length k. For example try:
#IISF([1,4,9,16,32],3);
IISF:=proc(L,k) local gu,gu1,lu,lu1,mu:
if k<0 then
 RETURN({}):
fi:

if k>nops(L) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if k=0 then
 RETURN({[]}):
fi:
lu:=divisors(L[k]):
if k=1 then
 RETURN({seq([lu1],lu1 in lu)}):
fi:

gu:=IISF(L,k-1):


mu:={}:
for gu1 in gu do
 for lu1 in lu do
  if gu1[nops(gu1)]<lu1 then
    mu:=mu union {[op(gu1),lu1]}:
  fi:
 od:
od:
mu:

end:


#FactorizeI1(C,k): Given a C finite integer sequence,C
#tries to factorize it into a C finite sequence of
#order k and order nops(C[1])-k. For example, try:
#FactorizeI1(RtoC(TDB(t)[1],t),2);
FactorizeI1:=proc(C,k) local gu,mu,mu1,lu,gu1,gu11,C1,i1:
gu:=SeqFromRec(C,4*k):
mu:=ISF(gu,2*k) minus {[1$(2*k)]}:

for mu1 in mu do
 lu:=GuessRec1mini(mu1,k):

  if lu<>FAIL then
   gu1:=SeqFromRec(lu,4*k):
     if {seq(type(gu1[i1],integer),i1=1..nops(gu1))}={true} then
      gu11:=[seq(gu[i1]/gu1[i1],i1=1..4*k)]:
         if {seq(type(gu11[i1],integer),i1=1..nops(gu11))}={true} then
          gu:=SeqFromRec(C,10*k):
          gu1:=SeqFromRec(lu,10*k):
          gu11:=[seq(gu[i1]/gu1[i1],i1=1..10*k)]:
           C1:=GuessRec1(gu11,nops(C[2])/k):
           if C1<>FAIL then
              RETURN(lu,C1):
           fi:
         fi:
     fi:
   fi:

  od:
FAIL:
end:




#GuessPR1(C1,C2,d,x,y): Inputs two C-finite sequences and a pos. integer
#d and outputs a polynomial P of degree d such that P(C1[n],C2[n])=0
#for all n. If none found it returns FAIL. For example to solve
#the Hugh Montgomery proposed (but not selected) Putnam problem,
#(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,
#where F(n) is the Fibonacci sequence type
#GuessPR1([[0,1],[1,1]],[[1,1],[1,1]],4,x,y);
GuessPR1:=proc(C1,C2,d,x,y) local eq,var,gu,a,S1,S2,N,var1,pol,n1,mu,v1:
 gu:=GenPol([x,y],a,d):
 var:=gu[2]:
 gu:=gu[1]:
 N:=nops(var):
 S1:=SeqFromRec(C1,N+10):
 S2:=SeqFromRec(C2,N+10):
 eq:={seq(subs({x=S1[n1],y=S2[n1]},gu),n1=1..N+10)}:
 var1:=solve(eq,var):
 pol:=subs(var1,gu):
  if pol=0 then
    RETURN(FAIL):
  fi:

 S1:=SeqFromRec(C1,N+30):
 S2:=SeqFromRec(C2,N+30):

if {seq(normal(subs({x=S1[n1],y=S2[n1]},pol)),n1=N+11..N+30)}<>{0} then
print(`The conjecture`, pol , `did not materialize`):
RETURN(FAIL):
fi:

mu:={}:

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

if nops(mu)>1 then
 print(`There are than one solution`):
fi:

factor(coeff(pol,mu[1],1)):

end:



#GuessPR(C1,C2,d,x,y): Inputs two C-finite sequences and a pos. integer
#d and outputs a polynomial P of degree <=d such that P(C1[n],C2[n])=0
#for all n. If none found it returns FAIL. For example to solve
#the Hugh Montgomery proposed (but not selected) Putnam problem,
#(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,
#where F(n) is the Fibonacci sequence type
#GuessPR([[0,1],[1,1]],[[1,1],[1,1]],6,x,y);
GuessPR:=proc(C1,C2,d,x,y) local gu,d1:
for d1 from 1 to d do
 gu:=GuessPR1(C1,C2,d1,x,y):
 if gu<>FAIL then
  RETURN(gu):
 fi:
od:
FAIL:
end:



#GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a,
#and a non-negative integer deg, outputs a generic polynomial
#of total degree deg,in the variables Resh, indexed by the letter a of
#followed by the set of coeffs. For example, try:
#GenPol([x,y],a,1);
GenPol:=proc(Resh,a,deg)
local var,POL1,i,x,Resh1,POL:

if not type(Resh,list) or nops(Resh)=0 then
 ERROR(`Bad Input`):
fi:

x:=Resh[1]:
if nops(Resh)=1 then
RETURN(convert([seq(a[i]*x^i,i=0..deg)],`+`),{seq(a[i],i=0..deg)}):
fi:

Resh1:=[op(2..nops(Resh),Resh)]:
var:={}:
POL:=0:

for i from 0 to deg do
POL1:=GenPol(Resh1,a[i],deg-i):
var:=var union POL1[2]:
POL:=expand(POL+POL1[1]*x^i):
od:
POL,var:
end:


#GuessNLR1(C1,r,d,x): Inputs a C-finite sequence C and a pos. integer
#d and outputs a polynomial P(x[1],..., x[r]) of degree d such that 
#P(C1[n],C1[n-1],C1[n-r])=0
#for all n. If none found it returns FAIL. For example to solve
#the Hugh Montgomery proposed (but not selected) Putnam problem,
#(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,
#where F(n) is the Fibonacci sequence type
#GuessNLR1([[0,1],[1,1]],1,4,x);
GuessNLR1:=proc(C1,r,d,x) local eq,var,gu,a,S1,N,var1,pol,n1,mu,v1,i,i1:
 gu:=GenPol([seq(x[i],i=0..r)],a,d):
 var:=gu[2]:
 gu:=gu[1]:
 N:=nops(var)+r:
 S1:=SeqFromRec(C1,N+10):
 eq:={seq(subs({seq(x[i1]=S1[n1+i1],i1=0..r)},gu),n1=1..N+10-r)}:
 var1:=solve(eq,var):
 pol:=subs(var1,gu):
  if pol=0 then
    RETURN(FAIL):
  fi:

 S1:=SeqFromRec(C1,N+30):

if 
{seq(
normal(subs({seq(x[i1]=S1[n1+i1-1],i1=0..r)},pol)),n1=N+11..N+30-r-1)}<>{0} then
print(`The conjecture`, pol , `did not materialize`):
RETURN(FAIL):
fi:

mu:={}:

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

if nops(mu)>1 then
 print(`There are than one solution`):
fi:

factor(coeff(pol,mu[1],1)):

end:



#GuessNLR(C1,r,d,x): Inputs a C-finite sequence C and a pos. integer
#d and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that 
#P(C1[n],C1[n-1],C1[n-r])=0
#for all n. If none found it returns FAIL. For example to solve
#the Hugh Montgomery proposed (but not selected) Putnam problem,
#(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,
#where F(n) is the Fibonacci sequence type
#GuessNLR([[0,1],[1,1]],1,4,x);
GuessNLR:=proc(C1,r,d,x) local gu,d1:

for d1 from 1 to d do
gu:=GuessNLR1(C1,r,d1,x):
if gu<>FAIL then
 RETURN(gu):
fi:
od:
FAIL:
end:




#GuessNLRv(C1,r,d): Verbose form of GuessNLR(C1,r,d,x), i.e.
#Inputs a C-finite sequence C and a pos. integer
#d and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that 
#P(C1[n],C1[n-1],C1[n-r])=0
#for all n. If none found it returns FAIL. For example to solve
#the Hugh Montgomery proposed (but not selected) Putnam problem,
#(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,
#where F(n) is the Fibonacci sequence type
#GuessNLRv([[0,1],[1,1]],1,4,x);
GuessNLRv:=proc(C1,r,d) local gu,x,F,n,i1:
gu:=GuessNLR(C1,r,d,x):

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

print(`Theorem: Let F(n) be the sequence defined by the recurrence`):
print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))):
print(`subject to the initial conditions`):
print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))):
print(`Then we have the following identity valid for every non-negative integer n`):
print(subs({seq(x[i1]=F(n+i1),i1=0..r)},gu)=0):
print(``):
print(`Proof: Routine! (since everything in sight is C-finite).`):
end:



#FactorizeI1v(C1,k): Verbose form of FactorizeI1(C1,k), i.e.
FactorizeI1v:=proc(C1,k) local gu,C2,C3,i1,n:
gu:=FactorizeI1(C1,k): 

C2:=gu[1]:
C3:=gu[2]:

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

print(`Theorem: Let F(n) be the sequence defined by the recurrence`):
print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))):
print(`subject to the initial conditions`):
print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))):

print(``):
print(`Let G(n) be the sequence defined on non-negative integers by the recurrence`):
print(G(n)=add(C2[2][i1]*G(n-i1),i1=1..nops(C2[2]))):
print(`subject to the initial conditions`):
print(seq(G(i1-1)=C2[1][i1],i1=1..nops(C2[1]))):

print(``):
print(`Let H(n) be the sequence defined by the recurrence`):
print(H(n)=add(C3[2][i1]*H(n-i1),i1=1..nops(C3[2]))):
print(`subject to the initial conditions`):
print(seq(H(i1-1)=C3[1][i1],i1=1..nops(C3[1]))):

print(`Then the following is true for every non-negative integer n`):
print(F(n)=G(n)*H(n)):
print(``):
print(`Proof: Routine! (since everything in sight is C-finite) .`):
end:


#mSect(C1,m,i): Given a C-finite sequence C1(n), finds
#the C-finite description of C(m*n+i) where m is a pos. integer
#and 0<=i<m. For example, try:
#mSect([[0,1],[1,1]],2,0);
mSect:=proc(C1,m,i) local S,n1:
S:=SeqFromRec(C1,2*nops(C1[1])*(m+1)+3*m+10):
S:=[seq(S[i+1+m*n1],n1=0..2*nops(C1[1])+3)]:
GuessRec(S):
end:



#mSectV(C1,m,i): Verbose form of mSect(C1,m,i)
#For example, try:
#mSectV([[0,1],[1,1]],2,0);
mSectV:=proc(C1,m,i) local C2,F,G,n,i1:
C2:=mSect(C1,m,i):

print(`Theorem: Let F(n) be the sequence defined by the recurrence`):
print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))):
print(`subject to the initial conditions`):
print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))):
print(``):
print(`Define a new sequence`):
print(G(n)=F(i+m*n)):
print(`then G(n) satisfies the recurrence`):
print(G(n)=add(C2[2][i1]*G(n-i1),i1=1..nops(C2[2]))):
print(`subject to the initial conditions`):
print(seq(G(i1-1)=C2[1][i1],i1=1..nops(C2[1]))):
print(``):
end:



#BT(C1): The Binomial Transform of the C-finite sequence C,
#in other words, the sequence add(binomial(n,i)*C1[i],i=0..n):
#For example, try:
#BT([[0,1],[1,1]]);
BT:=proc(C1) local S,i1,n1:
S:=SeqFromRec(C1,nops(C1[1])*3+5):
S:=[seq(add(S[i1+1]*binomial(n1,i1),i1=0..n1),n1=0..nops(S)-1)]:
GuessRec(S):
end:


#BTv(C1): Verbose form of mSect(C1,m,i)
#For example, try:
#BTv([[0,1],[1,1]]);
BTv:=proc(C1) local C2,F,G,n,i1,i:
C2:=BT(C1):

print(`Theorem: Let F(n) be the sequence defined by the recurrence`):
print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))):
print(`subject to the initial conditions`):
print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))):
print(``):
print(`Let G(n) be the sequence that is its Binomial Transform `):
print(`In other words, the sequence defined by`):
print(G(n)=Sum(binomial(n,i)*F(i),i=0..n)):
print(`then G(n) satisfies the recurrence`):
print(G(n)=add(C2[2][i1]*G(n-i1),i1=1..nops(C2[2]))):
print(`subject to the initial conditions`):
print(seq(G(i1-1)=C2[1][i1],i1=1..nops(C2[1]))):
print(``):
end:

#Ux(x): The C-finite description of the Chebychev polynomials of
#the second kind. Try: Ux(x);
Ux:=proc(x):[[1,2*x],[2*x,-1]]:end:

#Tx(x): The C-finite description of the Chebychev polynomials of
#the first kind
Tx:=proc(x):[[1,x],[2*x,-1]]:end:

#Fn(): The C-finite description of the Fibonacci sequence
Fn:=proc():[[0,1],[1,1]]:end:
#Fn(): The C-finite description of fibonacci(n+1)
fn:=proc():[[1,1],[1,1]]:end:
#Ln(): The C-finite description of the Lucas Numbers
Ln:=proc():[[2,1],[1,1]]:end:
#Pn(): The C-finite description of the Pell Numbers
Pn:=proc():[[0,1],[2,1]]:end:

#Trn(): The C-finite description of the Triboanacci numbers Numbers
Trn:=proc():[[1,1,2],[1,1,1]]:end:




#KefelV(C1,C2): Verbose version of Kefel(S1,S2);
#For example try:
#KefelV([[1,2],[3,4]],[[2,5],[-1,6]]);
KefelV:=proc(C1,C2) local F,G,H,n,t,C,i1,S1,S2,S,m:
C:=Kefel(C1,C2):
print(`A proof in the Spirit of Zeilberger of A C-finite Identity`):

print(`Theorem: Let F(n) be the sequence defined by the recurrence`):
print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))):
print(`subject to the initial conditions`):
print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))):
print(``):
print(`Equivalently, the ordinary generating function of F(n) w.r.t. to t is`):
print(Sum(F(n)*t^n,n=0..infinity)=CtoR(C1,t)):
print(`Also Let G(n) be the sequence defined by the recurrence`):
print(G(n)=add(C2[2][i1]*G(n-i1),i1=1..nops(C2[2]))):
print(`subject to the initial conditions`):
print(seq(G(i1-1)=C2[1][i1],i1=1..nops(C2[1]))):
print(`Equivalently, the ordinary generating function of G(n) w.r.t. to t is`):
print(Sum(G(n)*t^n,n=0..infinity)=CtoR(C2,t)):
print(`Let H(n)=F(n)G(n), then`):
print(Sum(H(n)*t^n,n=0..infinity)=CtoR(C,t)):
print(`or equivalently, H(n) is the sequence defined by the recurrence`):
print(H(n)=add(C[2][i1]*H(n-i1),i1=1..nops(C[2]))):
print(`subject to the initial conditions`):
print(seq(H(i1-1)=C[1][i1],i1=1..nops(C[1]))):

print(``):
m:=nops(C1[1])*nops(C2[1])+nops(C[1]):
print(`Proof: F(n) is a C-finite sequence of order`, nops(C1[2])):
print(`and G(n) is a C-finite sequence of order`, nops(C2[2])):
print(`Hence their product is a C-finite sequence of order`):
print( nops(C1[1])*nops(C2[1])):
print(`Since the claimed sequence for their product is a C-finite sequence`):
print(`of order`, nops(C[1]),` it suffices to check the equality for`):
print(`the first`, m, ` values `):

S1:=expand(SeqFromRec(C1,m)):
S2:=expand(SeqFromRec(C2,m)):
S:=expand(SeqFromRec(C,m)):

print(`The first`, m , `members of the sequence F(n) are`):
print(S1):

print(`The first`, m,  `members of the sequence G(n) are`):
print(S2):
print(`The first`, m, `members of the proposed   product are`):
print(S):

for i1 from 1 to m do
if expand(S1[i1]*S2[i1]-S[i1])<>0 then
 print(`Oops something is wrong `):
 RETURN(FAIL):
else
 print(S1[i1], `times `, S2[i1], `is indeed `, S[i1]):
fi:
od:
print(`QED!`):

end:




#CHseq(L,n,j,K): inputs  a list of pairs [C,e], where C is a C-finite
#sequence, and e is an affine-linear expression in the symbols
#n and j, of the form a*n+b*j+c, where a and a+b are non-negative
#integers. Outputs the first K terms of
#Sum(L[1][1](L[1][2])*L[2][1](L[2][2]),j=0..n-1)
#using up to K values. For example, try:
#CHseq([[Fn(),j],[Fn(),j],[Fn(),2*n-j],20);
CHseq:=proc(L,n,j,K) local i1,S,j1,n1,Max1:

Max1:=0:
for i1 from 1 to nops(L) do
 for n1 from 0 to K do
   for j1 from 0 to n1-1 do
   Max1:=max(Max1,subs({j=j1,n=n1},L[i1][2])):
    if subs({j=j1,n=n1},L[i1][2])<0 then
      RETURN(FAIL):
    fi:
   od:
 od:
od:

for i1 from 1 to nops(L) do
S[i1]:=SeqFromRec(L[i1][1],Max1+1):
od:

[
seq(
add(mul(S[i1][subs({j=j1,n=n1},L[i1][2])+1],i1=1..nops(L)),j1=0..n1-1),
n1=0..K)]:

end:


#CH(L,n,j): 
#Find a sequence of the style of Curtis Greene and
#Herb Wilf ("Closed form summation of C-finite sequences",
#Trans. AMS)
#inputs  a list of pairs [C,e], where C is a C-finite
#sequence, and e is an affine-linear expression in the symbols
#n and j, of the form a*n+b*j+c, where a and a+b are non-negative
#integers. Outputs the C-finite description of
#Sum(L[1][1](L[1][2])*L[2][1](L[2][2]),j=0..n-1)
#using up to K values. For example, try:
#CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]);
CH:=proc(L,n,j) local K, C:
C:=GuessRec(CHseq(L,n,j,20)):

if C<>FAIL then
 RETURN(C):
fi:
for K from 25 by 5 while C=FAIL do
C:=GuessRec(CHseq(L,n,j,K)):
od:

C:

end:






#####New stuff



#GuessCH1(C1,r,d,x,C2): Inputs a C-finite sequence C1 and a pos. integer
#d and outputs a polynomial P(x[1],..., x[r]) of degree d such that 
#P(C1[n],C1[n-1],C1[n-r])=C2[n]
#for all n. If none found it returns FAIL. For example to reproduce
#Eq. (12) of the Greene-Wilf article
# type
#GuessCH1([[0,1],[1,1]],1,3,x,CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j));
GuessCH1:=proc(C1,r,d,x,C2) local eq,var,gu,a,S1,N,var1,pol,n1,mu,v1,i,i1,S2,mu1:
 gu:=GenPol([seq(x[i],i=0..r)],a,d):
 var:=gu[2]:
 gu:=gu[1]:
 N:=nops(var)+r:
 S1:=SeqFromRec(C1,N+11):
 S2:=SeqFromRec(C2,N+11):
 eq:={seq(subs({seq(x[i1]=S1[n1+i1+1],i1=0..r)},gu)-S2[n1+1],n1=1..N+10-r)}:
 var1:=solve(eq,var):
 pol:=subs(var1,gu):
  if pol=0 then
    RETURN(FAIL):
  fi:


 S1:=SeqFromRec(C1,N+30):
 S2:=SeqFromRec(C2,N+30):

if 
{seq(
normal(subs({seq(x[i1]=S1[n1+i1+1],i1=0..r)},pol)-S2[n1+1]),n1=N+11..N+30-r-2)
}<>{0} then
RETURN(FAIL):
fi:

mu:={}:

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

if mu<>{} then
factor(subs({seq(mu1=0,mu1 in mu)},pol)):
else factor(pol):
fi:

end:




#FindCHg(C1,r,d, x,C2): Inputs a C-finite sequence C1 and a pos. integer
#d and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that 
#P(C1[n],C1[n-1],C1[n-r])=C2[n]
#for all n. If none found it returns FAIL. For example to reproduce
#Eq. (12) of the Greene-Wilf article.
#If there is none of degree <=d then it returns FAIL
# type
#FindCHg([[0,1],[1,1]],2,4,x,CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j));
FindCHg:=proc(C1,r,d,x,C2) local gu,d1:

for d1 from 1 to d do
gu:=GuessCH1(C1,r,d1,x,C2):
if gu<>FAIL then
 RETURN(gu):
fi:
od:
FAIL:
end:

###end new stuff



#FindCH(C1,r,d, x,L,n,j): Inputs a C-finite sequence C1 and pos. integers
#r and d, and a letter x, and a list L of affine-linear expressions
#in n,j, and
#and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that 
#P(C1[n],C1[n-1],C1[n-r])=C2[n]
#where C2(n) is the sum(C1[L[1]]*C1[L[2]*...C1[L[nops(L)],j=0..n-1):
#If none found it returns FAIL. For example to reproduce
#Eq. (12) of the Greene-Wilf article.
#If there is none of degree <=d then it returns FAIL
# type
#FindCH(Fn(),1,3,x,[j,j,2*n-j],n,j);
FindCH:=proc(C1,r,d,x,L,n,j) local C2,i1,L1:
L1:=[seq([C1,L[i1]],i1=1..nops(L))]:
FindCHg([[0,1],[1,1]],r,d,x,CH(L1,n,j)):
end:



#FindCHv(C1,r,d, L,n,j): verbose version of FindCH(C1,r,d,x, L,n,j): 
# type
#FindCHv(Fn(),1,3,[j,j,2*n-j],n,j);
FindCHv:=proc(C1,r,d,L,n,j) local gu,x,F,i1,t:
gu:=FindCH(C1,r,d,x,L,n,j):

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

print(`Theorem: Let F(n) be the sequence defined by the recurrence`):
print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))):
print(`subject to the initial conditions`):
print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))):
print(``):
print(`Equivalently, the ordinary generating function of F(n) w.r.t. to t is`):
print(Sum(F(n)*t^n,n=0..infinity)=CtoR(C1,t)):
print(`Also Let G(n) be the sequence defined by `):
print(G(n)=Sum(mul(F(L[i1]),i1=1..nops(L)),j=0..n-1)):
print(`then we have the following closed-form expression for G(n)`):
print(G(n)=subs({seq(x[i1]=F(n+i1),i1=0..r)},gu)):
print(``):
print(`Proof: Both sides are C-finite, so it is enough to check the`):
print(`first few terms, and this is left to the dear reader.`):

end:


#FindCHvTN(C1,r,d, L,n,j,mispar): verbose version of FindCH(C1,r,d,x, L,n,j): 
# type
#FindCHvTN(Fn(),1,3,[j,j,2*n-j],n,j,8);
FindCHvTN:=proc(C1,r,d,L,n,j,mispar) local gu,x,F,i1,t:
gu:=FindCH(C1,r,d,x,L,n,j):

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

print(`Theorem Number`, mispar):
print( `Let F(n) be the sequence defined by the recurrence`):
print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))):
print(`subject to the initial conditions`):
print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))):
print(``):
print(`Equivalently, the ordinary generating function of F(n) w.r.t. to t is`):
print(Sum(F(n)*t^n,n=0..infinity)=CtoR(C1,t)):
print(`Also Let G(n) be the sequence defined by `):
print(G(n)=Sum(mul(F(L[i1]),i1=1..nops(L)),j=0..n-1)):
print(`then we have the following closed-form expression for G(n)`):
print(G(n)=subs({seq(x[i1]=F(n+i1),i1=0..r)},gu)):
print(``):
print(`Proof: Both sides are C-finite, so it is enough to check the`):
print(`first few terms, and this is left to the dear reader.`):

end:




#Curt1(j,n,M): outputs 
#all {a1n+b1j} with 0<=a1<=M and 0<=a1+b1<=M
#Type: Curt1(j,n,4);
Curt1:=proc(j,n,M)  local a1,b1:
[seq(seq(a1*n+b1*j, b1=-a1..-1),a1=0..M),
seq(seq(a1*n+b1*j, b1=1..M-a1 ),a1=0..M)]:
end:

#IV1(m,K): all the lists of weakly-increasing positive integers of size m
#with largest element K. Do IV(3,5);
IV1:=proc(m,K) local gu,i,gu1,gu11:
option remember:
if m=1 then
 RETURN([seq([i],i=1..K)]):
fi:
gu:={}:
for i from 1 to K do
gu1:=IV1(m-1,i):
gu:=gu union {seq([op(gu11),K],gu11 in gu1)}:
od:
gu:
end:



#IV(m,K): all the lists of weakly-increasing positive integers of size m
#with largest element <=K. Do IV(3,5);
IV:=proc(m,K) local i:
{seq(op(IV1(m,i)),i=1..K)}:
end:


Khevre:=proc(m,j,n,M) local gu,mu,mu1,i1:
gu:=Curt1(j,n,M):
mu:=IV(m-1,nops(gu)):
[seq([j,seq(gu[mu1[i1]],i1=1..m-1)],mu1 in mu)]:
end:

#SeferGW(m,d,M): outputs a book of Fibonacci identities with
#right hand sides of degree at most d in F(n),F(n+1)
#for Greene-Wilf type sums where the summand is
#of the form F(j)*F(a1n+b1j)*... with m terms, for example, try:
#SeferGW(2,5,2);
SeferGW:=proc(m,d,M) local n,j,gu,mu,mu1,x,C1,co:
co:=0:
print(`A Book of Definite Summation Fibonacci Identities`):
print(`in the style of Curtis Greene and Herbert Wilf`):
print(``):
print(`By Shalosh B. Ekhad `):

C1:=Fn():
mu:=Khevre(m,j,n,M):
for mu1 in mu do
gu:=FindCH(C1,1,d,x,mu1,n,j):

 if gu<>FAIL then
co:=co+1:
print(`------------------------------------------------`):
   FindCHvTN(C1,1,d,mu1,n,j,co):
 fi:
od:

print(`------------------------------------------------`):
print(`I hope that you enjoyed, dear reader, these`, co, `beautiful and`):
print(`deep theorems. `):
end:



#GenHomPol(Resh,a,deg): Given a list of variables Resh, a symbol a,
#and a non-negative integer deg, outputs a generic 
#homogeneous polynomial
#of degree deg,in the variables Resh, indexed by the letter a of
#followed by the set of coeffs. For example, try:
#GenHomPol([x,y],a,1);
GenHomPol:=proc(Resh,a,deg)
local var,POL1,i,x,Resh1,POL:

if not type(Resh,list) or nops(Resh)=0 then
 ERROR(`Bad Input`):
fi:

x:=Resh[1]:
if nops(Resh)=1 then
RETURN(a[1]*Resh[1]^deg,{a[1]}):
fi:

Resh1:=[op(2..nops(Resh),Resh)]:
var:={}:
POL:=0:

for i from 0 to deg do
POL1:=GenHomPol(Resh1,a[i],deg-i):
var:=var union POL1[2]:
POL:=expand(POL+POL1[1]*x^i):
od:
POL,var:
end:


#GuessInNLR1(C1,r,d,x): Inputs a C-finite sequence C and a pos. integer
#d and outputs an
#interesting polynomial P(x[0],..., x[r]) of degree d such that 
#of the form x[0]*P1(x[1], ..., x[r])+P2(x[1], ...,x[r])+c
#where P1,P2 are homog. of degree d, d-1 resp. such that
#P(C1[n],C1[n-1],C1[n-r])=0
#for all n. If none found it returns FAIL. For example to solve
#the Hugh Montgomery proposed (but not selected) Putnam problem,
#(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,
#where F(n) is the Fibonacci sequence type
#GuessInNLR1([[0,1],[3,-1]],2,2,x);
GuessInNLR1:=proc(C1,r,d,x) local eq,var,
gu,a,b,S1,N,var1,pol,n1,mu,v1,i,i1,gu1,gu2,c,lu,mu1:

 gu1:=GenHomPol([seq(x[i],i=1..r)],a,d-1):
 gu2:=GenHomPol([seq(x[i],i=1..r)],b,d):
 var:=gu1[2] union gu2[2] union {c}:
 gu:=expand(gu1[1]*x[0]+gu2[1]+c):
 N:=nops(var)+r:
 S1:=SeqFromRec(C1,N+10):
 eq:={seq(subs({seq(x[i1]=S1[n1+i1],i1=0..r)},gu),n1=1..N+10-r)}:
 var1:=solve(eq,var):
 pol:=subs(var1,gu):
  if pol=0 then
    RETURN(FAIL):
  fi:

 S1:=SeqFromRec(C1,N+40):

if 
{seq(
normal(subs({seq(x[i1]=S1[n1+i1-1],i1=0..r)},pol)),n1=N+11..N+40-r-1)}<>{0} then
print(`The conjecture`, pol , `did not materialize`):
RETURN(FAIL):
fi:

mu:={}:

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

lu:={}:

for i from 1 to nops(mu) do

 
mu1:=factor(coeff(pol,mu[i],1)):

if not type(mu1,`*`) then
 lu:=lu union {mu1}:
fi:

od:

lu:

end:








#KefelSk(a,b): The recurrence operator
#satisfied by the product of any solution of
#the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1]
#and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1]
#where a and b are lists
#For example, try:
#KefelSk([a1,a2],[b1,b2]);
KefelSk:=proc(a,b) local d1,d2,S1,S2,i,x,y,L,C,var,eq,L1,L2,gu,j:

d1:=nops(a):
d2:=nops(b):

S1:=[[seq(x[-d1+i],i=0..d1-1)],a]:

S2:=[[seq(y[-d2+i],i=0..d2-1)],b]:

L1:=SeqFromRec(S1,d1*d2+1):
L2:=SeqFromRec(S2,d1*d2+1):

L:=expand([seq(L1[i]*L2[i],i=1..d1*d2+1)]):

gu:=expand(L[d1*d2+1]-add(C[i]*L[d1*d2+1-i],i=1..d1*d2)):

var:={seq(C[i],i=1..d1*d2)}:
eq:={seq(seq(coeff(coeff(gu,x[-i],1),y[-j],1),j=1..d2),i=1..d1)}:
var:=solve(eq,var):

[seq(subs(var,C[i]),i=1..d1*d2)]:


end:


#Agel(p,z): inputs a polynomial p in floating-point, rounds it. Try:
#Agel(1.00001+3.00001*z,z);
Agel:=proc(p,z) local i,gu,gu1:
gu:=add(round(coeff(p,z,i))*z^i,i=ldegree(p,z)..degree(p,z)):
gu1:=gu-p:
if max(seq(abs(coeff(gu1,z,i)),i=ldegree(p,z)..degree(p,z)))>10^(-6) then
 RETURN(FAIL):
else
 RETURN(gu):
fi:
end:




#DimerData(z): the list of length 10 whose i-th item is the C-finite representation (only the recurrence part)
#of the sequence of weight-enumerators of tilings of an i by m rectangle, with weight z^(NumberOfVerticalTiles).
#Try:
#DimerData(z);
DimerData:=proc(z)
[[0, 1], [z, 1], [0, 2*z^2+2, 0, -1], [z^2, 3*z^2+2, z^2, -1], [0, 3*z^4+8*z^2+
4, 0, -10*z^4-16*z^2-6, 0, 3*z^4+8*z^2+4, 0, -1], [z^3, 6*z^4+10*z^2+4, 5*z^5+5
*z^3, z^6-13*z^4-20*z^2-6, -5*z^5-5*z^3, 6*z^4+10*z^2+4, -z^3, -1], [0, 4*z^6+
20*z^4+24*z^2+8, 0, -52*z^8-192*z^6-256*z^4-144*z^2-28, 0, 64*z^10+416*z^8+892*
z^6+844*z^4+360*z^2+56, 0, -16*z^12-160*z^10-744*z^8-1408*z^6-1216*z^4-480*z^2-\
70, 0, 64*z^10+416*z^8+892*z^6+844*z^4+360*z^2+56, 0, -52*z^8-192*z^6-256*z^4-\
144*z^2-28, 0, 4*z^6+20*z^4+24*z^2+8, 0, -1], [z^4, 10*z^6+30*z^4+28*z^2+8, 15*
z^8+35*z^6+19*z^4, 7*z^10-78*z^8-300*z^6-354*z^4-168*z^2-28, z^12-75*z^10-230*z
^8-217*z^6-63*z^4, -15*z^12+148*z^10+822*z^8+1442*z^6+1146*z^4+420*z^2+56, 69*z
^12+232*z^10+303*z^8+182*z^6+43*z^4, -119*z^12-642*z^10-1673*z^8-2304*z^6-1644*
z^4-560*z^2-70, 69*z^12+232*z^10+303*z^8+182*z^6+43*z^4, -15*z^12+148*z^10+822*
z^8+1442*z^6+1146*z^4+420*z^2+56, z^12-75*z^10-230*z^8-217*z^6-63*z^4, 7*z^10-\
78*z^8-300*z^6-354*z^4-168*z^2-28, 15*z^8+35*z^6+19*z^4, 10*z^6+30*z^4+28*z^2+8
, z^4, -1], [0, 5*z^8+40*z^6+84*z^4+64*z^2+16, 0, -190*z^12-1200*z^10-3026*z^8-\
3872*z^6-2632*z^4-896*z^2-120, 0, 655*z^16+7400*z^14+30734*z^12+65392*z^10+
80039*z^8+58488*z^6+25116*z^4+5824*z^2+560, 0, -550*z^20-9600*z^18-75506*z^16-\
291824*z^14-635896*z^12-847104*z^10-715572*z^8-384192*z^6-126672*z^4-23296*z^2-\
1820, 0, 125*z^24+3000*z^22+41860*z^20+312560*z^18+1269923*z^16+3077832*z^14+
4746678*z^12+4822192*z^10+3260653*z^8+1448360*z^6+404404*z^4+64064*z^2+4368, 0,
-3275*z^24-63200*z^22-546320*z^20-2580000*z^18-7423128*z^16-13887552*z^14-\
17518258*z^12-15134672*z^10-8937678*z^8-3532000*z^6-888888*z^4-128128*z^2-8008,
0, 20775*z^24+326600*z^22+2165340*z^20+8191760*z^18+19832526*z^16+32436304*z^14
+36765756*z^12+29101024*z^10+15961703*z^8+5915064*z^6+1405404*z^4+192192*z^2+
11440, 0, -41650*z^24-558400*z^22-3359060*z^20-11855040*z^18-27215340*z^16-\
42684320*z^14-46777648*z^12-36011264*z^10-19292248*z^8-7003776*z^6-1633632*z^4-\
219648*z^2-12870, 0, 20775*z^24+326600*z^22+2165340*z^20+8191760*z^18+19832526*
z^16+32436304*z^14+36765756*z^12+29101024*z^10+15961703*z^8+5915064*z^6+1405404
*z^4+192192*z^2+11440, 0, -3275*z^24-63200*z^22-546320*z^20-2580000*z^18-\
7423128*z^16-13887552*z^14-17518258*z^12-15134672*z^10-8937678*z^8-3532000*z^6-\
888888*z^4-128128*z^2-8008, 0, 125*z^24+3000*z^22+41860*z^20+312560*z^18+
1269923*z^16+3077832*z^14+4746678*z^12+4822192*z^10+3260653*z^8+1448360*z^6+
404404*z^4+64064*z^2+4368, 0, -550*z^20-9600*z^18-75506*z^16-291824*z^14-635896
*z^12-847104*z^10-715572*z^8-384192*z^6-126672*z^4-23296*z^2-1820, 0, 655*z^16+
7400*z^14+30734*z^12+65392*z^10+80039*z^8+58488*z^6+25116*z^4+5824*z^2+560, 0,
-190*z^12-1200*z^10-3026*z^8-3872*z^6-2632*z^4-896*z^2-120, 0, 5*z^8+40*z^6+84*
z^4+64*z^2+16, 0, -1], [z^5, 15*z^8+70*z^6+112*z^4+72*z^2+16, 35*z^11+140*z^9+
171*z^7+65*z^5, 28*z^14-322*z^12-2265*z^10-5184*z^8-5768*z^6-3388*z^4-1008*z^2-\
120, 9*z^17-566*z^15-3215*z^13-6692*z^11-6587*z^9-3087*z^7-551*z^5, z^20-255*z^
18+1353*z^16+19012*z^14+68971*z^12+125890*z^10+133221*z^8+84938*z^6+32032*z^4+
6552*z^2+560, -35*z^21+1984*z^19+15158*z^17+46708*z^15+77312*z^13+74348*z^11+
41517*z^9+12474*z^7+1561*z^5, 411*z^22-3406*z^20-55574*z^18-291977*z^16-829849*
z^14-1442558*z^12-1611098*z^10-1174278*z^8-552608*z^6-160888*z^4-26208*z^2-1820
, -2164*z^23-15597*z^21-52968*z^19-105361*z^17-129494*z^15-99521*z^13-47920*z^
11-14785*z^9-3114*z^7-395*z^5, 5527*z^24+61547*z^22+386154*z^20+1608919*z^18+
4462986*z^16+8297046*z^14+10472463*z^12+9038642*z^10+5307687*z^8+2073470*z^6+
512512*z^4+72072*z^2+4368, -6907*z^25-62987*z^23-295050*z^21-909216*z^19-\
1898966*z^17-2657435*z^15-2455558*z^13-1463330*z^11-535745*z^9-108423*z^7-9163*
z^5, 4508*z^26-4591*z^24-378657*z^22-2676224*z^20-9843835*z^18-22728088*z^16-\
35204576*z^14-37654642*z^12-28066567*z^10-14480232*z^8-5043640*z^6-1125124*z^4-\
144144*z^2-8008, -1591*z^27+41809*z^25+525554*z^23+2582552*z^21+7213271*z^19+
12799308*z^17+15047773*z^15+11830501*z^13+6121965*z^11+1989750*z^9+365715*z^7+
28765*z^5, 300*z^28-25078*z^26-95157*z^24+783883*z^22+7101270*z^20+25736222*z^
18+55416381*z^16+78964798*z^14+77734646*z^12+53630348*z^10+25795557*z^8+8435826
*z^6+1777776*z^4+216216*z^2+11440, -28*z^29+6237*z^27-90066*z^25-1297784*z^23-\
6358401*z^21-17032260*z^19-28498360*z^17-31413624*z^15-23177509*z^13-11309780*z
^11-3490046*z^9-613764*z^7-46563*z^5, z^30-694*z^28+42087*z^26+192164*z^24-\
927407*z^22-9606494*z^20-34884094*z^18-73719263*z^16-102492918*z^14-98357116*z^
12-66229900*z^10-31153572*z^8-9984576*z^6-2066064*z^4-247104*z^2-12870, 28*z^29
-6237*z^27+90066*z^25+1297784*z^23+6358401*z^21+17032260*z^19+28498360*z^17+
31413624*z^15+23177509*z^13+11309780*z^11+3490046*z^9+613764*z^7+46563*z^5, 300
*z^28-25078*z^26-95157*z^24+783883*z^22+7101270*z^20+25736222*z^18+55416381*z^
16+78964798*z^14+77734646*z^12+53630348*z^10+25795557*z^8+8435826*z^6+1777776*z
^4+216216*z^2+11440, 1591*z^27-41809*z^25-525554*z^23-2582552*z^21-7213271*z^19
-12799308*z^17-15047773*z^15-11830501*z^13-6121965*z^11-1989750*z^9-365715*z^7-\
28765*z^5, 4508*z^26-4591*z^24-378657*z^22-2676224*z^20-9843835*z^18-22728088*z
^16-35204576*z^14-37654642*z^12-28066567*z^10-14480232*z^8-5043640*z^6-1125124*
z^4-144144*z^2-8008, 6907*z^25+62987*z^23+295050*z^21+909216*z^19+1898966*z^17+
2657435*z^15+2455558*z^13+1463330*z^11+535745*z^9+108423*z^7+9163*z^5, 5527*z^
24+61547*z^22+386154*z^20+1608919*z^18+4462986*z^16+8297046*z^14+10472463*z^12+
9038642*z^10+5307687*z^8+2073470*z^6+512512*z^4+72072*z^2+4368, 2164*z^23+15597
*z^21+52968*z^19+105361*z^17+129494*z^15+99521*z^13+47920*z^11+14785*z^9+3114*z
^7+395*z^5, 411*z^22-3406*z^20-55574*z^18-291977*z^16-829849*z^14-1442558*z^12-\
1611098*z^10-1174278*z^8-552608*z^6-160888*z^4-26208*z^2-1820, 35*z^21-1984*z^
19-15158*z^17-46708*z^15-77312*z^13-74348*z^11-41517*z^9-12474*z^7-1561*z^5, z^
20-255*z^18+1353*z^16+19012*z^14+68971*z^12+125890*z^10+133221*z^8+84938*z^6+
32032*z^4+6552*z^2+560, -9*z^17+566*z^15+3215*z^13+6692*z^11+6587*z^9+3087*z^7+
551*z^5, 28*z^14-322*z^12-2265*z^10-5184*z^8-5768*z^6-3388*z^4-1008*z^2-120, -\
35*z^11-140*z^9-171*z^7-65*z^5, 15*z^8+70*z^6+112*z^4+72*z^2+16, -z^5, -1]]:
end:

#KastF(n,z): The C-finite representation (using numeric-floating point)
#of the Kastelian solution to the Dimer problem for n by m rectangles
#where z->1, z'->z. Try:
#KastF(6,z):
KastF:=proc(n,z) local i,gu:

if not type(n/2,integer) then
 print(`n must be an even integer`):
fi:

gu:=[2*z*cos(1*Pi/(n+1)),1];

for i from 2 to n/2 do
 gu:=KefelSk(gu, [2*z*cos(i*Pi/(n+1)),1]);
od:

[seq(Agel(evalf(gu[i]),z),i=1..nops(gu))]:

end:


#KastF1(n,z0): Like KastF(n,z) but for numeric z0 (to make it faster)
#The C-finite representation (using numeric-floating point)
#of the Kastelian solution to the Dimer problem for n by m rectangles
#where z->1, z'->z0. z0 is a number.Try:
#KastF1(6,3):
KastF1:=proc(n,z0) local i,gu:


if not type(n/2,integer) then
 print(`n must be an even integer`):
fi:

gu:=evalf([2*z0*cos(1*Pi/(n+1)),1]);

for i from 2 to trunc(n/2) do
 gu:=KefelSk(gu, [2*z0*evalf(cos(i*Pi/(n+1))),1]);
od:

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

end:


#start stuff added July 22, 2022

Coe1:=proc(f, x)
    local i;
    {seq(expand(coeff(f, x, i)), i = ldegree(f, x) .. degree(f, x))};
end:

Coe:=proc(f, x)
    local gu,  x1, gu1;
    if nops(x) = 1 then RETURN(Coe1(f, x[1])); end if;
    x1 := [op(2 .. nops(x), x)];
    gu := Coe(f, x1);
    {seq(op(Coe1(gu1, x[1])), gu1 in gu)};
end:



#GuessNLRg1(INI,REC,d,x): Inputs a C-finite recurrence REC and a list of length nops(REC), INI and a positive integer d
#d and outputs a polynomial P(x[1],..., x[r]) of degree d such that 
#if x[1](n)=SeqFromRec([INI[1],REC]), ...,  x[r](n)=SeqFromRec([INI[r],REC]),  then
#P(x[1](n),.., x[r](n))=0
#for all n. If none found it returns FAIL. For example try:
#GuessNLRg1([[1,k],[0,b]],[2*k,-1],2,x);
#GuessNLRg1([[1,1,2], [0,1,1],[0,0,1]],[1,1,1],2,x);
GuessNLRg1:=proc(INI,REC,d,x) local eq,var,gu,a,S,N,var1,pol,n1,mu,v1,i,i1:


 gu:=GenPol([seq(x[i],i=1..nops(INI))],a,d):
 var:=gu[2]:
 gu:=gu[1]:
 N:=nops(var)+10:

for i from 1 to nops(INI) do
 S[i]:=SeqFromRec([INI[i],REC],N):
od:

 eq:={seq(subs({seq(x[i1]=S[i1][n1],i1=1..nops(INI))},gu),n1=1..N)}:
 var1:=solve(eq,var):
 pol:=subs(var1,gu):
  if pol=0 then
    RETURN(FAIL):
  fi:

mu:={}:

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

if nops(mu)>1 then
 print(`There are than one solution`):
fi:

factor(coeff(pol,mu[1],1)):

end:



#GuessNLRg(INI,REC,d,x): Inputs a C-finite recurrence REC and a list of length nops(REC), INI and a positive integer d
#d and outputs a polynomial P(x[1],..., x[r]) of degree <=d such that 
#if x[1](n)=SeqFromRec([INI[1],REC]), ...,  x[r](n)=SeqFromRec([INI[r],REC]),  then
#P(x[1](n),.., x[r](n))=0
#for all n. If none found it returns FAIL. For example try:
#GuessNLRg([[1,k],[0,b]],[2*k,-1],2,x);
#GuessNLRg([[1,1,2], [0,1,1],[0,0,1]],[1,1,1],2,x);
GuessNLRg:=proc(INI,REC,d,x) local d1,gu:

for d1 from 1 to d do

gu:=GuessNLRg1(INI,REC,d1,x):

if gu<>FAIL then
 RETURN(numer(gu)):
fi:
od:
FAIL:
end:


####New stuff




#GuessNLR1g(INI, REC,d,x): Inputs a a list of lists of length r (say) each list of lentgh k
#and a list REC of length r, say
#d and outputs a polynomial P(x[1],..., x[r]) of degree d such that 
#P(C1[n],C2[n],Cr[n])=0
#for all n. 
#where C1[n]=SeqFromRec([INI[1],REC]) etc.
#If none found it returns FAIL. For example to solve
#the Hugh Montgomery proposed (but not selected) Putnam problem,
#(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0,
#where F(n) is the Fibonacci sequence type
#GuessNLR1g([[0,1],[1,1]],[3,-1]],2,x);
GuessNLR1g:=proc(INI,REC,d,x) local eq,var,gu,a,S,N,var1,pol,n1,mu,v1,i,i1,r:
r:=nops(REC):
if not (nops(INI)=r and {seq(nops(INI[i1]),i1=1..nops(INI))}={r}) then
 RETURN(FAIL):
fi:

 gu:=GenPol([seq(x[i],i=1..r)],a,d):
 var:=gu[2]:
 gu:=gu[1]:
 N:=nops(var)+r:

S:=[]:

for i from 1 to nops(INI) do

 S[i]:=SeqFromRec([INI[i],REC],N+10):

od:





 eq:={seq(subs({seq(x[i1]=S[i1][n1],i1=1..r)},gu),n1=1..N+10-r)}:
 var1:=solve(eq,var):
 pol:=subs(var1,gu):
  if pol=0 then
    RETURN(FAIL):
  fi:


mu:={}:

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

if nops(mu)>1 then
 print(`There are than one solution`):
fi:

factor(coeff(pol,mu[1],1)):

end:



#GuessNLRg(INI,REC,d,x): 
GuessNLRg:=proc(INI,REC,d,x) local gu,d2:

for d2 from 1 to d do
gu:=GuessNLRg1(INI,REC,d2,x):

if gu<>FAIL then
 RETURN(gu):
fi:
od:
FAIL:
end:





###Start Diophantine equations

#SolG(a0,a1,b0,b1,t): The pair of rational functions [f1(t),f2(t)] such that its respective Maclaurin coefficients solve the generalized Pell equation. Try:
#SolPG(3,5,11,17,t);
SolPG:=proc(a0,a1,b0,b1,t)  local f1,f2,c,c0:
c0:=2*(a0*b0+a1*b1)/(a0*b1+a1*b0):

f1:=CtoR([[a0,a1],[c,-1]],t):
f2:=CtoR([[b0,b1],[c,-1]],t):
f1:=normal(subs(c=c0,f1)):
f2:=normal(subs(c=c0,f2)):
[f1,f2]:
end:


P:=CtoR([[a0,a1],[c0,-1]],x);

#PellG(a0,a1,b0,b1,x): The Pell-like Diophantine equation with parameters a0,a1,b0,b1
PellG:=proc(a0,a1,b0,b1,x,y)
local A0, B0, C0, g:

A0 := b0^2 - b1^2:
B0 := a0^2 - a1^2:
C0 := b0^2 * a1^2 - a0^2 * b1^2:
g := igcd(A0, B0, C0):
(A0 / g) * x^2 - (B0 / g) * y^2 - (C0 / g):
end:


#Theo2(a0,a1,b0,b1,x,y): outputs a theorem about the generalized Pell equation coming from C-finite sequenes of order 2. Try:
#Theore(5,6,6,7,x,y);
#Theo2(a0,a1,b0,b1,x,y):
Theo2:=proc(a0,a1,b0,b1,x,y)  local P,f1,f2,gu,t,a,b,n,i:

P:=PellG(-a0,a1,b0,b1,x,y):

gu:=SolPG(-a0,a1,b0,b1,t):

print(`Infinitely many solutions to the Diophantine Equation`, P=0):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem `):
print(``):
print(`Let `):
print(``):
print(Sum(a(n)*t^n,n=0..infinity)=gu[1]):
print(``):
print(Sum(b(n)*t^n,n=0..infinity)=gu[2]):
print(``):

f1:=taylor(gu[1],t=0,8):
f2:=taylor(gu[2],t=0,8):

print(`Here are first 8 pairs  starting at n=0`):
print(``):
lprint([seq([coeff(f1,t,i),coeff(f2,t,i)],i=0..7)]):
print(``):
print(`Then the infinite set of pairs, x=a[n],y=b[n] satisfy the Diophantine equation `):
print(``):
print(P=0):
print(``):
print(`Proof:  It suffices to check it for the n=1..7, and indeed plugging it in gives`):
f1:=taylor(gu[1],t=0,8):
f2:=taylor(gu[2],t=0,8):

print({seq(normal(subs({x=coeff(f1,t,i),y=coeff(f2,t,i)},P)),i=0..7)}):

end:







#GoodQ2(K): finds all quadruples [a0,a1,b0,b1] such that a0*b1-a1*b0=1 and 1<=a0,a1,b0<=K. Try:
#GoodQ2(K)
GoodQ2:=proc(K) local gu,a0,a1,b0,b1:

gu:={}:

for a0 from 1 to K do
for a1 from 1 to K do
for b0 from 1 to K do
 b1:=(1+a1*b0)/a0:
 if type(b1,integer) then
  gu:=gu union {[a0,a1,b0,b1]}:
 fi:
od:
od:
od:
gu:
end:


#GoodQ3(K): finds all quadruples [a0,a1,b0,b1] such that
#   c = 2(a0 b0 + a1 b1) / (a0 b1 + a1 b0)
# is an integer, b0^2 != b1^2, a0^2 != a1^2, and 1 <= a0, a1, b0 <= K.
GoodQ3:=proc(K) local gu,a0,a1,b0,b1, c0:

gu:={}:

for a0 from 0 to K do
for a1 from 1 to K do
    if a1 = a0 then
        next:
    fi:
for b0 from 0 to K do
for b1 from 1 to K do

    if b0 = b1 then
        next:
    fi:

    if a0 * b1 + a1 * b0 = 0 then
        next:
    fi:
 c0 := 2 * (a0 * b0 + a1 * b1) / (a0 * b1 + a1 * b0):
 if type(c0,integer) then
  gu:=gu union {[a0,a1,b0,b1]}:
 fi:
od:
od:
od:
od:
gu:
end:

PellC := proc(a0, a1, b0, b1)
    2*(a0*b0+a1*b1)/(a0*b1+a1*b0):
end:

findPells := proc(K, x, y)
    local P, a, b, c, g, goods, L:
    goods := {}:
    for L in GoodQ3(K) do
        P := PellG(op(L), x, y):
        a, b, c := coeff(P, x, 2), coeff(P, y, 2), coeff(coeff(P, x, 0), y, 0):
        # ax^2 + by^2 + c = 0

        if a = 1 and b < 0 then
            goods := goods union {[L, [PellC(op(L)), -1], P]}:
        fi:
    od:

    goods:
end: