###################################################################### ## SternCF.txt Save this file as StanleyCF.txt to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read StanlyCF.txt # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## DoronZeil at gmail dot com # ###################################################################### print(`First Written: Feb. 2021: tested for Maple 2018 `): print(`Version : Feb. 2021 `): print(): print(`This is SternCF.txt, one of the Maple packages`): print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's article: `): print(`Automated Generation of Generating Functions Related to Generalized Stern's Diatomic Arrays in the footsteps of Richard Stanley `): print(): print(`The most current version is available on WWW at:`): print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/StanleyStern.txt .`): print(`Please report all bugs to: DoronZeil at gmail dot com .`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(`For a list of the supporting functions type: ezra1();`): print(): print(`For a list of the STORY functions type: ezraS();`): print(): print(`For a list of Cfinite functions, type ezraCfinite();`): print(): ezraCfinite:=proc() if args=NULL then print(`The Cfinite procedures are: CtoR, GuessRec, SeqFromRec `): print(``): else ezra(args): fi: end: ezraRS:=proc() if args=NULL then print(`The procedures implementing the functions in Richard Stanley's "Theorems and Conjectures on Some Rational Generating functions" are: `): print(` Jrk, JrkS, J2kC, J3kC `): else ezra(args): fi: end: ezraS:=proc() if args=NULL then print(`The STORY procedures are: `): print(` Paper1, Paper1mat, Paper2, Paper2mat, RSmatAno `): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(`The SUPPORTING procedures are: `): print(` Atom1, Eq, Gn, ROp, vMseq, vn, vnL, vnSeq `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` SternCF.txt: A Maple package for computing generating functions dear to Richard Stanley based on C-finite sequences`): print(`The MAIN procedures are: `): print(` RS, RSe, RSmat, RSmatV,RSmatVt, RSv`): elif nargs=1 and args[1]=Atom1 then print(`Atom1(C,A,x,i): Inputs a C-finite sequence C=[INI,REC] (in our convention), `): print(`let d be the order of C (i.e. nops(C[2])). It also inputs a list of pairs`): print(`each of the form [e[0],[e[1],e[2],...e[d]]`): print(`Outputs the polynmial`): print(`Sum(e[0]*x^(e[1]*C[i]+e[2]*C[i+1]+...e[d]*C[i+d-1]), over all the members of A. Try:`): print(` Atom1([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],x,4); `): print(` Atom1([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]] ],x,4); `): print(` Atom1([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]], [1,[0,1]] ],x,4); `): elif nargs=1 and args[1]=Eq then print(`Eq(C,A,P,x,L,Z,t): An equation for the generating function of vnL(C,A,P,x,L,n) denoted by Z[L] in terms`): print(` of other Z[L']-s, followed by its children, Try: `): print(` Eq([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[[0,[0]], [0,[0]]],Z,t); `): print(`Eq([[1,1],[1,1]], [[1,[0,0]],[1,[1,0]], [1,[0,1]]] ,1,x,[ [0,[0,0]], [0,[0,0]] ],Z,t);`): elif nargs=1 and args[1]=Gn then print(`Gn(C,A,P,x,n):`): print(` Inputs a C-finite sequence C=[INI,REC] (in our convention), `): print(`let d be the order of C (i.e. nops(C[2])). It also inputs a list of pairs`): print(`each of the form [e[0],[e[1],e[2],...e[d]]`): print(`Outputs the product of P*Atom1(C,A,x,i) i=0..n-1`): print(` Gn([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,4); `): print(` Gn([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]] ],1,x,4); `): print(` Gn([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]], [1,[0,1]] ],1,x,4); `): elif nargs=1 and args[1]=Jrk then print(`Jrk(r,k,t,x,LIMIT1): The rational generating function J_r^{(k)}(t,x) defined on p.19 in Richard Stanley's paper`): print(` "Theorems and Conjectures on Some Rational Generating functions" THIS IS FOR NUMERIC k`): print(`LIMIT1 is the cap on the number of equations. Try:`): print(`Jrk(2,2,t,x,100);`): elif nargs=1 and args[1]=JrkS then print(`JrkS(r,k,x,LIMIT1,C): The CONJECTURED rational generating function J_r^{(k)}(1,x) defined on p.13 in Richard Stanley's paper`): print(`in terms of a SYMBOLIC k, assuming that indeed the coefficients do do not depend on k. `): print(`It checks that it true for k=2,3,..,C So if C=2 then it does not check it`): print(` "Theorems and Conjectures on Some Rational Generating functions" `): print(`N is the cap on the number of equations. Try:`): print(`JrkS(2,k,x,100,3);`): elif nargs=1 and args[1]=J2kC then print(`J2kC(k,t,x): The PROVED rational generating function J_2^{(k)}(t,x) defined on p.19. Theorem 5.3. Try:`): print(`J2kC(3,t,x);`): elif nargs=1 and args[1]=J3kC then print(`J3kC(k,t,x): The CONECTURED rational generating function J_3^{(k)}(t,x) defined on p.19 in Richard Stanley's paper`): print(`CORRECTING his conjecture. Try:`): print(`J3kC(3,t,x);`): elif nargs=1 and args[1]=Paper1 then print(`Paper1(C,A,P,x,t,LIMIT1 N): inputs C,A,P,x,B,t,LIMIT1`): print(`outputs the paper with N theorems RSv(p,q,b,[i],x,t) for i from 1 to N. Try:`): print(` Paper1([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,t,200,4);`); print(` Paper1([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,t,200,2);`); print(` Paper1([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,t,1000,2);`); print(` Paper1([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,t,2000,2);`); elif nargs=1 and args[1]=Paper1mat then print(`Paper1mat(C,A,P,x,LIMIT1, N,K1): inputs C,A,P,x,B,t,LIMIT1`): print(`outputs the paper with N theorems RSmatVt(p,q,b,[i],x,K1) for i from 1 to N. Try:`): print(` Paper1mat([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,200,4,20);`); print(` Paper1mat([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,200,2,20);`); print(` Paper1mat([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,1000,2,20);`); print(` Paper1mat([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,2000,2,20);`); elif nargs=1 and args[1]=Paper2 then print(`Paper2(C,A,P,x,t,LIMIT1 N): inputs C,A,P,x,B,t,LIMIT1`): print(`outputs the paper with N theorems RSv(p,q,b,[1$i],x,t) for i from 1 to N. Try:`): print(` Paper2([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,t,200,4);`); print(` Paper2([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,t,200,4);`); print(` Paper2([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,t,200,2);`); print(` Paper2([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,t,1000,2);`); print(` Paper2([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,t,2000,2);`); elif nargs=1 and args[1]=Paper2mat then print(` Paper2mat([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,200,2,20);`); print(` Paper2mat([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,200,2,20);`); print(` Paper2mat([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,1000,2,20);`); print(` Paper2mat([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,2000,2,20);`); elif nargs=1 and args[1]=ROp then print(`ROp(L,C): Given a list L and a C-finite sequence C=[INI,REC], (let d=nops(REC)) where nops(L)=d finds the vector L1 (of the same length, d)`): print(`such that L[1]*C[n+1]+ ...+ L[d]*C[n+d]=L'[1]*C[n]+...+ L'[d]*C[n+d-1]. Try:`): print(`ROp([2,3],[[1,1],[1,1]]);`): elif nargs=1 and args[1]=RS then print(`RS(C,A,P,x,B,t,LIMIT1):`): print(` Inputs a C-finite sequence C=[INI,REC] (in our convention), `): print(`let d be the order of C (i.e. nops(C[2])). It also inputs a list of pairs`): print(`each of the form [e[0],[e[1],e[2],...e[d]]. It also inputs a pos. integer LIMIT1 such that it returns FAIL if the system of equations`): print(`is larger than it.`): print(`outputs the generating function for the sequence vn(C,A,P,x,B,n) (q.v.). `): print(` It uses the NON-EXPERIMENTAL approach, by solving a system of linear equations. Try`): print(` RS([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],t,50); `): print(` RS([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],t,50);`); print(` RS([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],t,50);`); print(` RS([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,[2],t,100);`); print(` RS([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,[2],t,1000);`); print(` RS([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,[2],t,2000);`); elif nargs=1 and args[1]=RSmat then print(`RSmat(C,A,P,x,B,LIMIT1):`): print(` Inputs a C-finite sequence C=[INI,REC] (in our convention), `): print(`let d be the order of C (i.e. nops(C[2])). It also inputs a list of pairs`): print(`each of the form [e[0],[e[1],e[2],...e[d]]. It also inputs a pos. integer LIMIT1 such that it returns FAIL if the system of equations`): print(`is larger than it.`): print(`outputs a pair [v,A] where v is a vector of integers and A is a matrix of integers such that`): print(`the generating function for the sequence vn(C,A,P,x,B,n) (q.v.). `): print(`is (I-A*t)^(-1)v[1].`): print(` RSmat([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],50); `): print(` RSmat([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],50);`); print(` RSmat([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],50);`); print(` RSmat([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,[2],100);`); print(` RSmat([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,[2],1000);`); print(` RSmat([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,[2],2000);`); elif nargs=1 and args[1]=RSmatAno then print(`RSmatAno(C,A,P,x,B,LIMIT1):`): print(`An annotated version of RSmatAno(C,A,P,x,B,LIMIT1) (q.v.). Try:`): print(` RSmatAno([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],50); `): print(` RSmatAno([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],50);`); print(` RSmatAno([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],50);`); print(` RSmatAno([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,[2],100);`); print(` RSmatAno([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,[2],1000);`); print(` RSmatAno([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,[2],2000);`); elif nargs=1 and args[1]=RSmatV then print(`RSmatV(C,A,P,x,B,LIMIT1):`): print(`verbose form of RSmat(C,A,P,x,B,LIMIT1): (q.v.)`): print(` RSmatV([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],50); `): print(` RSmatV([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],50);`); print(` RSmatV([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],50);`); elif nargs=1 and args[1]=RSmatVt then print(`RSmatVt(C,A,P,x,B,LIMIT1,K1):`): print(`terse form of RSmatV(C,A,P,x,B,LIMIT1): (q.v.), without showing the matrix and vector, and computing the frst K1 terms`): print(` RSmatV([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],50,30); `): print(` RSmatV([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],50,30);`); print(` RSmatV([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],50,30);`); elif nargs=1 and args[1]=RSviaMat then print(`RSviaMat(C,A,P,x,B,t,LIMIT1):`): print(`Does the same thing as RS(C,A,P,x,B,t,LIMIT1) (q.v.) buy via RSmat, using matrix inverse. For checking purposes only.`): print(` RSviaMat([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],t,50); `): print(` RSviaMat([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],t,50);`); print(` RSviaMat([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],t,50);`); elif nargs=1 and args[1]=RSe then print(`RSe(C,A,P,x,B,t,MAX):`): print(` Inputs a C-finite sequence C=[INI,REC] (in our convention), `): print(`let d be the order of C (i.e. nops(C[2])). It also inputs a list of pairs`): print(`each of the form [e[0],[e[1],e[2],...e[d]]`): print(`outputs the generating function for the sequence vn(C,A,P,x,B,n) (q.v.). `): print(` It uses the EXPERIMENTAL approach. Try`): print(` RSe([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],t,20); `): print(` RSe([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],t,20);`); print(` RSe([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],t,25);`); elif nargs=1 and args[1]=RSv then print(`RSv(C,A,P,x,B,t,LIMIT1):`): print(`Verbose form of RS(C,A,P,x,B,t,LIMIT1) (q.v.). Try:`): print(` RSv([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],t,50); `): print(` RSv([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],t,50);`); print(` RSv([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],t,50);`); elif nargs=1 and args[1]=vMseq then print(`vMseq(v,M,K): Given a vector of rational functions satisfying f(t)=v+t M f(t), for some vector v and matrix M`): print(`returns the first K+1 terms of the first component of f. Try:`): print(`vMseq([1,2],[[1,2],[2,3]],4);`): elif nargs=1 and args[1]=vn then print(`vn(C,A,P,x,B,n): Same inputs as Gn(C,A,P,x,n) (q.v.) as well as a list of non-negative integers B, outputs`): print(`the number add(mul(coeff(Gn(C,A,P,x,n),x,p+i-1)^B[i],i=1..nops(B)),p=0..infinity). Try:`): print(` vn([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],4); `): print(` vn([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]] ],1,x,[2],4); `): print(` vn([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]], [1,[0,1]] ],1,x,[2],4); `): elif nargs=1 and args[1]=vnL then print(`vnL(C,A,P,x,L,n): Same inputs as Gn(C,A,P,x,n) (q.v.) as well as a list L of pairs (let d=nops(C[1])=nops(C[2]) be the order `): print(`of the underlying C-finite sequence`): print(`[[L[1][1],L[1][2]], [[L[2][1],L[2][2]], ... ..., [where L[i][1] is a non-negative integer and L[i][2] is a list of integers of length d`): print(`It outputs`): print(`the number add(mul(coeff(Gn(C,A,P,x,n),x,p+L[i][1]-add(L[i][2][j]*C[n+j]),j=1..d),i=1..nops(L)),p=0..infinity). Try:`): print(` vnL([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[ [0,[0]], [0,[0]] ], 4); `): print(` vnL([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]] ],1,x,[[0,[0,0]],[0,[0,0]]], 4); `): print(` vnL([[0,1],[1,1]],[[1,[0,0]], [1,[1,0]], [1,[0,1]] ],1,x, [[0,[0,0]],[0,[0,0]] ], 4); `): elif nargs=1 and args[1]=vnSeq then print(`vnSeq(C,A,P,x,B,n): outputs the first n+1 terms of the sequence vn(C,A,P,x,n) (q.v.) starting at n=0. Try:`): print(` vnSeq([[1],[2]],[[1,[0]],[1,[1]],[1,[2]] ],1,x,[2],4); `): print(` vnSeq([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]] ],1,x,[2],4);`); print(` vnSeq([[1,1],[1,1]],[[1,[0,0]], [1,[0,1]],[1,[1,1]] ],1,x,[2],4);`); print(` vnSeq([[1,1,1],[1,1,1]],[[1,[0,0,0]], [1,[0,0,1]],[1,[0,1,1]], [1,[1,1,1]] ],1,x,[2],4);`); print(` vnSeq([[1,1,1,1],[1,1,1,1]],[[1,[0,0,0,0]], [1,[0,0,0,1]],[1,[0,0,1,1]], [1,[0,1,1,1]] ,[1,[1,1,1,1]] ],1,x,[2],4);`); print(` vnSeq([[1,1,1,1,1],[1,1,1,1,1]],[[1,[0,0,0,0,0]], [1,[0,0,0,0,1]],[1,[0,0,0,1,1]], [1,[0,0,1,1,1]] ,[1,[0,1,1,1,1]],[1,[1,1,1,1,1]] ],1,x,[2],4);`); else print(`There is no such thing as`, args): fi: end: ###From Cfinite.txt #SeqFromRec(S,N): Inputs S=[INI,ope] #where INI is the list of initial conditions, a ope a list of #size L, say, and a recurrence operator ope, codes a list of #size L, finds the first N0 terms of the sequence satisfying #the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L). #For example, for the first 20 Fibonacci numbers, #try: SeqFromRec([[1,1],[1,1]],20); SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope: INI:=S[1]:ope:=S[2]: if not type(INI,list) or not type(ope,list) then print(`The first two arguments must be lists `): RETURN(FAIL): fi: L:=nops(INI): if nops(ope)<>L then print(`The first two arguments must be lists of the same size`): RETURN(FAIL): fi: if not type(N,integer) then print(`The third argument must be an integer`, L): RETURN(FAIL): fi: if N=`, 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: #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: ###End From Cfinite.txt #Atom1(C,A,x,i): Inputs a C-finite sequence C=[INI,REC] (in our convention), a sequence A=[a1,...,ah] (of numbers or symbols), Atom1:=proc(C,A,x,i) local gu,j,d,i1,gadol: d:=nops(C[1]): if nops(C[2])<> d then RETURN(FAIL): fi: if not (type(A,list) and {seq(nops(A[i1]),i1=1..nops(A))}={2} and {seq(nops(A[i1][2]),i1=1..nops(A))}={d}) then RETURN(FAIL): fi: gu:=SeqFromRec(C,i+d+2): add(A[i1][1]*x^add(A[i1][2][j]*gu[i+j],j=1..d),i1=1..nops(A)): end: #Gn(C,A,P,x,n): Inputs a C-finite sequence C=[INI,REC] (in our convention), a sequence A=[a1,...,ah] (of numbers or symbols), #a variable x, and a non-negative integer n, and a polynomial P of x, outputs #P(x)*Prod(Atom1(C,A,x,i),i=0..n); Gn:=proc(C,A,P,x,n) local i: expand(P*mul(Atom1(C,A,x,i),i=0..n-1)): end: #vn(C,A,P,x,B,n): Same inputs as Gn(C,A,P,x,n) (q.v.) as well as a list of non-negative integers B, outputs #the number add(mul(coeff(Gn(C,A,P,x,n),x,p+i-1)^B[i],i=1..nops(B)),p=0..infinity). Try: vn:=proc(C,A,P,x,B,n) local gu,k,m,i: option remember: m:=nops(B): gu:=Gn(C,A,P,x,n): add(mul(coeff(gu,x,k+i-1)^B[i],i=1..m),k=0..degree(gu,x)): end: #vnSeq(C,A,P,x,B,N): Same inputs as vn(C,A,P,x,n) (q.v.) but outputs the first N+1 terms of the sequence vn(C,A,P,x,n) starting at n=0. Try: #vnSeq([[1],[2]],[1],1,x,[2],4); #vnSeq([[1],[2]],[1,1],1,x,[2],4); #vnSeq([[1,2],[1,1]],[0,1],1,x,[2],4); vnSeq:=proc(C,A,P,x,B,N) local n: [seq(vn(C,A,P,x,B,n),n=0..N)]: end: #RSe(C,A,P,x,B,t,MAX): inputs a C-finite seqeunce C=[INI,REC], a list A of numbers (or symbols), a polynomial P in the variable x, #a list of non-negative integers B and a variable t, and outputs the (ordinary) generating function, in t, of the sequence #vn(C,A,P,x,B,n) (q.v.). Try RSe:=proc(C,A,P,x,B,t,MAX) local gu,mu,K: if MAX<10 then print(MAX, `should have been at least 10`): RETURN(FAIL): fi: for K from 10 by 5 to MAX do gu:=vnSeq(C,A,P,x,B,K): mu:=GuessRec(gu): if mu<>FAIL then RETURN(factor(CtoR(mu,t))): fi: od: print(MAX, `is not large enough, make it bigger. `): FAIL: end: #vnL(C,A,P,x,L,n): Same inputs as Gn(C,A,P,x,n) (q.v.) as well as a list L of pairs (let d=nops(C[1])=nops(C[2]) be the order #of the underlying C-finite sequence #[[L[1][1],L[1][2]], [[L[2][1],L[2][2]], ... ..., [where L[i][1] is a positive integer and L[i][2] is a list of integers of length d #It outputs #the number add(mul(coeff(Gn(C,A,P,x,n),x,p+L[i][1]-add(L[i][2][j]*C[n+j]),j=1..d),i=1..nops(L)),p=0..infinity). Try: vnL:=proc(C,A,P,x,L,n) local d,gu,p,i,S,gadol,j,su,pro: option remember: d:=nops(C[2]): gu:=Gn(C,A,P,x,n): S:=SeqFromRec(C,n+d+2): gadol:=max(seq(add(L[i][2][j]*S[n+j+1],j=1..d),i=1..nops(L))): su:=0: for p from 0 to degree(gu,x)+gadol do pro:=1: for i from 1 to nops(L) do pro:=pro*coeff(gu,x,p+L[i][1]- add( L[i][2][j]*S[n+j], j=1..d)): od: su:=su+pro: od: su: end: #ROp(L,C): Given a list L and a C-finite sequence C=[INI,REC], (let d=nops(REC)) where nops(L)=d finds the vector L1 (of the same length, d) #such that L[1]*C[n+1]+ ...+ L[d]*C[n+d]=L'[1]*C[n]+...+ L'[d]*C[n+d-1]. Try: #ROp([2,3],[[1,1],[1,1]]); ROp:=proc(L,C) local i,d,L1,S,n1: if not (nops(C)=2 and nops(C[1])=nops(C[2]) and nops(C[2])=nops(L)) then print(`bad input`): RETURN(FAIL): fi: d:=nops(L): L1:=[C[2][d]*L[d],seq(L[i]+C[2][d-i]*L[d],i=1..d-1)]: S:=SeqFromRec(C,30): if {seq(add(L[i]*S[n1+i+1],i=1..d)-add(L1[i]*S[n1+i],i=1..d),n1=1..nops(S)-d-1)}<>{0} then print(`Something is wrong`): RETURN(FAIL): fi: L1: end: #CanF(S): The canonical form of the list of lists L. Try #CanF([[1,1,1],[2,1,1],[3,1,1]]); CanF:=proc(L) local L1, d,i,j,katan: d:=nops(L[1][2]): for j from 1 to d do katan[j]:=min(seq(L[i][2][j],i=1..nops(L))): od: L1:=[seq([L[i][1],[seq(L[i][2][j]-katan[j],j=1..d)]],i=1..nops(L))]: sort(L1): end: #Eq(C,A,P,x,L,Z,t): An equation for the generating function of vnL(C,A,P,x,L,n) denoted by Z[L] in terms #of other Z[L']-s, followed by its children, Try: #Eq([[1,1],[1,1]], [[1,[0,0]],[1,[1,0]], [1,[0,1]] ,1,x,[ [0,[0,0]], [0,[0,0]] ],Z,t); Eq:=proc(C,A,P,x,L,Z,t) local d,gu,Y,i,j,eq,Lnew, gu1,coe,L1,var,i1,n1,kak: d:=nops(C[2]): Lnew:=[seq([L[i][1],ROp(L[i][2],C )],i=1..nops(L))]: gu:=mul(x[i]^Lnew[i][1],i=1..nops(Lnew)): gu:=gu*mul(mul(Y[i1,j]^Lnew[i1][2][j],j=1..d),i1=1..nops(Lnew)): for i from 1 to nops(L) do kak:=add(A[i1][1]*mul(Y[i,j]^A[i1][2][j],j=1..d),i1=1..nops(A)): gu:=expand(gu*kak): od: eq:=Z[L]-vnL(C,A,P,x,L,0): var:={}: for i from 1 to nops(gu) do gu1:=op(i,gu): coe:=gu1/mul(x[i1]^degree(gu1,x[i1]),i1=1..nops(L)): coe:=normal(coe/mul(mul(Y[i1,j]^degree(gu1,Y[i1,j]),j=1..d),i1=1..nops(L))): L1:=[ seq( [degree(gu1,x[i1]),[seq(degree(gu1,Y[i1,j]),j=1..d) ]] , i1=1..nops(L)) ]: L1:=CanF(L1): if {seq(vnL(C,A,P,x,L1,n1),n1=0..8)}<>{0} then eq:=eq-coe*t*Z[L1]: var:=var union {L1}: fi: od: eq,var: end: #RS(C,A,P,x,B,t,LIMIT1): inputs a C-finite seqeunce C=[INI,REC], a list A of numbers (or symbols), a polynomial P in the variable x, #a list of non-negative integers B and a variable t, and outputs the (ordinary) generating function, in t, of the sequence #vn(C,A,P,x,B,n) (q.v.). Try #RS([[1],[2]],[1],1,x,[2],t,400); #RS([[1],[2]],[1,1],1,x,[2],t,400); #RS([[1,2],[1,1]],[1],1,x,[2],t,400); #RS([[1,2],[1,1]],[1,1],1,x,[2],t,400); RS:=proc(C,A,P,x,B,t,LIMIT1) local i,d,eq,var,L,L1,StillToDo, AlreadyDone,eq1,Z,L2,lu: option remember: d:=nops(C[2]): L:=[seq([i-1,[0$d]]$B[i],i=1..nops(B))]: StillToDo:={L}: eq:={}: var:={}: AlreadyDone:={}: while StillToDo<>{} do L1:=StillToDo[1]: eq1:=Eq(C,A,P,x,L1,Z,t) : eq:=eq union {eq1[1]}: if nops(eq)>LIMIT1 then print(`Exceeded the limit of number of equations`): RETURN(FAIL): fi: AlreadyDone:=AlreadyDone union {L1}: StillToDo:=StillToDo union eq1[2] minus AlreadyDone: var:=var union {seq(Z[L2],L2 in eq1[2])}: od: var:=solve(eq,var): lu:=factor(subs(var,Z[L])): #if vnSeq(C,A,P,x,B,15)<>[seq(coeff(taylor(lu,t=0,16),t,i1),i1=0..15)] then # print(`Something is wrong`): # RETURN(FAIL): #fi: lu: end: #RSv(C,A,P,x,B,t,LIMIT1): inputs a C-finite seqeunce C=[INI,REC], a list A of numbers (or symbols), a polynomial P in the variable x, #a list of non-negative integers B and a variable t, and outputs the (ordinary) generating function, in t, of the sequence #vn(C,A,P,x,B,n) (q.v.). Try #RSv([[1],[2]],[1],1,x,[2],t,400); #RSv([[1],[2]],[1,1],1,x,[2],t,400); #RSv([[1,2],[1,1]],[1],1,x,[2],t,400); #RSv([[1,2],[1,1]],[1,1],1,x,[2],t,400); RSv:=proc(C,A,P,x,B,t,LIMIT1) local gu,f1,Z,i,j,d,n,t0,H: t0:=time(): gu:=RS(C,A,P,x,B,t,LIMIT1): if gu=FAIL then RETURN(FAIL): fi: d:=nops(C[2]): f1:=CtoR(C,t): print(``): print(`---------------------------------`): print(``): print(` Let Z[n] be the integer sequence whose generating function is`): print(``): print(Sum(Z[j]*t^j,j=0..infinity)=f1): print(``): print(`Let`): print(``): print(F[n](x)=P*Product(add(A[i1][1]*x^add(A[i1][2][j]*Z[i+j-1],j=1..d),i1=1..nops(A)),i=0..n-1)): print(``): print(`Write: `): print(``): print(F[n](x)=Sum(a(n,i)*x^i,i=0..infinity)): print(``): print(`Let :`): print(``): print(H(n)=Sum(mul(a(n,k+i-1)^B[i],i=1..nops(B)),k=0..infinity)): print(``): print(`Then `): print(``): print(Sum(H(n)*t^n,n=0..infinity)=gu): print(``): print(`and in Maple notation`): print(``): lprint(gu): print(``): print(`For the sake of the OEIS, here are the first 30 terms, starting at n=0`): print(``): gu:=taylor(gu,t=0,31): print(seq(coeff(gu,t,i),i=0..30)): print(``): print(`-----------------------------`): print(``): print(`This took`, time()-t0, `seconds. `): end: #Paper1(C,A,P,x,t,LIMIT1,N): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t #outputs the paper with N theorems RSv((C,A,P,x,[i],t,LIMIT1) for i from 1 to N. Try: Paper1:=proc(C,A,P,x,t,LIMIT1,N) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanley-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N do print(`Theorem Number`, i): print(``): if RS (C,A,P,x,[i],t,LIMIT1)<>FAIL then RSv(C,A,P,x,[i],t,LIMIT1,N): fi: print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: #Paper1mat(C,A,P,x,LIMIT1,N,K1): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t #outputs the paper with N theorems RSv((C,A,P,x,[i],t,LIMIT1) for i from 1 to N. Try: Paper1mat:=proc(C,A,P,x,LIMIT1,N,K1) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanley-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N do print(`Theorem Number`, i): print(``): if RSmat(C,A,P,x,[i],LIMIT1)<>FAIL then RSmatVt(C,A,P,x,[i],LIMIT1,K1): fi: print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: #Paper2(C,A,P,x,t,LIMIT1,N): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t #outputs the paper with N theorems RSv((C,A,P,x,[i],t,LIMIT1) for i from 1 to N. Try: Paper2:=proc(C,A,P,x,t,LIMIT1,N) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanley-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N do print(`Theorem Number`, i): print(``): if RS(C,A,P,x,[i],t,LIMIT1)<>FAIL then RSv(C,A,P,x,[1$i],t,LIMIT1,N): fi: print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: #Paper2mat(C,A,P,x,LIMIT1,N,K1): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t #outputs the paper with N theorems RSv((C,A,P,x,[i],t,LIMIT1) for i from 1 to N. Try: Paper2mat:=proc(C,A,P,x,LIMIT1,N,K1) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanley-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N do print(`Theorem Number`, i): print(``): if RSmat(C,A,P,x,[i],LIMIT1)<>FAIL then RSmatVt(C,A,P,x,[1$i],LIMIT1,K1): fi: print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: #EqM(C,A,P,x,L,Z): An equation for the generating function of vnL(C,A,P,x,L,n) denoted by Z[L] in terms #of other Z[L']-s, followed by its children, Try: #EqM([[1,1],[1,1]], [[1,[0,0]],[1,[1,0]], [1,[0,1]] ,1,x,[ [0,[0,0]], [0,[0,0]] ],Z,t); EqM:=proc(C,A,P,x,L,Z) local d,gu,Y,i,j,eq,Lnew, gu1,coe,L1,var,i1,n1,kak,eqM: d:=nops(C[2]): Lnew:=[seq([L[i][1],ROp(L[i][2],C )],i=1..nops(L))]: gu:=mul(x[i]^Lnew[i][1],i=1..nops(Lnew)): gu:=gu*mul(mul(Y[i1,j]^Lnew[i1][2][j],j=1..d),i1=1..nops(Lnew)): for i from 1 to nops(L) do kak:=add(A[i1][1]*mul(Y[i,j]^A[i1][2][j],j=1..d),i1=1..nops(A)): gu:=expand(gu*kak): od: eqM:=[Z[L],vnL(C,A,P,x,L,0)]: eq:=0: var:={}: for i from 1 to nops(gu) do gu1:=op(i,gu): coe:=gu1/mul(x[i1]^degree(gu1,x[i1]),i1=1..nops(L)): coe:=normal(coe/mul(mul(Y[i1,j]^degree(gu1,Y[i1,j]),j=1..d),i1=1..nops(L))): if not type(coe,numeric) then RETURN(FAIL): fi: L1:=[ seq( [degree(gu1,x[i1]),[seq(degree(gu1,Y[i1,j]),j=1..d) ]] , i1=1..nops(L)) ]: L1:=CanF(L1): if {seq(vnL(C,A,P,x,L1,n1),n1=0..8)}<>{0} then eq:=eq+coe*Z[L1]: var:=var union {L1}: fi: od: [op(eqM),eq] ,var: end: #RSmat(C,A,P,x,B,LIMIT1): inputs a C-finite seqeunce C=[INI,REC], a list A of numbers (or symbols), a polynomial P in the variable x, #a list of non-negative integers B and a variable t, and outputs the (ordinary) generating function, in t, of the sequence #vn(C,A,P,x,B,n) (q.v.). outputs a square matrix M, and a colum with intger entries, of such that the desired output is (I-t*M)^(-1)[1,0,0,0....]. #LIMIT1 is the max size of the matrix #Try #RSmat([[1],[2]],[1],1,x,[2],400); #RSmat([[1],[2]],[1,1],1,x,[2],400); #RSmat([[1,2],[1,1]],[1],1,x,[2],400); #RSmat([[1,2],[1,1]],[1,1],1,x,[2],400); RSmat:=proc(C,A,P,x,B,LIMIT1) local i,d,eq,var,L,L1,StillToDo, AlreadyDone,eq1,Z,L2,lu,v,TA, MAT,j,i1: option remember: d:=nops(C[2]): L:=[seq([i-1,[0$d]]$B[i],i=1..nops(B))]: StillToDo:={L}: eq:={}: var:={}: AlreadyDone:={}: while StillToDo<>{} do L1:=StillToDo[1]: eq1:=EqM(C,A,P,x,L1,Z) : eq:=eq union {eq1[1]}: if nops(eq)>LIMIT1 then print(`Exceeded the limit of number of equations`): RETURN(FAIL): fi: AlreadyDone:=AlreadyDone union {L1}: StillToDo:=StillToDo union eq1[2] minus AlreadyDone: var:=var union {seq(Z[L2],L2 in eq1[2])}: od: var:=[Z[L],op(convert(var minus {Z[L]},list))]: for i from 1 to nops(var) do TA[var[i]]:=i: od: for eq1 in eq do v[TA[eq1[1]]]:=eq1[2]: for j from 1 to nops(var) do MAT[TA[eq1[1]],TA[var[j]]]:=coeff(eq1[3],var[j],1): od: od: [[seq(v[i1],i1=1..nops(var))], [seq([seq(MAT[i,j],j=1..nops(var))],i=1..nops(var))]]: end: with(linalg): RSviaMat:=proc(C,A,P,x,B,t,LIMIT1) local gu,v,M,M1,i,j,v1: gu:=RSmat(C,A,P,x,B,LIMIT1): v:=gu[1]: M:=gu[2]: for i from 1 to nops(v) do for j from 1 to nops(v) do if i=j then M1[i,j]:=1-t*M[i,j]: else M1[i,j]:=-t*M[i,j]: fi: od: od: M1:=matrix([seq([seq(M1[i,j],j=1..nops(v))],i=1..nops(v))]): M1:=inverse(M1): v1:=matrix([seq([v[i1]],i1=1..nops(v))]): factor(multiply(M1,v1)[1,1]): end: #vMseq(v,M,K): Given a vector of rational functions satisfying f(t)=v+t M f(t), for some vector v and matrix M #returns the first K+1 terms of the first component of f. Try: #vMseq([1,2],[[1,2],[2,3]],4); vMseq:=proc(v,M,K) local f,i,i1: if not (type(v,list) and type(M,list) and nops(v)=nops(M) and {seq(nops(M[i]),i=1..nops(M))}={nops(M)} ) then print(`bad input`): RETURN(FAIL): fi: for i from 1 to nops(v) do f[i]:=v[1]: od: for i1 from 1 to K do for i from 1 to nops(v) do f[i]:=expand(v[i]+t*add(M[i][j]*f[j],j=1..nops(v))): od: od: [seq(coeff(f[1],t,i),i=0..K)]: end: #RSmatV((C,A,P,x,B,LIMIT1): Verbose form of RSmat(C,A,P,x,B,LIMIT1) (q.v.) RSmatV:=proc(C,A,P,x,B,LIMIT1) local f1,t,gu,Z,F,H,n,i1,t0,lu,i,j,d,k: t0:=time(): gu:=RSmat(C,A,P,x,B,LIMIT1): if gu=FAIL then RETURN(FAIL): fi: d:=nops(C[2]): f1:=CtoR(C,t): print(``): print(`---------------------------------`): print(``): print(` Let Z[n] be the integer sequence whose generating function is`): print(``): print(Sum(Z[j]*t^j,j=0..infinity)=f1): print(``): print(` Let `): print(``): print(F[n](x)=P*Product(add(A[i1][1]*x^add(A[i1][2][j]*Z[i+j-1],j=1..d),i1=1..nops(A)),i=0..n-1)): print(``): print(`Write: `): print(``): print(F[n](x)=Sum(a(n,i)*x^i,i=0..infinity)): print(``): print(` Let : `): print(``): print(H(n)=Sum(mul(a(n,k+i-1)^B[i],i=1..nops(B)),k=0..infinity)): print(``): print(` Then `): print(``): print(Sum(H(n)*t^n,n=0..infinity)): print(`equals (I-Mt)^(-1) v [1]`): print(``): print(`where M is the matrix`): print(``): print(matrix(gu[2])): print(``): print(`and v is the vector `): print(``): print(gu[1]): print(``): lu:=vMseq(op(gu),30): print(`For the sake of the OEIS, here are the first 31 terms, starting at n=0`): print(``): print(lu): print(``): print(`-----------------------------`): print(``): print(`This took`, time()-t0, `seconds. `): end: #RSmatVt((C,A,P,x,B,LIMIT1): Verbose form of RSmat(C,A,P,x,B,LIMIT1) (q.v.) RSmatVt:=proc(C,A,P,x,B,LIMIT1,K1) local f1,t,gu,Z,F,H,n,i1,t0,lu,i,j,d,k: t0:=time(): gu:=RSmat(C,A,P,x,B,LIMIT1): if gu=FAIL then RETURN(FAIL): fi: d:=nops(C[2]): f1:=CtoR(C,t): print(``): print(`---------------------------------`): print(``): print(` Let Z[n] be the integer sequence whose generating function is`): print(``): print(Sum(Z[j]*t^j,j=0..infinity)=f1): print(``): print(` Let `): print(``): print(F[n](x)=P*Product(add(A[i1][1]*x^add(A[i1][2][j]*Z[i+j-1],j=1..d),i1=1..nops(A)),i=0..n-1)): print(``): print(`Write: `): print(``): print(F[n](x)=Sum(a(n,i)*x^i,i=0..infinity)): print(``): print(` Let : `): print(``): print(H(n)=Sum(mul(a(n,k+i-1)^B[i],i=1..nops(B)),k=0..infinity)): print(``): print(` Then `): print(``): print(Sum(H(n)*t^n,n=0..infinity)): print(`equals (I-Mt)^(-1) v [1]`): print(``): print(`where M is a certain square matrix of dimension`, nops(gu[1])): print(``): print(`and v is a certain vector of length`, nops(gu[2])): print(``): print(`that are too big to display. At any rate we can use them to find the first`, K1+1, `terms starting at n=0, for the sake of the OEIS.`): print(``): lu:=vMseq(op(gu),K1): print(``): print(lu): print(``): print(`-----------------------------`): print(``): print(`This took`, time()-t0, `seconds. `): end: #EqMv(C,A,P,x,L,Z): An equation for the generating function of vnL(C,A,P,x,L,n) denoted by Z[L] in terms #of other Z[L']-s, followed by its children, Try: #EqM([[1,1],[1,1]], [[1,[0,0]],[1,[1,0]], [1,[0,1]] ,1,x,[ [0,[0,0]], [0,[0,0]] ],Z,t); EqMv:=proc(C,A,P,x,L,Z) local d,gu,Y,i,j,eq,Lnew, gu1,coe,L1,var,i1,n1,kak,eqM,t0: d:=nops(C[2]): print(`the time is`, t0): print(`We are finding an equation for state`, L): Lnew:=[seq([L[i][1],ROp(L[i][2],C )],i=1..nops(L))]: print(`adjusting, since we are going down we have sate`, Lnew): gu:=mul(x[i]^Lnew[i][1],i=1..nops(Lnew)): gu:=gu*mul(mul(Y[i1,j]^Lnew[i1][2][j],j=1..d),i1=1..nops(Lnew)): print(`My are trying to find all the children by finding the childern polynomial`): print(gu): for i from 1 to nops(L) do kak:=add(A[i1][1]*mul(Y[i,j]^A[i1][2][j],j=1..d),i1=1..nops(A)): gu:=expand(gu*kak): od: print(`The initial symbolic evolution expression is`): print(``): print(gu): print(``): eqM:=[Z[L],vnL(C,A,P,x,L,0)]: print(`We are forming the equation`): eq:=0: var:={}: for i from 1 to nops(gu) do gu1:=op(i,gu): coe:=gu1/mul(x[i1]^degree(gu1,x[i1]),i1=1..nops(L)): coe:=normal(coe/mul(mul(Y[i1,j]^degree(gu1,Y[i1,j]),j=1..d),i1=1..nops(L))): if not type(coe,numeric) then RETURN(FAIL): fi: L1:=[ seq( [degree(gu1,x[i1]),[seq(degree(gu1,Y[i1,j]),j=1..d) ]] , i1=1..nops(L)) ]: print(`One of the chidren is`): print(L1): L1:=CanF(L1): print(`whose canonical form is`): print(``): print(L1): if {seq(vnL(C,A,P,x,L1,n1),n1=0..8)}<>{0} then print(`It is not 0 `): eq:=eq+coe*Z[L1]: var:=var union {L1}: else print(`It gives a 0 , so we discard it`): fi: od: [op(eqM),eq] ,var: end: #RSmatAno(C,A,P,x,B,LIMIT1): An annotated version of RSmat #Try #RSmatAno([[1],[2]],[1],1,x,[2],400); #RSmatAno([[1],[2]],[1,1],1,x,[2],400); #RSmatAno([[1,2],[1,1]],[1],1,x,[2],400); #RSmatAno([[1,2],[1,1]],[1,1],1,x,[2],400); RSmatAno:=proc(C,A,P,x,B,LIMIT1) local i,d,eq,var,L,L1,StillToDo, AlreadyDone,eq1,Z,L2,v,TA, MAT,j: d:=nops(C[2]): L:=[seq([i-1,[0$d]]$B[i],i=1..nops(B))]: StillToDo:={L}: eq:={}: var:={}: AlreadyDone:={}: while StillToDo<>{} do L1:=StillToDo[1]: print(`We are doing state`, L1): eq1:=EqMv(C,A,P,x,L1,Z) : eq:=eq union {eq1[1]}: print(`So far we already did`, nops(AlreadyDone), `states `): print(``): print(`We still need to do`, nops(StillToDo), `of them (and more may come up) `): print(``): print(`Right now there are`, nops(eq), `equations `): print(``): if nops(eq)>LIMIT1 then print(`Exceeded the limit of number of equations`): print(`The set of states that we did is`): print(``): print(AlreadyDone): print(``): print(`We still need to do`, nops(StillToDo), `states `): RETURN(FAIL): fi: AlreadyDone:=AlreadyDone union {L1}: StillToDo:=StillToDo union eq1[2] minus AlreadyDone: var:=var union {seq(Z[L2],L2 in eq1[2])}: od: var:=[Z[L],op(convert(var minus {Z[L]},list))]: for i from 1 to nops(var) do TA[var[i]]:=i: od: for eq1 in eq do v[TA[eq1[1]]]:=eq1[2]: for j from 1 to nops(var) do MAT[TA[eq1[1]],TA[var[j]]]:=coeff(eq1[3],var[j],1): od: od: print(`There are`, nops(AlreadyDone), `states at the end, here there are`): print(``): print(AlreadyDone): print(``): print(``): #print(`The final output is`): #ka:=[[seq(v[i1],i1=1..nops(var))], [seq([seq(MAT[i,j],j=1..nops(var))],i=1..nops(var))]]: #print(ka): #ka: end: #Jrk(r,k,t,x,LIMIT1): The rational generating function J_r^{(k)}(t,x) defined on p.13 in Richard Stanley's paper # "Theorems and Conjectures on Some Rational Generating functions" #N is the cap on the number of equations. Try: #Jrk(2,2,t,x,100); Jrk:=proc(r,k,t,x,LIMIT1) local X,lu,lu1,lu2,i: option remember: lu:= RS([[1$k],[1$k]],[[1,[0$k]], [t,[1$k]] ],1,X,[r],x,LIMIT1); if lu=FAIL then RETURN(FAIL): fi: lu1:=numer(lu): lu2:=denom(lu): lu1:=add(factor(coeff(lu1,x,i))*x^i,i=0..degree(lu1,x)): lu2:=add(factor(coeff(lu2,x,i))*x^i,i=0..degree(lu2,x)): lu1/lu2: end: #JrkS(r,k,x,LIMIT1,C): The CONJECTURED rational generating function J_r^{(k)}(1,x) defined on p.13 in Richard Stanley's paper #in terms of a SYMBOLIC k, assuming that indeed the coefficients do do not depend on k. #It checks that it true for k=2,3,..,C So if C=2 then it does not check it # "Theorems and Conjectures on Some Rational Generating functions" #N is the cap on the number of equations. Try: #JrkS(2,k,x,100,3); JrkS:=proc(r,k,x,LIMIT1,C) local lu,lu1,lu2,lu1S,lu2S,i,luS,k1: if not type(k, symbol) then print(k, `must be a symbol`): RETURN(FAIL): fi: lu:=Jrk(r,2,1,x,LIMIT1): if lu=FAIL then RETURN(FAIL): fi: lu1:=numer(lu): lu2:=denom(lu): lu1S:=0: for i from 0 to trunc(degree(lu1,x)/2)+1 do lu1S:=lu1S+coeff(lu1,x,2*i)*x^(k*i)+coeff(lu1,x,2*i+1)*x^(k*i+1): od: lu2S:=0: for i from 0 to trunc(degree(lu2,x)/2)+1 do lu2S:=lu2S+coeff(lu2,x,2*i)*x^(k*i)+coeff(lu2,x,2*i+1)*x^(k*i+1): od: luS:=lu1S/lu2S: luS: for k1 from 2 to C do if normal(subs(k=k1,luS))<>normal(Jrk(r,k1,1,x,LIMIT1)) then print(`when k1 is`, k1): print(``): lprint([normal(subs(k=k1,luS)),normal(Jrk(r,k1,1,x,LIMIT1))]): print(luS, ` did not work out `): RETURN(FAIL): fi: od: luS: end: #J2kC(k,t,x): Stanley's CONECTURED rational generating function J_2^{(k)}(t,x) defined on p.19 in Richard Stanley's paper #CORRECTING his conjecture. Try: #J2kC(3,t,x); J2kC:=proc(k,t,x) :(1-t^(k-1)*(1+t^2)*x^k)/(1-(1+t^2)*x-t^(k-1)*(1+t^2)*x^k+t^(k-1)*(1+t^4)*x^(k+1)):end: #J3kC(k,t,x): The CONECTURED rational generating function J_3^{(k)}(t,x) defined on p.13 in Richard Stanley's paper #CORRECTING his conjecture. Try: #J3kC(3,t,x); J3kC:=proc(k,t,x) : (1-t^(k-1)*(1+t^(k-1))*(1+t^3)*x^k+t^(3*(k-1))*(1-t^3)^2*x^(2*k))/ ( 1-(1+t^3)*x-t^(k-1)*(1+t^(k-1))*(1+t^3)*x^k+t^(k-1)*(1+t^(k-1))*(1-t^3+t^6)*x^(k+1)+t^(3*(k-1))*(1-t^3)^2*x^(2*k) -t^(3*(k-1))*(1-t^3)*(1-t^6)*x^(2*k+1) ): end: