### Homework 8.  Patrick Devlin.  February 14, 2012  ###
### I have no preferrences about keeping my work private ###

### Programs from homework 7 and class 6 at the bottom ###

Help:=proc():
  print(`From homework 8:  findPolySpecialProblem(), GuessMPR(L, X, ORDER), RecToPR(C, DEGREE, ORDER, f, n)`):
  print(`From class 8:  GuessPR(L, X, ORDER, DEGREE)`):
  print(`From homework 7:  GenHomPol(X,d,a,co:=1), GenNonHomPol(X,d,a,co:=1)`):
  print(`From class 6:  GuessRec1(L,r),  GuessRec(L) `):
  print(`From homework 6:  RecToSeq(C,n)`):
end proc:



###############################
# Problem 0 (Special Problem) #
###############################

findPolySpecialProblem:=proc() local L, g, tempAns, x, y, A, var, solvedStuff:
 tempAns:=GenNonHomPol([x, y], 4, A):
 g:=tempAns[1]:
 var:=tempAns[2]:
 L:=[[8, 29, 104, 293, 680], [61, 112, 237, 496, 973], [216, 309, 504, 861, 1464], [557, 704, 989, 1472, 2237], [1192, 1405, 1800, 2437, 3400]]: 
 g:=eval(g, solve( {seq(seq( eval(g,[x=i, y=j])=L[i][j], j=1..nops(L[i]))  ,i=1..nops(L))}, var)):
 with(plots, textplot, display):
 display(plots[implicitplot](g,x=-6..6,y=-5..5,axes=none, style = POINT, symbol=point, numpoints=10000),
  textplot([-3.2, 1.5, `proc():`]),
  textplot([-2.02, 1.15, `print(``I love Experimental Math!``)`]),
  textplot([-3.2, .8, `g:=`]), textplot([-1.5,.8,g]),
  textplot([-1.59, .475, `display(plots[implicitplot](g, x=-6..6, y=-5..5, ...`]),
  textplot([-3.1,.17,`end proc:`])):
end proc:




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

GuessMPR:=proc(L, X, ORDER) local DegreeCounter, P:
 for DegreeCounter from 0 to (nops(L)-6) do
  P:=GuessPR(L, X, ORDER, DegreeCounter):
  if(P <> FAIL) then return P: fi:
 od:
 return FAIL:
end proc:


#############
# Problem 2 #
#############

RecToPR:=proc(C, maxDeg, maxOrd, f, n) local degreeOrderSum, tempDeg, tempOrd, L, maxNumberOfCoeffs, varList, answer:
 maxNumberOfCoeffs:=binomial(maxDeg+maxOrd+1, maxOrd+1): # see note below for this count, we double it and add 10 for good measure
 L:=RecToSeq(C,(maxNumberOfCoeffs * 2 + 10)): # we won't need a list bigger than this
 if(L=FAIL) then return FAIL: fi:

 for degreeOrderSum from 0 to (maxDeg+maxOrd) do
  for tempDeg from max(0, degreeOrderSum-maxOrd) to min(degreeOrderSum, maxDeg) do
   tempOrd:=degreeOrderSum-tempDeg:
   varList:=[seq(f(n+i), i=0..tempOrd)]:
   answer:=GuessPR(L, varList, tempOrd, tempDeg, 1):
   if(answer<> FAIL) then return answer: fi:
  od:
 od:

 return FAIL:
end proc:

 ## Note!
 #
 # Let C(D,R) be the number of coefficients in a general polynomial of total degree at most D in at most R variables.
 # Then C(D,R) is the same as the number of terms  x1 ^d1 * x2 ^d2 * ... *
 # xr^dR, where each di is a nonnegative integer, and d1 + d2 + d3 + ... + dR <= D.
 # This is the same as how many ways there are to put AT MOST D
 # indistinguishable objects into R distinguishable containers,
 # which (if you add a new container called "I'm not using this object right now") is the same as
 # the number of ways to be EXACTLY D indistinguishable objects into R+1 distinguishable containers.
 #
 #  By the usual "stars and bars" argument, this gives us  C(D,R) = binomial(D+R, R).
 ###  That is to say, let D stars:  ****** be the objects, and use R bars | |... | to divide up the objects into the R+1 containers,
 ###  then every ordering looks something like   ***|*|**||****| ,  which would be for instance  3 in box 1, 1 in box 2, 2 in box 3, 0 in box 4, 4 in box 5, and 0 in box 6.
 ###  There are binomial(D+R, R) ways of arranging these stars and bars
 # 
 # Therefore!  Since a recurrence of order  ORD needs ORD+1 variables, we know
 # that the largest possible polynomial will have  binomial(maxDeg+maxOrd+1, maxOrd+1)  coefficients



#############
# Problem 3 #
#############

# Running  RecToPR on C:=[[-1,-1,-1], [0,0,1]] with various values of degree and order reveals:
#
# The recurrence with smallest order [then smallest degree] has order 2  and  degree 3.
# The recurrence with smallest degree [then smallest order] has degree 1  and  order 3.
#
# The recurrence with (degree + order) minimized has degree 1 and order 3
#
# The "actual" relation of the tribonacci numbers is most naturally:
# [[-1,-1,-1], [1,1,1]], which is to say,  f(1) = f(2) = f(3) = 1, and  f(n) + f(n+1) + f(n+2) = f(n+3)


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

