### 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.