###################################################################### ## StanleyStern.txt Save this file as StanleyStern.txt to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read StanlyStern.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 StanleyStern.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: ezraS:=proc() if args=NULL then print(`The STORY procedures are: `): print(` Paper1, Paper1e, Paper2, Paper2e, RSano `): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(`The SUPPORTING procedures are: CheckRS, CF, Eq, Fpqbnx , SeqUpqbAnx, SeqUnL, UpqbAnx, UnL `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` StanleyStern.txt: A Maple package for computing generating functions dear to Richard Stanely`): print(`The MAIN procedures are: `): print(` RS, RSe, RSv`): elif nargs=1 and args[1]=CF then print(`CF(L): The canonical form of the list of pairs L. Try`): print(`CF([[1,3],[1,4],[2,5]]);`): elif nargs=1 and args[1]=CheckRS then print(`CheckRS(p,q,b,A,x,N): checks procedure RS(p,q,b,A,x,t) up to N terms. Try:`): print(`CheckRS(1+x+x^2,1,2,[3],x,15);`): elif nargs=1 and args[1]=Eq then print(`Eq(p,q,b,x,L,t,Z): The eqation expressing the generating function of UnL(p,q,b,n,L) called Z[L], in terms of other`): print(`UnL(p,q,b,n,L1)'s called Z[L1]. It also returns the set of L1's that show up. Try:`): print(`Eq(1+x+x^2,1,2,x,[[0,0],[0,0]],t,Z);`): elif nargs=1 and args[1]=Fpqbnx then print(`Fpqbnx(p,q,b,n,x): inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x)`): print(`outputs expand(q*mul(subs(x=x^(b^i),i=0..n-1)), defined on p.7 of Richard Stanley's paper`): print(`"Some Linear Recurrences Motivated by Stern's Diatomic arrays arXiv:1901.04647v1 [math.CO] 15 Jan 2019. Try (for the original case (with n=5) of Stern)`): print(` Fpqbnx(1+x+x^2,1,2,5,x);`): elif nargs=1 and args[1]=Paper1 then print(`Paper1(p,q,b,x,t,N): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t`): print(`outputs the paper with N theorems RSv(p,q,b,[i],x,t) for i from 1 to N. Try:`): print(`Paper1(1+x+x^2,1,2,x,t,5):`): elif nargs=1 and args[1]=Paper1E then print(`Paper1e(p,q,b,x,t,N1,N2): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t`): print(`outputs the paper with N1 theorems RSvE(p,q,b,[i],N2,x,t) for i from 1 to N1. Try:`): print(`Paper1E(1+x+x^2,1,2,x,t,5,20):`): elif nargs=1 and args[1]=Paper2 then print(`Paper2(p,q,b,x,t,N): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t`): print(`outputs the paper with N theorems RSv(p,q,b,[1$i],x,t) for i from 1 to N1. Try:`): print(`Paper2(1+x+x^2,1,2,x,t,5):`): elif nargs=1 and args[1]=Paper2e then print(`Paper2e(p,q,b,x,t,N1,N2): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t`): print(`outputs the paper with N theorems RSv(p,q,b,[1$i],N2,x,t) for i from 1 to N. Try:`): print(`Paper2e(1+x+x^2,1,2,x,t,5,20):`): elif nargs=1 and args[1]=RS then print(`RS(p,q,b,A,x,t): `): print(`Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x)`): print(`and a list of non-neg. integers A , and a variable name t, outputs (if possible) the rational function, in t,`): print(`whose coefficient of t^n (in its Taylor expansion) is UpqbAnx(p,q,b,A,n,x) (q.v.) It does it`): print(`by using symbolic dynamical programming. Try:`): print(`RS(1+x+x^2,1,2,[3],x,t);`): elif nargs=1 and args[1]=RSano then print(`RSano(p,q,b,A,x,t): `): print(` an annotated version of RS(p,q,b,A,x,t) (q.v.). Try: `): print(`RSano(1+x+x^2,1,2,[3],x,t);`): elif nargs=1 and args[1]=RSe then print(`RSe(p,q,b,A,N,x,t): `): print(`Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x)`): print(`and a list of non-neg. integers A , and a variable name t, outputs (if possible) the rational function, in t,`): print(`whose coefficient of t^n (in its Taylor expansion) is UpqbAnx(p,q,b,A,n,x) (q.v.) It does it`): print(`by guessing up to N terms. If N terms do not suffice, it return FAIL. Try:`): print(`RSe(1+x+x^2,1,2,[3],10,x,t);`): elif nargs=1 and args[1]=SeqUnL then print(`SeqUnL(p,q,b,N,x,L): The list consisting of the first N+1 terms, starting at n=0 of UnL(p,q,b,n,x,L) (q.v.) Try:`): print(`SeqUnL(1+x+x^2,1,2,10,x,[[0,0],[0,0], [0,0]]);`): elif nargs=1 and args[1]=SeqUpqbAnx then print(`SeqUpqbAnx(p,q,b,A,N,x): Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer N (and a variable name x)`): print(`and a list of non-neg. integers A (of length m, say) outputs the first N+1 terms (starting at n=0) of`): print(`UpqbAnx(1+x+x^2,1,2,[3],n,x); (q.v.) Try:`): print(`SeqUpqbAnx(1+x+x^2,1,2,[3],10,x);`): elif nargs=1 and args[1]=UnL then print(`UnL(p,q,b,n,x,L): `): print(`Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x)`): print(`and a list of pairs of non-neg. integers (of length m, say), L, let `): print(`c_k(n) be the coefficient of Fpqbnx(p,q,b,n,x) (q.v.). It outputs`): print(`Sum(c_{k+A[1][2]-A[1][1]*b^n}(n)*c_{k+A[2][2]-A[2][1]*b^n}(n)*...c_{k+A[m][2]-A[m][1]*b^n}(n),k=0..infinity). Try:`): print(`UnL(1+x+x^2,1,2,10,x,[[0,0],[0,0], [0,0]]);`): elif nargs=1 and args[1]=UpqbAnx then print(`UpqbAnx(p,q,b,A,n,x): Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x)`): print(`and a list of non-neg. integers A (of length m, say), let`): print(`c_k(n) be the coefficient of Fpqbnx(p,q,b,n,x) (q.v.). It outputs`): print(`Sum(c_k(n)^A[1]*c_{k+1}(n)^A[2]*...*c_{k+m-1}^A[m],k=0..infinity): Try:`): print(`UpqbAnx(1+x+x^2,1,2,[3],10,x);`): 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 #Fpqbnx(p,q,b,n,x): Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x) #outputs expand(q*mul(subs(x=x^(b^i),i=0..n-1)), defined on p.7 of Richard Stanley's paper #"Some Linear Recurrences Motivated by Stern's Diatomic arrays arXiv:1901.04647v1 [math.CO] 15 Jan 2019 Fpqbnx:=proc(p,q,b,n,x) local i: option remember: sort(expand(q*mul(subs(x=x^(b^i),p),i=0..n-1))): end: #UpqbAnx(p,q,b,A,n,x): #Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x) #and a list of non-neg. integers A (of length m, say) , let #c_k(n) be the coefficient of Fpqbnx(p,q,b,n,x) (q.v.). It outputs #Sum(c_k(n)^A[1]*c_{k+1}(n)^A[2]*...*c_{k+m-1}^A[m],k=0..infinity): Try: #UpqbAnx(1+x+x^2,1,2,[3],10,x); UpqbAnx:=proc(p,q,b,A,n,x) local gu,k,m,i: option remember: m:=nops(A): gu:=Fpqbnx(p,q,b,n,x): add(mul(coeff(gu,x,k+i-1)^A[i],i=1..m),k=0..degree(gu,x)): end: #SeqUpqbAnx(p,q,b,A,N,x): The list consisting of the first N+1 terms, starting at n=0 of UpqbAnx(p,q,b,A,n,x) (q.v.) Try: #SeqUpqbAnx(1+x+x^2,1,2,[3],10,x); SeqUpqbAnx:=proc(p,q,b,A,N,x) local n: [seq(UpqbAnx(p,q,b,A,n,x),n=0..N)]:end: #RSe(p,q,b,A,N,x,t): #Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x) #and a list of non-neg. integers A , and a variable name t, outputs (if possible) the rational function, in t, #whose coefficient of t^n (in its Taylor expansion) is UpqbAnx(p,q,b,A,n,x) (q.v.) It does it #by guessing up to N terms. If N terms do not suffice, it return FAIL. Try: #RSe(1+x+x^2,1,2,[3],10,x,t); RSe:=proc(p,q,b,A,N,x,t) local gu,mu,K: if N<10 then print(N, `should have been at least 10`): RETURN(FAIL): fi: for K from 10 by 5 to N do gu:=SeqUpqbAnx(p,q,b,A,K,x): mu:=GuessRec(gu): if mu<>FAIL then RETURN(factor(CtoR(mu,t))): fi: od: print(N, `is not large enough, make it bigger. `): FAIL: end: #UnL(p,q,b,n,x,L): #Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x) #and a list of pairs of non-neg. integers (of length m, say), L, let #c_k(n) be the coefficient of Fpqbnx(p,q,b,n,x) (q.v.). It outputs #Sum(c_{k+A[1][2]-A[1][1]*b^n}(n)*c_{k+A[2][2]-A[2][1]*b^n}(n)*...c_{k+A[m][2]-A[m][1]*b^n}(n),k=0..infinity). Try #UnL(1+x+x^2,1,2,10,x,[[0,0],[0,0], [0,0]]); UnL:=proc(p,q,b,n,x,L) local i,gu,m,k,gadol: option remember: gadol:=max(seq(L[i][2],i=1..nops(L))): m:=nops(L): gu:=Fpqbnx(p,q,b,n,x): add(mul(coeff(gu,x,k+L[i][2]-b^n*L[i][1]),i=1..m),k=0..degree(gu,x)+b^n*gadol): end: #SeqUnL(p,q,b,N,x,L): The list consisting of the first N+1 terms, starting at n=0 of UnL(p,q,b,n,x,L) (q.v.) Try: #SeqUnL(1+x+x^2,1,2,10,x,[[0,0],[0,0], [0,0]]); SeqUnL:=proc(p,q,b,N,x,L) local n: [seq(UnL(p,q,b,n,x,L),n=0..N)]:end: #CF(L): The canonical form of the list of pairs L. Try #CF([[1,3],[1,4],[2,5]]); CF:=proc(L) local katan,i: katan:=min(seq(L[i][1],i=1..nops(L))): sort([seq([L[i][1]-katan,L[i][2]],i=1..nops(L))]): end: #Eq(p,q,b,x,L,t,Z): The eqation expressing the generating function of UnL(p,q,b,n,L) called Z[L], in terms of other #UnL(p,q,b,n,L1)'s called Z[L1]. It also returns the set of L1's that show up. Try: #Eq(1+x+x^2,1,2,x,[[0,0],[0,0]],t,Z); Eq:=proc(p,q,b,x,L,t,Z) local gu,eq,y,Y,gu1,L1,coe,yeladim,i,i1: eq:=Z[L]-UnL(p,q,b,0,x,L): gu:=mul(Y[i1]^(b*L[i1][1]),i1=1..nops(L))*mul(y[i1]^L[i1][2],i1=1..nops(L)): gu:=expand(gu*mul(subs(x=Y[i1],p),i1=1..nops(L))): yeladim:={}: for i from 1 to nops(gu) do gu1:=op(i,gu): L1:=[seq([degree(gu1,Y[i1]),degree(gu1,y[i1])],i1=1..nops(L))]: coe:=gu1/mul(Y[i1]^L1[i1][1],i1=1..nops(L1))/mul(y[i1]^L1[i1][2],i1=1..nops(L1)): L1:=CF(L1): if SeqUnL(p,q,b,5,x,L1)<>[0$6] then eq:=eq-coe*t*Z[L1]: yeladim:=yeladim union {L1}: fi: od: eq,yeladim: end: #RS(p,q,b,A,N,x,t): #Inputs polynomials p, q in the variable x, a positive integer b, a non-neg. integer n (and a variable name x) #and a list of non-neg. integers A , and a variable name t, outputs the rational function, in t, #whose coefficient of t^n (in its Taylor expansion) is UpqbAnx(p,q,b,A,n,x) (q.v.) It does it #RIGOROUSLY by setting up, dynamically, a system of linear equations, and then #solving it. Try: #RS(1+x+x^2,1,2,[3],10,x,t); RS:=proc(p,q,b,A,x,t) local L,L1,eq,var,StillToDo,Z,AlreadyDone,gu,i,k1: L:=[seq([0,i-1]$A[i],i=1..nops(A))]: eq:={}: var:={}: StillToDo:={L}: AlreadyDone:={}: while StillToDo<>{} do L1:=StillToDo[1]: StillToDo:=StillToDo minus {L1}: AlreadyDone:=AlreadyDone union {L1}: gu:=Eq(p,q,b,x,L1,t,Z): eq:=eq union {gu[1]}: StillToDo:=StillToDo union (gu[2] minus AlreadyDone): var:=var union {seq(Z[k1],k1 in gu[2])}: od: var:=solve(eq,var): factor(subs(var,Z[L])): end: #CheckRS(p,q,b,A,x,N): checks procedure RS(p,q,b,A,x,t) up to N terms. Try: #CheckRS(1+x+x^2,1,2,[3],x,15); CheckRS:=proc(p,q,b,A,x,N) local gu,t,lu,i: gu:=RS(p,q,b,A,x,t): lu:=SeqUpqbAnx(p,q,b,A,N,x): gu:=taylor(gu,t=0,N+1): evalb([seq(coeff(gu,t,i),i=0..N)]=lu): end: #RSv(p,q,b,A,N,x,t): Verbose version of RS(p,q,b,A,N,x,t) (q.v.). Try: #RSv(1+x+x^2,1,2,[3],10,x,t); RSv:=proc(p,q,b,A,x,t) local gu,i,n,B,a,k,t0: t0:=time(): gu:=RS(p,q,b,A,x,t): print(``): print(`---------------------------------`): print(``): print(` Let `): print(``): print(F[n](x)=q*Product(subs(x=x^(b^i),p),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(B(n)=Sum(mul(a(n,k+i-1)^A[i],i=1..nops(A)),k=0..infinity)): print(``): print(`Then `): print(``): print(Sum(B(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(p,q,b,x,t,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(p,q,b,[i],x,t) for i from 1 to N. Try: #Paper1(1+x+x^2,1,2,x,t,5): Paper1:=proc(p,q,b,x,t,N) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanely-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N do print(`Theorem Number`, i): print(``): RSv(p,q,b,[i],x,t): print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: #Paper2(p,q,b,x,t,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(p,q,b,[1$i],x,t) for i from 1 to N. Try: #Paper2(1+x+x^2,1,2,x,t,5): Paper2:=proc(p,q,b,x,t,N) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanely-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N do print(`Theorem Number`, i): print(``): RSv(p,q,b,[1$i],x,t): print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: #RSvE(p,q,b,A,N,x,t): Verbose version of RSe(p,q,b,A,N,x,t) (q.v.). Try: #RSvE(1+x+x^2,1,2,[3],10,x,t); RSvE:=proc(p,q,b,A,N,x,t) local gu,i,n,B,a,k,t0: t0:=time(): gu:=RSe(p,q,b,A,N,x,t): print(``): print(`---------------------------------`): print(``): print(` Let `): print(``): print(F[n](x)=q*Product(subs(x=x^(b^i),p),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(B(n)=Sum(mul(a(n,k+i-1)^A[i],i=1..nops(A)),k=0..infinity)): print(``): print(`Then `): print(``): print(Sum(B(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: #Paper1e(p,q,b,x,t,N1,N2): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t #outputs the paper with N1 theorems RSvE(p,q,b,[i],N2,x,t) for i from 1 to N. Try: #Paper1e(1+x+x^2,1,2,x,t,5,20): Paper1e:=proc(p,q,b,x,t,N1,N2) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanely-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N1 do print(`Theorem Number`, i): print(``): RSvE(p,q,b,[i],N2,x,t): print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: #Paper2e(p,q,b,x,t,N1,N2): inputs polynomials p,q, in the variable x and a positive integer b and a variable name t #outputs the paper with N1 theorems RSvE(p,q,b,[1$i],N2,x,t) for i from 1 to N. Try: #Paper2e(1+x+x^2,1,2,x,t,5,20): Paper2e:=proc(p,q,b,x,t,N1,N2) local t0,i: t0:=time(): print(`Rational generating functions for the Certain Stanely-Stern Sums`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i from 1 to N1 do print(`Theorem Number`, i): print(``): RSvE(p,q,b,[1$i],N2,x,t): print(``): od: print(``): print(`-----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to produce. `): print(``): end: ####Start Annotated part #RSano(p,q,b,A,N,x,t): An annotated version of RS(p,q,b,A,N,x,t) (q.v.). Try: #RSano(1+x+x^2,1,2,[3],10,x,t); RSano:=proc(p,q,b,A,x,t) local L,L1,eq,var,StillToDo,Z,AlreadyDone,gu,i,k1,yofee: L:=[seq([0,i-1]$A[i],i=1..nops(A))]: eq:={}: var:={}: StillToDo:={L}: AlreadyDone:={}: while StillToDo<>{} do L1:=StillToDo[1]: StillToDo:=StillToDo minus {L1}: AlreadyDone:=AlreadyDone union {L1}: gu:=Eq(p,q,b,x,L1,t,Z): eq:=eq union {gu[1]}: StillToDo:=StillToDo union (gu[2] minus AlreadyDone): var:=var union {seq(Z[k1],k1 in gu[2])}: od: print(`We have `, nops(eq), ` equations and `, nops(var), `variables `): print(``): print(`The set of states has`, nops(AlreadyDone), `members, here there are `): print(``): print(AlreadyDone): print(``): var:=solve(eq,var): print(`The solution is`): print(``): print(var): print(``): yofee:=factor(subs(var,Z[L])): print(`The degree of the denominator of the generating function (alias order of the recurrence of the sequence of interest) is`, degree(denom(yofee),t)): print(``): print(`Finally, the generating function itself is`): print(``): print(yofee): print(``): print(`and in Maple notation `): print(``): lprint(yofee): end: