### Homework 13 - Patrick Devlin - March 4 2012 ###
## I have no preference about keeping my work private ##

Help:=proc():
 print(`GuessH(L, n, N), OTN(S,N), AllOTN(S,N)`):
 print(`From class 13:`):
 print(`LIF(p,u,N)`):
 print(`GuessHGivenDegAndOrd(L,DEG,ORD,n,N,printErrors:=false)`):
 print(`-- Type  HelpAll()  to see all functions --`):
end:

HelpAll:=proc():
 print(`From class 12:`):
 print(`AlgToSeq(F, P, x, K)\n`):
 print(`From class 10:`):
 print(`GuessAlg(L, P, x, MaxDegP:=(nops(L)/2 - 10)),  GuessAlgHelp(L, P, x, degP, maxDegX:=infinity)`):
 print(`GuessAlgExactDegrees(L, P, x, degP, degX)\n`):
 print(`From class 9:`):
 print(`GenHomPolG(X, d, a, DegList:=[infinity], co:=0),  GPg(X, d, a, DegList:=[infinity], co:=0)`):
 print(`GuessPR1(L, x, ORDER, DEGREE)`):
 print(`From previous classes:`):
 print(`GenPol(X, d, a, co), GenHOMPol(X, d, a, co)`):
 print(`GP(X, d, a), GuessPR(L, x, ORDER, DEGREE)`):
end:



#############
# Problem 1 #
#############

GuessH:=proc(L, n, N) local orderPlusDeg, order, temp:
  for orderPlusDeg from 0 to nops(L) do
   for order from 0 to orderPlusDeg do
     temp:=GuessHGivenDegAndOrd(L, order, orderPlusDeg-order, n, N):
     if(temp<>FAIL) then return temp: fi:
   od:
  od:
  return FAIL:
end proc:


#############
# Problem 2 #
#############
OTN:=proc(S,N) local u:
  return LIF(add(u^s, s in S), u, N):
end proc:

#############
# Problem 3 #
#############
with(combinat):
AllOTN:=proc(S,N) local K, L, S1:
  K:=subsets(S):
  L:={}:
  while not K[finished] do
    S1:=K[nextvalue]():
    L:=L union {[S1, OTN( {0} union S1, N)]}:
  od:
  return L:
end proc:



#############
# Problem 4 #
#############

# The subsets not in Sloane are:
#
# {1, 5}, {2, 4}, {2, 5}, {3, 4}, {3, 5}, {4, 5}, {1, 2, 4}, {1, 2, 5},
# {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}.
#
# Of these, {1, 5} is the MOST interesting, and even still, I personally don't
# think it is "interesting enough" to be in OEIS.



####################
# From other stuff #
####################

LIF:=proc(p, u, N) local n, s:
  s:=taylor(p,u,N+2):
  return [seq(coeftayl(s^n, u=0, n-1)/n,n=1..N)]:   # not even checking if p is in fact a function in u holomorphic in the neighborhood u=0.
end proc:

## This takes a list, L, and tries to output a linear recurrence operator with polynomial coefficients
#    of order exactly ORD and maximum degree exactly DEG that annihilates L.
#  It's in terms of the shift operator N and polynomials n.
# 
## It will output something like  (a_0 (n) + a_1 (n) N + a_2 (n) N^2 + ... + a_ORD (n) N^ORD) L = 0
#  or FAIL if it fails.
#  If the polynomials  a_0 (n) = a_1 (n) = ... = a_ORD (n) = 0  are the only solution, then it will return FAIL


GuessHGivenDegAndOrd:=proc(L, ORD, DEG, n, N, printErrors:=false) local ope, a, i, j, vars, eqns, n1, solvedVars:
  if(nops(L) < (1+ORD)*(1+DEG)+ORD+6) then if(printErrors) then print(`List is too short!  Need at least `, (1+ORD)*(1+DEG)+ORD+6, ` terms in L`): fi: return FAIL: fi:
  
  ope:=add(add(a[i,j]*n^i * N^j, i=0..DEG),j=0..ORD):
  vars:={seq(seq(a[i,j], i=0..DEG),j=0..ORD)}:

  eqns:=[seq(
          add(eval(coeff(ope,N,i)*L[n1+i],{n=n1,N=1}), i=0..ORD) = 0,
         n1=1..min((1+DEG)*(1+ORD)+6, (nops(L)-ORD)))]: #this "min" is actually just for good measure.  We already checked this at first.
       # We don't solve for ALL the equations cause that would be solving perhaps a LOT of stuff if L is large.

  solvedVars:=solve(eqns, vars):
  ope:=subs(solvedVars, ope):

  if(testeq(ope,0)) then return FAIL: fi:

  eqns:=[seq( # after passing first "fast" test, we check ALL the data to make sure it worked
          add(eval(coeff(ope,N,i)*L[n1+i],{n=n1,N=1}), i=0..ORD) = 0,
           n1=1..(nops(L)-ORD))]:
  solvedVars:=solve(eqns, vars):
  ope:=subs(solvedVars, ope):

  if(testeq(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)):
  return ope:
