### Kristen Lew, Homework 13 ## ExpMath, Spring 2012 # Post what you wish. Help:=proc(): print(`Previous processes: GuessH1(L,ORD,DEG,n,N), LIF(P,u,N) `): print(`New processes: GuessH(L,n,N), OTN(S,N), AllOTN(S,N) `): end: ### Problem 1 ## Using GuessH1(L,ORD,DEG,n,N) as a subroutine, write a procedure # GuessH(L,n,N) that tries to guess a linear recurrence operator with # polynomial coefficients, ope, phrased in terms of the discrete variable n and # the corresponding shift operator, N (where Nf(n):=f(n+1)), of order+degree as # small as possible, and within the same order+degree, with order as small as possible, # annihilating the sequence L. # If there is no such operator with (1+order)*(1+degree)+6 ≤ nops(L), return FAIL. # If the conjectured operator, ope, has (1+order)*(1+degree)+6 < nops(L) then it # should be checked that ope annihilated the whole list # L (and if it does not, it should keep trying until (1+order)*(1+degree)+6=nops(L) ) GuessH:=proc(L,n,N) local m, ORD, DEG, g : m:=nops(L): for ORD from 1 to m-7 do for DEG from 0 to (m-6)/(1+order)-1 do g:=GuessH1(L,ORD, DEG, n, N): if g<>FAIL then RETURN(g): fi: od: od: FAIL: end: GuessH1:=proc(L,ORD,DEG,n,N) local ope,a,i,j,var,eq,n1,var1: if nops(L)<(1+ORD)*(1+DEG)+ORD+6 then print(`We need at least`, (1+ORD)*(1+DEG)+ORD+6, `terms in L`): RETURN(FAIL): fi: ope:=add(add(a[i,j]*n^i*N^j,i=0..DEG),j=0..ORD): var:={seq(seq(a[i,j],i=0..DEG),j=0..ORD)}: eq:={seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD), n1=1..(1+DEG)*(1+ORD)+6)}: var1:=solve(eq,var): ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:=subs({seq(v=1, v in var)}, ope): ope:=add(factor(coeff(ope,N,j))*N^j, j=0..degree(ope,N)): end: ### Problem 2 ## By using LIF(P,u,N), write a procedure OTN(S,N) that inputs a set of # non-negative integers, S, and outputs the list of non-neg. integers, # let's call it L, such that L[n] is the number of ordered trees with n vertices # where every vertex must have a number of children that is a member of S. # For example, OTN({0,2},7); should give [1,0,1,0,2,0,5] # and OTN({0,1,2},7); # should give the first seven terms of the Motzkin numbers sequence. OTN:=proc(S,N) local s, U, k, u, L: s:=nops(S): U:=0: for k from 1 to s do U:=U+u^(S[k]): od: L:=LIF(U,u,N): end: LIF:=proc(P,u,N) local n: [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]: end: ### Problem 3 ## Write a program AllOTN(S,N) that inputs a set S and # outputs the set {[S1,OTN( {0} union S1,N)]} for all subsets S1 of S. AllOTN:=proc(S,N) local output, PSet, p, k, otn: with(combinat): PSet:=powerset(S): output:={}: p:=nops(PSet): for k from 1 to p do otn:=OTN({0} union PSet[k], N): output:=output union {[PSet[k],otn]}: od: output: end: ### Problem 4 ## By looking at the output of AllOTN({1,2,3,4,5},20); # find out which of the 32 subsets, S1, of {1,2,3,4,5}, have the property # that the number of ordered trees where the number of children of every # vertex must belong to 0 union S1 is currently in the OEIS. # {} : A000007 The characteristic function of 0: a(n) = 0^n. # {1} : A000012 The simplest sequence of positive numbers: the all 1's sequence. # {2} : A126120 A000108 (Catalan numbers) interpolated with 0's. # {1,2} : A001006 Motzkin numbers # {1,3} : A071879 Coefficients of power series solution to g(x) = 1 + x*g(x) + (x*g(x))^3. # {1,4} : A127902 Series reversion of x/(1+x+x^4) # {2,3} : A001005 Number of ways of partitioning n points on a circle into # subsets only of sizes 2 and 3. # {1,2,3} : A036765 Number of rooted trees with a degree constraint. # {1,3,4} : A198951 G.f. satisfies: A(x) = (1 + x*A(x))*(1 + x^3*A(x)^3). # {1,2,3,4} : A036766 Number of rooted trees with a degree constraint. # {1,2,3,4,5} : A036767 Number of rooted trees with a degree constraint.