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

Digits:=100:
 
print(`Created:  Dec. 21, 2015`):
print(``):
print(` This is CfiniteIntegral.txt, a Maple package that `):
print(` accompanyies the article `):
print(` "The C-finite Ansatz Meets the Holonomic Ansatz" `):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
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 procedures from the package Findrec , type ezraCStory();, for help with`):
 print(`a specific procedure, type ezraC(procedure_name);   .`):
 print(``):
print(`-------------------------------`):

print(`-----------------------`):
 print(`For a list of the Story procedures from the package Cfinie, type ezraCStory();, for help with`):
 print(`a specific procedure, type ezraC(procedure_name);   .`):
 print(``):
print(`-------------------------------`):

print(`-----------------------`):
 print(`For a list of the procedures coding famous sequences from the C-finite , type`):
 print(`ezraCFamous();`):
 print(` for help with`):
 print(`a specific procedure, type ezraC(procedure_name);   .`):
 print(``):
print(`-------------------------------`):

print(`-----------------------`):
 print(`For a list of the supporting procedures for the C-finite package, type ezraC1();, for help with`):
 print(`a specific procedure, type ezraC(procedure_name);   .`):
 print(``):
print(`-------------------------------`):

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



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



with(combinat):
with(numtheory):

###start stuff from the package Cfinite
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, CtoR, Curt1,`):
 print(` etop, GenHomPol, `):
 print(` GenPol,  FindCH1, FindCHg, GuessInNLR1, GuessNLR1, GuessPR1`):
 print(` GuessRec1, GuessRec1mini,IISF, ISF, Krovim, ptoe `):
 print(`  ProdIndicator, ProdIndicatorG, ProdProfile, `):
 print(` ProdProfileE, RtoC, TDB`):

else
ezra(args):
fi:

end:


ezraC:=proc()

if args=NULL then
 print(`The main procedures are: BT,  CH, DimerData, Factorize, FactorizeI1`):
 print(` FindCH,GuessNLR, GuessPR, GuessRec, IsProd, IsProdG `):
 print(` Kast, Kefel, KefelR, KefelS, KefelSk, KefelSm, Khibur, KhiburS`):
 print(` mSect, SeqFromRecC `):
 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])=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])=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])=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])=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])=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])=SeqFromRecC then
print(`SeqFromRecC(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: SeqFromRecC([[1,1],[1,1]],20);`):

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])=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:


###from EKHAD
gzor:=proc(f,x,n)
local i,gu:
gu:=f:
for i from 1 to n do
gu:=diff(gu,x):
od:
RETURN(gu):
end:

goremc:=proc(D,ORDER)
local i,gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*D^i:
od:
RETURN(gu):
end:

pashetc:=proc(p,D)
local gu,p1,i,mekh:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,D) do
gu:=gu+factor(coeff(p1,D,i))*D^i:
od:
RETURN(gu,mekh):
end:


duisc:= proc(A,ORDER,y,x,D)
local KAMA,gorem,conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,i,shad:
 
KAMA:=10:
gorem:=goremc(D,ORDER):
 
conj:=gorem:
yakhas:=0:
 
for i from 0 to degree(conj,D) do
yakhas:=yakhas+normal(coeff(conj,D,i)*gzor(A,x,i)/A):
yakhas:=normal(yakhas):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=SolveTools[Linear](eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
shad:=pashetc(gorem,D):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
shad[1],S:
end:

AZc:=proc(A,y,x,D)
local ORDER,gu,KAMA:
KAMA:=12:
 
for ORDER from 1 to KAMA do
gu:=duisc(A,ORDER,y,x,D):
if gu<>0 and gu[1]<>0 then
 RETURN(gu):
fi:
od:
0:
end:

#AZcI(A,y,x,D,a,b): An Inhomog. differential equation, in x, and D=d/dx, for Int(A(x,y),y=a..b);
#It returns the pair [DiffEq,RightSide]. Try:
#AZcI(1/(1-x*t-x^2*t),x,t,Dt,0,1);
AZcI:=proc(A,y,x,D,a,b)
local ORDER,gu,KAMA,ope,rhs1,cert:
KAMA:=12:
 
for ORDER from 1 to KAMA do
gu:=duisc(A,ORDER,y,x,D):
if gu<>0 and gu[1]<>0 then
ope:=gu[1]:
cert:=A*gu[2]:

rhs1:=normal(normal(subs(y=b,cert))- normal(subs(y=a,cert))):
RETURN([ope,rhs1]):
fi:
od:
FAIL:
end:
 

###end from EKHAD

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:

#SeqFromRecC(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: SeqFromRecC([[1,1],[1,1]],20);
SeqFromRecC:=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:=SeqFromRecC(S1,K):
L2:=SeqFromRecC(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:=SeqFromRecC(S1,K):
L2:=SeqFromRecC(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:=SeqFromRecC(S1,d1+d2+1):
L2:=SeqFromRecC(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:=SeqFromRecC(S1,d1*d2+1):
L2:=SeqFromRecC(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(SeqFromRecC(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRecC(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:=SeqFromRecC(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:=SeqFromRecC(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:=SeqFromRecC(C,10*k):
          gu1:=SeqFromRecC(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:=SeqFromRecC(C1,N+10):
 S2:=SeqFromRecC(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:=SeqFromRecC(C1,N+30):
 S2:=SeqFromRecC(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:=SeqFromRecC(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:=SeqFromRecC(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:=SeqFromRecC(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:=SeqFromRecC(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(SeqFromRecC(C1,m)):
S2:=expand(SeqFromRecC(C2,m)):
S:=expand(SeqFromRecC(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]:=SeqFromRecC(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:
 gu:=GenPol([seq(x[i],i=0..r)],a,d):
 var:=gu[2]:
 gu:=gu[1]:
 N:=nops(var)+r:
 S1:=SeqFromRecC(C1,N+11):
 S2:=SeqFromRecC(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:=SeqFromRecC(C1,N+30):
 S2:=SeqFromRecC(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:=SeqFromRecC(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:=SeqFromRecC(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:=SeqFromRecC(S1,d1*d2+1):
L2:=SeqFromRecC(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:


###End stuff from the package Cfinite
##start stuff from the package FindRec
ezraF:=proc()
if args=NULL then

print(` FindRec: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF`):
print():

elif nargs=1 and args[1]=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 nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRecC then
print(`SeqFromRecC(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRecC(N-n-1,n,N,[1],10);`):


fi:
 
end:


###Findrec 
#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);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+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)+2 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
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

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:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L,i1:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:

 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):


     if ope<>FAIL then
      if normal({seq(add(subs(n=n1,coeff(ope,N,i1))*f[n1+i1],i1=0..degree(ope,N)),n1=1..nops(f)-degree(ope,N))})={0} then
        RETURN(ope):
      else
       RETURN(FAIL):
      fi:
    fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(H,n,N,K): Inputs H=[Ini,ope]
#Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(H,n,N,K)
local ope1,gu,L,n1,j1,ope,Ini:

Ini:=H[1]:
ope:=H[2]:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
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,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
##end stuff from the package FindRec


ezra:=proc()

if args=NULL then
 print(`The procedures are: ApplyU, IntC, IntCaz, IntCazV, IntCS, IntCSaz, IntCSv, IntCv, SeqFromRec, SHK1a, SHK1b, Star `):

elif nops([args])=1 and op(1,[args])=ApplyU then
print(`ApplyU(P,x,n,U): inputs a polynomial P in the variable x, a symbol n, and an Umbra U, an expression in n`):
print(`applies the Umbra defined by x^n->U(n)`):
print(` Try: `):
print(` ApplyU(x+x^5,x,n,1/(n+1)); `):


elif nops([args])=1 and op(1,[args])=AZcI  then
print(`AZcI(A,y,x,D,a,b): An Inhomog. differential equation, in x, and D=d/dx, for Int(A(x,y),y=a..b);`):
print(`It returns the pair [DiffEq,RightSide]. Try:`):
print(`AZcI(1/(1-x*t-x^2*t),x,t,Dt,0,1);`):

elif nops([args])=1 and op(1,[args])=IntC  then
print(`IntC(C,x,p,n,N,U,MaxC): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)`):
print(`in the variable x, a positive integer p, a symbol n, and an Umbra U, and a positive integer MaxC,`):
print(` guesses a (rigorously-provable) linear recurrence operator with polynomial coefficients, using the `):
print(` shift operator N, `):
print(` for the sequence U(C_n(x)^p). It returns the pair [INI,ope]. Try: `):
print(` IntC([[1,1],[1,x]],x,1,n,N,1/(n+1),10); `):
print(`For the orthogonality of the Chebychev polynomials of the First kind, try:`):
print(` IntC([[1,x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n,10); `):
print(`For the orthogonality of the Chebychev polynomials of the Second kind, try:`):
print(` IntC([[1,2*x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n/(n+2),10); `):
print(`Here are some random examples: `):
print(` IntC([[1,2*x,x],[1,x,x^2]],x,1,n,N, 1/(n+1),10); `):
print(` IntC([[1,2*x,x,x^2],[1,x,x^2,x]],x,1,n,N, 1/(n+1)/(n+2),18); `):


elif nops([args])=1 and op(1,[args])=IntCaz  then
print(`IntCaz(C,x,K,p,t,Dt,a,b): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)`):
print(`in the variable x, a nice function K of x, a positive integer p, a symbol t, a symbol Dt for diffferentiatin with respect to t,`):
print(`and numbers a and b`):
print(`Outputs a linear differential equation ope(t,Dt)F(t)=Rhs(t)`):
print(`for the generating function of the sequence  int(K(x)*C_n(x)^p,x=a..b):`):
print(`It returns the pair [ope,Rhs]`):
print(` Try: `):
print(` IntCaz([[1,1],[1,x]],x,1,1,t,Dt,0,1); `):
print(` IntCaz([[1,1],[1,x]],x,1,2,t,Dt,0,2); `):
print(` IntCaz([[1,1,1],[1,x,x^2]],x,1,1,t,Dt,0,1); `):

elif nops([args])=1 and op(1,[args])=IntCazV  then
print(`IntCazV(C,x,K,p,t,a,b):  Verbose version of IntCaz(C,x,K,p,t,Dt,a,b) (q.v)`):
print(` Try: `):
print(` IntCazV([[1,1],[1,x]],x,1,1,t,0,1); `):
print(` IntCazV([[1,1],[1,x]],x,1,2,t,0,2); `):
print(` IntCazV([[1,1,1],[1,x,x^2]],x,1,1,t,0,1); `):



elif nops([args])=1 and op(1,[args])=IntCS  then
print(`IntCS(C,x,r,n,N,U,MaxC): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)`):
print(`in the variable x, a positive integer p, a symbol n, and an Umbra U, and a positive integer MaxC,`):
print(` guesses a (rigorously-provable) linear recurrence operator with polynomial coefficients, using the `):
print(` shift operator N`):
print(` for the sequence U(C_n(x)*Star(C_n(x))), where Star(P(x))=P(1/x)*x^(degree(P)) . It returns the pair [INI,ope]. Try: `):
print(` It returns the pair [INI,ope]`):
print(``):
print(`For example to test Theorem 1 of Seon-Hong Kim's article:`):
print(` On some integrals involving Chebyshev polynomials", Ramanujan J. v. 38 (2015),`):
print(`(629-639) `):
print(` First part: i.e. Int(U_(2n)(x)*U*x_(2n)(x),x=-1..1).. Try: `):
print(` H1:=IntCS([[1,2*x],[2*x,-1]],x,2,n,N,(1+(-1)^n)/(n+1),10); `):
print(` evalb(SeqFromRec(H1,n,N,50)=SHK1a(50));`):
print(` Second part: i.e. Int(T_(2n)(x)*Star(T_(2n))(x),x=-1..1).. Try: `):
print(` H2:=IntCS([[1,x],[2*x,-1]],x,2,n,N,(1+(-1)^n)/(n+1),10); `):
print(` evalb(SeqFromRec(H2,n,N,50)=SHK1b(50));`):
print(``):
print(`For Theorem 2 of that paper, First part: `):
print(` IntCS([[1,x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n,10); `):
print(`Second part: `):
print(` IntCS([[1,2*x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n/(n+2),10); `):

elif nops([args])=1 and op(1,[args])=IntCSaz  then
print(`IntCSaz(C,x,r,t,Dt,a,b): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)`):
print(`in the variable x, a positive integer r,  symbols t and Dt (Dt denotes differentiation with repect to t), and numbers a and b such that a<b,`):
print(`finds an inhomogeneous linear  differential equation of the form ope(t,Dt)F(t)=Rhs(t) for`):
print(`for the generating function in t, for  Int(C_(rn)(x)Star(C_rn(x),x=a..b).`):
print(`It returns the pair [ope,Rhs], where INI are the initial conditions and ope is the operator`):
print(` Try: `):
print(` IntCSaz([[1,1],[1,x]],x,2,t,Dt,0,1); `):


elif nops([args])=1 and op(1,[args])=IntCSv  then
print(`IntCSv(C,x,r,n,N,U,MaxC,GADOL): Verbose version of IntCS(C,x,r,n,N,U,MaxC) (q.v) `):
print(`For example to test Theorem 1 of Seon-Hong Kim's article:`):
print(` On some integrals involving Chebyshev polynomials", Ramanujan J. v. 38 (2015),`):
print(`(629-639) `):
print(` First part: i.e. Int(U_(2n)(x)*U*x_(2n)(x),x=-1..1).. Try: `):
print(` IntCSv([[1,2*x],[2*x,-1]],x,2,n,N,(1+(-1)^n)/(n+1),10,1000); `):
print(` Second part: i.e. Int(T_(2n)(x)*Star(T_(2n))(x),x=-1..1).. Try: `):
print(` IntCSv([[1,x],[2*x,-1]],x,2,n,N,(1+(-1)^n)/(n+1),10,1000); `):
print(``):
print(`For Theorem 2 of that paper, First part: `):
print(` IntCSv([[1,x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n,10,1000); `):
print(`Second part: `):
print(` IntCSv([[1,2*x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n/(n+2),10,1000); `):

elif nops([args])=1 and op(1,[args])=IntCv  then
print(`IntCv(C,x,p,n,N,U,MaxC,GADOL): Verbose version of IntC(C,x,p,n,N,U,MaxC) (q.v.) `):
print(`The inputs are the same as for IntC, except that there is an additional input, GADOL, a large positive integer`):
print(`that the article computes the GADOL's term of the sequence`):
print(`Try: `):
print(` IntCv([[1,1],[1,x]],x,1,n,N,1/(n+1),10,500); `):
print(`For the orthogonality of the Chebychev polynomials of the First kind, try:`):
print(` IntCv([[1,x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n,10,500); `):
print(`For the orthogonality of the Chebychev polynomials of the Second kind, try:`):
print(` IntCv([[1,2*x],[2*x,-1]],x,2,n,N, (1+(-1)^n)/2*binomial(n,n/2)/2^n/(n+2),10,500); `):



elif nops([args])=1 and op(1,[args])=SeqFromRec  then
print(`SeqFromRec(H,n,N,K): Inputs H=[Ini,ope]`):
print(`Given the first L`):
print(`terms of the sequence Ini=[f(1), ..., f(L)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):

elif nops([args])=1 and op(1,[args])=Star  then
print(`Star(P,x): P(1/x)*x^n, when n is the degree of x, Try`):
print(`Star(1+3*x+2*x^2,x);`):

elif nops([args])=1 and op(1,[args])=SHK1a  then
print(`SHK1a(N): The first N terms of the left hand-side of Theorem 1, first part`):
print(`Seon-Hong Kim's paper (Rmanujan j. 38(2015), 629-639, Eq. (1.3)`):
print(`Try:`):
print(`SHK1a(10);`):

elif nops([args])=1 and op(1,[args])=SHK1b  then
print(`SHK1a(b): The first N terms of the left hand-side of Theorem 1, second part`):
print(`Seon-Hong Kim's paper (Rmanujan j. 38(2015), 629-639, Eq. (1.4)`):
print(`SHK1b(10);`):

else

print(`There is no help for`, args):

fi:

end:


#ApplyU(P,x,n,U): inputs a polynomial P in the variable x, a symbol n, and an Umbra U, an expression in n
#applies the Umbra defined by x^n->U(n)
#Try:
#ApplyU(x+x^5,x,n,1/(n+1));
ApplyU:=proc(P,x,n,U) local i1:
add(coeff(P,x,i1)*subs(n=i1,U),i1=0..degree(P,x)):
end:


#Star(P,x): P(1/x)*x^n, when n is the degree of x, Try
#Star(1+3*x+2*x^2,x);
Star:=proc(P,x) local n:
n:=degree(P,x):
expand(subs(x=1/x,P)*x^n):
end:


#IntC(C,x,p,n,N,U,MaxC): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)
#in the variable x, a positive integer p, a symbol n, and an Umbra U, and a positive integer MaxC,
#guesses a (rigorously-provable) linear recurrence operator with polynomial coefficients, using the
#shift operator N,
#for the sequence U(C_n(x)^p). 
#It returns the pair [INI,ope]
#Try:
#IntC([[1,1],[1,x]],x,1,n,N,1/(n+1),10);

IntC:=proc(C,x,p,n,N,U,MaxC) local gu,ORDER,L,K,i,ope,INI:
K:=max(seq(seq((2+L-ORDER)*(1+ORDER),ORDER=0..L),L=0..MaxC))+11:

gu:=SeqFromRecC(C,K):
gu:=[op(2..nops(gu),gu)]:


gu:=eval([seq(ApplyU(gu[i]^p,x,n,U),i=1..nops(gu))]):

if {seq(gu[2*i],i=1..trunc(nops(gu)/2))}={0} then
   print(`The sequence of values is 0 in even indices we will take let n be (n-1)/2 `):
  gu:=[seq(gu[2*i+1],i=0..nops(gu)/2)]:
fi:

if {seq(gu[2*i+1],i=0..trunc(nops(gu)/2)-1)}={0} then
   print(`The sequence of values is 0 in odd indices we will take let n be n/2 `):
  gu:=[seq(gu[2*i],i=1..trunc(nops(gu)/2))]:
fi:



ope:=Findrec(gu,n,N,MaxC):


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

INI:=eval([op(1..degree(ope,N),gu)]):
[INI,collect(ope,N)]:
end:


#IntCS(C,x,r,n,N,U,MaxC): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)
#in the variable x, a positive integer r, a symbol n, and an Umbra U, and a positive integer MaxC,
#guesses a (rigorously-provable) linear recurrence operator with polynomial coefficients, using the
#shift operator N,
#for the sequence U(C_(rn)(x)Star(C_rn(x),x)). 
#It returns the pair [INI,ope], where INI are the initial conditions and ope is the operator
#Try:
#IntCS([[1,1],[1,x]],x,2,n,N,1/(n+1),10);

IntCS:=proc(C,x,r,n,N,U,MaxC) local gu,ORDER,L,K,i,ope,INI:
K:=max(seq(seq((2+L-ORDER)*(1+ORDER),ORDER=0..L),L=0..MaxC))+11:

gu:=SeqFromRecC(C,r*K+1):
gu:=[op(2..nops(gu),gu)]:
gu:=[seq(gu[r*i],i=1..K)]:
gu:=[seq(expand(gu[i]*Star(gu[i],x)),i=1..nops(gu))]:
gu:=eval([seq(ApplyU(gu[i],x,n,U),i=1..nops(gu))]):


ope:=Findrec(gu,n,N,MaxC):


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

INI:=eval([op(1..degree(ope,N),gu)]):
[INI,collect(ope,N)]:
end:


#IntCv(C,x,p,n,N,U,MaxC,GADOL):  Verbose version of IntC(C,x,p,n,N,U,MaxC) (q.v.)
#Try
#IntCv([[1,1],[1,x]],x,1,n,N,1/(n+1),10);

IntCv:=proc(C,x,p,n,N,U,MaxC,GADOL) local gu,lu,t,P,rec,ini,i1,s,INI,ope,G,t0:

t0:=time():

gu:=IntC(C,x,p,n,N,U,MaxC):

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

INI:=gu[1]:
ope:=gu[2]:

lu:=CtoR(C,t):

ini:=C[1]:
rec:=C[2]:

print(``):
print(`--------------------------------------------------------------------------------------------`):
print(``):
print(`A Linear Recurrence With Polynomial Coefficients Satisfied by The Inegral of the`, p, `-th power of  Certain C-finite Polynomial Sequence`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Proposition: Let `, P[n](x), `be the sequence of polynomial definied by the following recurrence `):
print(``):
print(P[n](x)=add(rec[i1]*P[n-i1](x),i1=1..nops(rec))):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(P[i1]=ini[i1+1],i1=0..nops(ini)-1)):
print(``):
print(`Equivalently, in terms of generating function`):
print(``):
print(Sum(P[n](x)*t^n,n=0..infinity)=lu ):
print(``):
print(`Define the Umbra, U, to be equal to`, U, `when applied to `, x^n, `and extended linearly `):
print(``):
print(`Let's define the sequence of numbers`):
print(``):
print(` s(n), to be the action of U applied to`, P[n](x)^p):
print(``):
print(`Then s(n) satisfies the following linear recurrence equation with polynomial coefficients `):
print(``):
print(add(coeff(ope,N,i1)*s(n+i1),i1=0..degree(ope,N))=0):
print(``):
print(`Subject to the initial condition`):
print(``):
print(seq(s(i1)=INI[i1],i1=1..nops(INI))):
print(``):

G:=SeqFromRec(gu,n,N,GADOL)[GADOL]:

print(`Just for fun, using this recurrence we get that `):
print(``):
print(s(GADOL)=G):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate. `):
print(``):
print(`----------------------------------------------------`):
print(``):
end:



#IntCSv(C,x,p,n,N,U,MaxC,GADOL):  Verbose version of IntCS(C,x,r,n,N,U,MaxC) (q.v.)
#Try
#IntCSv([[1,1],[1,x]],x,1,n,N,1/(n+1),10);

IntCSv:=proc(C,x,r,n,N,U,MaxC,GADOL) local gu,lu,t,P,rec,ini,i1,s,INI,ope,G,t0,StarP:

t0:=time():

gu:=IntCS(C,x,r,n,N,U,MaxC):

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

INI:=gu[1]:
ope:=gu[2]:

lu:=CtoR(C,t):

ini:=C[1]:
rec:=C[2]:

print(``):
print(`--------------------------------------------------------------------------------------------`):
print(``):
print(`A Linear Recurrence With Polynomial Coefficients Satisfied by The Inegral of the`, r, `-th Part of a  Certain C-finite Polynomial Sequence`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Proposition: Let `, P[n](x), `be the sequence of polynomial definied by the following recurrence `):
print(``):
print(P[n](x)=add(rec[i1]*P[n-i1](x),i1=1..nops(rec))):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(P[i1]=ini[i1+1],i1=0..nops(ini)-1)):
print(``):
print(`Equivalently, in terms of generating function`):
print(``):
print(Sum(P[n](x)*t^n,n=0..infinity)=lu ):
print(``):
print(`Define the Umbra, U, to be equal to`, U, `when applied to `, x^n, `and extended linearly `):
print(``):
print(`Also define Star(P)(x)=P(1/x)*x^n, where n is the degree of the polynomaial P(x)`):
print(``):
print(`Let's define the sequence of numbers`):
print(``):
print(` s(n), to be the action of U applied to`, P[n*r](x)*StarP[n*r](x)) :
print(``):
print(`Then s(n) satisfies the following linear recurrence equation with polynomial coefficients `):
print(``):
print(add(coeff(ope,N,i1)*s(n+i1),i1=0..degree(ope,N))=0):
print(``):
print(`Subject to the initial condition`):
print(``):
print(seq(s(i1)=INI[i1],i1=1..nops(INI))):
print(``):

G:=SeqFromRec(gu,n,N,GADOL)[GADOL]:

print(`Just for fun, using this recurrence we get that `):
print(``):
print(s(GADOL)=G):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate. `):
print(``):
print(`----------------------------------------------------`):
print(``):
end:



#SHK1a(N): The first N terms of the left hand-side of Theorem 1, first part
#Seon-Hong Kim's paper (Rmanujan j. 38(2015), 629-639, Eq. (1.3)
SHK1a:=proc(N) local n,j:
[seq(2*add((-1)^j/(2*j+1),j=0..2*n),n=1..N)]:
end:


#SHK1b(N): The first N terms of the left hand-side of Theorem 1, second part
#Seon-Hong Kim's paper (Rmanujan j. 38(2015), 629-639, Eq. (1.4)
SHK1b:=proc(N) local n,j:
[seq( (-1)^n*n/4^(n-1)+2/(4*n+1)+(-1)^(n+1)*n/4^n*
(add((-1)^j*4^j/(2*j-1),j=1..n)+
add((-1)^trunc((j+2)/2)*2^(j+2)/(2*j+1),j=1..2*n)
),n=1..N)]:
end:




#IntCaz(C,x,K,p,t,Dt,a,b): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)
#in the variable x, K, a nice function of x, a positive integer p, a symbol t, a symbol Dt for diffferentiatin with respect to t,
#and numbers a and b
#Outputs a linear differential equation ope(t,Dt)F(t)=Rhs(t)
#for the generating function of the sequence  int(C_n(x)*K(x),x=a..b):
#It returns the pair [ope,Rhs]
#Try:
#IntCaz([[1,1],[1,x]],x,1,t,Dt,0,1);

IntCaz:=proc(C,x,K,p,t,Dt,a,b) local gu,i1:
gu:=C:
for i1 from 2 to p do
gu:=Kefel(gu,C):
od:

gu:=CtoR(gu,t):


AZcI(gu*K,x,t,Dt,a,b);


end:




#IntCSaz(C,x,r,t,Dt,a,b): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)
#in the variable x, a positive integer r,  symbols t and Dt (Dt denotes differentiation with repect to t), and numbers a and b such that a<b,
#finds an inhomogeneous linear  differential equation of the form ope(t,Dt)F(t)=Rhs(t) for
#for the generating function in t, for  Int(C_(rn)(x)Star(C_rn(x),x=a..b).
#It returns the pair [ope,Rhs], where INI are the initial conditions and ope is the operator
#Try:
#IntCSaz([[1,1],[1,x]],x,2,t,Dt,0,1);

IntCSaz:=proc(C,x,r,t,Dt,a,b) local gu,i,K,gu1:

for K from 5 by 5 to 100 do 
gu:=SeqFromRecC(C,K):
gu:=[seq(expand(gu[i]*Star(gu[i],x)),i=1..nops(gu))]:
gu1:=subs(x=7,gu):
if GuessRec(gu1)<>FAIL then
gu:=GuessRec(gu):
if gu=FAIL then
 RETURN(FAIL):
fi:

gu:=CtoR(gu,t):

RETURN(AZcI(gu,x,t,Dt,a,b)):
fi:

od:

FAIL:

end:



#IntCS(C,x,r,n,N,U,MaxC): Inputs a C finite sequence C=[INI,REC], defining a polynomial sequence C_n(x)
#in the variable x, a positive integer r, a symbol n, and an Umbra U, and a positive integer MaxC,
#guesses a (rigorously-provable) linear recurrence operator with polynomial coefficients, using the
#shift operator N,
#for the sequence U(C_(rn)(x)Star(C_rn(x),x)). 
#It returns the pair [INI,ope], where INI are the initial conditions and ope is the operator
#Try:
#IntCS([[1,1],[1,x]],x,2,n,N,1/(n+1),10);

IntCS:=proc(C,x,r,n,N,U,MaxC) local gu,ORDER,L,K,i,ope,INI:
K:=max(seq(seq((2+L-ORDER)*(1+ORDER),ORDER=0..L),L=0..MaxC))+11:

gu:=SeqFromRecC(C,r*K+1):
gu:=[op(2..nops(gu),gu)]:
gu:=[seq(gu[r*i],i=1..K)]:
gu:=[seq(expand(gu[i]*Star(gu[i],x)),i=1..nops(gu))]:
gu:=eval([seq(ApplyU(gu[i],x,n,U),i=1..nops(gu))]):


ope:=Findrec(gu,n,N,MaxC):


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

INI:=eval([op(1..degree(ope,N),gu)]):
[INI,collect(ope,N)]:
end:


#IntCazV(C,x,K,p,t,a,b):  Verbose version of IntCaz(C,x,K,p,t,Dt,a,b) (q.v)
#Try
#IntCazV([[1,1],[1,x]],x,1,t,0,1);

IntCazV:=proc(C,x,K,p,t,a,b) local gu, Dt,ini,lu,rec,ope,yemin,s,n,f,i1,P,t0:

t0:=time():

lu:=CtoR(C,t):

gu:=IntCaz(C,x,K,p,t,Dt,a,b):

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

ope:=gu[1]:
yemin:=gu[2]:

ini:=C[1]:
rec:=C[2]:

print(``):
print(`--------------------------------------------------------------------------------------------`):
print(``):
print(`A Linear Differential Equation With Polynomial Coefficients Satisfied by The Generating Function`):
print(`of the Inegral of the`, p, `-th power of  Certain C-finite Polynomial Sequence times`, K, `From x=`, a, `to x=`, b):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Proposition: Let `, P[n](x), `be the sequence of polynomial definied by the following recurrence `):
print(``):
print(P[n](x)=add(rec[i1]*P[n-i1](x),i1=1..nops(rec))):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(P[i1]=ini[i1+1],i1=0..nops(ini)-1)):
print(``):
print(`Equivalently, in terms of generating function`):
print(``):
print(Sum(P[n](x)*t^n,n=0..infinity)=lu ):
print(``):
print(``):
print(`Let's define the sequence of numbers`):
print(``):
print(s(n)=Int(P_n(x)^p*K,x=a..b)):
print(``):
print(`and let`, f(t), `be the ordinary generating function`):
print(``):
print(f(t)=Sum(s(n)*t^n,n=0..infinity)):
print(``):
print(`Then f(t) satisfies the following inhomogeneous linear differntial equation with polynomial coefficients `):
print(``):
print(coeff(ope,Dt,0)*f(t)+ add(coeff(ope,Dt,i1)*diff(f(t),t$i1),i1=1..degree(ope,Dt))=yemin):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(coeff(ope,Dt,i1)*diff(f(t),t$i1),i1=1..degree(ope,Dt))=yemin):
print(``):

print(`The proof, using the beautiful continuous Almkvist-Zeilberger algorithm is left to the readers. `):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate. `):
print(``):
print(`----------------------------------------------------`):
print(``):
end: