#John Kim
#Use whatever you like

#read "C:/Users/John Y. Kim/Documents/Maple/hw13.txt":

 Help:=proc(): print(` LIF(P,u,N), GuessH1(L,ORD,DEG,n,N)`): end:
 #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:


#GuessH(L,n,N): Inputs a sequence of numbers L and letters 
#n and N and outputs an linear recurrence operator with 
#poly. coeffs. (conj.) annihilating the sequence L 
#with minimum DEGREE+ORDER (with min ORDER).

GuessH:=proc(L,n,N) local ORD, DEG, ope:

for ORD from 0 to (nops(L)-7)/2 do
for DEG from 0 to (nops(L)-ORD-6)/(1+ORD)-1 do

ope:=GuessH1(L,ORD,DEG,n,N):
if ope<>FAIL then
 RETURN(ope):
fi:

od:
od:

FAIL:

end:



#OTN(S,N): 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 P,i,u:

P:=add(u^S[i],i=1..nops(S)):
LIF(P,u,N):

end:



#sets(S): Inputs a set S and outputs the set of subsets of S.

sets:=proc(S) option remember: local P,P1,i:

if S={} then
 {{}}:
else
 P:=sets({op(1..nops(S)-1,S)}):
 P1:={seq(P[i] union {S[nops(S)]},i=1..nops(P))}:
 P union P1:
fi:
end:



#AllOTN(S,N): 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 P,i,All:

P:=sets(S):
All:={}:

for i from 1 to nops(P) do

All:=All union {[P[i],op(OTN(P[i] union {0},N))]}:

od:

All:

end:


#AllOTN({1,2,3,4,5},20);
#All subsets w/ one element are in OEIS
#Similar to this form: Binomial(5n,n)/(4n+1)
#{1,2} yields Motzkin numbers
#{1,3} yields Coefficients of power series solution to g(x) = 1 + x*g(x) + (x*g(x))^3.
#{1,4} yields Series reversion of x/(1+x+x^4). 
#{2,3} yields Number of ways of partitioning n points on a circle into subsets only of sizes 2 and 3. 
#{2,4} yields 	Number of modes of connections of 2n points. 
#{1,2,3} yields Number of rooted trees with a degree constraint. 
#{1,2,3,4} yields above.
#{1,2,3,4,5} yields above.