Help16:=proc(): print(`LIFtoS(F,P,n), NRope(F,P,n,N,K)`): print(`Rope(F,P,n,N) `): end: #LIFtoS(F,P,n): #Inputs an expression F in the variable P. # and a pos. integern #outputs the list of the first n coefficients in the formal power #series that is the UNIQUE solution of the FUNCTIONAL equatation #P(x)=x*F(P(x)) #For example, to get the the enumerating sequence for the #number of complete binary trees with i vertices, up to 10 #do LIFtoS(1+P^2,P,x,10); #P(x)=x*(1+P(x)^2) LIFtoS:=proc(F,P,n) local x,Y,i,p: Y:=coeff(F,P,0)*x: for i from 1 to n do Y:=x*subs(P=Y,F): Y:=add(coeff(Y,x,p)*x^p,p=0..n): od: [seq(coeff(Y,x,i),i=1..n)]: end: ############stuff from AZ.txt############## #AZ.txt: The Almkvist-Zeilberger algorithm #taken from #http://www.math.rutgers.edu/~zeilberg/EKHAD HelpAZ:=proc(): print(` AZd(A,y,n,N) `): end: pashetd:=proc(p,N) local gu,p1,mekh,i: p1:=normal(p): mekh:=denom(p1): p1:=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,mekh): end: goremd:=proc(N,ORDER) local i, gu: gu:=0: for i from 0 to ORDER do gu:=gu+(bpc[i])*N^i: od: RETURN(gu): end: duisd:= proc(A,ORDER,y,n,N) local 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,shad,KAMA,i1,apc,apd: KAMA:=10: gorem:=goremd(N,ORDER): conj:=gorem: yakhas:=0: for i1 from 0 to degree(conj,N) do yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A): 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:=pashetd(gorem,N): 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: RETURN(shad[1],S): end: #AZd(A,y,n,N): inputs a hyperexponential expression A of the #contintuous variable y and the discrete variable n and #outputs a linear recurrence operator annihilating the #integral of A(y,n) along an interval where #A vanishes at the endpoints (including infinite intervals) #or a closed circle around the orign. For example try #AZd((1+y+y^2)^n/y^(n+1),y,n,N); AZd:= proc(A,y,n,N) local ORDER,gu,KAMA: KAMA:=200: for ORDER from 0 to KAMA do gu:=duisd(A,ORDER,y,n,N): if gu<>0 then RETURN(gu): fi: od: FAIL: end: #####End of stuff from AZ.txt #NRope(F,P,n,N,K): Inputs an expression F in variable P and #symbols n and N and a pos. integer K for #guessing, outputs a non-rogorously guessed #linear recurrence operator in the shift operator N # #annihilating the sequence of coeffs. of the unique f.p.s. #solution of #P(x)=x*F(P(x)).For example, #NRope(1+P+P^2,P,n,N,30); NRope:=proc(F,P,n,N,K) local L: L:=LIFtoS(F,P,K): GuessH(L,n,N): end: #Rope(F,P,n,N): Inputs an expression F in variable P and #symbols n and N and a pos. integer K for #guessing, outputs a rigorously derived (using the #not-as-famous-as-it-should-be Almkvist-Zeilberger algorithm #AND Lagrange-Inversion, a #linear recurrence operator in the shift operator N # #annihilating the sequence of coeffs. of the unique f.p.s. #solution of #P(x)=x*F(P(x)).For example, #Rope(1+P+P^2,P,n,N,30); Rope:=proc(F,P,n,N): #LIF: the n-th term is the coeff. of P^(n-1) in #F(P)^n divided by n AZd((F/P)^n/n,P,n,N)[1]: end: #Asy(ope,n,N,Ini): inputs a linear recurrence operator ope #in n and the shift operator N, and the initial values #(of length the order of ope, i.e. degree(ope,N) and #outputs (parially non-rig.) the asymptic expression #for the sequence in the form #C*mu^n*n^theta Asy:=proc(ope,n,N,Ini) local mu,ope1,sol,max: ope1:=lcoeff(ope,n): sol:=[solve(ope1,N)]: end: