#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 13 # Experimental Mathematics # It is okay to link to this assignment on the course webpage. Help := proc(): print(`Help(), GuessH(L, n, N), CheckOper(L, oper, n, N), OTN(S, N), AllOTN(S,N)`): end: # Stuff from Class 13: #Lagrange Inversion applied to tree-enumeration #of ordered #Starting the holonomic ansatz (a.k.a P-recursive) #LIF: Inputs an expression P of the variable #u that possesses #a Maclaurin expansion in u, and outputs the #first N terms of the series expansion of the #unique f.p.s u(t) satisfying the functional eq. #u(t)=t*P(u(t)). For example, #LIF(1+u^2,u,10); should return something with #Catalan numbers LIF:=proc(P,u,N) local n: [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]: end: #GuessH1(L,ORD,DEG,n,N): Inputs a sequence of numbers L #and pos. integer ORD and non-neg. integer DEG and letters #n and N and outputs an linear recurrence operator with #poly. coeffs. (conj.) annihilating the sequence L #nops(L)>=(1+ORD)*(1+DEG)+ORD+6 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 1 # ############# GuessH := proc(L, n, N) local oper, imax, i, ord, deg: imax := floor((nops(L)-7)/2): for i from 1 to imax do: for ord from 1 to i do: for deg from 0 to i-ord do: if nops(L) >= (1+ord)*(1+deg)+ord+6 then: oper := GuessH1(L,ord,deg,n,N): if CheckOper(L, oper, n, N) then: return oper: fi: fi: od: od: od: return FAIL: end: # CheckOper(L, oper, n, N) # Returns true if L is annihliated by oper, where n is the independent # variable and N is the shift operator. CheckOper := proc(L, oper, n, N) local ord,deg, i, j, k, L2: ord := degree(oper, N): deg := degree(oper, n): L2 := add( add( coeff( coeff(oper, n, j), N, i )*[seq(L[k+i]*k^j, k=1..nops(L)-i)][1..nops(L)-ord], j=0..deg ), i=0..ord ); if L2 = [0\$nops(L2)] then: return true: else: return false: fi: end: ############# # Problem 2 # ############# OTN := proc(S, N) local i: return LIF(add(x^i, i=S), x, N): end: ############# # Problem 3 # ############# with(combinat): AllOTN := proc(S,N) local T: return {seq([T, OTN({0} union T, N)], T=powerset(S))}: end: ############# # Problem 4 # ############# # In OEIS: # (empty set) # 1 # 2 # 1,2 # 1,3 # 1,4 # 2,3 # 1,2,3 # 1,3,4 # 1,2,3,4 # 1,2,3,4,5 # In OEIS, but only if zeroes are disregarded: # 3 # 4 # 5 # 2,4 # Not in OEIS: # 1,5 # 2,5 # 3,4 # 3,5 # 4,5 # 1,2,4 # 1,3,5 # 1,4,5 # 2,3,4 # 2,3,5 # 2,4,5 # 3,4,5 # 1,2,3,5 # 1,2,4,5 # 1,3,4,5 # 2,3,4,5