#SCHUTZENBERGER: A PACKAGE FOR HANDLING AND CONJECTURING AND PROVING #ALGEBRAIC EQUATIONS SATISFIED BY INTEGER SEQUENCES Written by Doron #Zeilberger, formerly Temple University #now (2002-) Rutgers University (zeilberg@math.rutgers.edu) #Last Update: June 27, 2017 (in collaboration with Mingjia Yang) #adding two new procedures:radtoalg,radtorecV #PREVIOUS UPDATES: Jan. 14, 1997 #Last UPDATE: Nov. 4, 2009 (to add qfindrec, by request of Stavros Garoufalidis) #The most current version is available from #http://www.math.temple.edu/~zeilberg #To use it get it first, call it `SCHUTZENBERGER` then go into Maple #by typing `maple`. Then, in Maple , type `read SCHUTZENBERGER;` #and see the help files by doing ezra() and ezra(function_name); #(Lots of overlap with Salvy and Zimmerman's Maple package #gfun available by anon. ftp to inria.fr) print(`SCHUTZENBERGER: A PACKAGE FOR HANDLING AND CONJECTURING AND PROVING`): print(`ALGEBRAIC EQUATIONS SATISFIED BY INTEGER SEQUENCES, and RELATED STUFF`): print(`This version: June 21, 2017, by Mingjia Yang and Doron Zeilberger`): print(``): print(`Previous version: Nov. 4, 2009 (to add procedure qfindrec, requested by Stravros Garoufalidis `): print(`[Previous Version of Jan. 14, 1997]`): print(`Named after our great master Marco Schutzenberger(1920-1996)`): print(`written by Doron Zeilberger(zeilberg@math.temple.edu).`): print(`The most current version is always available from`): print(`http://www.math.temple.edu/~zeilberg`): print(`(Lots of overlap with Salvy and Zimmerman's Maple package`): print(`gfun available by anon. ftp to inria.fr).`): print(``): print(`For a list of the ORIGINAL procedures type ezra(), for help with`): print(`a specific procedure, type ezra(procedure_name)`): print(``): print(`For a list of the NEW procedures, for the OEIS paper by Shalosh B. Ekhad and Mingjia Yang, type ezraOEIS(), for help with`): print(`a specific procedure, still type ezra(procedure_name)`): print(``): ezraOEIS:=proc() if args=NULL then print(`The new OEIS procedures are:`): print(`algtorecV, algtorecVwp, MatharConjs, OEISpaper, radtoalg,radtorec, radtorecV,radtorecVwp, `): print(``): print(`a specific procedure, type ezra(procedure_name)`): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`Contains the following procedures:`): print(`algprod, algtodiff,algtorec, algtorecE, algtorecV, algtorecVwp, radtoalg,radtorecV,radtorecVwp`): print(` decomp,decompx, difftorec, empir,Empir ,,findrec2, Findrec, genseq,genseqnew,,logalg, MatharConjs, mfindrec`): print(`newts,newt,OperDiv, pol_in_terms_of_p`): print(` OperToRecEq, qfac1,qbin1, qfindrec, radtoalg, radtorec, radtorecV, radtorecVwp`): print(``): print(`a specific procedure, type ezra(procedure_name)`): print(`Warning: q is a global variable, do not use it!`): fi: if nops([args])=1 and op(1,[args])=OEISpaper then print(`OEISpaper(L,P,x,n,a): inputs a list of pairs [OEISAnumber, GeneratingFunctionInTermsOfRadicals], outputs an article`): print(`with rigorously-proved theorems of linear recurrences satisfied by the sequence of coefficients, that`): print(`were formerly (June 2017) stated as "conjectured". For example, try:`): print(`OEISpaper(MatharConjs(),P,x,n,a);`): fi: if nops([args])=1 and op(1,[args])=MatharConjs then print(`MatharConjs(): the set of R.J. Mathar's conjectures from the OEIS. Try:`): print(`MatharConjs();`): fi: if nops([args])=1 and op(1,[args])=radtorecV then print(`radtorecV(F,P,x,n,N,a) takes a generating function F`): print(`and rigorously find and prove the reccurence satisfied by F`): print(`Try: radtorecV(sqrt((1-2*x - sqrt(1-4*x-32*x^2))/2)/(3*x),P,x,n,N,a)--the recurrence we get corresponds to the recurrence conjectured in A200375`): print(`Also try: radtorecV(1+(2*x^2) / (-1+4*x-2*x^2+sqrt(1-4*x)),P,x,n,N,a)--in this case we get a more complicated recurrence than the one conjectured in A129775`): fi: if nops([args])=1 and op(1,[args])=radtorecVwp then print(`Similar to radtorecV, but with proof. radtorecV(F,P,x,n,a) takes a generating function F`): print(`and rigorously find and prove the reccurence satisfied by F`): print(`Try: radtorecVwp(sqrt((1-2*x - sqrt(1-4*x-32*x^2))/2)/(3*x),P,x,n,a)--the recurrence we get corresponds to the recurrence conjectured in A200375`): print(`Also try: radtorecVwp(1+(2*x^2) / (-1+4*x-2*x^2+sqrt(1-4*x)),P,x,n,a)--in this case we get a more complicated recurrence than the one conjectured in A129775`): fi: if nops([args])=1 and op(1,[args])=radtoalg then print(`radtoalg(LHS,RHS) takes two sides of a radical equation (LHS,RHS):in our case, LHS=P, RHS=the generating function F`): print(`where F is usually in the form (poly(x)+Rad(x))^k/poly(x) or poly(x)/(poly(x)+Rad(x))`): print(`and returns the left hand side of the algebraic equation (without radicals) in x and P. (P symbolizes F)`): print(`Try:radtoalg(P,sqrt((1-2*x - sqrt(1-4*x-32*x^2))/2)/(3*x))--output:324*P^4*x^4+36*x^2*(2*x-1)*P^2+36*x^2`): fi: if nops([args])=1 and op(1,[args])=radtorec then print(`radtorec(F,x,n,N): inputs a generating function given in terms of radicals and outputs a linear recurrence operator`): print(`annihilating the sequence of coefficients. Try:`): print(` radtorec((1-sqrt(1-4*x))/(2*x),x,n,N); `): fi: if nops([args])=1 and op(1,[args])=OperToRecEq then print(`OperToRecEq(ope,n,N,a): inputs a linear recurrence operator ope in n and the forward shift operator N,`): print(`and a symbol a, outputs the linear recurrence equation that asserts that the sequence is annihilated by`): print(`the operator, in terms of decreasing arguments of n. Try`): print(`OperToRecEq(N-(n+1),n,N,a);`): fi: if nops([args])=1 and op(1,[args])=OperDiv then print(`OperDiv(ope1,ope2,n,N): given linear recurrence operators ope1`): print(`and (monin in N) ope2 in n and N (the forward shift operator in N),`): print(`outputs operators Q(N,n) and R(N,n) such that`): print(`ope1(N)=Q(N,n)ope2(n,N)+R(N,n) where the order of Q(N,n)(i.e. degree in N)`): print(`is the difference in the oders of ope1 and ope2 and the order of`): print(`R(N,n) is less than the order of ope2 `): print(` For example, try: `): print(` OperDiv(n*N^2+n/(n+1)*N,N+1/n,n,N); `): fi: if nops([args])=1 and op(1,[args])=qfac1 then print(`qfac1(q,n): (1-q)*...*(1-q^n) `): fi: if nops([args])=1 and op(1,[args])=qbin1 then print(`qbin1(q,n,k): qfac1(q,n)/qfac1(q,k1)/qfac1(q,n-k) `): fi: if nops([args])=1 and op(1,[args])=`decompx` then print(`decompx(pol,N,x): given a polynomial in the variables N and x`): print(`investigates whether, when viewed as a polynomial in N`): print(` it can be written as the alg. product (i.e. its roots are`): print(`the direct product of the sets of the roots of its " factors")`): print(`of a quadratic in N`): print(`of the form N^2+2*a*x*N-1 where a is an (algebraic) pure number`): print(`and another polynomial in N. with coeffs. that are polynomials in x`): fi: if nops([args])=1 and op(1,[args])=`findrec2` then print(` findrec2(f,DEGREEn,DEGREEk,ORDER,n,k,N)`): print(` finds empirically an ordi. linear recurrence in n`): print(` with polynomial coeffs. The input is a double sequence`): print(`f[n,k] given as a list of lists.`): print(`STARTING at f[1,1],i.e. f[0] is not considered`): print(` where DEGREEn:=the maximal degree of the coefficients in n`): print(` where DEGREEk:=the maximal degree of the coefficients in k`): print(`and ORDER:=the order of the recurrence in n.`): print(` The output is the operator`): print(` in n, k and N, where N is the forward unit shift: Nf(n):=f(n+1).`): print(`For example findrec2([[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]],0,0,1,n,k,N)`): print(` should yield N-1`): fi: if nops([args])=1 and op(1,[args])=`Findrec` then print(`Findrec(f,n,N,MaxC) finds empirically an ordi. linear recurrence`): print(` with polynomial coeffs. The input is a sequence f given as a list`): print(`STARTING at f[1],i.e. f[0] is not considered`): print(` in n and N, where N is the forward unit shift: Nf(n):=f(n+1).`): print(`Of course, if a sequence is not P-recursive there would never`): print(`enough data`): print(`For example, try:`): print(` Findrec([1,1,2,3,5,8,13,21,34],n,N,2); `): print( `Findrec([1,2,5,14,42,132,429],m,M,2);`): fi: if nops([args])=1 and op(1,[args])=`qfindrec` then print(`qfindrec(f,q,DEGREE,ORDER,n,N) tries to finds empirically an ordi. linear recurrence`): print(` with polynomial coeffs. in (q,q^n) a recurrence of degree DEGREE and order ORDER.`): print(` The input is a sequence f given as a list`): print(`of rational functions in q`): print(`STARTING at f[1],i.e. f[0] is not considered`): print(` in n and N, where N is the forward unit shift: Nf(n):=f(n+1).`): print(`For example:`): print(` qfindrec([seq(q^i,i=1..20)],q,0,1,n,N); `): print(` should yield N-q .`): print(`If there is no recurrence of the specified order and degree, it returns FAIL`): fi: if nops([args])=1 and op(1,[args])=`qFindrec` then print(`qFindrec(f,q,n,N, MaxC): guesses a linear recurrence operator (in N, and q^n) with polynomial coefficients`): print(` annihilating the sequence of complexity <=MaxC. Try: `): print(` qFindrec([seq(qfac1(q,i),i=1..10)],q,n,N,2); `): fi: if nops([args])=1 and op(1,[args])=`decomp` then print(`decomp(pol,x): given a polynomial in x, pol, decides whether it is`): print(`the algebraic product of a quadratic polynomial pol1, and another pol`): print(`pol2. If yes, it returns pol1 and pol2, if not, it returns 0`): fi: if nops([args])=1 and op(1,[args])=`pol_in_terms_of_p` then print(`pol_in_terms_of_p(resh,x) :given the list of power sums [p_1,..,p_deg]`): print(`finds the monic polynomial, in x, whose roots are the power sums of`): fi: if nops([args])=1 and op(1,[args])=`newts` then print(`newts(pol,x): gives a list of the power sums of the roots pol(x)`): print(`from p_1 to p_(degree of pol(x)), using newtson's equations`): fi: if nops([args])=1 and op(1,[args])=`newt` then print(`newt(pol,x,R): gives a list of the power sums of the roots pol(x)`): print(`from p_1 to p_R, using newtson's equations`): fi: if nops([args])=1 and op(1,[args])=`algprod` then print(`algprod(P,Q,x): inputs polynomials P and Q in x , and outputs`): print(`a poly R, in x, of degree deg(P)deg(Q) whose roots are all the`): print(`possible products of the roots of P and Q`): fi: if nops([args])=1 and op(1,[args])=`logalg` then print(`logalg(F,P,x,ORDER): Using Comtet's method, given a polynomial`): print(`F(P,x), finds the linear diff. eq. of order ORDER satisfied by the`): print(`f.p.s. log(f(x)), where f(x) satisfies the alg.eq. F(f(x),x)=0 `): print(`E.g. try logalg(P-1-x*P^2,P,x,3) to get the differential `): print(`operator annihilating the log of the g.f for the Catalan numbers`): print(`ORDER=degree(F,P)+1 should always work, but you may be able to`): print(`to get away with lower order`): fi: if nops([args])=1 and op(1,[args])=`difftorec` then print(`difftorec(ope1,D,x,n,N): Converts the differential eq. (linear)`): print( `with poly coeff. satisfied by the g.f f(x) given by the opertor`): print(`ope (D,x),to the recurrence, in n and N satisfied `): print(`by the coeffs. of f(x). For example`): print(`difftorec(x*D-1,D,x,n,N) will yield n-1`): fi: if nops([args])=1 and op(1,[args])=`algtorec` then print(` algtorec(F,P,x,n,N): Inputs a polynomial F in the two variables P and x, where P is the ordinary generating function `): print(` of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator `): print(` for the Catalan numbers when they are defined by P=1+x*P^2, try `): print(` algtorec(P-1-x*P^2,P,x,n,N); `): print(`To rigorously prove Mathar's conjecture in OEIS A(x) = (1 - sqrt(1 - 4*(x-x^2)/(1-x+x^2) ))/2, https://oeis.org/A090826 `): print(` Type: `): print(` ope:=algtorec(4*P^2-8*P^2*x-4*P^2*x^2-4*P+8*P^2*x^3+4*P*x+4*P^2*x^4+4*P*x^2+4*x,P,x,n,N); `): print(` OperToRecEq(ope,n,N,a); `): fi: if nops([args])=1 and op(1,[args])=`algtorecE` then print(`algtorecE(F,P,x,n,N,MaxC): Like algtorecE(F,P,x,n,N) but by pure guessing, MaxC is the maximal complexity`): print(`of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator `): print(` for the Catalan numbers when they are defined by P=1+x*P^2, try: `): print(` algtorecE(P-1-x*P^2,P,x,n,N,2); `): print(`To empirically (in fact semi-rigorously) discover Mathar's conjecture in OEIS A(x) = (1 - sqrt(1 - 4*(x-x^2)/(1-x+x^2) ))/2, https://oeis.org/A090826 `): print(` Type: `): print(` ope:=algtorecE(4*P^2-8*P^2*x-4*P^2*x^2-4*P+8*P^2*x^3+4*P*x+4*P^2*x^4+4*P*x^2+4*x,P,x,n,N,4); `): print(` OperToRecEq(ope,n,N,a); `): fi: if nops([args])=1 and op(1,[args])=`algtorecV` then print(`algtorecV(F,P,x,n,N,a): Verbose form of algtorec(F,P,x,n,N) and algtorecE(F,P,x,n,N,MaxC) (q.c) `): print(`To rigorously discover and prove Mathar's conjecture in OEIS A(x) = (1 - sqrt(1 - 4*(x-x^2)/(1-x+x^2) ))/2, https://oeis.org/A090826 `): print(`and have a human-readable statement `): print(` Type: `): print(` algtorecV(4*P^2-8*P^2*x-4*P^2*x^2-4*P+8*P^2*x^3+4*P*x+4*P^2*x^4+4*P*x^2+4*x,P,x,n,N,a); `): print(`Also try: `): print(` algtorecV(P-1-(x+3*x^2)*P^2-x*P^3,P,x,n,N,a); `): fi: if nops([args])=1 and op(1,[args])=`algtorecVwp` then print(`algtorecVwp(F,P,x,n,N,a): Like algtorecV(F,P,x,n,N,a) (q.v.) but with proof. `): print(`To rigorously discover and prove Mathar's conjecture in OEIS A(x) = (1 - sqrt(1 - 4*(x-x^2)/(1-x+x^2) ))/2, https://oeis.org/A090826 `): print(`and have a human-readable statement `): print(` Type: `): print(` algtorecVwp(4*P^2-8*P^2*x-4*P^2*x^2-4*P+8*P^2*x^3+4*P*x+4*P^2*x^4+4*P*x^2+4*x,P,x,n,N,a); `): print(`Also try: `): print(` algtorecVwp(P-1-(x+3*x^2)*P^2-x*P^3,P,x,n,N,a); `): fi: if nops([args])=1 and op(1,[args])=`algtodiff` then print(`algtodiff(F,P,x,ORDER): Using Comtet's method, given a polynomial`): print(`F(P,x), finds the linear diff. eq. of order ORDER satisfied by the`): print(`formal power series f(x) that satisfies the alg.eq. F(f(x),x)=0 `): print(`Given an algebraic equation F(P(x),x) and desired`): print(`order, finds, if one exists, a linear diff. eq. of order ORDER`): print(`if you then do difftorec(",D,x,n,N), you would get a new proof`): print(`that the Catalan numbers, defined by their combinatorial meaning`): print(`indeed equals (2n)!/(n!(n+1)!)`): print(`ORDER=degree(F,P) should always work, but you may be able to`): print(`to get away with lower order`): print(`Warning: it does not always gives the MINIMAL operator`): print(``): print(`For example try: `): print(` algtodiff(P-1-x*P^2,P,x,2); `): fi: if nops([args])=1 and op(1,[args])=`genseq` then print(`genseq(eq,P,x,K) Generates the coefficients from x^0`): print(`through x^K of the formal`): print(` powers series that is a solution of the eq P=eq(x,P)`): print(`For example genseq(1/(1-x-x^2)-1,P,x,4) should yield`): print(`[1,1,2,3,5] and genseq(1+x*P^2,P,x,5) should yield [1,1,2,5,14,42]`): fi: if nops([args])=1 and op(1,[args])=`genseqnew` then print(`genseqnew(eq,P,x,K,PHI) outputs, the coeffs. from x^0 to x^K, as`): print(`a list of polynomials in t (a global variable) of the f.p.s. `): print(`PHI(x,t) that is a solution of the algebraic functional equation`): print(`PHI(x,t)=eq(PHI(x,t),P(x),x,t), where P(x)=PHI(x,1)`): print(`Such a functional equation came up in Doron Zeilbeger's proof`): print(`of Julian West's conj. about 2-Stack-sortable permutations.`): print(`See Discrete Math 102(1992) 85-93. For example, try `): print(` genseqnew(PHI-1/(1-x*t)-x*t*(P-t*PHI)*(P-PHI)/(1-t)^2,P,x,5,PHI)`): print(`to get the first 6 terms in the f.p.s. of interest in the`): print(`above paper`): print(`by taking K=40 and plugging in t=1, you can get a semi-rigorous`): print(`proof of Julian West's conj. Of course the above paper has a `): print(`proof`): fi: if nops([args])=1 and op(1,[args])=`empir` then print(`empir(gu,degx,degP,x,P) empirically finds an algebraic equation`): print(` F(P(x),x)=0 of degree degP in P(x) and degx in x for `): print(`formal power series P(x):=sum_i gu[i]*x^i, where gu is a list`): print(`for example empir([seq(1,i=1..20)],1,1,x,P) should yield`): print(` P-xP-1=0`): fi: if nops([args])=1 and op(1,[args])=`Empir` then print(`Empir(gu,x,P) empirically finds an algebraic equation`): print(` F(P(x),x)=0 for the `): print(`formal power series P(x):=sum_{i=0} gu[i]*x^i, where gu is a list`): print(`for example Empir([seq(1,i=1..20)],x,P) should yield`): print(` P-xP-1=0`): print(`If there is not enough data, it reurns 0. You should then`): print(`try it again with a longer sequence`): fi: if nops([args])=1 and op(1,[args])=`mfindrec` then print(`mfindrec(DEG,ORDER,f,p): like findrec(DEG,ORDER,f), but mod p`): fi: end: pashet:=proc(p,n,N) local i,gu1,gu,p1,ra: p1:=normal(p): gu1:=denom(p1): ra:=degree(gu1,N): p1:=subs(n=n+ra,numer(p1)): p1:=expand(p1): gu:=0: for i from 0 to degree(p1,N) do gu:=gu+factor(coeff(p1,N,i))*N^i: od: RETURN(gu): end: difftorec:=proc(ope1,D,x,n,N) local i,j,gu,ord,ope,r,mu: ope:=expand(ope1): gu:=0: ord:=degree(ope,D): for i from 0 to ord do mu:=coeff(ope,D,i): for j from 0 to degree(mu,x) do gu:=gu+coeff(mu,x,j)*product(n+r-j,r=1..i)*N^(i-j): od: od: pashet(gu,n,N): end: simp:=proc(pa,deg,bitui,P,x) local mone,mekh,gu: gu:=normal(bitui): mone:=numer(gu): mekh:=denom(gu): mone:=simp1(pa,deg,mone,P,x): mekh:=simp1(pa,deg,mekh,P,x): expand(mone)/expand(mekh): end: simp1:=proc(pa,deg,bitui,P,x) local i,gu: gu:=expand(bitui): for i from deg to 3*deg+1 do gu:=subs(P(x)^i=pa[i],gu): od: gu: end: algtodiff:=proc(F,P,x,ORDER) local i,gu,a,lu,pa,deg,eq1,lu1,y,var,eq,ope,pip: deg:=degree(F,P): pa:=array(deg..3*deg+1): eq1:=subs(P^deg=y,F): pa[deg]:=solve(eq1=0,y): pa[deg]:=subs(P=P(x),pa[deg]): for i from 1 to 2*deg+1 do lu:=pa[deg+i-1]*P(x): lu:=expand(lu): lu:=subs(P(x)^deg=pa[deg],lu): pa[deg+i]:=lu: od: gu:=a[0]*P(x): lu:=subs(P=P(x),F): lu1:=diff(lu,x): lu1:=subs(diff(P(x),x)=y,lu1): lu1:=solve(lu1=0,y): lu1:=simp(pa,deg,lu1,P,x): lu1:=normal(lu1): gu:=gu+a[1]*lu1: gu:=normal(gu): gu:=simp(pa,deg,gu,P,x): var:={a[0],a[1]}: ope:=a[0]+a[1]*D: lu:=lu1: for i from 2 to ORDER do lu:=diff(lu,x): lu:=subs(diff(P(x),x)=lu1,lu): lu:=normal(lu): lu:=simp(pa,deg,lu,P,x): gu:=gu+a[i]*lu: gu:=normal(gu): gu:=simp(pa,deg,gu,P,x): var:=var union {a[i]}: ope:=ope+a[i]*D^i: od: gu:=normal(gu): gu:=numer(gu): gu:=expand(gu): eq:={}: for i from 0 to deg-1 do eq:=eq union {coeff( gu,P(x),i)=0}: od: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: ope:=subs(var,ope): if ope=0 then RETURN(FAIL): fi: ope:=normal(ope): ope:=numer(ope): ope:=expand(ope): pip:=coeff(ope,D,degree(ope,D)): pip:=coeff(pip,x,degree(pip,x)): normal(ope/pip): end: genseq:=proc(eq,P,x,K) local i,gu,P1: gu:=[coeff(coeff(eq,P,0),x,0)]: P1:=coeff(eq,P,0): for i from 1 to K do P1:=subs(P=P1,eq): P1:=taylor(P1,x=0,i+1): gu:=[op(gu),coeff(P1,x,i)]: od: gu: end: #genseqnew:=proc(eq,P,x,K,PHI) #PHI(x,t) that is a solution of the algebraic functional equation #PHI(x,t)=eq(PHI(x,t),P(x),x,t), where P(x)=PHI(x,1) genseqnew:=proc(eq,P,x,K,PHI) local i,gu,P1,eq1,eq2,PHI1: eq1:=normal(eq): gu:=[1]: PHI1:=1: P1:=1: for i from 1 to K do eq2:=subs(PHI=PHI1,eq1): eq2:=normal(eq2): eq2:=subs(P=P1,eq2): eq2:=normal(eq2): PHI1:=taylor(eq2,x=0,i+1): PHI1:=normal(PHI1): P1:=subs(t=1,PHI1): gu:=[op(gu),coeff(PHI1,x,i)]: od: gu: end: #logalg(F,P,x,ORDER): Using #an adaptation of Comtet's method, given a polynomial #F(P,x), finds the linear diff. eq. of order ORDER satisfied by the #formal power series log (f(x)), where f(x) satisfies the alg.eq. F(f(x),x)=0 logalg:=proc(F,P,x,ORDER) local i,gu,a,lu,pa,deg,eq1,lu1,lu1a,y,var,eq,ope: deg:=degree(F,P): pa:=array(deg..3*deg+1): eq1:=subs(P^deg=y,F): pa[deg]:=solve(eq1=0,y): pa[deg]:=subs(P=P(x),pa[deg]): for i from 1 to 2*deg+1 do lu:=pa[deg+i-1]*P(x): lu:=expand(lu): lu:=subs(P(x)^deg=pa[deg],lu): pa[deg+i]:=lu: od: gu:=0: lu:=subs(P=P(x),F): lu1:=diff(lu,x): lu1:=subs(diff(P(x),x)=y,lu1): lu1:=solve(lu1=0,y): lu1a:=lu1/P(x): lu1a:=simp(pa,deg,lu1a,P,x): lu1a:=normal(lu1a): gu:=gu+a[1]*lu1a: gu:=normal(gu): gu:=simp(pa,deg,gu,P,x): var:={a[1]}: ope:=a[1]*D: lu:=lu1a: for i from 2 to ORDER do lu:=diff(lu,x): lu:=subs(diff(P(x),x)=lu1,lu): lu:=normal(lu): lu:=simp(pa,deg,lu,P,x): gu:=gu+a[i]*lu: gu:=normal(gu): gu:=simp(pa,deg,gu,P,x): var:=var union {a[i]}: ope:=ope+a[i]*D^i: od: gu:=normal(gu): gu:=numer(gu): gu:=expand(gu): eq:={}: for i from 0 to deg-1 do eq:=eq union {coeff( gu,P(x),i)=0}: od: var:=solve(eq,var): ope:=subs(var,ope): ope:=normal(ope): ope:=numer(ope): ope:=expand(ope): #pip:=coeff(ope,D,degree(ope,D)): #pip:=coeff(pip,x,degree(pip,x)): #normal(ope/pip): end: #algprod(P,Q,x): inputs polynomials P and Q in x , and outputs #a poly R, in x, of degree deg(P)deg(Q) whose roots are all the #possible products of the roots of P and Q algprod:=proc(P,Q,x) local degp,degq,deg,pol,var,eq,a,i,mu,lu,reshp,reshq,xp,xq,pol1,i1,i2: degp:=degree(P,x): degq:=degree(Q,x): deg:=degp*degq: reshp:=[]: for i from 0 to degp-1 do reshp:=[op(reshp),x^i]: od: mu:=P/coeff(P,x,degp): mu:=x^degp-mu: reshp:=[op(reshp),mu]: lu:=mu: for i from 1 to deg-degp do lu:=expand(x*lu): lu:=subs(x^(degp)=mu,lu): lu:=expand(lu): reshp:=[op(reshp),lu]: od: reshq:=[]: for i from 0 to degq-1 do reshq:=[op(reshq),x^i]: od: mu:=Q/coeff(Q,x,degq): mu:=x^degq-mu: reshq:=[op(reshq),mu]: lu:=mu: for i from 1 to deg-degq do lu:=expand(x*lu): lu:=subs(x^(degq)=mu,lu): lu:=expand(lu): reshq:=[op(reshq),lu]: od: reshp:=subs(x=xp,reshp): reshq:=subs(x=xq,reshq): pol:=0: pol1:=0: eq:={}: var:={}: for i from 0 to deg do pol:=pol+a[i]*x^i: var:=var union {a[i]}: pol1:=expand(pol1+a[i]*op(i+1,reshp)*op(i+1,reshq)): od: for i1 from 0 to degp-1 do for i2 from 0 to degq-1 do eq:=eq union {coeff(coeff(pol1,xp,i1),xq,i2)}: od: od: var:=solve(eq,var): pol:=subs(var,pol): normal(pol/coeff(pol,x,degree(pol,x))): end: newts:=proc(pol,x) local mu,deg,resh,pol1,r,i: deg:=degree(pol,x): pol1:=pol/coeff(pol,x,deg): pol1:=expand(pol1): resh:=[-coeff(pol1,x,deg-1)]: for r from 2 to deg do mu:=0: for i from 1 to r-1 do mu:=mu-coeff(pol1,x,deg-i)*op(r-i,resh): od: mu:=mu-r*coeff(pol1,x,deg-r): mu:=expand(mu): resh:=[op(resh),mu]: od: resh: end: newt:=proc(pol,x,R) local mu,deg,resh,pol1,r,i: if R<=degree(pol,x) then RETURN([op(1..R,newts(pol,x))]): fi: deg:=degree(pol,x): pol1:=pol/coeff(pol,x,deg): pol1:=expand(pol1): resh:=newts(pol,x): for r from deg+1 to R do mu:=0: for i from 1 to deg do mu:=mu-coeff(pol1,x,deg-i)*op(r-i,resh): od: mu:=expand(mu): resh:=[op(resh),mu]: od: resh: end: #pol_in_terms_of_p(resh) :given the list of power sums [p_1,..,p_deg] #finds the monic polynomial whose roots are the power sums of pol_in_terms_of_p:=proc(resh,x) local pol,mu,deg,r,i: deg:=nops(resh): pol:=x^deg: for r from 1 to deg do mu:=-op(r,resh): for i from 1 to r-1 do mu:=mu-coeff(pol,x,deg-i)*op(r-i,resh): od: mu:=mu/r: mu:=expand(mu): pol:=pol+mu*x^(deg-r): od: pol: end: decomp:=proc(pol,x) local i,deg,resh1,resh2,resh,a,b,pol1,pol2,eq,var,resh2a,eq1: deg:=degree(pol,x): deg:=deg/2: if not type(deg, integer) then RETURN(0): fi: pol1:=x^2+a*x+b: var:={a,b}: resh1:=newt(pol1,x,2*deg): resh:=newts(pol,x): resh2:=[]: for i from 1 to 2*deg do resh2:=[op(resh2),op(i,resh)/op(i,resh1)]: od: pol2:=pol_in_terms_of_p([op(1..deg,resh2)],x): resh2a:=newt(pol2,x,2*deg): eq:={}: for i from 1 to deg do eq1:=normal(op(i,resh2)-op(i,resh2a)): if eq1<>0 then ERROR(`Something is wrong`): fi: od: for i from deg+1 to 2*deg do eq1:=normal(op(i,resh2)-op(i,resh2a)): eq:=eq union {eq1}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: pol1:=subs(var,pol1): pol2:=subs(var,pol2): pol1,pol2: end: decompx:=proc(pol,N,x) local i1,i,mu,deg,resh1,resh2,resh,a,pol1,eq,var: deg:=degree(pol,N): deg:=deg/2: if not type(deg, integer) then RETURN(0): fi: pol1:=N^2+2*a*x*N-1: var:={a}: resh1:=newt(pol1,N,2*deg): resh:=newts(pol,N): resh2:=[]: for i from 1 to 2*deg do resh2:=[op(resh2), rem(op(i,resh),op(i,resh1),x)]: od: eq:={}: for i from 1 to 3 do mu:=op(i,resh2): mu:=expand(numer(normal(mu))): for i1 from 0 to degree(mu,x) do eq:=eq union {numer(normal(coeff(mu,x,i1)))=0}: od: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: var: end: findrec2:=proc(f,DEGREEn,DEGREEk,ORDER,n,k,N) local ope,var,eq,a,i,i1,j,n0,k0,kv,var1,eq1,mu: if (1+DEGREEn)*(1+DEGREEk)*(2+ORDER)+1+ORDER>nops(f)*nops(op(1,f)) then ERROR(` Not enough date`): fi: ope:=0: var:={}: for i from 0 to ORDER do for i1 from 0 to DEGREEn do for j from 0 to DEGREEk do ope:=ope+a[i,i1,j]*n^i1*k^j*N^i: var:=var union {a[i,i1,j]}: od: od: od: eq:={}: mu:=trunc(evalf(sqrt((1+DEGREEn)*(1+ORDER)*(1+DEGREEk)))): #if mu+2+ORDER for n0 from 1 to mu+2 do for k0 from 1 to mu+2 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs({k=k0,n=n0},coeff(ope,N,i))*op(k0,op(n0+i,f)): od: eq:= eq union {eq1}: od: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if nops(kv)>1 then print(` either DEGREE or ORDER are too high`): print(`The output is not the minimal possible operator`): fi: for i from 1 to nops(kv) do ope:=subs(op(i,kv)=1,ope): od: ope: end: findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): RETURN(FAIL): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if nops(kv)>1 then RETURN(0): fi: for i from 1 to nops(kv) do ope:=subs(op(i,kv)=1,ope): od: ope: end: Findrec:=proc(f,n,N,MaxC) local i1,kal,DEGREE, ORDER,L,gu,gu1,i: for L from 1 to MaxC do for ORDER from 1 to L do for DEGREE from 0 to L-ORDER do gu:=findrec(f,DEGREE,ORDER,n,N): if gu<>0 then gu1:=0: for i from 0 to degree(gu,N) do gu1:=gu1+factor(coeff(gu,N,i))*N^i: od: kal:=checkrec(gu1,n,N,f): if kal=[seq(0,i1=1..nops(kal))] then RETURN(gu1): fi: fi: od: od: od: FAIL: end: #checkrec(ope,n,N,Seq1) applies the operator ope(N,n) to the sequence gu checkrec:=proc(ope,n,N,Seq1) local i,j,lu,ru: lu:=[]: for i from 1 to nops(Seq1)-degree(ope,N) do ru:=0: for j from 0 to degree(ope,N) do ru:=ru+subs(n=i,coeff(ope,N,j))*op(i+j,Seq1): od: lu:=[op(lu),ru]: od: normal(lu): end: #empir(gu,degx,degP,x,P) #to "fit" an algebraic equation F(P(x),x)=0 of degree #degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i empir:=proc(gu,degx,degP,x,P) local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip: if (1+degx)*(1+degP) > nops(gu)-3 then RETURN(`sequence too small`): fi: F:=0: var:={}: for i1 from 0 to degx do for i2 from 0 to degP do F:=F+a[i1,i2]*x^i1*P^i2: var:=var union {a[i1,i2]}: od: od: cand:=0: for i1 from 0 to nops(gu)-1 do cand:=cand+op(i1+1,gu)*x^i1: od: lu:=subs(P=cand,F): lu:=taylor(lu,x=0,nops(gu)-1): eq:={}: for i1 from 0 to nops(gu)-2 do eq:=eq union {coeff(lu,x,i1)=0} od: mu:=solve(eq,var): F:=subs(mu,F): if F=0 then RETURN(0): fi: flo:=degree(F,P): pip:=coeff(F,P,flo): flo:=degree(pip,x): pip:=coeff(pip,x,flo): F:=F/pip: normal(F): end: Empir:=proc(gu,x,P) local degx,degP,L,lu: for L from 1 to (nops(gu)-3)/3 do for degP from 1 to L do for degx from 0 to min(trunc((nops(gu)-3)/(1+degP))-1,L-degP) do lu:=empir(gu,degx,degP,x,P): if lu<>0 then RETURN(lu): fi: od: od: od: 0: end: ###qfindrec qfindrec:=proc(f,q,DEGREE,ORDER,n,N) local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu,var1T,eqT,opeT: if (1+DEGREE)*(1+ORDER)+4+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*q^(n*j)*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: eqT:=subs(q=2,eq): var1T:=solve(eqT,var): opeT:=subs(var1T,ope): if opeT=0 then RETURN(FAIL): fi: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if nops(kv)>1 then print(`Too much slack`): fi: if ope<>0 then RETURN(yafe(ope,N)[2]): else RETURN(FAIL): fi: end: ###End qfindrec #qFindrec(f,q,n,N,MaxC): guesses a linear recurrence operator (in N, and q^n) with polynomial coefficients #annihilating the sequence f. Try: #qFindrec([seq(qfac1(q,i),i=1..10)],q,n,N,2); qFindrec:=proc(f,q,n,N,MaxC) local i1,kal,DEG, ORD,L,gu,gu1,i: for L from 1 to MaxC do for ORD from 1 to L do for DEG from 0 to L-ORD do gu:=qfindrec(f,q,DEG,ORD,n,N): if gu<>FAIL then gu1:=0: for i from 0 to degree(gu,N) do gu1:=gu1+factor(coeff(gu,N,i))*N^i: od: kal:=checkrec(gu1,n,N,f): if kal=[seq(0,i1=1..nops(kal))] then RETURN(gu1): else print(gu1, `did not work out`): fi: fi: od: od: od: FAIL: end: yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: qfac1:=proc(q,n1) local i: mul(1-q^i,i=1..n1): end: qbin1:=proc(q,n1,k1): normal(qfac1(q,n1)/qfac1(q,k1)/qfac1(q,n1-k1)):end: #MyGcd(L): inputs a list of polynomials and finds their gcd. Try #MyGcd([n,n^2,n^3]); MyGcd:=proc(L) local gu: if L=[] then RETURN(FAIL): elif nops(L)=1 then RETURN(L[1]): else gu:=MyGcd([op(1..nops(L)-1,L)]): RETURN(gcd(gu,L[nops(L)])): fi: end: #algtorec(F,P,x,n,N): Inputs a polynomial F in the two variables P and x, where P is the ordinary generating function #of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator #for the Catalan numbers when they are defined by P=1+x*P^2, try #algtorec(P-1-x*P^2,P,x,n,N); algtorec:=proc(F,P,x,n,N) local ope,d,i, mu: for d from 1 to degree(F,P) do ope:=algtodiff(F,P,x,d): if ope<>FAIL then ope:=difftorec(ope,D,x,n,N): mu:=MyGcd([seq(coeff(ope,N,i),i=0..degree(ope,N))]): ope:=normal(ope/mu): ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)): RETURN(ope): fi: od: FAIL: end: #algtorecE(F,P,x,n,N,MaxC): Like algtorecE(F,P,x,n,N) but by pure guessing, MaxC is the maximal comlexity #of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator #for the Catalan numbers when they are defined by P=1+x*P^2, try #algtodiffE(P-1-x*P^2,P,x,n,N); algtorecE:=proc(F,P,x,n,N,MaxC) local deg, lu,mu,F1,kama: lu:=coeff(coeff(F,P,1),x,0): if lu=0 then RETURN(FAIL): fi: F1:=expand(F-lu*P): F1:=-F1/lu: if coeff(coeff(F1,P,0),x,0)<>0 then RETURN(FAIL): fi: kama:= max(seq((1+deg)*(1+MaxC-deg)+7+MaxC-deg,deg=0..MaxC)): mu:=genseq(F1,P,x,kama): mu:=[op(2..nops(mu),mu)]: Findrec(mu,n,N,MaxC): end: #OperDiv(ope1,ope2,n,N): given linear recurrence operators ope1 #and (monin in N) ope2 in n and N (the forward shift operator in N), #outputs operators Q(N,n) and R(N,n) such that #ope1(N)=Q(N,n)ope2(n,N)+R(N,n) where the order of Q(N,n)(i.e. degree in N) #is the difference in the oders of ope1 and ope2 and the order of #R(N,n) is less than the order of ope2 #For example, try #OperDiv(n*N^2+n/(n+1)*N,N+1/n,n,N); OperDiv:=proc(ope1,ope2,n,N) local gu,lu,ORDER1,ORDER2,ope1A,i1: ORDER1:=degree(ope1,N): ORDER2:=degree(ope2,N): if coeff(ope2,N,ORDER2)<>1 then RETURN(FAIL): fi: if degree(ope1,N)FAIL then MaxC1:=degree(ope2,N)+degree(ope2,n): else MaxC1:=MaxC: fi: if MaxC1FAIL then MaxC1:=degree(ope2,N)+degree(ope2,n): else MaxC1:=MaxC: fi: if MaxC1FAIL then MaxC1:=degree(ope2,N)+degree(ope2,n): else MaxC1:=MaxC: fi: MaxC1:=degree(ope2,N)+degree(ope2,n): #if ope2=FAIL then # INI:=[seq(coeff(taylor(normal(F),x=0,degree(ope1,N)+12),x,i),i=1..degree(ope1,N))]: #else # INI:=[seq(coeff(taylor(normal(F),x=0,degree(ope2,N)+12),x,i),i=1..degree(ope2,N))]: #fi: if MaxC1FAIL then MaxC1:=degree(ope2,N)+degree(ope2,n): else MaxC1:=MaxC: fi: MaxC1:=degree(ope2,N)+degree(ope2,n): if ope2=FAIL then INI:=[seq(coeff(taylor(F,x=0,degree(ope1,N)+12),x,i),i=1..degree(ope1,N))]: else INI:=[seq(coeff(taylor(F,x=0,degree(ope2,N)+12),x,i),i=1..degree(ope2,N))]: fi: if MaxC1