# Honestly, I don't see a pattern other than the (perhaps "trivially obvious"?) pattern that
# [[-1, -1, -1, ... , -1], [0, 0, 0, ... , 0, 1]] defines a shift of the sequence
# [[-1, -1, -1, ... , -1], [1, 1, 1, ... , 1, 1]]
#
# This seems to be the most natural recurrence for the k-bonacci numbers
## (It's possible there is some simple-looking quasipolynomial that fits well, but I didn't check that)


######################
# Stuff from Class 8 #
######################

# GuessPR(L,X, ORDER, DEGREE)

## Inputs a sequence of numbers L, a formal symbol X (for indexed variables X[0], X[1], ..., X[ORDER]), and nonnegative integers, ORDER and DEGREE
## Outputs a polynomial P of total degree at most DEGREE in the variables x[0], x[1], ..., x[ORDER]
#  such that P(L[n], L[n+1], L[n+2], ... , L[n+order]) = 0 for all relevant n
#  
#  If it fails to do so, then it returns FAIL

GuessPR:=proc(L,X, ORDER, DEGREE, offset:=0) local ans, P, var, x, eqns, iter, k, a:
  x:=[seq(X[iter+offset], iter=0..ORDER)]:
  ans:=[GenNonHomPol(x,DEGREE,a,0)]:
  P:=ans[1]:
  var:=ans[2]:
  # The following line tries to solve a linear system that is as big as possible
  eqns:=solve([seq(0 = eval(P, [seq(X[k+offset]=L[iter+k], k=0..ORDER)] ),iter=1..(nops(L)-ORDER))], var):
  if(eqns <> NULL) then
    try:
      P:=eval(P,eqns):
    catch:
      return FAIL:
    end try:
    if(P=0) then return FAIL: fi:    
    return factor(P):
  fi:
  return FAIL:
end:



#########################
# Stuff from Homework 7 #
#########################

# Takes a list of variables X, the total degree d, a list of coefficients a,
# and an optional parameter co, which keeps track of depth in recursive calls.
# 
# A modified version of GenPol.  It returns a polynomial in the variables
# X[1], ... , X[nops(X)],  where each term is a monomial of total degree d.

GenHomPol:=proc(X,d,a,co:=1) local P, var, co1,m,i,x,X1,P1: option remember:
 m:=nops(X):
 x:=X[m]:

 if m=1 then
  RETURN(x^d * a[co], {a[co]},co+1):
 fi:

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

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

 return (P,var,co1):
end:

# Takes a list of variables X, a maximum total degree d, a variable name for
# the list of coefficients, a,  and an option argument co, indicating which
# index of the list of variables should be the first used for the coefficients.
# 
# Outputs a general polynomial in X, whose maximum total degree is exactly d

GenNonHomPol:= proc(X,d,a, co:=1) local P, var, i, tempAns, co1:
  P:=0:  var:={}:  co1:=co:
  for i from 0 to d do
    tempAns:=GenHomPol(X,i,a,co1):
    P:=expand(P+tempAns[1]):
    var:=var union tempAns[2]:
    co1:=tempAns[3]:
  od:
  return (P, var, co):
end proc:



######################
# Stuff from Class 6 #
######################

# RecToSeq(C,n) takes a C-finite sequence C and returns a list of length n
# generated according to the rules dictated by C
#
# Example: RecToSeq([[-1,-1], [1,1]], 10) --outputs-->
# [1,1,2,3,5,8,13,21,34,55]

RecToSeq:=proc(C,n) local i, i2, L:
 if (nops(C) <> 2) then return FAIL: fi:
 if (nops(C[2]) < nops(C[1])) then return FAIL: fi: # two checks that the input is good

 L:=[seq(C[2][i], i=1..min(n, nops(C[2])))]:

 for i from (1+nops(C[2])) to n do
  L:=[op(L), add(-1*(C[1][j])*(L[i-j])  , j=1..nops(C[2]))]:
 od:

 return L:
end proc:


# GuessRec1(L,r):  given a list of numbers and a positive integer, r,
#  it tries to guess a linear recurrence relation of order r
#  that fits the list

# Looking for a recurrence of the form   L[n] + c_1 L[n-1] + ... + c_r L[n-r] = 0   for all relevant n
# We need r data points to determine this system of (c_1, c_2, ... , c_{r})  unknowns

# So we'll output [ [c_1, c_2, c_3, ... , c_r], [initial_conditions]].
# (i.e., Fibonacci would be  F[n] - F[n-1] - F[n-2] = 0   -> [[-1, -1], [1,1]]

GuessRec1:=proc(L,r) local var, c, i, n, eqns:
 if (nops(L) <= 2*r + 6) then return FAIL; fi: # The list isn't long enough, so return fail
 var:={seq(c[i], i=1..r)}:
 eqns:=[ seq( (L[n] + sum(c[i] * L[n-i],i=1..r))=0, n=(r+1)..nops(L))]:
 var:=solve(eqns,var);
 if(var = NULL) then return FAIL: fi:
 return [[seq(eval(c[i],var),i=1..r)],L[1..r]];
end proc:

# GuessRec(L) returns a linear recurrence relation of minimal order that fits L

GuessRec:=proc(L) local r, ans:
 for r from 1 to (nops(L)-6)/2 do
  ans:=GuessRec1(L,r):
  if (ans <> FAIL) then return ans: fi:
 od:
 return FAIL:
end proc: