#Charles Wolf
#hw13

#Feel free to share

#Programs from C13:

 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:

###################################################################################################################################################################
#1)

#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 d,ORD,ope:

for d from 1 to nops(L) do
for ORD from 0 to d do
ope:=GuessH1(L,ORD,d-ORD,n,N):
if (1+ORD)*(1+d-ORD)+6 > nops(L) then RETURN(FAIL):fi:
if ope<>FAIL then RETURN(ope):fi:

od:
od:

FAIL:
end:

####################################################################################################################################################################
#2)

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

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

#####################################################################################################################################################################
#3)

#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 i,P,L:
with(combinat):
P:=powerset(S):
L:={}:
for i from 1 to 2^nops(S) do
L:={op(L),{{op(P[i]),0},OTN({op(P[i]),0},N)}}:
od:
L:
end:

#####################################################################################################################################################################
#4)
#The following 0 union S1 have the desired property, namely that OTN(0 union S1,20) is in OEIS, for subsets S1 of {1,2,3,4,5}:  
#0
#0,5
#0,4
#0,3
#0,2
#0,2,4
#0,1
#0,2,3
#0,1,4
#0,1,3
#0,1,3,4
#0,1,2
#0,1,2,3
#0,1,2,3,4
#0,1,2,3,4,5