end proc:

############
# Class 12 #
############

### AlgToSeq(F,P,x,K)
## Inputs:
# (i)    an algebraic expression F in variables P and x, where P is a the gen func. of a mystery sequence in the variable x
# (ii)   P, the generating function
# (iii)  an indeterminant x that is the dummy variable in the gen function P
# (iv)   a positive integer K
## Outputs:
#  A sequence of length K+1 denoting the coefficient of x^0 , x^1 , ... , x^K  in the Maclaurin expansion of P
#  - If this fails, it outputs FAIL. -
#
## Examples:
#  AlgToSeq(P*(1-x)-1, P, x, 4) --outputs--> [1,1,1,1,1]
#  AlgToSeq(P-1-x*P^2, P, x, 4) --outputs--> [1,1,2,5,14]

AlgToSeq:=proc(F,P,x,K,start:=1) local H, G, i, F1, f:
  H:=coeff(F,P,1):

  if(testeq(0, eval(H,x=0))) then return FAIL: fi:  # Need to check this or the method will fail

  # Think...!  The expression F says...  P = (P*H-F)/H.  So set this

  F1:=(P*H-F)/H:  # Then F1 = -(a0 (x) + a2(x)*P^2 + ...) / H(x).
                  # So if we iterate this, the P's ought to go further and further away, and we'll be left with F1(x) = STUFF(x) + P^BIG * OTHERSTUFF(P,x)
  f:=start:
  
  for i from 1 to K do
    f:=expand(eval(F1, P=f)):
    f:=series(f,x=0,2*K+10): # truncate it for speed
  od:
  f:=series(f,x=0,K+1):
  return [seq(coeff(f,x,i), i=0..K)]:
end proc:


############
# Class 10 #
############

# GuessAlg(L, P, x, MaxDegP:=(nops(L)/2 - 10)):  Takes an input sequence of numbers L and
# tries to guess and algebraic equation of the form:
#  F(P,x),  where P is the symbol for the generating function of L, [as best as we can tell!],
#  and x is an indeterminant (corresponding to generating function for L).

GuessAlg:=proc(L, P, x, MaxDegP:=(nops(L)/2-10), MaxDegX:=infinity) local g, DegP:

  for DegP from 1 to MaxDegP do
    g:=GuessAlgHelp(L, P, x, DegP, MaxDegX):
    if(g <> FAIL) then return g: fi:
  od:

  return FAIL:
end proc:




GuessAlgHelp:=proc(L, P, x, DegP, MaxDegX:=infinity) local g, DegX:

  for DegX from 0 while( ((1+DegX)*(1+DegP) < nops(L) - 6) and (DegX <= MaxDegX)) do
    g:=GuessAlgExactDegrees(L, P, x, DegP, DegX):
    if(g <> FAIL) then return g: fi:
  od:

  return FAIL:
end proc:


###
## Given a list L, formal symbols, P, and x, and two integers, DegP, DegX,
## it tries to return a polynomial of the form
## F(P, x) [not equal 0], with max degree of P equal to DegP, and max degree of x equal to DegX
## such that if genP(x) is the generating function of L, then F(genP, x) is close enough to 0.

GuessAlgExactDegrees:=proc(L, P, x, DegP, DegX) local F, a, i, j, vars, genP, FofP, eqns, solvedVars, v:
  F:=add(add(a[i,j]*x^i * P^j, i=0..DegX), j=0..DegP):
  vars:={seq(seq(a[i,j],i=0..DegX), j=0..DegP)}:
  genP:=add(L[i] * x^(i-1), i=1..nops(L)):
  
  FofP:=series(eval(F, P=genP), x, nops(L)+10):

  eqns:={ seq( coeff(FofP, x, i) = 0, i=0..nops(L)-1)}:

  solvedVars:=solve(eqns, vars):
  
  F:=eval(F, solvedVars):

  if (testeq(F,0)) then return FAIL: fi:
  # This next line is "optional".  Dr. Z. put it in to get rid of "stupid looking" free parameters
  F:=eval(F, {seq(v = 1, v in vars)}):

  return add(factor(coeff(F,P, i))*P^i, i=0..degree(F,P)): # this is just F cleaned up a little
end proc:



###########
# Class 9 #
###########


# GenHomPolG(X,d,a,DegList, co): 
#
## Inputs:
# (i)    a list of variables X e.g. [x,y,z] (let m=nops(X), the number of variables)
# (ii)   a pos. integer d, for the degree
# (iii)  a symbol a by which the coeff. (undetermined) of the pol. are expressed
# (iv)   a list of degrees for the individual variables (corresponding to X)
# (v)    co, starting counter
#
## Outputs: 
# (i)    A generic HOMOG. polynomial in the variables X of degree d
#         with coeff. expressed as a[co],a[co+1], ...,  a[co+AsNeeded]
#         and with the restriction that the degree of X[i] <=DegList[i] 
# (ii)   the set of coefficients
# (iii)  the new value of the counter
#
# For example, GenHomPolG([x,y],3,a,[2,2],0) --outputs--> a[0]*x^2*y + a[1]*x*y^2, {a[0],a[1]}, 2

GenHomPolG:=proc(X,d,a,DegList:=[infinity],co:=0) local P, var, co1, m, i, x, X1, P1, newDegList: option remember:
 m:=nops(X):
 newDegList:=DegList:
 if(m>nops(DegList)) then newDegList:=[op(DegList), infinity$m]: fi: # or maybe we wanna return FAIL instead?

 if(convert(newDegList,`+`) < d) then return (0, {}, co): fi:

 x:=X[1]:

 if m=1 then
  RETURN(a[co]*x^d, {a[co]},co+1):  # We don't need to wonder if d > DegList[1] by that above check
 fi:

 co1:=co:
 X1:=[op(2..m,X)]:
 newDegList:=[op(2..m, newDegList)]:
 P:=0:
 var:={}:

 for i from 0 to min(d,DegList[1]) do
  P1:=GenHomPolG(X1, d-i, a, newDegList, co1):
  P:=expand(P+P1[1]*x^i):
  var:=var union P1[2]:
  co1:=P1[3]:
 od:

 return (P,var,co1):
end proc:


#GPg(X,d,a, DegList:=[infinity]):
#
## Inputs:
# (i)    a list of variables X e.g. [x,y,z] (let m=nops(X), the number of variables)
# (ii)   a pos. integer d, for the degree
# (iii)  a symbol a by which the coeff. (undetermined) of the pol. are expressed
# (iv)   an optional parameter for DegList, which is a list of maximum degrees for the variables
# (v)    an optional parameter for where the index of the coefficients should start
#
## Outputs: 
# (i)    A generic polynomial in the variables X of total degree d with coeff. expressed as a[co],ac[co+1], ...,  a[co+AsNeeded]
# (ii)   the set of coefficients
# (iii)  the new value of the counter
#
# For example, GenPol([x],1,a,0) --outputs--> a[0]+a[1]*x,{a[0],a[1]},2    

GPg:=proc(X,d,a, DegList:=[infinity], co:=0) local P, var,co1,d1,P1:
 co1:=co:
 P:=0:
 var:={}:

 for d1 from 0 to d do
  P1:=GenHomPolG(X,d1,a,DegList,co1):
  var:=var union P1[2]:
  co1:=P1[3]:
  P:=P+P1[1]:
 od:

 return (P,var):
end proc:


# GuessPR1(L,x,ORDER,DEGREE)
## Inputs 
# (i)    a sequene L of numbers
# (ii)   a SYMBOL x (for the indexed variables x[1],x[2], ...)
# (iii)  a pos. integer, ORDER
# (iv)   a pos. integer Degree
## Output
# A polynomial, let's call it, P(x[0],x[1], ..., x[ORDER]) such
#   that for all n P(L[n],L[n+1], ..., L[n+ORDER])=0
#   and the degree of X[ORDER] is 1

GuessPR1:=proc(L,x,ORDER,DEGREE) local X,i,P,var,eq,a:
  X:=[seq(x[i],i=0..ORDER)]:

  P:=GPg(X,DEGREE,a,[DEGREE$ORDER, 1]): # Note, this is the only line that was changed from GuessPR [below]

  var:=P[2]:

  if nops(var)+6>nops(L) then
   print(`Please make the list longer, we need`, nops(var)+6 , `terms`):
   RETURN(FAIL):
  fi:

  P:=P[1]:

  eq:={ seq( subs( {seq(x[i]=L[n+i],i=0..ORDER)},  P)=0, n=1.. nops(var)+6 ) }:
  var:=solve(eq,var):
  factor(subs(var,P)):

end:




#GuessPR(L,x,ORDER,DEGREE)
#Inputs 
#(i)a sequene L of numbers
#(ii) a SYMBOL x (for the indexed variables x[1],x[2], ...)
#(iii) a pos. integer, ORDER
#(iv) a pos. integer Degree
#Output
#A polynomial, let's call it, P(x[0],x[1], ..., x[ORDER]) such
#that for all n P(L[n],L[n+1], ..., L[n+ORDER])=0
GuessPR:=proc(L,x,ORDER,DEGREE) local X,i,P,var,eq,a:

X:=[seq(x[i],i=0..ORDER)]:

P:=GP(X,DEGREE,a):

var:=P[2]:

if nops(var)+6>nops(L) then
 print(`Please make the list longer, we need`, nops(var)+6 , `terms`):
 RETURN(FAIL):
fi:

P:=P[1]:

#eq:={ seq( subs( {x[0]=L[n],x[1]=L[n+1], ..., x[ORDER]=L[n+ORDER]},
#  P)=0, n=1..min(nops(L)-ORDER, nops(var)+6)  }:


eq:={ seq( subs( {seq(x[i]=L[n+i],i=0..ORDER)},  P)=0, n=1.. nops(var)+6 ) }:

var:=solve(eq,var):

factor(subs(var,P)):

end: