###  Homework 6.  Pat Devlin.  February 11, 2012.  ###
###  I have no preferrence about keeping my work private ###



## Program from class 6 at the end ##

Help:=proc(): print(`For Homework 6: RecToSeq(C,n),  Add1(C1,C2),  Mult1(C1,C2),  BT(C),  Tsect(C,t,r),  `): print(`From class 6: GuessRec1(L,r),  GuessRec(L) `): end:



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

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




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

# Add1(C1,C2) takes two C-finite sequences, C1 and C2, and returns a C-finite sequence
# for their term-wise sum of smallest order possible.
#
# Example: Add1([[-1],[1]],[[-2],[1]]) --outputs-->  [[-3,2], [2,3]]

Add1:=proc(C1,C2) local L1, L2, termsNeeded, i:
 try:
  termsNeeded := 2*nops(C1[1])+2*nops(C2[1])+10:
 catch:
  print(`Invalid input for Add1(C1,C2)`):
  return FAIL:
 end try:
 
 
 L1:=RecToSeq(C1,termsNeeded):
 L2:=RecToSeq(C2,termsNeeded):
 if((L1 = FAIL) or (L2 = FAIL)) then
  print(`Insufficient initial conditions or invalid input form for Add1(C1,C2)`):
  return FAIL:
 fi:
 
 return GuessRec([seq(L1[i]+L2[i],i=1..min(nops(L1),nops(L2)))]):
end proc:



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


# Mult1(C1,C2) takes two C-finite sequences, C1 and C2, and returns a C-finite sequence
# for their term-wise product of smallest order possible.
#
# Example:  Mult1([[-2],[1]],[[-3],[1]]) --outputs-->  [[-6],[1]].

Mult1:=proc(C1,C2) local L1, L2, termsNeeded, i:
 try:
  termsNeeded := 2*nops(C1[1])*nops(C2[1])+10:
 catch:
  print(`Invalid input for Mult1(C1,C2)`):
  return FAIL:
 end try:


 L1:=RecToSeq(C1,termsNeeded):
 L2:=RecToSeq(C2,termsNeeded):
 if((L1 = FAIL) or (L2 = FAIL)) then
  print(`Insufficient initial conditions or invalid input form for Mult1(C1,C2)`):
  return FAIL:
 fi:

 return GuessRec([seq(L1[i]*L2[i],i=1..min(nops(L1),nops(L2)))]):
end proc:


#############
# Problem 5 #
#############


# BT(C) takes a C-finite sequence, C, and returns a C-finite
# sequence for its binomial transformation, which is given by...
# 
# If a[n]  (n=0..infinity)  is the sequence satisfying the rule C, then the sequence
# b[n] := sum( binomial(n,k) * a[k], k=0..n)  is the binomial transform for a.
# 
# We then return b[n] as a C-finite sequence as in our notation
# 
# Example: Add1([[-1],[1]],[[-2],[1]]) --outputs-->  [[-3,2], [2,3]]

BT:=proc(C) local L, termsNeeded, i,k:
 try:
  termsNeeded := 2*nops(C[1])+10:
 catch:
  print(`Invalid input for BT(C)`):
  return FAIL:
 end try:

 L:=RecToSeq(C,termsNeeded):

 if(L = FAIL) then
  print(`Insufficient initial conditions or invalid input form for BT(C)`):
  return FAIL:
 fi:

 return GuessRec([seq(sum( binomial(i-1, k-1) * L[k], k=1..i) ,i=1..nops(L))]):
end proc:




#############
# Problem 6 #
#############

# Tsect(C,t,r) takes a rule for a C-finite sequence, C, and integers
# 0 <= r < t.  Now let a[n] be the sequence satisfying rule C.
# Then Tsect outputs the rule for the sequence a[t*n + r].

Tsect:=proc(C,t,r) local L, termsNeeded, i, t2, r2:

 t2:=round(t):  r2:=round(r):

 try:
  termsNeeded := (2*nops(C[1])+10)*t2 + r2+3:
 catch:
  print(`Invalid input for Tsect(C,t,r)`):
  return FAIL:
 end try:

 L:=RecToSeq(C,termsNeeded):

 if(L = FAIL) then
  print(`Insufficient initial conditions or invalid input form for Tsect(C,t,r)`):
  return FAIL:
 fi:
 
 if ( (r2<0) or (r2 >= t2)) then
  print(`Invalid values of r or t for Tsect(C,t,r)`):
  return FAIL:
 fi:

 return GuessRec([seq(L[t2*i+r2+1],i=0..(nops(L)-r2-1)/t2)]):  # t*i+r+1 > nops(L)  <=>  i > (nops(L)-r-1)/t
end proc:



########################
# Program from class 6 #
########################

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