###################################################################### ##Prodinger.txt: Save this file as Prodinger.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read `Prodinger.txt` # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Aug. 5, 2021 Digits:=100: ##start Cfinite.txt print(`Created: Aug. 2021`): print(`version of Feb. 2013`): print(` This is Prodinger.txt `): print(`It accompanyies the article `): print(`Automatic Generation of Convolution Identities for C-finite sequences`): 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 type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`------------------------------`): ezra:=proc() if args=NULL then print(`The main procedures are: HelmutT, ProdingerT, Sipur1, Sipur2`): elif nargs=1 and args[1]=HelmutT then print(`HelmutT(R1,R2, x): Given rational functions R1 and R2 (possibly the same) that are the generating function of a sequence a(k), b(k) finds the rational function whose coefficiets are Sum(binomial(n,k)*a[k]*b[n-k],k=0..n).`): print(`Try: HelmutT(x/(1-x-x^2),x/(1-x-x^2), x);`): elif nargs=1 and args[1]=ProdingerT then print(`ProdingerT(R, x): Given rational functions R that are the generating function of a sequence a(k), finds the rational function whose coefficiets are Sum(binomial(n,k)*a[k]*a[n-k],k=0..n).`): print(`Try: ProdingerTT(x/(1-x-x^2), x);`): elif nargs=1 and args[1]=Sipur1 then print(`Sipur1(K): inputs a positive integer N and outputs an article about the convolution identities for k-bonacci number from k=2 (the Fibonacci case) to k=N. Try:`): print(`Sipur1(10);`): elif nargs=1 and args[1]=Sipur2 then print(`Sipur2(K): inputs a positive integer N and outputs an article about the convolution identities for k1-bonacci numbers and k2-bonacci from 2<=k1<=k2<=K. Try:`): print(`Sipur2(5);`): else print(`There is no ezra for`, args): fi: end: ##start Cfinite.txt print(`-----------------------`): print(`For a list of the Cfinite procedures type ezraC();, for help with`): print(`a specific procedure, type ezraC(procedure_name); .`): print(``): print(`------------------------------`): with(combinat): with(numtheory): 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: ezra1C:=proc() if args=NULL then print(` The supporting procedures are: 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, Factorize, FactorizeI1`): print(` FindCH, HelmutT, GuessNLR, GuessPR, GuessRec, IsProd, IsProdG `): print(` Kefel, KefelR, KefelS, KefelSm, Khibur, KhiburS`): print(` mSect, SeqFromRec `): print(` `): Cfinite.txt 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])=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])=HelmutT then print(`HelmutT(R,x): Given a rational function R that is the generating function of a sequence a(k), finds the rational function whose coefficiets are Sum(binomial(n,k)*a[k]*a[n-k],k=0..n).`): print(`Try: HelmutT(x/(1-x-x^2),x);`): 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])=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])=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])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=`, 2*d+3 ): RETURN(FAIL): fi: var:={seq(a[i],i=1..d)}: eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): else RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]): fi: end: #GuessRec(L): inputs a sequence L and tries to guess #a recurrence operator with constant cofficients #satisfying it. It returns the initial values and the operator # as a list. For example try: #GuessRec([1,1,1,1,1,1]); GuessRec:=proc(L) local gu,d: for d from 1 to trunc(nops(L)/2)-2 do gu:=GuessRec1(L,d): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: #Khibur(S1,S2): inputs two C-finite sequences S1,S2, and outputs #their sum, For example, try: #Khibur([[1,2],[3,4]],[[2,5],[-1,6]]); Khibur:=proc(S1,S2) local d1,d2,L1,L2,L,i,K: if nops(S1[1])<>nops(S1[2]) or nops(S2[1])<>nops(S2[2]) then print(`Bad input`): RETURN(FAIL): fi: d1:=nops(S1[1]): d2:=nops(S2[1]): K:=2*(d1+d2)+4: L1:=SeqFromRec(S1,K): L2:=SeqFromRec(S2,K): L:=[seq(L1[i]+L2[i],i=1..K)]: GuessRec(L): end: #Kefel(S1,S2): inputs two C-finite sequences S1,S2, and outputs #their product, For example, try: #Kefel([[1,2],[3,4]],[[2,5],[-1,6]]); Kefel:=proc(S1,S2) local d1,d2,L1,L2,L,i,K: if nops(S1[1])<>nops(S1[2]) or nops(S2[1])<>nops(S2[2]) then print(`Bad input`): RETURN(FAIL): fi: d1:=nops(S1[1]): d2:=nops(S2[1]): K:=2*(d1*d2)+4: L1:=SeqFromRec(S1,K): L2:=SeqFromRec(S2,K): L:=[seq(L1[i]*L2[i],i=1..K)]: GuessRec(L): end: #KhiburSslow(a,b,d1,d2): The recurrence operator #satisfied by the sum of any solution of #the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1] #and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1] #For example, try: #KhiburSslow(a,b,2,2); KhiburSslow:=proc(a,b,d1,d2) local S1,S2,i: S1:=[[seq(ithprime(5+i),i=1..d1)],[seq(a[i],i=1..d1)]]: S2:=[[seq(ithprime(9+d1+i),i=1..d2)],[seq(b[i],i=1..d2)]]: Khibur(S1,S2)[2]: end: #KefelSslow(a,b,d1,d2): The recurrence operator #satisfied by the product of any solution of #the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1] #and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1] #For example, try: #KefelSslow(a,b,2,2); KefelSslow:=proc(a,b,d1,d2) local S1,S2,i: S1:=[[seq(ithprime(5+i),i=1..d1)],[seq(a[i],i=1..d1)]]: S2:=[[seq(ithprime(9+d1+i),i=1..d2)],[seq(b[i],i=1..d2)]]: Kefel(S1,S2)[2]: end: #KhiburS(a,b,d1,d2): The recurrence operator #satisfied by the sum of any solution of #the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1] #and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1] #For example, try: #KhiburS(a,b,2,2); KhiburS:=proc(a,b,d1,d2) local S1,S2,i,x,y,L,C,var,eq,L1,L2,gu: S1:=[[seq(x[-d1+i],i=0..d1-1)],[seq(a[i],i=1..d1)]]: S2:=[[seq(y[-d2+i],i=0..d2-1)],[seq(b[i],i=1..d2)]]: L1:=SeqFromRec(S1,d1+d2+1): L2:=SeqFromRec(S2,d1+d2+1): L:=expand([seq(L1[i]+L2[i],i=1..d1+d2+1)]): gu:=expand(L[d1+d2+1]-add(C[i]*L[d1+d2+1-i],i=1..d1+d2)): var:={seq(C[i],i=1..d1+d2)}: eq:={seq(coeff(gu,x[-i],1),i=1..d1),seq(coeff(gu,y[-i],1),i=1..d2)}: var:=solve(eq,var): [seq(subs(var,C[i]),i=1..d1+d2)]: end: #KefelS(a,b,d1,d2): The recurrence operator #satisfied by the product of any solution of #the recurrence f[n]=a[1]*f[n-1]+...+a[d1]*f[n-d1] #and g[n]=b[1]*g[n-1]+...+b[d1]*g[n-d1] #For example, try: #KefelS(a,b,2,2); KefelS:=proc(a,b,d1,d2) local S1,S2,i,x,y,L,C,var,eq,L1,L2,gu,j: S1:=[[seq(x[-d1+i],i=0..d1-1)],[seq(a[i],i=1..d1)]]: S2:=[[seq(y[-d2+i],i=0..d2-1)],[seq(b[i],i=1..d2)]]: L1:=SeqFromRec(S1,d1*d2+1): L2:=SeqFromRec(S2,d1*d2+1): L:=expand([seq(L1[i]*L2[i],i=1..d1*d2+1)]): gu:=expand(L[d1*d2+1]-add(C[i]*L[d1*d2+1-i],i=1..d1*d2)): var:={seq(C[i],i=1..d1*d2)}: eq:={seq(seq(coeff(coeff(gu,x[-i],1),y[-j],1),j=1..d2),i=1..d1)}: var:=solve(eq,var): [seq(subs(var,C[i]),i=1..d1*d2)]: end: #ProdIndicatorN(m,n): The profile of multiplicity of #the (m*n)^2 ratios of a list of m*n numbers to #be the Cartesian product of a set of m numbers with #a set of n numbers. For example, try: #ProdIndicatorN(2,2); ProdIndicatorN:=proc(m,n) local a,b,i,j,gu,x,lu,mu: option remember: gu:=[seq(seq(a[i]*b[j],j=1..n),i=1..m)]: ProdProfileE(gu): end: #ProdIndicator(m,n): The profile of multiplicity of #the (m*n)^2 ratios of a list of m*n numbers to #be the Cartesian product of a set of m numbers with #a set of n numbers. For example, try: #ProdIndicator(2,2); ProdIndicator:=proc(m,n) local a,b,i,j,gu,x,lu,mu: option remember: gu:=[seq(seq(a[i]*b[j],j=1..n),i=1..m)]: lu:=[seq(seq(gu[i]/gu[j],j=1..m*n),i=1..m*n)]: lu:=add(x[lu[i]],i=1..nops(lu)): mu:=[]: for i from 1 to nops(lu) do if type(op(1,op(i,lu)),integer) then mu:=[op(mu),op(1,op(i,lu))]: else: mu:=[op(mu),1]: fi: od: sort(mu): end: #ProdProfile(L,eps): Given a list of complex numbers finds its #(approximate profile). For example, try: #ProdProfile([1,2,3,6],10^(-4)); ProdProfile:=proc(L,eps) local i,j,lu,mu,S,gu: gu:=[seq(seq(L[i]/L[j],j=1..nops(L)),i=1..nops(L))]: mu:=[]: S:={seq(i,i=1..nops(gu))}: while S<>{} do lu:=Krovim(gu,S[1],eps): mu:=[op(mu),nops(lu)]: S:=S minus lu: od: sort(mu): end: #IsProd(f,t,m,n,eps): Given a rational function f of t #decides whether the C-finite sequence of its coeeficients #is a product of a C-finite sequences of order m and n. #For example, try: #IsProd(1/((1-t)*(1-2*t)*(1-3*t)*(1-6*t)),t,2,2); IsProd:=proc(f,t,m,n,eps) local lu,gu: lu:=denom(f): if degree(lu,t)<>m*n then print(`The degree of the denominator of`, f, `should be`, m*n): RETURN(FAIL): fi: gu:=evalf([solve(lu,t)]): if ProdProfile(gu,eps)=ProdIndicator(m,n) then RETURN(true): else RETURN(false,gu): fi: end: #CtoR(S,t), the rational function, in t, whose coefficients #are the members of the C-finite sequence S. For example, try: #CtoR([[1,1],[1,1]],t); CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L: if not (type(S,list) and nops(S)=2 and type(S[1],list) and type(S[2],list) and nops(S[1])=nops(S[2]) and type(t, symbol) ) then print(`Bad input`): RETURN(FAIL): fi: D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])): N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])): L1:=expand(D1*N1): L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1): f:=L1/D1: L:=degree(D1,t)+10: f1:=taylor(f,t=0,L+1): if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)): RETURN(FAIL): else RETURN(f): fi: end: #RtoC(R,t): Given a rational function R(t) finds the #C-finite description of its sequence of coefficients #For example, try: #RtoC(1/(1-t-t^2),t); RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a: R1:=normal(R): D1:=denom(R1): d:=degree(D1,t): if degree(numer(R1),t)>=d then print(`The degree of the top must be less than the degree of the bottom`): RETURN(FAIL): fi: a:=coeff(D1,t,0): S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]: f1:=taylor(R,t=0,d+1): S1:=[seq(coeff(f1,t,i),i=0..d-1)]: [S1,S2]: end: #KefelR(R1,R2,t): inputs rational functions of t, #R1 and R2, #and outputs the rational function #whose coeffs. is their Hadamard product, for example try: #KefelR(1/(1-t-t^2),1/(1-t-t^3),t); KefelR:=proc(R1,R2,t) local S1,S2,S: S1:=RtoC(R1,t): S2:=RtoC(R2,t): S:=Kefel(S1,S2): normal(CtoR(S,t)): end: #Krovim(L,i,eps): Given a list L and an index between #1 and nops(L), finds the set of indices (including i) #such that abs(L[j]-L[i]){} 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)]FAIL then gu1:=SeqFromRec(lu,4*k): if {seq(type(gu1[i1],integer),i1=1..nops(gu1))}={true} then gu11:=[seq(gu[i1]/gu1[i1],i1=1..4*k)]: if {seq(type(gu11[i1],integer),i1=1..nops(gu11))}={true} then gu:=SeqFromRec(C,10*k): gu1:=SeqFromRec(lu,10*k): gu11:=[seq(gu[i1]/gu1[i1],i1=1..10*k)]: C1:=GuessRec1(gu11,nops(C[2])/k): if C1<>FAIL then RETURN(lu,C1): fi: fi: fi: fi: od: FAIL: end: #GuessPR1(C1,C2,d,x,y): Inputs two C-finite sequences and a pos. integer #d and outputs a polynomial P of degree d such that P(C1[n],C2[n])=0 #for all n. If none found it returns FAIL. For example to solve #the Hugh Montgomery proposed (but not selected) Putnam problem, #(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0, #where F(n) is the Fibonacci sequence type #GuessPR1([[0,1],[1,1]],[[1,1],[1,1]],4,x,y); GuessPR1:=proc(C1,C2,d,x,y) local eq,var,gu,a,S1,S2,N,var1,pol,n1,mu,v1: gu:=GenPol([x,y],a,d): var:=gu[2]: gu:=gu[1]: N:=nops(var): S1:=SeqFromRec(C1,N+10): S2:=SeqFromRec(C2,N+10): eq:={seq(subs({x=S1[n1],y=S2[n1]},gu),n1=1..N+10)}: var1:=solve(eq,var): pol:=subs(var1,gu): if pol=0 then RETURN(FAIL): fi: S1:=SeqFromRec(C1,N+30): S2:=SeqFromRec(C2,N+30): if {seq(normal(subs({x=S1[n1],y=S2[n1]},pol)),n1=N+11..N+30)}<>{0} then print(`The conjecture`, pol , `did not materialize`): RETURN(FAIL): fi: mu:={}: for v1 in var1 do if op(1,v1)=op(2,v1) then mu:=mu union {op(1,v1)}: fi: od: if nops(mu)>1 then print(`There are than one solution`): fi: factor(coeff(pol,mu[1],1)): end: #GuessPR(C1,C2,d,x,y): Inputs two C-finite sequences and a pos. integer #d and outputs a polynomial P of degree <=d such that P(C1[n],C2[n])=0 #for all n. If none found it returns FAIL. For example to solve #the Hugh Montgomery proposed (but not selected) Putnam problem, #(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0, #where F(n) is the Fibonacci sequence type #GuessPR([[0,1],[1,1]],[[1,1],[1,1]],6,x,y); GuessPR:=proc(C1,C2,d,x,y) local gu,d1: for d1 from 1 to d do gu:=GuessPR1(C1,C2,d1,x,y): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: #GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a, #and a non-negative integer deg, outputs a generic polynomial #of total degree deg,in the variables Resh, indexed by the letter a of #followed by the set of coeffs. For example, try: #GenPol([x,y],a,1); GenPol:=proc(Resh,a,deg) local var,POL1,i,x,Resh1,POL: if not type(Resh,list) or nops(Resh)=0 then ERROR(`Bad Input`): fi: x:=Resh[1]: if nops(Resh)=1 then RETURN(convert([seq(a[i]*x^i,i=0..deg)],`+`),{seq(a[i],i=0..deg)}): fi: Resh1:=[op(2..nops(Resh),Resh)]: var:={}: POL:=0: for i from 0 to deg do POL1:=GenPol(Resh1,a[i],deg-i): var:=var union POL1[2]: POL:=expand(POL+POL1[1]*x^i): od: POL,var: end: #GuessNLR1(C1,r,d,x): Inputs a C-finite sequence C and a pos. integer #d and outputs a polynomial P(x[1],..., x[r]) of degree d such that #P(C1[n],C1[n-1],C1[n-r])=0 #for all n. If none found it returns FAIL. For example to solve #the Hugh Montgomery proposed (but not selected) Putnam problem, #(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0, #where F(n) is the Fibonacci sequence type #GuessNLR1([[0,1],[1,1]],1,4,x); GuessNLR1:=proc(C1,r,d,x) local eq,var,gu,a,S1,N,var1,pol,n1,mu,v1,i,i1: gu:=GenPol([seq(x[i],i=0..r)],a,d): var:=gu[2]: gu:=gu[1]: N:=nops(var)+r: S1:=SeqFromRec(C1,N+10): eq:={seq(subs({seq(x[i1]=S1[n1+i1],i1=0..r)},gu),n1=1..N+10-r)}: var1:=solve(eq,var): pol:=subs(var1,gu): if pol=0 then RETURN(FAIL): fi: S1:=SeqFromRec(C1,N+30): if {seq( normal(subs({seq(x[i1]=S1[n1+i1-1],i1=0..r)},pol)),n1=N+11..N+30-r-1)}<>{0} then print(`The conjecture`, pol , `did not materialize`): RETURN(FAIL): fi: mu:={}: for v1 in var1 do if op(1,v1)=op(2,v1) then mu:=mu union {op(1,v1)}: fi: od: if nops(mu)>1 then print(`There are than one solution`): fi: factor(coeff(pol,mu[1],1)): end: #GuessNLR(C1,r,d,x): Inputs a C-finite sequence C and a pos. integer #d and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that #P(C1[n],C1[n-1],C1[n-r])=0 #for all n. If none found it returns FAIL. For example to solve #the Hugh Montgomery proposed (but not selected) Putnam problem, #(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0, #where F(n) is the Fibonacci sequence type #GuessNLR([[0,1],[1,1]],1,4,x); GuessNLR:=proc(C1,r,d,x) local gu,d1: for d1 from 1 to d do gu:=GuessNLR1(C1,r,d1,x): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: #GuessNLRv(C1,r,d): Verbose form of GuessNLR(C1,r,d,x), i.e. #Inputs a C-finite sequence C and a pos. integer #d and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that #P(C1[n],C1[n-1],C1[n-r])=0 #for all n. If none found it returns FAIL. For example to solve #the Hugh Montgomery proposed (but not selected) Putnam problem, #(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0, #where F(n) is the Fibonacci sequence type #GuessNLRv([[0,1],[1,1]],1,4,x); GuessNLRv:=proc(C1,r,d) local gu,x,F,n,i1: gu:=GuessNLR(C1,r,d,x): if gu=FAIL then RETURN(FAIL): fi: print(`Theorem: Let F(n) be the sequence defined by the recurrence`): print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))): print(`subject to the initial conditions`): print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))): print(`Then we have the following identity valid for every non-negative integer n`): print(subs({seq(x[i1]=F(n+i1),i1=0..r)},gu)=0): print(``): print(`Proof: Routine! (since everything in sight is C-finite).`): end: #FactorizeI1v(C1,k): Verbose form of FactorizeI1(C1,k), i.e. FactorizeI1v:=proc(C1,k) local gu,C2,C3,i1,n: gu:=FactorizeI1(C1,k): C2:=gu[1]: C3:=gu[2]: if gu=FAIL then RETURN(FAIL): fi: print(`Theorem: Let F(n) be the sequence defined by the recurrence`): print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))): print(`subject to the initial conditions`): print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))): print(``): print(`Let G(n) be the sequence defined on non-negative integers by the recurrence`): print(G(n)=add(C2[2][i1]*G(n-i1),i1=1..nops(C2[2]))): print(`subject to the initial conditions`): print(seq(G(i1-1)=C2[1][i1],i1=1..nops(C2[1]))): print(``): print(`Let H(n) be the sequence defined by the recurrence`): print(H(n)=add(C3[2][i1]*H(n-i1),i1=1..nops(C3[2]))): print(`subject to the initial conditions`): print(seq(H(i1-1)=C3[1][i1],i1=1..nops(C3[1]))): print(`Then the following is true for every non-negative integer n`): print(F(n)=G(n)*H(n)): print(``): print(`Proof: Routine! (since everything in sight is C-finite) .`): end: #mSect(C1,m,i): Given a C-finite sequence C1(n), finds #the C-finite description of C(m*n+i) where m is a pos. integer #and 0<=i0 then print(`Oops something is wrong `): RETURN(FAIL): else print(S1[i1], `times `, S2[i1], `is indeed `, S[i1]): fi: od: print(`QED!`): end: #CHseq(L,n,j,K): inputs a list of pairs [C,e], where C is a C-finite #sequence, and e is an affine-linear expression in the symbols #n and j, of the form a*n+b*j+c, where a and a+b are non-negative #integers. Outputs the first K terms of #Sum(L[1][1](L[1][2])*L[2][1](L[2][2]),j=0..n-1) #using up to K values. For example, try: #CHseq([[Fn(),j],[Fn(),j],[Fn(),2*n-j],20); CHseq:=proc(L,n,j,K) local i1,S,j1,n1,Max1: Max1:=0: for i1 from 1 to nops(L) do for n1 from 0 to K do for j1 from 0 to n1-1 do Max1:=max(Max1,subs({j=j1,n=n1},L[i1][2])): if subs({j=j1,n=n1},L[i1][2])<0 then RETURN(FAIL): fi: od: od: od: for i1 from 1 to nops(L) do S[i1]:=SeqFromRec(L[i1][1],Max1+1): od: [ seq( add(mul(S[i1][subs({j=j1,n=n1},L[i1][2])+1],i1=1..nops(L)),j1=0..n1-1), n1=0..K)]: end: #CH(L,n,j): #Find a sequence of the style of Curtis Greene and #Herb Wilf ("Closed form summation of C-finite sequences", #Trans. AMS) #inputs a list of pairs [C,e], where C is a C-finite #sequence, and e is an affine-linear expression in the symbols #n and j, of the form a*n+b*j+c, where a and a+b are non-negative #integers. Outputs the C-finite description of #Sum(L[1][1](L[1][2])*L[2][1](L[2][2]),j=0..n-1) #using up to K values. For example, try: #CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]); CH:=proc(L,n,j) local K, C: C:=GuessRec(CHseq(L,n,j,20)): if C<>FAIL then RETURN(C): fi: for K from 25 by 5 while C=FAIL do C:=GuessRec(CHseq(L,n,j,K)): od: C: end: #####New stuff #GuessCH1(C1,r,d,x,C2): Inputs a C-finite sequence C1 and a pos. integer #d and outputs a polynomial P(x[1],..., x[r]) of degree d such that #P(C1[n],C1[n-1],C1[n-r])=C2[n] #for all n. If none found it returns FAIL. For example to reproduce #Eq. (12) of the Greene-Wilf article # type #GuessCH1([[0,1],[1,1]],1,3,x,CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j)); GuessCH1:=proc(C1,r,d,x,C2) local eq,var,gu,a,S1,N,var1,pol,n1,mu,v1,i,i1,S2: gu:=GenPol([seq(x[i],i=0..r)],a,d): var:=gu[2]: gu:=gu[1]: N:=nops(var)+r: S1:=SeqFromRec(C1,N+11): S2:=SeqFromRec(C2,N+11): eq:={seq(subs({seq(x[i1]=S1[n1+i1+1],i1=0..r)},gu)-S2[n1+1],n1=1..N+10-r)}: var1:=solve(eq,var): pol:=subs(var1,gu): if pol=0 then RETURN(FAIL): fi: S1:=SeqFromRec(C1,N+30): S2:=SeqFromRec(C2,N+30): if {seq( normal(subs({seq(x[i1]=S1[n1+i1+1],i1=0..r)},pol)-S2[n1+1]),n1=N+11..N+30-r-2) }<>{0} then RETURN(FAIL): fi: mu:={}: for v1 in var1 do if op(1,v1)=op(2,v1) then mu:=mu union {op(1,v1)}: fi: od: if mu<>{} then factor(subs({seq(mu1=0,mu1 in mu)},pol)): else factor(pol): fi: end: #FindCHg(C1,r,d, x,C2): Inputs a C-finite sequence C1 and a pos. integer #d and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that #P(C1[n],C1[n-1],C1[n-r])=C2[n] #for all n. If none found it returns FAIL. For example to reproduce #Eq. (12) of the Greene-Wilf article. #If there is none of degree <=d then it returns FAIL # type #FindCHg([[0,1],[1,1]],2,4,x,CH([[Fn(),j],[Fn(),j],[Fn(),2*n-j]],n,j)); FindCHg:=proc(C1,r,d,x,C2) local gu,d1: for d1 from 1 to d do gu:=GuessCH1(C1,r,d1,x,C2): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: ###end new stuff #FindCH(C1,r,d, x,L,n,j): Inputs a C-finite sequence C1 and pos. integers #r and d, and a letter x, and a list L of affine-linear expressions #in n,j, and #and outputs a polynomial P(x[1],..., x[r]) of degree<= d such that #P(C1[n],C1[n-1],C1[n-r])=C2[n] #where C2(n) is the sum(C1[L[1]]*C1[L[2]*...C1[L[nops(L)],j=0..n-1): #If none found it returns FAIL. For example to reproduce #Eq. (12) of the Greene-Wilf article. #If there is none of degree <=d then it returns FAIL # type #FindCH(Fn(),1,3,x,[j,j,2*n-j],n,j); FindCH:=proc(C1,r,d,x,L,n,j) local C2,i1,L1: L1:=[seq([C1,L[i1]],i1=1..nops(L))]: FindCHg([[0,1],[1,1]],r,d,x,CH(L1,n,j)): end: #FindCHv(C1,r,d, L,n,j): verbose version of FindCH(C1,r,d,x, L,n,j): # type #FindCHv(Fn(),1,3,[j,j,2*n-j],n,j); FindCHv:=proc(C1,r,d,L,n,j) local gu,x,F,i1,t: gu:=FindCH(C1,r,d,x,L,n,j): if gu=FAIL then RETURN(gu): fi: print(`Theorem: Let F(n) be the sequence defined by the recurrence`): print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))): print(`subject to the initial conditions`): print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))): print(``): print(`Equivalently, the ordinary generating function of F(n) w.r.t. to t is`): print(Sum(F(n)*t^n,n=0..infinity)=CtoR(C1,t)): print(`Also Let G(n) be the sequence defined by `): print(G(n)=Sum(mul(F(L[i1]),i1=1..nops(L)),j=0..n-1)): print(`then we have the following closed-form expression for G(n)`): print(G(n)=subs({seq(x[i1]=F(n+i1),i1=0..r)},gu)): print(``): print(`Proof: Both sides are C-finite, so it is enough to check the`): print(`first few terms, and this is left to the dear reader.`): end: #FindCHvTN(C1,r,d, L,n,j,mispar): verbose version of FindCH(C1,r,d,x, L,n,j): # type #FindCHvTN(Fn(),1,3,[j,j,2*n-j],n,j,8); FindCHvTN:=proc(C1,r,d,L,n,j,mispar) local gu,x,F,i1,t: gu:=FindCH(C1,r,d,x,L,n,j): if gu=FAIL then RETURN(gu): fi: print(`Theorem Number`, mispar): print( `Let F(n) be the sequence defined by the recurrence`): print(F(n)=add(C1[2][i1]*F(n-i1),i1=1..nops(C1[2]))): print(`subject to the initial conditions`): print(seq(F(i1-1)=C1[1][i1],i1=1..nops(C1[1]))): print(``): print(`Equivalently, the ordinary generating function of F(n) w.r.t. to t is`): print(Sum(F(n)*t^n,n=0..infinity)=CtoR(C1,t)): print(`Also Let G(n) be the sequence defined by `): print(G(n)=Sum(mul(F(L[i1]),i1=1..nops(L)),j=0..n-1)): print(`then we have the following closed-form expression for G(n)`): print(G(n)=subs({seq(x[i1]=F(n+i1),i1=0..r)},gu)): print(``): print(`Proof: Both sides are C-finite, so it is enough to check the`): print(`first few terms, and this is left to the dear reader.`): end: #Curt1(j,n,M): outputs #all {a1n+b1j} with 0<=a1<=M and 0<=a1+b1<=M #Type: Curt1(j,n,4); Curt1:=proc(j,n,M) local a1,b1: [seq(seq(a1*n+b1*j, b1=-a1..-1),a1=0..M), seq(seq(a1*n+b1*j, b1=1..M-a1 ),a1=0..M)]: end: #IV1(m,K): all the lists of weakly-increasing positive integers of size m #with largest element K. Do IV(3,5); IV1:=proc(m,K) local gu,i,gu1,gu11: option remember: if m=1 then RETURN([seq([i],i=1..K)]): fi: gu:={}: for i from 1 to K do gu1:=IV1(m-1,i): gu:=gu union {seq([op(gu11),K],gu11 in gu1)}: od: gu: end: #IV(m,K): all the lists of weakly-increasing positive integers of size m #with largest element <=K. Do IV(3,5); IV:=proc(m,K) local i: {seq(op(IV1(m,i)),i=1..K)}: end: Khevre:=proc(m,j,n,M) local gu,mu,mu1,i1: gu:=Curt1(j,n,M): mu:=IV(m-1,nops(gu)): [seq([j,seq(gu[mu1[i1]],i1=1..m-1)],mu1 in mu)]: end: #SeferGW(m,d,M): outputs a book of Fibonacci identities with #right hand sides of degree at most d in F(n),F(n+1) #for Greene-Wilf type sums where the summand is #of the form F(j)*F(a1n+b1j)*... with m terms, for example, try: #SeferGW(2,5,2); SeferGW:=proc(m,d,M) local n,j,gu,mu,mu1,x,C1,co: co:=0: print(`A Book of Definite Summation Fibonacci Identities`): print(`in the style of Curtis Greene and Herbert Wilf`): print(``): print(`By Shalosh B. Ekhad `): C1:=Fn(): mu:=Khevre(m,j,n,M): for mu1 in mu do gu:=FindCH(C1,1,d,x,mu1,n,j): if gu<>FAIL then co:=co+1: print(`------------------------------------------------`): FindCHvTN(C1,1,d,mu1,n,j,co): fi: od: print(`------------------------------------------------`): print(`I hope that you enjoyed, dear reader, these`, co, `beautiful and`): print(`deep theorems. `): end: #GenHomPol(Resh,a,deg): Given a list of variables Resh, a symbol a, #and a non-negative integer deg, outputs a generic #homogeneous polynomial #of degree deg,in the variables Resh, indexed by the letter a of #followed by the set of coeffs. For example, try: #GenHomPol([x,y],a,1); GenHomPol:=proc(Resh,a,deg) local var,POL1,i,x,Resh1,POL: if not type(Resh,list) or nops(Resh)=0 then ERROR(`Bad Input`): fi: x:=Resh[1]: if nops(Resh)=1 then RETURN(a[1]*Resh[1]^deg,{a[1]}): fi: Resh1:=[op(2..nops(Resh),Resh)]: var:={}: POL:=0: for i from 0 to deg do POL1:=GenHomPol(Resh1,a[i],deg-i): var:=var union POL1[2]: POL:=expand(POL+POL1[1]*x^i): od: POL,var: end: #GuessInNLR1(C1,r,d,x): Inputs a C-finite sequence C and a pos. integer #d and outputs an #interesting polynomial P(x[0],..., x[r]) of degree d such that #of the form x[0]*P1(x[1], ..., x[r])+P2(x[1], ...,x[r])+c #where P1,P2 are homog. of degree d, d-1 resp. such that #P(C1[n],C1[n-1],C1[n-r])=0 #for all n. If none found it returns FAIL. For example to solve #the Hugh Montgomery proposed (but not selected) Putnam problem, #(finding a polynomial of degree 4 such that P(F(n),F(n+1))=0, #where F(n) is the Fibonacci sequence type #GuessInNLR1([[0,1],[3,-1]],2,2,x); GuessInNLR1:=proc(C1,r,d,x) local eq,var, gu,a,b,S1,N,var1,pol,n1,mu,v1,i,i1,gu1,gu2,c,lu,mu1: gu1:=GenHomPol([seq(x[i],i=1..r)],a,d-1): gu2:=GenHomPol([seq(x[i],i=1..r)],b,d): var:=gu1[2] union gu2[2] union {c}: gu:=expand(gu1[1]*x[0]+gu2[1]+c): N:=nops(var)+r: S1:=SeqFromRec(C1,N+10): eq:={seq(subs({seq(x[i1]=S1[n1+i1],i1=0..r)},gu),n1=1..N+10-r)}: var1:=solve(eq,var): pol:=subs(var1,gu): if pol=0 then RETURN(FAIL): fi: S1:=SeqFromRec(C1,N+40): if {seq( normal(subs({seq(x[i1]=S1[n1+i1-1],i1=0..r)},pol)),n1=N+11..N+40-r-1)}<>{0} then print(`The conjecture`, pol , `did not materialize`): RETURN(FAIL): fi: mu:={}: for v1 in var1 do if op(1,v1)=op(2,v1) then mu:=mu union {op(1,v1)}: fi: od: lu:={}: for i from 1 to nops(mu) do mu1:=factor(coeff(pol,mu[i],1)): if not type(mu1,`*`) then lu:=lu union {mu1}: fi: od: lu: end: ##End Cfinite #ProdingerT(R,x): Given a rational function R that is the generating function of a sequence a(k), finds the rational function whose coefficiets are Sum(binomial(n,k)*a[k]*a[n-k],k=0..n). #Try: ProdingerT(,x/(1-x-x^2),x); ProdingerT:=proc(R,x) local L,f,N,k,C,n: N:=2*binomial(degree(denom(R),x)+1,2)+3: f:=taylor(R,x=0,N+1): L:=[seq(add(binomial(n,k)*coeff(f,x,k)*coeff(f,x,n-k),k=0..n),n=0..N)]: C:=GuessRec(L): if C=FAIL then RETURN(FAIL): fi: factor(CtoR(C,x)): end: #HelmutT(R1,R2,x): Given a rational function R1 that is the generating function of a sequence a(k), #and a rational function R2 that is the generating function of a sequence b(k), #finds the rational function whose coefficiets are Sum(binomial(n,k)*a[k]*b[n-k],k=0..n). #Try: HelmutT(x/(1-x-x^2),x/(1-x-x^2),x); HelmutT:=proc(R1,R2,x) local L,f1,f2,N,k,C,n: N:=2*(degree(denom(R1),x)*degree(denom(R2),x)+2): f1:=taylor(R1,x=0,N+1): f2:=taylor(R2,x=0,N+1): L:=[seq(add(binomial(n,k)*coeff(f1,x,k)*coeff(f2,x,n-k),k=0..n),n=0..N)]: C:=GuessRec(L): if C=FAIL then RETURN(FAIL): fi: factor(CtoR(C,x)): end: #Sipur1(K): inputs a positive integer N and outputs an article about the convolution identities for k-bonacci number from k=2 (the Fibonacci case) to k=N. Try: #Sipur1(10); Sipur1:=proc(K) local k,co,lu,R,t0,a,b,n,r,i,x: t0:=time(): co:=0: print(`Convolution Identities for the k-Bonacci numbers from k=2 (the Fibonacci case) all the way to k=`,K): print(``): print(`By Shalosh B. Ekhad`): print(``): for k from 2 to K do co:=co+1: R:=x/(1-add(x^i,i=1..k)): lu:=ProdingerT(R,x): print(``): print(`Theorem number`, co, `Let a(n) be the`, k, `-bonacci numbers that are defined via the generating function `): print(``): print(Sum(a(n)*x^n,n=0..infinity)= R): print(``): print(`Let b(n) be the binomial convolution`): print(``): print(b(n)=Sum(binomial(n,r)*a(r)*a(n-r),r=0..n)): print(``): print(`Then the generating function of the sequence b(n) is`): print(``): print(Sum(b(n)*x^n,n=0..infinity)= lu): print(``): print(`and in Maple notation`): print(``): lprint(lu): print(``): od: print(``): print(`-------------------------`): print(``): print(`This ends this article that took`, time()-t0, `seconds to produce. `): end: #Sipur2(K): inputs a positive integer N and outputs an article about the convolution identities for k-bonacci number and s-bonacci numbers for 2<=k