###################################################################### ##RPS: Save this file as RPS # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read RPS # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Dec. 9, 2011 print(`Created: Dec. 9, 2011`): print(` This is RPS `): print(`Named after Richard Peter Stanley`): print(`It is one of the packages that accompany the article `): print(`Automatic Solution of Richard Stanley's Amer. Math. Monthly Problem #11610`): print(`and ANY Problem of That Type`): print(`by Shalosh B. Ekhad and Doron Zeilberger`): print(`available from http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/rps.html`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): print(`For a list of the procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): with(combinat): ezra11:=proc() if args=NULL then print(` The secondarysupporting procedures are: HSP `): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(` The supporting procedures are: algtodiff, algtodiff1, AsyC `): print(`AZc, CheckAsy, CheckGWtE, CheckRec,`): print(` DiffToRec1, Empir, EqsVars, findrec, Findrec,`): print(` Images1, Khalek, Milim, SeqFromRec, Wt, WtEd, Zugot`): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: `): print(` ALGg, DOaz, GWtE, Info1, Info2, Mishpat`): print(` ProveGF, RECaz, RECg, RSct, RSseq, TerseBook, VerboseBook`): print(` WtE `): elif nops([args])=1 and op(1,[args])=ALGg then print(`ALGg(m,w1,w2,K,x,P): Given a positive integer`): print(`m and two words (lists), w1, w2 of the same length in the`): print(`alphabet {1, ...,m}, `): print(`and a positive integer K, and symbols x, P`): print(`uses guessing with K data points`): print(`to conjecture a polynomial in x and P`): print(`such that you get zero when`): print(`you plug-in the ordinary generating function of`): print(` the sequence`): print(`"number of words in {1, ...m} of length n that has`): print(`as many factors (substrings) that are w1 and`): print(`factors that are w2 equals `): print(`For example, to get a recurrence operator for`): print(`Richard Stanley's`): print(`Amer. Math. Monthly problem #11610 (page 937`): dprint(` of v. 118 (No. 10), [Dec. 2011]), please type:`): print(`ALGg(2,[1,1],[1,2],40,x,P);`): elif nops([args])=1 and op(1,[args])=algtodiff then print(`algtodiff(F,P,x,D1): Inputs a polynomial F,`): print(`in P and x, and outputs`): print(`a differential operator annihilating`): print(`P(x) where P(x) is the alg. formal power series`): print(`satisfying the algebraic equation`): print(`F(P(x),x))=0 . For example, try:`): print(`algtodiff(P-1-x*P^2,P,x,D1);`): elif nops([args])=1 and op(1,[args])=algtodiff1 then print(`algtodiff1(F,P,x,D1,ORDER): Inputs a polynomial F,`): print(`in P and x, and an ORDER, and outputs`): print(`a differential operator annihilating`): print(`P(x) where P(x) is the alg. formal power series`): print(`satisfying the algebraic equation`): print(`F(P(x),x))=0 . For example, try:`): print(`algtodiff1(P-1-x*P^2,P,x,D1,2);`): elif nops([args])=1 and op(1,[args])=AsyC then print(`AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions `): print(`to ope(n,N)f(n)=0, with the given initial conditions`): print(`where ope(n,N) is a recurrence operator`): print(`up to the M's term,`): print(`and complete with an empirically derived constant in front`): print(`using K terms`): elif nops([args])=1 and op(1,[args])=AZc then print(` AZc(A,y,x,D): The Almkvist-Zeilberger Continuous algorithm`): print(`Given an hyperexponential expression A, of y and x`): print(`returns a differential operator ope(x,D) that annihilates`): print(` the function x->Integral( A(x,y),y);`): print(`For example, try:`): print(`AZc(1/(1-x*(y+1/y))/y,y,x,Dx);`): elif nops([args])=1 and op(1,[args])=CheckAsy then print(`CheckAsy(a,n,L,eps): checks whether the asympotic formula a, an`): print(`expression in n, is a good one for the sequence L. With tolerance`): print(`eps. For example, try:`): print(`CheckAsy(2^n,n,[seq(2^i,i=1..100)],1/100);`): elif nops([args])=1 and op(1,[args])=CheckGWtE then print(`CheckGWtE(m,k,N): checks the first N terms of`): print(`the output of CheckGWtE(m,k,z,t); For example, try:`): print(`CheckGWtE(2,2,5);`): elif nops([args])=1 and op(1,[args])=CheckRec then print(`CheckRec(L,ope,n,N): Given a sequence of numbers L and a`): print(`recurrence operator ope(n,N) checks whether indeed`): print(`L is annihilated by ope(n,N). For example, try:`): print(`CheckRec([1,2,4,8,16],N-2,n,N);`): elif nops([args])=1 and op(1,[args])=DiffToRec1 then print(`DiffToRec1(ope,t,Dt,n,N): Given a linear differential operator`): print(`ope(t,Dt) and symbols n,N, outputs the linear recurrence`): print(`operator with polynomial coefficients satisfied by`): print(`the sequence of coeff. of any f.p.s. annihilating ope(t,Dt)`): print(`where N is the shift in n. For example, try:`): print(`DiffToRec1(Dt+t, t,Dt,n,N);`): elif nops([args])=1 and op(1,[args])=DOaz then print(`DOaz(m,w1,w2,z,Dz): Given a positive integer`): print(`m and two words (lists) of the same length, w1, w2 in the`): print(`alphabet {1, ...,m}, `): print(`uses the continuous Almkvist-Zeilberger algorithm,`): print(`to find a linear differential operator with`): print(`polynomial coefficients, in terms of z and the`): print(`differentiation w.r.t. z, Dz, that annihilates`): print(`the (ordinary) generating function`): print(`of the sequence`): print(`"number of words in {1, ...m} of length n that has`): print(`as many factors (substrings) that are w1 as it has`): print(`factors that are w2. `): print(``): print(`For example, to get the differential operator for`): print(`the generating function for the sequence described in`): print(`Richard Stanley's `): print(`Amer. Math. Monthly problem #11610 (page 937`): print(` of v. 118 (No. 10), [Dec. 2011]), please type:`): print(`DOaz(2,[1,1],[1,2],t,Dt);`): elif nops([args])=1 and op(1,[args])=Empir then print(`Empir(L,x,P): inputs a sequence L and `): print(`and conjectures an polynomial equation`): print(`F(x,P)=0 satisfies by the formal power`): print(`series L[1]+L[2]*x+L[3]*x^2+...`): print(`For example, try:`): print(`Empir([seq(i1,i1=1..40)],x,P);`): elif nops([args])=1 and op(1,[args])=EqsVars then print(`EqsVars(m,k,z,t,A):The set of equations, followed by the`): print(`set of variables that enables to compute the generating`): print(`function in z and t[w1] (w1 in [1,...,m]^k)`): print(`that is the weight-enumerator of all words in`): print(`the alphabet {1, ...,m} with weight z^length(w) times`): print(`product of t[w_i,...,w_i+k-1] for i from 1 to nops(w)-k+1`): print(`A is a symbol for the indexing of the weight enumerator`): print(`of words,`): print(`For example, try`): print(`EqsVars(2,2,z,t,A);`): elif nargs=1 and args[1]=findrec then print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`): print(`the sequence f of degree DEGREE and order ORDER.`): print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`): elif nargs=1 and args[1]=Findrec then print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`): print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `): print(`of maximum DEGREE+ORDER<=MaxC`): print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`): elif nops([args])=1 and op(1,[args])=HSP then print(`HSP(ope,n,N): the highest singularity of the recurrence`): print(`operator ope(n,N). For example, try:`): print(`HSP((n-3)*(n-4)*N+1,n,N);`): elif nops([args])=1 and op(1,[args])=Images1 then print(`Images1(m,S): Given a set of words S in the alphabet`): print(`{1, ...,m}, finds the set of images under permuting`): print(`the letters and reversing the words.`): print(`For example, try:`): print(`Images1(2,{[1,1],[1,2]});`): elif nops([args])=1 and op(1,[args])=Info1 then print(`Info1(m,w1,w2,n,N,P,t,K): Given a pos. integer m,`): print(`and two words w1 and w2 of the same length in the`): print(`alphabet {1, ..., m}, and symbols n and N`): print(`P and t. It also inputs a positive integer`): print(`K (the number of terms for guessing the alg. relation, see below)`): print(`We are interested in the integer sequence`): print(`a(n):=number of n-lettered words in the alphabet {1, ...,m}`): print(`that has as many substrings (alias factors, `): print(`alias consecutive subwords) that are w1 as there are w2.`): print(`The output this procedure is a list consisting`): print(`(i) a (conjectured) recurrence operator ope(n,N) annihilating`): print(`the sequence a(n) together with the initial values `): print(`(starting at n=1) providing`): print(`a full holonomic (P-recursive) description of the sequence`): print(`(ii) The minimal polynomial F(P,x) such that`): print(`F(P(x),x)=0 where P(x) is the ordinary generating function of`): print(`the sequence: i.e.,`): print(`P(x)=Sum(a[n]*x^n,n=0..infinity);`): print(`For example, for Richard Stanley's Amer. Math. Monthly Problem`): print(`11610 (118(10)[Dec. 2011], p. 937)`): print(`type`): print(`Info1(2,[1,1],[1,2],n,N,P,t,40);`): elif nops([args])=1 and op(1,[args])=Info2 then print(`Info2(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m,`): print(`and two words w1 and w2 of the same length in the`): print(`alphabet {1, ..., m}, and symbols n and N`): print(`P and t. It also inputs a positive integer`): print(`K (the number of terms for guessing the alg. relation, see below)`): print(`We are interested in the integer sequence`): print(`a(n):=number of n-lettered words in the alphabet {1, ...,m}`): print(`that has as many substrings (alias factors, `): print(`alias consecutive subwords) that are w1 as there are w2.`): print(`The output this procedure is a list consisting`): print(`(i) a (conjectured) recurrence operator ope(n,N) annihilating`): print(`the sequence a(n) together with the initial values `): print(`(starting at n=1) providing`): print(`a full holonomic (P-recursive) description of the sequence`): print(`(ii) The minimal polynomial F(P,x) such that`): print(`F(P(x),x)=0 where P(x) is the ordinary generating function of`): print(`the sequence: i.e.,`): print(`P(x)=Sum(a[n]*x^n,n=0..infinity);`): print(`(iii) an asymptotic expression for a(n) to order M`): print(`(iv): the ratios of the sequence to the asymptotic expression`): print(`for the terms 1990 to 2000 (just for checking)`): print(`The sequene of ratios of the members of the sequence to the asymptotics`): print(`only for checking purposes`): print(`(vi): a conjectured value for the constant, if avaiable, or else 0`): print(`For example, for Richard Stanley's Amer. Math. Monthly Problem`): print(`11610 (118(10)[Dec. 2011], p. 937)`): print(`type`): print(`Info2(2,[1,1],[1,2],n,N,P,t,40,3);`): elif nops([args])=1 and op(1,[args])=Khalek then print(`Khalek(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) `): print(` 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(` Khalek(n*N^2+n/(n+1)*N,N+1/n,n,N); `): elif nops([args])=1 and op(1,[args])=GWtE then print(`GWtE(m,k,z,t): The global weight-enumerator of words in`): print(`the alphabet {1,..,m} according to the weight z^length`): print(`times the product of t[FactorsOfLength_k]. For`): print(`example, try:`): print(`GWtE(2,2,z,t);`): elif nops([args])=1 and op(1,[args])=Milim then print(`Milim(m,k): all the words of length k in the alphabet {1, ...,m}`): print(`For example, try:`): print(`Milim(2,4);`): elif nops([args])=1 and op(1,[args])=Mishpat then print(`Mishpat(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m,`): print(`and two words w1 and w2 of the same length in the`): print(`alphabet {1, ..., m}, and symbols n and N`): print(`P and t. It also inputs a positive integer`): print(`K (the number of terms for guessing the alg. relation, see below).`): print(`We are interested in the integer sequence`): print(`a(n):=number of n-lettered words in the alphabet {1, ...,m}`): print(`that has as many substrings (alias factors, `): print(`alias consecutive subwords) that are w1 as there are w2.`): print(`The output this procedure is a verbose theorem and sketch of proof`): print(`about the sequence a(n)`): print(`For example, for Richard Stanley's Amer. Math. Monthly Problem`): print(`11610 (118(10)[Dec. 2011], p. 937)`): print(`type: `): print(`Mishpat(2,[1,1],[1,2],n,N,P,t,40,3);`): elif nops([args])=1 and op(1,[args])=ProveGF then print(`ProveGF(m,w1,w2,z,F): Given a positive integer`): print(`m and two words (lists),w1,w2, of the same length in the`): print(`alphabet {1, ...,m}, and an expression F that`); print(`depends on z,`): print(`proves (rigorously!), using`): print(`the continuous Almkvist-Zeilberger algorithm,`): print(`that F is the ordinary generating function (using `): print(`the variable z) of the sequence`): print(`"number of words in {1, ...m} of length n that has`): print(`as many factors (substrings) that are w1 and`): print(`factors that are w2 equals `): print(`such that the constant term of s is the generating`): print(`function for words in {1,..,m}, according to length that`); print(`have the same number of occurences of the factor w1 and w2.`): print(`For example, to completely solve Richard Stanley's`): print(`Amer. Math. Monthly problem #11610 (page 937`): print(` of v. 118 (No. 10), [Dec. 2011]), please type:`): print(`ProveGF(2,[1,1],[1,2],t,1/(1-t)/2+(1+2*t)/sqrt((1-t)*(1-2*t)*(1+t+2*t^2))/2);`): elif nops([args])=1 and op(1,[args])=RECaz then print(`RECaz(m,w1,w2,n,N): Given a positive integer`): print(`m and two words (lists),w1,w2, of the same length in the`): print(`alphabet {1, ...,m}, `): print(`uses the continuous Almkvist-Zeilberger algorithm,`): print(`to find a linear recurrence operator with`): print(`polynomial coefficients, ope(n,N) in terms of n and the`): print(`shift operator in n, N, that annihilates`): print(`the sequence`): print(`"number of words in {1, ...m} of length n that has`): print(`as many factors (substrings) that are w1 and`): print(`factors that are w2 equals `): print(`For example, to get a recurrence operator for`): print(`Richard Stanley's`): print(`Amer. Math. Monthly problem #11610 (page 937`): dprint(` of v. 118 (No. 10), [Dec. 2011]), please type:`): print(`RECaz(2,[1,1],[1,2],n,N);`): elif nops([args])=1 and op(1,[args])=RECg then print(`RECg(m,w1,w2,n,N,MaxC): Given a positive integer`): print(`m and two words (lists),w1,w2, of the same length in the`): print(`alphabet {1, ...,m}, `): print(`uses guessing`): print(`to find a linear recurrence operator with`): print(`polynomial coefficients, ope(n,N) in terms of n and the`): print(`shift operator in n, N`): print(`such that the degree plus order is <=MaxC ,`): print(` that annihilates the sequence`): print(`"number of words in {1, ...m} of length n that has`): print(`as many factors (substrings) that are w1 and`): print(`factors that are w2 equals `): print(`It also gives the initial values.`): print(`For example, to get a recurrence operator for`): print(`Richard Stanley's`): print(`Amer. Math. Monthly problem #11610 (page 937`): dprint(` of v. 118 (No. 10), [Dec. 2011]), please type:`): print(`RECg(2,[1,1],[1,2],n,N,10);`): elif nops([args])=1 and op(1,[args])=RSct then print(`RSct(m,w1,w2,z,s): The rational function in z and s`): print(`such that the constant term of s is the generating`): print(`function for words in {1,..,m}, according to length that`): print(`have the same number of occurences of the factor w1 and w2.`): print(`For example, try:`): print(`RSct(2,[1,1],[1,2],z,s);`): elif nops([args])=1 and op(1,[args])=RSseq then print(`RSseq(m,w1,w2,N): Given a pos. integer m and and`): print(`two words of the same length w1,w2 in the alphabet `): print(`{1,...,m}, finds the first N terms of the sequence`): print(`enumerating the number of n-letter words in the`): print(`alphabet {1,...,m} that have the same number of `): print(`factors that are w1 and that are factors that are w2`): print(`For example, try:`): print(`RSseq(2,[1,1],[1,2],10);`): elif nops([args])=1 and op(1,[args])=SeqFromRec then print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`): print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`extends it to the first K values`): print(`For example, try:`): print(`SeqFromRec(N-2,n,N,[1],10);`): elif nops([args])=1 and op(1,[args])=TerseBook then print(`TerseBook(m,k,n,N,P,t,K): inputs positive integers m and k`): print(`and outputs a set of lists consiting of `): print(`the input and output of Info1(m,k,n,N,P,t,K)`): print(`(qv) for all members of Zugot(m,k) (qv).`): print(`For example, try:`): print(`TerseBook(2,2,n,N,P,t,40);`): elif nops([args])=1 and op(1,[args])=VerboseBook then print(`VerboseBook(m,k,n,P,t,K,M): inputs positive integers m and k`): print(`and outputs a web-book with generating functions`): print(`and recurrences for sequences enumerating`): print(`analogous quantities to Richard Stanley's Monthly`): print(`Problem 11610 (118(10)[Dec. 2011], p. 937)`): print(`concerning all pairs of words in {1, ...,m} `): print(`of length k (up to equivalence). `): print(`m is the size of the alphabet, k is the length of the`): print(`distinguished words, n,P,t is the symbol for displaying`): print(`K and M are as in Info2 `): print(`For example try:`): print(`VerboseBook(2,2,n,P,t,40,4);`): elif nops([args])=1 and op(1,[args])=Wt then print(`Wt(w,k,z,t): The weight of the word w in terms of z and`): print(`the indexed variables t by factors of length k. For`): print(`example, try:`): print(`Wt([1,2,3,2,1],2,z,t);`): elif nops([args])=1 and op(1,[args])=WtE then print(`WtE(m,k,z,t,S): The weight-enumerator of words in`): print(`the alphabet {1,..,m} according to the weight z^length`): print(`times the product of t[FactorsOfLength_kThatAreInTheSetS]. For`): print(`example, try:`): print(`WtE(2,2,z,t,{[1,1],[1,2]});`): elif nops([args])=1 and op(1,[args])=WtEd then print(`WtEd(m,k,z,t): Like WtE(m,k,z,t,S) `): print(`but done directly. For checking purposes only. For example, try:`): print(`WtEd(2,2,z,t,{[1,1],[1,2]});`): elif nops([args])=1 and op(1,[args])=Zugot then print(`Zugot(m,k): Representatives of sets of pairs of words of`): print(`length k in the alphabet {1, ..., m}. For example, try:`): print(`Zugot(2,2);`): else print(`There is no ezra for`,args): fi: end: #########Begin AsyRec######### ##From AsyRec with(SolveTools): with(numtheory): #MonoN(ope,n,N,k): Given a linear recurrence operator with #poly coeffs. ope(n,N), and a pos. integer k, outputs the expression of #N^k as a linear combination of 1, N, ..., N^(ORDER-1) #For example, try #MonoN(N-n-1,n,N,3); MonoN:=proc(ope,n,N,k) local ORDER,coe0,i,lu1,lu2: ORDER:=degree(ope,N): if k=1 then RETURN(FAIL): fi: if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #HafelOper(ope,n,N,L): applies the operator ope(n,N) to #the sequence L, for example, try: #HafelOper(N-(n+1),n,N,[1,2,6,24,120]); HafelOper:=proc(ope,n,N,L) local i,n1: [seq(add(subs(n=n1,coeff(ope,N,i))*L[n1+i],i=0..degree(ope,N)), n1=1..nops(L)-degree(ope,N))]: end: #TestKsect(ope,n,N,K,r,M,Ini): tests procedure Ksect #with initial conditions Ini up to M terms TestKsect:=proc(ope,n,N,K,r,M,Ini) local gu,Ope,n1: gu:=SeqFromRec(ope,n,N,Ini,M*K+r): if gu=FAIL then RETURN(FAIL): fi: Ope:=Ksect(ope,n,N,K,r): gu:=[seq(gu[K*n1+r],n1=1..M-1)]: evalb({op(HafelOper(Ope,n,N,gu))}={0}): end: rf:=proc(a,n) local i:mul(a+i,i=0..n-1):end: #CODV1(ope,n,N,k): Given a linear recurrence operator ope(n,N) #annihilating a[n], say, outputs the operator #annihilating b[n]:=a[n]/n!^k #For example, try: CODV1(N-(n+1),n,N,1): CODV1:=proc(ope,n,N,k) local i: add(coeff(ope,N,i)*(rf(n+1,i))^k*N^i,i=0..degree(ope,N)): end: #CODV(ope,n,N,k,K): Given a linear recurrence operator ope(n,N) #annihilating a[n], say, outputs the operator #annihilating b[n]:=a[n*K]/n!^k #For example, try: CODV(N-(n+1),n,N,1,1): CODV:=proc(ope,n,N,k,K) local Ope: Ope:=Ksect(ope,n,N,K,0): CODV1(Ope,n,N,k): end: #FindKk(ope,n,N): Given a linear recurrence operator with polynomial #coefficients, ope(n,N), finds the integers K and k such #that if a(n) is a solution of ope(n,N)a(n)=0, then #b(n):=a(K*n)/n!^k is annihilated by a standard operator #the output is the pair [k,K] and the transformed operator #For exanple, try: FindKk(N^2-n,n,N); FindKk:=proc(ope,n,N) local ope1,Ld,deg,sp,yakhas,K,k,pu: pu:=Aope1(ope,n,N): if {solve(pu,N)}<>{} and {solve(pu,N)}<>{0} then RETURN([1,0],ope): fi: ope1:=expand(ope): Ld:=ldegree(ope,N): ope1:=expand(subs(n=n-Ld,ope1)/N^Ld ): deg:=degree(ope,N): sp:=degree(coeff(ope,N,0),n)-degree(coeff(ope,N,deg),n): yakhas:=sp/deg: k:=numer(yakhas): K:=denom(yakhas): [K,k],CODV(ope,n,N,k,K): end: #Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N) Aope1:=proc(ope,n,N) local gu: gu:=expand(numer(normal(ope))): gu:=coeff(gu,n,degree(gu,n)): gu:=expand(gu/coeff(gu,N,degree(gu,N))): factor(gu): end: C3:=proc(a,b,c) option remember: if [a,b,c]=[0,0,0] then RETURN(1): fi: if a<0 or b<0 or c<0 then RETURN(0): fi: C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1): end: #Aluf1(pit): is it positive dominant? Aluf1:=proc(pit) local aluf,shi,i,makom: shi:=evalf(abs(pit[1])): aluf:={evalf(evalc(pit[1]))}: makom:={1}: for i from 2 to nops(pit) do if evalf(abs(pit[i]))>evalf(shi) then shi:=abs(evalf(evalc(pit[i]))): aluf:={pit[i]}: makom:={i}: elif evalf(abs(pit[i]))=shi then aluf:=aluf union {evalf(evalc(pit[i]))}: makom:=makom union {i}: fi: od: aluf,shi,makom: end: #OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepAS1:=proc(ope1,n,N,alpha,f,S1) local x,f1,L,F,A,mu,i,A1: f1:=subs(n=1/x,f): L:=degree(f1,x): F:=f1+A*x^(L+S1): mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F) ,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): mu:=taylor(mu,x=0,L+11): mu:=simplify(mu): for i from L+1 to L+9 while coeff(mu,x,i)=0 do od: if i=L+10 then RETURN(FAIL): fi: mu:=coeff(mu,x,i): A1:=subs(Linear({mu},{A}),A): if A1=0 then RETURN(FAIL): fi: subs({x=1/n,A=A1},F): end: #OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepA:=proc(ope1,n,N,alpha,f) local S1,pu: for S1 from 1 to 5 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od: pu:=OneStepAS1(ope1,n,N,alpha,f,S1): if pu=FAIL then RETURN(f): else RETURN(pu): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #Nor(ope,n,N): the Normalizer of the linear recurrence #operator with polynomial coefficients ope(n,N) #followed by its exponential growth constant #For example, try: #Nor(N^2-3*N+2,n,N); Nor:=proc(ope,n,N) local gu,pit,lu,makom,ope1: gu:=Aope1(ope,n,N): if degree(gu,N)=0 then print(`Not of exponential growth, case not implemented`): RETURN(FAIL): fi: if not type(expand(normal(ope/gu)),`+`) then pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: lu:=pit[makom[1]]: ope1:=Yafe(subs(N=lu*N,ope),N)[2]: RETURN(ope1,lu): fi: pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: if nops(makom)<>1 then lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then RETURN(FAIL): elif pit[makom[1]]+pit[makom[2]]=0 then RETURN(FAIL): fi: fi: lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then RETURN(FAIL): fi: lu:=pit[makom[1]]: ope1:=Yafe(subs(N=lu*N,ope),N)[2]: ope1,lu: end: #Atom(s,r,i,x,K): Expanding (n+i)^(s/r)-n^(s/r) in terms of #x=n^(-1/r) up to K terms #For example try: #Atom(2,3,2,x,3); Atom:=proc(s,r,i,x,K) local gu,y,i1: gu:=(1+i*y)^(s/r)-1: gu:=taylor(gu,y=0,K+1): gu:=add(coeff(gu,y,i1)*y^i1,i1=0..K): expand(subs(y=x^r,gu)/x^s): end: #FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics #in terms of x=n^(1/r) #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x); FindExpP1old:=proc(ope,n,N,r,x) local c,eq,var,s,sof,gu,ope1,deg,i,i1,v: sof:=add(c[s]*x^s,s=1..r-1): var:={seq(c[i],i=1..r-1)}: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=expand(subs(n=1/x^r,ope1)): gu:=0: for i from 0 to degree(ope1,N) do gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..r-1) ): od: gu:=taylor(gu,x=0,2*r-1): eq:={seq(coeff(gu,x,i1),i1=r..2*r-2)}: var:=[solve(eq,var)][1]: for v in var do if op(1,v)=op(2,v) then RETURN(FAIL): fi: od: subs(var,sof): end: #Findr(ope,n,N): Given an operator ope(n,N) #such that the leading terms in n is a multiple of #(N-1), finds the largest r such that it is # a multiple of (N-1)^r. For example, try: #Findr((N-1)^4,n,N); Findr:=proc(ope,n,N) local ope1,r: ope1:=coeff(expand(ope),n,degree(ope,n)): if expand(subs(N=1,ope1))<>0 then RETURN(FAIL): fi: for r from 1 while expand(subs(N=1,diff(ope1,N$r)))=0 do od: r: end: FindExpP:=proc(ope,n,N) local r,x,lu: r:=Findr(ope,n,N): if r=FAIL then RETURN(FAIL): fi: lu:=FindExpP1(ope,n,N,r,x) : if lu=FAIL then RETURN(FAIL): fi: subs(x=n^(1/r),lu): end: #NewOpe(ope,n,N,K): the exponential part + the transformed operator #(up to degree K asymp. in the coefficients) #For example, try: #NewOpe((n+1)*N^2-2*(n+5)*N+n+3,n,N,4); NewOpe:=proc(ope,n,N,K,x) local r,lu,lu1,ope1,deg,s,i,gu,mu: r:=Findr(ope,n,N): if r=FAIL then RETURN(FAIL): fi: lu:=FindExpP1(ope,n,N,r,x) : lu1:=FindExpP(ope,n,N) : if lu=FAIL then RETURN(FAIL): fi: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=subs(n=1/x^r,ope1): ope1:=expand(ope1): gu:=0: for i from 0 to degree(ope1,N) do mu:=0: for s from 1 to r-1 do mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3): od: gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i: od: gu:=taylor(gu,x=0,K+3): lu1,add(coeff(gu,x,i)*x^i,i=0..K): end: #Bdok(ope,n,N,K): the exponential part + the transformed operator #(up to degree K asymp. in the coefficients) #For example, try: #Bdok((n+1)*N^2-2*(n+5)*N+n+3,n,N,4); Bdok:=proc(ope,n,N,K) local r,x,lu,lu1,ope1,deg,s,i,gu,mu: r:=Findr(ope,n,N): if r=FAIL then RETURN(FAIL): fi: lu:=FindExpP1(ope,n,N,r,x) : lu1:=FindExpP(ope,n,N) : if lu=FAIL then RETURN(FAIL): fi: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=subs(n=1/x^r,ope1): ope1:=expand(ope1): gu:=0: for i from 0 to degree(ope1,N) do mu:=0: for s from 1 to r-1 do mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3): od: gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i: od: gu:=taylor(gu,x=0,K+3): lu1,add(coeff(gu,x,i)*x^i,i=0..K): end: OpePerN:=proc(N,n,r) local i,j,ope: ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r): ope:=expand(N^r*subs(n=n+r,ope)): ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r): ope:=expand(FindKk(ope,n,N)[2]): ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r): ope:=numer(Nor(ope,n,N)[1]): ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r): end: OpePer:=proc(N,n,r) local i,j,ope: ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r): ope:=expand(N^r*subs(n=n+r,ope)): end: OpePerG:=proc(N,n,S) local i,j,ope,r,gu,t: ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i in S): r:=max(op(S)): ope:=expand(N^r*subs(n=n+r,ope)): gu:=exp(add(t^i/i,i in S)): gu:=taylor(gu,t=0,degree(ope,N)+2): ope,[seq(coeff(gu,t,i)*i!,i=1..degree(ope,N))]: end: #Finda(ope,N,x,r): finds the first power x^a in the #asymptotic solution of ope(N,n)f(n)=0, where x=1/n^(1/r) #For example, try: #Finda((1+x)-(1+3*x)*N,N,x,1); Finda:=proc(ope,N,x,r) local gu,i,a,a1: gu:=add(coeff(ope,N,i)*(1+i*x^r)^(a/r),i=0..degree(ope,N)): gu:=taylor(gu,x=0,5*r+5): gu:=expand(add(coeff(gu,x,i)*x^i,i=0..5*r+4)): for i from 0 to 5*r+4 while expand(simplify(coeff(gu,x,i)))=0 do od: if i=5*r+5 then RETURN(FAIL): fi: a1:=[solve(coeff(gu,x,i),a)][1]: if a1=NULL then RETURN(FAIL): elif not type(a1,integer) then RETURN(FAIL): else RETURN(-a1): fi: end: #OneStepG(ope,N,x,r,Cu): Given a partial #asympt. expansion, in terms of x=1/n^(1/r), #for solutions of the recurrence equation #ope(n,N)a(n)=0, where ope is normalized of type #r (i.e. its leading coeff. is a multiple of (N-1)^r) #finds the next term #OneStepG((1+x)-(1+3*x)*N,N,x,1,x^2); OneStepG:=proc(ope,N,x,r,Cu) local gu,i,a,c,c1,Cu1,K,ka: K:=degree(Cu,x): Cu1:=Cu+c*x^(K+1): gu:=add(coeff(ope,N,i)* add(coeff(Cu1,x,a)*x^a*(1+i*x^r)^(-a/r),a=ldegree(Cu1,x)..degree(Cu1,x)), i=0..degree(ope,N)): gu:=expand(gu): gu:=taylor(gu,x=0,5*r+8+K): gu:= simplify(expand(add(coeff(gu,x,i)*x^i,i=0..5*r+7+K))): gu:=pashet(gu,x,5*r+7+K,c): ka:=expand(simplify(coeff(gu,x,ldegree(gu,x)))): ka:=simplify(ka): c1:=[solve(simplify(coeff(gu,x,ldegree(gu,x))),c)][1]: if c1=NULL then RETURN(FAIL): else RETURN(subs(c=c1,Cu1)): fi: end: #Asy1(ope,n,N,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the K's term Asy1:=proc(ope,n,N,K) local gu,lu,alpha,mu,ope1,ku,i,f,x,ka,vu,ope1A,r,Kk,pu,a: gu:=FindKk(ope,n,N): Kk:=gu[1]: ope1:=gu[2]: vu:=Nor(ope1,n,N): if vu=FAIL then RETURN(FAIL): fi: ope1:=vu[1]: lu:=vu[2]: ope1A:=Aope1(ope1,n,N): for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od: if r=1 then ###r=1 case mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,2),x,1)): ka:=simplify([solve(ku,alpha)]): if coeff(ka[1],I,1)<>0 then RETURN(FAIL): fi: alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN(lu^n*n^alpha,Kk): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: if f=FAIL then RETURN(FAIL): fi: RETURN(lu^n*n^alpha*f,Kk): #end r=1 case else ope1:=numer(ope1): if expand(diff(ope1,n))=0 then ka:=n^(r-1)*lu^n: RETURN(ka,Kk): fi: pu:=FindExpP(ope1,n,N): if pu=FAIL then RETURN(FAIL): fi: ope1:=NewOpe(ope1,n,N,K+5,x)[2]: a:=Finda(ope1,N,x,r): if a=FAIL then RETURN(FAIL): fi: ka:=x^a: for i from a to K do ka:=OneStepG(ope1,N,x,r,ka): od: ka:=exp(pu)*subs(x=1/n^(1/r),ka)*lu^n: RETURN(ka,Kk): fi: end: pashet:=proc(gu,x,K,c) local gu1,i,pu: Digits:=300: gu1:=0: for i from 0 to K do pu:=evalf(coeff(gu,x,i)): if not (abs(coeff(pu,c,1))<10^(-20) and abs(coeff(pu,c,0))<10^(-20) ) then gu1:=gu1+coeff(gu,x,i)*x^i: fi: od: gu1: end: #NakedStirling(n,K): the asymptotic expansion of #n!/((n/e)^n*sqrt(2*Pi*n). For example, try: #NakedStirling(n,5); NakedStirling:=proc(n,K) local gu,i: gu:=n!/(n/exp(1))^n/sqrt(2*Pi*n): gu:=expand(asympt(gu,n,K+1)): (n/exp(1))^n*sqrt(2*Pi*n), add(coeff(op(i,gu),n,-(i-1))/n^(i-1),i=1..K+1): end: #AsyF(ope,n,N,M): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the M's term, in terms of the factorial function AsyF:=proc(ope,n,N,M) local K,k,gu: gu:=Asy1(ope,n,N,M): if gu=FAIL then RETURN(FAIL): fi: K:=gu[2][1]: k:=gu[2][2]: gu:=gu[1]: (n/K)!^k*subs(n=n/K,gu): end: #AsyFC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0, with the given initial conditions #where ope(n,N) is a recurrence operator #up to the M's term, in terms of the factorial function #and complete with an empirically derived constant in front #using K terms AsyFC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1: Digits:=100: gu:=AsyF(ope,n,N,M): L:=SeqFromRec(ope,n,N,evalf(Ini),K): if L=FAIL then RETURN(FAIL): fi: mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]: er:=mu[nops(mu)]-mu[nops(mu)-1]: D1:=-trunc(log(abs(er))/log(10))-3: if D1<2 then print(`can't determine the constant`): RETURN(gu): fi: C:=evalf(mu[nops(mu)],D1): C,gu: end: #Asy1special(ope,n,N,K,x): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,given as a list #[pu,lu,expansion,r] where it is #exp(pu)*lu^n*expansion(x), where x=1/n^(1/r), and r is #a positive integer. It also returns [K,k] (see Asy1) #where ope(n,N) is a recurrence operator #up to the K's term Asy1special:=proc(ope,n,N,K,x) local gu,lu,alpha,mu,ope1,ku,i,f,ka,vu,ope1A,r,Kk,pu,a: gu:=FindKk(ope,n,N): Kk:=gu[1]: ope1:=gu[2]: vu:=Nor(ope1,n,N): if vu=FAIL then RETURN(FAIL): fi: ope1:=vu[1]: lu:=vu[2]: ope1A:=Aope1(ope1,n,N): for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od: if degree(ope1,n)=0 then RETURN([0,lu,x^(-r+1),1,1],Kk): fi: if r=1 then ###r=1 case mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,2),x,1)): ka:=simplify([solve(ku,alpha)]): if ka=[] then RETURN(FAIL): fi: if coeff(ka[1],I,1)<>0 then RETURN(FAIL): fi: alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN([0,lu,x^(-alpha),1,1],Kk): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: if f=FAIL then RETURN(FAIL): fi: f:=subs(n=1/x,f): RETURN([0,lu,x^(-alpha),f,1],Kk): #end r=1 case else ope1:=numer(ope1): if expand(diff(ope1,n))=0 then ka:=n^(r-1)*lu^n: RETURN(ka,Kk): fi: pu:=FindExpP(ope1,n,N): if pu=FAIL then RETURN(FAIL): fi: ope1:=NewOpe(ope1,n,N,K+5,x)[2]: a:=Finda(ope1,N,x,r): if a=FAIL then RETURN(FAIL): fi: ka:=x^a: for i from a to K do ka:=OneStepG(ope1,N,x,r,ka): od: RETURN([pu,lu,1,ka,r],Kk): fi: end: #IsMispar(a): is a (generalized) numeric type IsMispar:=proc(a) : if type(a,numeric) then true: elif type(a,`^`) and type(op(1,a),numeric) and type(op(2,a),numeric) then true: elif a=Pi then true: elif type(a,`^`) and op(1,a)=Pi and type(op(2,a),numeric) then true: elif type(a,function) and type(op(1,a),numeric) then true: else false: fi: end: #GetRidOfConst(P): gets rid of the constant in front GetRidOfConst:=proc(P) local i,P1: if not type(P,`*`) then RETURN(P): fi: P1:=1: for i from 1 to nops(P) do if not IsMispar(op(i,P)) then P1:=P1*op(i,P): fi: od: P1: end: #Asy(ope,n,N,M): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the M's term Asy:=proc(ope,n,N,M) local K,k,gu,vu,vu1,vu2,pu,lu,ka1,ka2,r,x,i,vu2a: gu:=Asy1special(ope,n,N,M,x): if gu=FAIL then RETURN(FAIL): fi: K:=gu[2][1]: k:=gu[2][2]: gu:=gu[1]: pu:=subs(n=n/K,gu[1]): lu:=gu[2]: ka1:=gu[3]: ka2:=gu[4]: r:=gu[5]: ka2:=subs(x=K^(1/r)*x,ka2): vu:=NakedStirling(n,M+1): vu1:=vu[1]: vu2:=vu[2]: vu1:=subs(n=n/K,vu1)^k: vu2:=subs(n=n/K,vu2)^k: vu2:=expand(subs(n=1/x^r,vu2)): vu2:=expand(ka2*vu2): vu2:=add(simplify(coeff(vu2,x,i))*x^i,i=0..M*r): vu2:=subs(x=1/n^(1/r),vu2): ka1:=subs(x=1/n^(1/r),ka1): ka1:=subs(n=n/K,ka1): vu2a:=op(1,vu2): vu2:=expand(normal(vu2/vu2a)): gu:=simplify(vu2a*vu1*lu^(n/K))*exp(pu)*ka1*vu2: GetRidOfConst(gu): end: #AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0, with the given initial conditions #where ope(n,N) is a recurrence operator #up to the M's term, #and complete with an empirically derived constant in front #using K terms AsyC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1,gu1: Digits:=100: gu:=Asy(ope,n,N,M): if gu=FAIL then RETURN(FAIL): fi: L:=SeqFromRec(ope,n,N,evalf(Ini),K): if L=FAIL then RETURN(FAIL): fi: mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]: er:=mu[nops(mu)]-mu[nops(mu)-1]: if er=0 then D1:=Digits: else D1:=-trunc(evalf(log(abs(er))/log(10)))-3: fi: if D1<2 then print(`can't determine the constant`): RETURN(FAIL): fi: C:=evalf(mu[nops(mu)],D1): if C<1/10000 then RETURN(FAIL): fi: gu1:=C*gu: if not CheckAsy(gu1,n,L,1/100) then RETURN(FAIL): fi: C,gu,D1: end: #FindExpP1g(ope,n,N,r,x,ds): finds the exponential part of the asymptotics #in terms of x=n^(1/r) as a poly. of degree ds in x #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1g((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1); FindExpP1g:=proc(ope,n,N,r,x,ds) local c,eq,var,s,sof,gu,ope1, deg,i,i1,v,ka,i11: sof:=add(c[s]*x^s,s=1..ds): var:={seq(c[i],i=1..ds)}: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=expand(subs(n=1/x^r,ope1)): gu:=0: for i from 0 to degree(ope1,N) do gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ): od: gu:=taylor(gu,x=0,3*r+10): for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od: eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1)}: eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds)}: lprint(eq,var): ka:=[solve(eq,var)]: if ka=[] then RETURN(FAIL): fi: var:=ka[1]: for v in var do if op(1,v)=op(2,v) then RETURN(FAIL): fi: od: subs(var,sof): end: #FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics #in terms of x=n^(1/r) #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x); FindExpP1:=proc(ope,n,N,r,x) local s,gu,nosaf: for s from 1 to r-1 do for nosaf from 0 to 1 do gu:=FindExpP1gExtra(ope,n,N,r,x,s,nosaf): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: #FindExpP1gExtra(ope,n,N,r,x,ds,nosaf) #: finds the exponential part of the asymptotics #in terms of x=n^(1/r) as a poly. of degree ds in x #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1gExtra((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1); FindExpP1gExtra:=proc(ope,n,N,r,x,ds,nosaf) local c,eq,var,s,sof,gu,ope1, deg,i,i1,v,ka,i11,ku,varf: sof:=add(c[s]*x^s,s=1..ds): var:={seq(c[i],i=1..ds)}: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=expand(subs(n=1/x^r,ope1)): gu:=0: for i from 0 to degree(ope1,N) do gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ): od: gu:=taylor(gu,x=0,3*r+10): for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od: eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1+nosaf)}: ka:=[solve(eq,var)]: if ka=[] then RETURN(FAIL): fi: var:=ka[1]: for v in var do if op(1,v)=op(2,v) then RETURN(FAIL): fi: od: #varf:=evalf(var): #ku:=[seq(subs(varf,c[s]),s=1..ds)]: #print(`ku is`, ku): #add(ku[s]*x^s,s=1..ds): subs(var,sof): end: #########End from AsyRec######### ###Begin from SCHUTZENBERGER##### 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(F,P,x,D1): Inputs a polynomial F, #in P and x, and outputs #a differential operator annihilating #P(x) where P(x) is the alg. formal power series #satisfying the algebraic equation #F(P(x),x))=0 . For example, try: #algtodiff(P-1-x*P^2,P,x); algtodiff:=proc(F,P,x,D1) local ORDER,gu: for ORDER from 1 do gu:=algtodiff1(F,P,x,D1,ORDER): if gu<>FAIL then gu:=add(factor(coeff(gu,D1,i))*D1^i,i=0..ORDER): RETURN(gu): fi: od: end: #algtodiff1(F,P,x,D1,ORDER): Inputs a polynomial F, #in P and x, and an ORDER, and outputs #a differential operator annihilating #P(x) where P(x) is the alg. formal power series #satisfying the algebraic equation #F(P(x),x))=0 . For example, try: #algtodiff1(P-1-x*P^2,P,x,2); algtodiff1:=proc(F,P,x,D1,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]*D1: 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]*D1^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): if ope=0 then RETURN(FAIL): fi: ope:=expand(ope): pip:=expand(coeff(ope,D1,degree(ope,D1))): pip:=coeff(pip,x,degree(pip,x)): normal(ope/pip): 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,ka,mu1,K,richard,ka1: K:=(1+degx)*(1+degP)+6: if K > nops(gu) then print(`sequence too small`): RETURN(FAIL): 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,K+1): eq:={}: for i1 from 0 to K do eq:=eq union {coeff(lu,x,i1)=0} od: mu:=solve(eq,var): ka:={}: for mu1 in mu do if op(1,mu1)=op(2,mu1) then ka:=ka union {op(1,mu1)}: fi: od: F:=subs(mu,F): if F=0 then RETURN(0): fi: richard:={seq(factor(coeff(F,ka1,1)) ,ka1 in ka)}: if nops(richard)>1 then RETURN(0): else RETURN(richard[1]): fi: end: #Empir(L,x,P): inputs a sequence L and #and conjectures an polynomial equation #F(x,P)=0 satisfies by the formal power #series L[1]+L[2]*x+L[3]*x^2+... #For example, try: #Empir([seq(i1,i1=1..40)],x,P); Empir:=proc(L,x,P) local degx,degP,lu,L1,i: for L1 from 1 to (nops(L)-3)/3 do for degP from 1 to L1 do for degx from 0 to min(trunc((nops(L)-6)/(1+degP))-1,L1-degP) do lu:=empir(L,degx,degP,x,P): if lu<>0 then RETURN(add(factor(coeff(lu,P,i))*P^i,i=0..degP)): fi: od: od: od: 0: end: ###End from SCHUTZENBERGER##### ###From Findrec######### #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecVerbose:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then print(`There is some slack, there are `, nops(ope)): print(ope): RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a,ope1: option remember: if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then RETURN(FAIL): fi: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then RETURN(FAIL): elif nops(ope)=1 then if HSP(ope[1],n,N)>=1 then RETURN(FAIL): else ope1:=Yafe(ope[1],N)[2]: if not CheckRec(f,ope1,n,N) then RETURN(FAIL): else RETURN(ope1): fi: fi: else RETURN(FAIL): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: ################# #Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with #poly coffs. #of maximum DEGREE+ORDER<=MaxC #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); Findrec:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 0 to L do DEGREE:=L-ORDER: if (2+DEGREE)*(1+ORDER)+4>=nops(f) then print(`Insufficient data for degree`, DEGREE, `and order `,ORDER): RETURN(FAIL): fi: ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #SeqFromRecOld(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRecOld:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #end Findrec with(linalg): #findrecEx(f,DEGREE,ORDER,m1): Explores whether thre #is a good chance that there is a recurrence of degree DEGREE #and order ORDER, using the prime m1 #For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003); findrecEx:=proc(f,DEGREE,ORDER,m1) local ope,var,eq,i,j,n0,eq1,a,A1, D1,E1,Eq,Var,f1,n,N: option remember: f1:=f mod m1: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1: od: eq:= eq union {eq1}: od: Eq:= convert(eq,list): Var:= convert(var,list): D1:=nops(Var): E1:=nops(Eq): if E10 then RETURN(false): fi: if E1-D1>=1 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+1],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: if E1-D1>=2 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+2],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: true: end: ###End From Findrec######### ###Begin From EKHAD pashetc:=proc(p,D) local gu,p1,i,mekh: p1:=normal(p): mekh:=denom(p1): p1:=numer(p1): p1:=expand(p1): gu:=0: for i from 0 to degree(p1,D) do gu:=gu+factor(coeff(p1,D,i))*D^i: od: RETURN(gu,mekh): end: gzor:=proc(f,x,n) local i,gu: gu:=f: for i from 1 to n do gu:=diff(gu,x): od: RETURN(gu): end: goremc:=proc(D,ORDER) local i,gu: gu:=0: for i from 0 to ORDER do gu:=gu+(bpc[i])*D^i: od: RETURN(gu): end: duisc:= proc(A,ORDER,y,x,D) local KAMA,gorem,conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu, mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,i,shad: KAMA:=50: gorem:=goremc(D,ORDER): conj:=gorem: yakhas:=0: for i from 0 to degree(conj,D) do yakhas:=yakhas+normal(coeff(conj,D,i)*gzor(A,x,i)/A): yakhas:=normal(yakhas): od: lu1a:=yakhas: LU1:=numer(yakhas): yakhas:=1/denom(yakhas): yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A): P:=LU1: Q:=numer(yakhas): R:=denom(yakhas): j1:=0: while j1 <= KAMA do g:=gcd(R,Q-j1*diff(R,y)): if g <> 1 then Q2:=(Q-j1*diff(R,y))/g: R:=normal(R/g): Q:=normal(Q2+j1*diff(R,y)): P:=P*g^j1: j1:=-1: fi: j1:=j1+1: od: P:=expand(P): R:=expand(R): Q:=expand(Q): Q1:=Q+diff(R,y): Q1:=expand(Q1): l1:=degree(R,y): l2:=degree(Q1,y): meka1:=coeff(R,y,l1): mekb1:=coeff(Q1,y,l2): if l1-1 <>l2 then k1:=degree(P,y)-max(l1-1,l2): else mumu:= -mekb1/meka1: if type(mumu,integer) and mumu > 0 then k1:=max(mumu, degree(P,y)-l2): else k1:=degree(P,y)-l2: fi: fi: fu:=0: if k1 < 0 then RETURN(0): fi: for ia1 from 0 to k1 do fu:=fu+apc[ia1]*y^ia1: od: gugu:=Q1*fu+R*diff(fu,y)-P: gugu:=expand(gugu): degg:=degree(gugu,y): for ia1 from 0 to degg do eq[ia1+1]:=coeff(gugu,y,ia1)=0: od: for ia1 from 0 to k1 do va1[ia1+1]:=apc[ia1]: od: for ia1 from 0 to ORDER do va2[ia1+1]:=bpc[ia1]: od: eq:=convert(eq,set): va1:=convert(va1,set): va2:=convert(va2,set): va1:=va1 union va2 : va1:=solve(eq,va1): fu:=subs(va1,fu): gorem:=subs(va1,gorem): if fu=0 and gorem=0 then RETURN(0): fi: for ia1 from 0 to k1 do gorem:=subs(apc[ia1]=1,gorem): od: for ia1 from 0 to ORDER do gorem:=subs(bpc[ia1]=1,gorem): od: fu:=normal(fu): shad:=pashetc(gorem,D): S:=lu1a*R*fu*shad[2]/P: S:=subs(va1,S): S:=normal(S): S:=factor(S): for ia1 from 0 to k1 do S:=subs(apc[ia1]=1,S): od: for ia1 from 0 to ORDER do S:=subs(bpc[ia1]=1,S): od: shad[1],S: end: Mygcd:=proc(L) local i,L1: if nops(L)=1 then RETURN(1): fi: if nops(L)=2 then RETURN(gcd(L[1],L[2])): fi: L1:=[op(1..nops(L)-1,L)]: gcd(Mygcd(L1),L[nops(L)]): end: AZc:=proc(A,y,x,D) local ORDER,gu,KAMA,ope,C,kak,i1: KAMA:=50: for ORDER from 1 to KAMA do gu:=duisc(A,ORDER,y,x,D): if gu[1]<>0 then ope:=gu[1]: C:=gu[2]: kak:=Mygcd([seq(coeff(ope,D,i1),i1=0..ORDER)]): if kak<>1 then ope:=add(normal(coeff(ope,D,i1)/kak)*D^i1,i1=0..ORDER): C:=normal(C/kak): fi: RETURN(ope,C): fi: od: 0: end: ###End From EKHAD ####Start DiffToREc1 #DiffToRec1(ope,t,Dt,n,N): Given a linear differential operator #ope(t,Dt) and symbols n,N, outputs the linear recurrence #operator with polynomial coefficients satisfied by #the sequence of coeff. of any f.p.s. annihilating ope(t,Dt) #where N is the shift in n. For example, try: #DiffToRec1(Dt+t, t,Dt,n,N); DiffToRec1:=proc(ope,t,Dt,n,N) local gu,i,j,ope1,d,j1: ope1:=expand(ope): gu:=0: for i from ldegree(ope1,t) to degree(ope1,t) do for j from 0 to degree(ope1,Dt) do gu:=gu+coeff(coeff(ope1,t,i),Dt,j)*mul(n-i+j1,j1=1..j)*N^(j-i): od: od: d:=ldegree(gu,N): gu:=N^(-d)*subs(n=n-d,gu): Yafe(gu,N)[2]: 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: ####End DiffToRec1 #Milim(m,k): all the words of length k in the alphabet {1, ...,m} #For example, try: #Milim(2,4); Milim:=proc(m,k) local mu,i,w: if k=0 then RETURN({[]}): fi: mu:=Milim(m,k-1): {seq(seq([op(w),i],i=1..m),w in mu)}: end: #EqsVars(m,k,z,t,A):The set of equations, followed by the #set of variables that enables to compute the generating #function in z and t[w1] (w1 in [1,...,m]^k) #that is the weight-enumerator of all words in #the alpbhabet {1, ...,m} with weight z^length(w) times #product of t[w_i,...,w_i+k-1] for i from 1 to nops(w)-k+1 #A is a symbol for the indexing of the weight enumerator #of words, #For example, try #EqsVars(2,2,z,t,A); EqsVars:=proc(m,k,z,t,A) local eq,var,gu,eq1,i,w: gu:=Milim(m,k-1): var:={}: eq:={}: for w in gu do var:=var union {A[op(w)]}: eq1:=A[op(w)]-z^(k-1)-z*add(t[op(w),i]*A[op(2..nops(w),w),i],i=1..m): eq:=eq union {eq1}: od: eq,var: end: #GWtE(m,k,z,t): The global weight-enumerator of words in #the alphabet {1,..,m} according to the weight z^length #times the product of t[FactorsOfLength_k]. For #example, try: #GWtE(2,2,z,t); GWtE:=proc(m,k,z,t) local gu,A,mu,lu,i1,mu1: mu:=EqsVars(m,k,z,t,A): gu:=solve(mu): lu:=add(subs(gu,mu1),mu1 in mu[2]): lu:=lu+add(m^i1*z^i1,i1=0..k-2): normal(lu): end: #Wt(w,k,z,t): The weight of the word w in terms of z and #the indexed variables t by factors of length k. For #example, try: #Wt([1,2,3,2,1],2,z,t); Wt:=proc(w,k,z,t) local i: z^nops(w)*mul(t[op(i..i+k-1,w)],i=1..nops(w)-k+1): end: #CheckGWtE(m,k,N): checks the first N terms of #the output of CheckGWtE(m,k,z,t); For example, try: #CheckGWtE(2,2,5); CheckGWtE:=proc(m,k,N) local z,t,gu,mu,lu,w,i: gu:=GWtE(m,k,z,t): gu:=taylor(gu,z=0,N+1): for i from 0 to N do lu:=expand(coeff(gu,z,i))*z^i: mu:=Milim(m,i): if expand(add(Wt(w,k,z,t),w in mu)-lu)<>0 then RETURN(false): fi: od: true: end: #WtE(m,k,z,t,S): The weight-enumerator of words in #the alphabet {1,..,m} according to the weight z^length #times the product of t[FactorsOfLength_kThatAreInTheSetS]. For #example, try: #WtE(2,2,z,t,{[1,1],[1,2]}); WtE:=proc(m,k,z,t,S) local gu,A,mu,lu,i1,mu1,S1,eq,var,s1: mu:=EqsVars(m,k,z,t,A): eq:=mu[1]: var:=mu[2]: S1:=Milim(m,k) minus S: eq:=subs({seq(t[op(s1)]=1,s1 in S1)},eq): gu:=solve(eq,var): lu:=add(subs(gu,mu1),mu1 in mu[2]): lu:=lu+add(m^i1*z^i1,i1=0..k-2): normal(lu): end: #WtEd(m,k,z,t,S): Like WtE(2,2,z,t,{[1,1],[1,2]}); #but doing it directly, for checking purposes only #Do #WtEd(2,2,z,t,{[1,1],[1,2]}); WtEd:=proc(m,k,z,t,S) local gu,S1,s1: S1:=Milim(m,k) minus S: gu:=GWtE(m,k,z,t): gu:=subs({seq(t[op(s1)]=1,s1 in S1)},gu): normal(gu): end: #RSseq(m,w1,w2,N): Given a pos. integer m and and #two words of the same length w1,w2 in the alphabet #{1,...,m}, finds the first N terms of the sequence #enumerating the number of n-letter words in the #alphabet {1,...,m} that have the same number of #factors that are w1 and that are factors that are w2 #For example, try: #RSseq(2,[1,1],[1,2],10); RSseq:=proc(m,w1,w2,N) local z,t,gu,i1,s,L,i,ku: if not type(m, integer) and m>=1 then print(`Bad input`): RETURN(FAIL): fi: if not type(N, integer) and N>=1 then print(`Bad input`): RETURN(FAIL): fi: if not (type(w1,list) and type(w2,list)) then print(`Bad input`): RETURN(FAIL): fi: if nops(w1)<>nops(w2) then print(`Bad input`): RETURN(FAIL): fi: if {op(w1)} minus {seq(i1,i1=1..m)}<>{} then print(`Bad input`): RETURN(FAIL): fi: if {op(w2)} minus {seq(i1,i1=1..m)}<>{} then print(`Bad input`): RETURN(FAIL): fi: gu:=WtE(m,nops(w1),z,t,{w1,w2}): gu:=subs({t[op(w1)]=s,t[op(w2)]=1/s},gu): gu:=taylor(gu,z=0,N+1): L:=[]: for i from 1 to N do ku:=expand(coeff(gu,z,i)): L:=[op(L),coeff(ku,s,0)]: od: L: end: #RSct(m,w1,w2,z,s): The rational function in z and s #such that the constant term of s is the generating #function for words in {1,..,m}, according to length that #have the same number of occurences of the factor w1 #and w2. #For example, try: #RSct(2,[1,1],[1,2],z,s); RSct:=proc(m,w1,w2,z,s) local gu,t,i1: if not type(m, integer) and m>=1 then print(`Bad input`): RETURN(FAIL): fi: if not (type(w1,list) and type(w2,list)) then print(`Bad input`): RETURN(FAIL): fi: if nops(w1)<>nops(w2) then print(`Bad input`): RETURN(FAIL): fi: if {op(w1)} minus {seq(i1,i1=1..m)}<>{} then print(`Bad input`): RETURN(FAIL): fi: if {op(w2)} minus {seq(i1,i1=1..m)}<>{} then print(`Bad input`): RETURN(FAIL): fi: if not (type(w1,list) and type(w2,list) and {op(w1)} minus {seq(i1,i1=1..m)}={} and {op(w2)} minus {seq(i1,i1=1..m)}={} and w1<>w2) then print(`bad input`): RETURN(FAIL): fi: gu:=WtE(m,nops(w1),z,t,{w1,w2}): normal(subs({t[op(w1)]=s,t[op(w2)]=1/s},gu)): end: #DOaz(m,w1,w2,z,Dz): Given a positive integer #m and two words (lists),w1,w2, of the same length in the #alphabet {1, ...,m}, #uses the continuous Almkvist-Zeilberger algorithm, #to find a linear differential operator with #polynomial coefficients, in terms of z and the #differentiation w.r.t. z, Dz, that annihilates #the (ordinary) generating function #of the sequence #"number of words in {1, ...m} of length n that has #as many factors (substrings) that are w1 and #factors that are w2 equals #such that the constant term of s is the generating #function for words in {1,..,m}, according to length that #have the same number of occurences of the factor w1 #and w2. #For example, to completely solve Richard Stanley's #Amer. Math. Monthly problem #11610 (page 937 # of v. 118 (No. 10), [Dec. 2011]), please type: #DOaz(2,[1,1],[1,2],t,Dt); DOaz:=proc(m,w1,w2,z,Dz) local mu,gu,s: option remember: if not (type(w1,list) and type(w2,list) and {op(w1)} minus {seq(i1,i1=1..m)}={} and {op(w2)} minus {seq(i1,i1=1..m)}={} and w1<>w2) then print(`bad input`): RETURN(FAIL): fi: mu:=RSct(m,w1,w2,z,s)/s: gu:=AZc(mu,s,z,Dz): if gu<>0 then RETURN(gu[1]): else RETURN(FAIL): fi: end: #ProveGF(m,w1,w2,z,F): Given a positive integer #m and two words (lists),w1,w2, of the same length in the #alphabet {1, ...,m}, and an expression F that #depends on z, #proves (rigorously!), using #the continuous Almkvist-Zeilberger algorithm, #that F is the ordinary generating function (using #the variable z) of the sequence #"number of words in {1, ...m} of length n that has #as many factors (substrings) that are w1 and #factors that are w2 equals #For example, to completely solve Richard Stanley's #Amer. Math. Monthly problem #11610 (page 937 # of v. 118 (No. 10), [Dec. 2011]), please type: #ProveGF(2,[1,1],[1,2],t,1/(1-t)/2+(1+2*t)/sqrt((1-t)*(1-2*t)*(1+t+2*t^2))/2); ProveGF:=proc(m,w1,w2,z,F) local mu,gu,s,richard,Dz,i,mu1,F1,Deg1,i1, ku1,ku2: mu:=RSct(m,w1,w2,z,s): gu:=AZc(mu/s,s,z,Dz): if gu=0 then RETURN(FAIL): fi: gu:=gu[1]: Deg1:=degree(gu,Dz)+degree(gu,z): F1:=taylor(F,z=0,Deg1+1): mu1:=taylor(mu,z=0,Deg1+1): ku1:=[seq(coeff(expand(coeff(mu1,z,i1)),s,0),i1=0..Deg1)]: ku2:=[seq(coeff(F1,z,i1),i1=0..Deg1)]: if ku1<>ku2 then print(`Boundary conditions do not match`): print(ku1,ku2): RETURN(false): fi: richard:=add(coeff(gu,Dz,i)*diff(F,z$i),i=1..degree(gu,Dz))+ coeff(gu,Dz,0)*F: richard:=normal(simplify(richard)): if richard<>0 then print(`Diff. eq. not satisfied`): RETURN(false): fi: true: end: #RECaz(m,w1,w2,n,N): Given a positive integer #m and two words (lists),w1,w2, of the same length in the #alphabet {1, ...,m}, #uses the continuous Almkvist-Zeilberger algorithm, #to find a linear recurrence operator with #polynomial coefficients, ope(n,N) in terms of n and the #shift operator in n, N, that annihilates #the sequence #"number of words in {1, ...m} of length n that has #as many factors (substrings) that are w1 and #factors that are w2 equals #such that the constant term of s is the generating #function for words in {1,..,m}, according to length that #have the same number of occurences of the factor w1 #and w2. It also gives the initial values. #For example, to get a recurrence operator for #Richard Stanley's #Amer. Math. Monthly problem #11610 (page 937 # of v. 118 (No. 10), [Dec. 2011]), please type: #RECaz(2,[1,1],[1,2],n,N); RECaz:=proc(m,w1,w2,n,N) local ope,z,Dz,ORDER,lu,i1: option remember: ope:=DOaz(m,w1,w2,z,Dz): ope:=DiffToRec1(ope,z,Dz,n,N): ORDER:=degree(ope,N): ope:=add(factor(coeff(ope,N,i1))*N^i1,i1=0..ORDER): lu:=RSseq(m,w1,w2,ORDER): [ope,lu]: end: #RECg(m,w1,w2,n,N,MaxC): Given a positive integer #m and two words (lists),w1,w2, of the same length in the #alphabet {1, ...,m}, #uses guessing #to find a linear recurrence operator with #polynomial coefficients, ope(n,N) in terms of n and the #shift operator in n, N, that annihilates #the sequence #"number of words in {1, ...m} of length n that has #as many factors (substrings) that are w1 and #factors that are w2 equals #such that the constant term of s is the generating #function for words in {1,..,m}, according to length that #have the same number of occurences of the factor w1 #and w2. It also gives the initial values. #For example, to get a recurrence operator for #Richard Stanley's #Amer. Math. Monthly problem #11610 (page 937 # of v. 118 (No. 10), [Dec. 2011]), please type: #RECg(2,[1,1],[1,2],n,N,10); RECg:=proc(m,w1,w2,n,N,MaxC) local gu,K,ope: option remember: K:=max(trunc(MaxC^2/4+12),7*MaxC+6): gu:=RSseq(m,w1,w2,K): ope:=Findrec(gu,n,N,MaxC): if ope=FAIL then RETURN(FAIL): fi: [ope,[op(1..degree(ope,N),gu)]]: end: #Khalek(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 #Khalek(n*N^2+n/(n+1)*N,N+1/n,n,N); Khalek:=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)1 and type(K,integer) and K>0 and type(x,symbolic) and type(P,symbolic) then print(`bad input`): RETURN(FAIL): fi: if not (type(w1,list) and type(w2,list) and {op(w1)} minus {seq(i1,i1=1..m)}={} and {op(w2)} minus {seq(i1,i1=1..m)}={} and w1<>w2) then print(`bad input`): RETURN(FAIL): fi: gu:=RSseq(m,w1,w2,K): gu:=[1,op(gu)]: Empir(gu,x,P): end: rev11:=proc(w) local i: [seq(w[nops(w)-i+1],i=1..nops(w))]: end: rev1:=proc(S) local w: {seq(rev11(w),w in S)}: end: #Images1(m,S): Given a set of words S in the alphabet #{1, ...,m}, finds the set of images under permuting #the letters and reversing the words #For example, try: #Images1(2,{[1,1],[1,2]}); Images1:=proc(m,S) local mu,gu1,pi,gu11,i,gu2: mu:=permute(m): gu1:={seq(subs({seq(i=pi[i],i=1..m)},S),pi in mu)}: gu2:={seq(rev1(gu11),gu11 in gu1)}: gu1 union gu2: end: #Zugot(m,k): Representatives of sets of pairs of words of #length k in the alphabet {1, ..., m}. For example, try: #Zugot(2,2); Zugot:=proc(m,k) local gu,mu: mu:=Milim(m,k): mu:=choose(mu,2); gu:={}: while mu<>{} do gu:=gu union {mu[1]}: mu:=mu minus Images1(m,mu[1]): od: gu: end: #Info1(m,w1,w2,n,N,P,t,K): Given a pos. integer m, #and two words w1 and w2 of the same length in the #alphabet {1, ..., m}, and symbols n and N #P and t. It also inputs a positive integer #K (the number of terms for guessing the alg. relation, see below) #We are interested in the integer sequence #a(n):=number of n-lettered words in the alphabet {1, ...,m} #that has as many substrings (alias factors, #alias consecutive subwords) that are w1 as there are w2. #The output this procedure is a list consisting #(i) a (conjectured) recurrence operator ope(n,N) annihilating #the sequence a(n) together with the initial values #(starting at n=1) providing #a full holonomic (P-recursive) description of the sequence #(ii) The minimal polynomial F(P,x) such that #F(P(x),x)=0 where P(x) is the ordinary generating function of #the sequence: i.e., #P(x)=Sum(a[n]*x^n,n=0..infinity); #(iii) the "explicit" expression for P(x) in terms of radicals. #For example, for Richard Stanley's Amer. Math. Monthly Problem #11610 (118(10)[Dec. 2011], p. 937) #type #Info1(2,[1,1],[1,2],n,N,P,t,40); Info1:=proc(m,w1,w2,n,N,P,t,K) local MaxC,opeD,REC1,ku: option remember: opeD:=RECaz(m,w1,w2,n,N): MaxC:=degree(opeD[1],N)+degree(numer(opeD[1]),n): REC1:=RECg(m,w1,w2,n,N,MaxC): ku:=ALGg(m,w1,w2,K,t,P): if ku=0 then print(`K =`, K, `is too small, make it bigger`): RETURN(FAIL): fi: [[m,{w1,w2}],REC1,ku]: end: #Info2(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m, #and two words w1 and w2 of the same length in the #alphabet {1, ..., m}, and symbols n and N #P and t. It also inputs a positive integer #K (the number of terms for guessing the alg. relation, see below) #We are interested in the integer sequence #a(n):=number of n-lettered words in the alphabet {1, ...,m} #that has as many substrings (alias factors, #alias consecutive subwords) that are w1 as there are w2. #The output this procedure is a list consisting #(i) a (conjectured) recurrence operator ope(n,N) annihilating #the sequence a(n) together with the initial values #(starting at n=1) providing #a full holonomic (P-recursive) description of the sequence #(ii) The minimal polynomial F(P,x) such that #F(P(x),x)=0 where P(x) is the ordinary generating function of #the sequence: i.e., #P(x)=Sum(a[n]*x^n,n=0..infinity); #(iii) The asymptotics of a(n) to order M #For example, for Richard Stanley's Amer. Math. Monthly Problem #11610 (118(10)[Dec. 2011], p. 937) #type #Info2(2,[1,1],[1,2],n,N,P,t,40,3); Info2:=proc(m,w1,w2,n,N,P,t,K,M) local MaxC,opeD,REC1,ku,lu,C,vu,i1,C1: option remember: opeD:=RECaz(m,w1,w2,n,N): MaxC:=degree(opeD[1],N)+degree(numer(opeD[1]),n): REC1:=RECg(m,w1,w2,n,N,MaxC): ku:=ALGg(m,w1,w2,K,t,P): if ku=0 then print(`K =`, K, `is too small, make it bigger`): RETURN(FAIL): fi: lu:=AsyC(REC1[1],n,N,M,REC1[2],2000); if lu<>FAIL then C:=lu[1]: if lu[3]>8 then C1:=identify(lu[1]): else C1:=0: fi: lu:=C*lu[2]: fi: vu:=SeqFromRec(REC1[1],n,N,evalf(REC1[2]),2000): vu:=evalf([seq(vu[i1]/subs(n=i1,lu),i1=2000-10..2000)]): [[m,{w1,w2}],REC1,ku,lu,vu,C1,C]: end: #TerseBook(m,k,n,N,P,t,K): inputs positive integers m and k #and outputs a set of lists consiting of #the input and output of Info1(m,k,n,N,P,t,K,M) #(qv) for all members of Zugot(m,k) (qv). #For example, try: #TerseBook(2,2,n,N,P,t,40); TerseBook:=proc(m,k,n,N,P,t,K) local gu,i: option remember: gu:=Zugot(m,k): {seq(Info2(m,op(gu[i]),n,N,P,t,K),i=1..nops(gu))}: end: #HSP(ope,n,N): the highest singularity of the recurrence #operator ope(n,N). For example, try: #HSP((n-3)*(n-4)*N+1,n,N); HSP:=proc(ope,n,N) local ope1,pol,gu,i,gu1: ope1:=numer(normal(ope)): pol:=lcoeff(ope1,N): gu:=[solve(pol,n)]: gu1:={}: for i from 1 to nops(gu) do if type(gu[i],integer) and gu[i]>=0 then gu1:=gu1 union {gu[i]}: fi: od: max(op(gu1)): end: #Mishpat(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m, #and two words w1 and w2 of the same length in the #alphabet {1, ..., m}, and symbols n and N #P and t. It also inputs a positive integer #K (the number of terms for guessing the alg. relation, see below) #We are interested in the integer sequence #a(n):=number of n-lettered words in the alphabet {1, ...,m} #that has as many substrings (alias factors, #alias consecutive subwords) that are w1 as there are w2. #The output this procedure is a verbose theorem and sketch of proof #about the sequence a(n) #For example, for Richard Stanley's Amer. Math. Monthly Problem #11610 (118(10)[Dec. 2011], p. 937) #type #Mishpat(2,[1,1],[1,2],n,N,P,t,40,3); Mishpat:=proc(m,w1,w2,n,N,P,t,K,M) local gu,i,mu,ope,i1,a,d1: gu:=Info2(m,w1,w2,n,N,P,t,K,M): print(`Proposition: Let a(n) be the number of n-lettered words`): print(`in the alphabet`, {seq(i,i=1..m)}): print(`With as many occurrences of the substring (consecutive subword)`, w1): print(`as those of`, w2): print(``): print(`Let P(t), be the ordinary generating function of a(n)`): print(P(t)=Sum(a(n)*t^n,n=0..infinity)): mu:=gu[3]: if degree(mu,P)=2 then print(`Then P=P(t) satisfies the quadratic equation`): print(mu=0): else print(`Then P=P(t) satisfies the algebaric equation`): print(mu=0): fi: ope:=gu[2]: print(`Not only that, a(n) safisfies the following linear recurrence`): d1:=degree(ope[1],N): print(a(n+degree(ope[1],N))= -add(coeff(ope[1],N,d1-i1)*a(n+d1-i1),i1=1..d1)): print(``): print(`subject to the initial conditions`): print(seq(a[i1]=ope[2][i1],i1=1..nops(ope[2]))): if gu[4]<>FAIL then print(`Finally, the asymptotics of a(n) to order`,M, `is : `): print(gu[4]): print(`Note that everything is rigorous except the constant in front`,gu[7]): print(`that is a mere non-rigorous estimate`): if gu[6]<>0 and gu[6]<>gu[7] then print(` Maple conjectures that its exact value`): print(`is equal to`, gu[6]): fi: fi: print(``): print(`Sketch of proof: The above has been proved completely rigorously`): print(`internally, and I will spare you the details. If you insist you`): print(` can get them yourself by tracing the code, or else donate `): print(` 50 dollars to the Society for the Liberation of Computerkind.`): print(`QED. `): print(``): print(`For the sake of people who want to use these formulas here`): print(`are the quardaric equation, recurrence operator, and asymptotics`): print(`in Maple input format`): lprint(mu,ope,gu[4]): print(``): print(`For the sake of Sloane, here are the first 50 terms of the sequece`): print([1,op(RSseq(m,w1,w2,50))]): end: #Mishpat1(m,w1,w2,n,N,P,t,K,M,k): #Like Mishpat(m,w1,w2,n,N,P,t,K,M,k) but with the Proposition #numbered k #For example, for Richard Stanley's Amer. Math. Monthly Problem #11610 (118(10)[Dec. 2011], p. 937) #type #Mishpat1(2,[1,1],[1,2],n,N,P,t,40,3,1); Mishpat1:=proc(m,w1,w2,n,N,P,t,K,M,k) local gu,i,mu,ope,i1,a,d1: gu:=Info2(m,w1,w2,n,N,P,t,K,M): print(`Proposition Number `, k, `: Let a(n) be the number of n-lettered words`): print(`in the alphabet`, {seq(i,i=1..m)}): print(`With as many occurrences of the substring (consecutive subword)`, w1): print(`as those of`, w2): print(``): print(`Let P(t), be the ordinary generating function of a(n)`): print(P(t)=Sum(a(n)*t^n,n=0..infinity)): mu:=gu[3]: if degree(mu,P)=2 then print(`Then P=P(t) satisfies the quadratic equation`): print(mu=0): else print(`Then P=P(t) satisfies the algebaric equation`): print(mu=0): fi: ope:=gu[2]: print(`Not only that, a(n) safisfies the following linear recurrence`): d1:=degree(ope[1],N): print(a(n+degree(ope[1],N))= -add(coeff(ope[1],N,d1-i1)*a(n+d1-i1),i1=1..d1)): print(``): print(`subject to the initial conditions`): print(seq(a[i1]=ope[2][i1],i1=1..nops(ope[2]))): if gu[4]<>FAIL then print(`Finally, the asymptotics of a(n) to order`,M, `is : `): print(gu[4]): print(`Note that everything is rigorous except the constant in front`,gu[7]): print(`that is a mere non-rigorous estimate`): if gu[6]<>0 and gu[6]<>gu[7] then print(` Maple conjectures that its exact value`): print(`is equal to`, gu[6]): fi: fi: print(``): print(`Sketch of proof: The above has been proved completely rigorously`): print(`internally, and I will spare you the details. If you insist you`): print(` can get them yourself by tracing the code, or else donate `): print(` 50 dollars to the Society for the Liberation of Computerkind.`): print(`QED. `): print(``): print(`For the sake of people who want to use these formulas here`): print(`are the quardaric equation, recurrence operator, and asymptotics`): print(`in Maple input format`): lprint(mu,ope,gu[4]): print(``): print(`For the sake of Sloane, here are the first 50 terms of the sequece`): print([1,op(RSseq(m,w1,w2,50))]): end: #VerboseBook(m,k,n,P,t,K,M): inputs positive integers m and k #and outputs a web-book with generating functions #and recurrences for sequences enumerating #analogous quantities to Richard Stanley's Monthly #Problem 11610 (118(10)[Dec. 2011], p. 937) #concerning all pairs of words in {1, ...,m} #of length k (up to equivalence). #m is the size of the alphabet, k is the length of the #distinguished words, n,P,t is the symbol for displaying #K and M are as in Info2 #For example try: #VerboseBook(2,2,n,P,t,40,4); VerboseBook:=proc(m,k,n,P,t,K,M) local gu,i1,t0,N: gu:=Zugot(m,k): t0:=time(): print(`Analogs of the Richard Stanley Amer. Math. Monthly Problem 11610`): print(`for ALL pairs of words of length`, k, `in an alphabet of`,m): print(`letters.`): print(``): print(`By Shalosh B. Ekhad `): print(``): for i1 from 1 to nops(gu) do Mishpat1(m,op(gu[i1]),n,N,P,t,K,M,i1): od: print(`It took`, time()-t0, `seconds to generate this book. `): end: #CheckRec(L,ope,n,N): Given a sequence of numbers L and a #recurrence operator ope(n,N) checks whether indeed #L is annihilated by ope(n,N). For example, try: #CheckRec([1,2,4,8,16],N-2,n,N); CheckRec:=proc(L,ope,n,N) local n1,i1: evalb({0}= {seq( add(subs(n=n1,coeff(ope,N,i1))*L[n1+i1],i1=0..degree(ope,N)), n1=1..nops(L)-degree(ope,N))}): end: #CheckAsy(a,n,L,eps): checks whether the asympotic formula a, an #expression in n, is a good one for the sequence L. With tolerance #eps. For example, try: #CheckAsy(2^n,n,[seq(2^i,i=1..100)],1/100); CheckAsy:=proc(a,n,L,eps) local n1,gu: gu:= min( seq(evalf(abs(L[n1]/subs(n=n1,a)-1.)),n1=trunc(nops(L)/2)..nops(L))): if gu