#John Kim
#Use whatever you like
#read "C:/Users/John Y. Kim/Documents/Maple/hw8.txt":


Help:=proc(): 
print(`GenPol(X,d,a,co), GenHOMPol(X,d,a,co)`):
print(` GP(X,d,a) , GuessPR(L,x,ORDER,DEGREE)`):
end:

#GenPol(X,d,a,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
#(iiii) co, starting counter
#Outputs: 
#(i)A generic polynomial in the variables X of 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); should output
#a[0]+a[1]*x,{a[0],a[1]},2    
GenPol:=proc(X,d,a,co) local P, var, co1,m,i,x,X1,P1:
option remember:
m:=nops(X):


x:=X[1]:

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


co1:=co:

X1:=[op(2..m,X)]:

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

P,var,co1:
end:


#GenHOMPol(X,d,a,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
#(iiii) co, starting counter
#Outputs: 
#(i)A generic HOMOG. polynomial in the variables X of 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); should output
#a[0]+a[1]*x,{a[0],a[1]},2    
GenHOMPol:=proc(X,d,a,co) local P, var, co1,m,i,x,X1,P1:
option remember:
m:=nops(X):


x:=X[1]:

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


co1:=co:

X1:=[op(2..m,X)]:

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

P,var,co1:
end:


#GP(X,d,a): 
#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

#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); should output
#a[0]+a[1]*x,{a[0],a[1]},2    
GP:=proc(X,d,a) local P, var,co1,d1,P1:

co1:=0:
P:=0:

var:={}:

for d1 from 0 to d do

P1:=GenHOMPol(X,d1,a,co1):

var:=var union P1[2] :

co1:=P1[3]:

P:=P+P1[1]:

od:

P,var:

end:

#GenHomPol(X,d,a,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
#(iiii) co, starting counter
#Outputs: 
#(i)A generic HOMOG. polynomial in the variables X of 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); should output
#a[0]+a[1]*x,{a[0],a[1]},2    
GenHOMPol:=proc(X,d,a,co) local P, var, co1,m,i,x,X1,P1:
option remember:
m:=nops(X):


x:=X[1]:

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


co1:=co:

X1:=[op(2..m,X)]:

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

P,var,co1:
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,var1,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):

var1:={}:

for i from 1 to nops(var) do
var1:=var1 union {lhs(var[i])=subs({seq(a[i]=1,i=0..nops(var))},rhs(var[i]))}:
od:

factor(subs(var1,P)):

end:


#GuessMPR(L,x,ORDER)
#Inputs
#(i)a sequene L of numbers
#(ii) a SYMBOL x (for the indexed variables x[1],x[2], ...)
#(iii) a pos. integer, ORDER
#Outputs
#A polynomial of smallest possible degree, 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.
#If there is no such polynomial, return FAIL.

GuessMPR:=proc(L,x,ORDER) local DEGREE,X,a,P,var,poly:

DEGREE:=1:
X:=[seq(x[i],i=0..ORDER)]:
P:=GP(X,DEGREE,a):
var:=P[2]:

while (nops(var)+6<=nops(L)) do
 poly:=GuessPR(L,x,ORDER,DEGREE):

 if poly<>FAIL and poly<>0 then
  RETURN(poly):
 fi:

 DEGREE:=DEGREE+1:
 P:=GP(X,DEGREE,a):
 var:=P[2]:
#var:=var+1:

od:

RETURN(FAIL):

end:


#RecToSeq(C,n): inputs a C-finite sequence, in our notation, and a pos. integer n, 
#and outputs the list of length n with the first n terms of the sequence. 
#For example, RecToSeq([[-1,-1],[1,1]],10); should output [1,1,2,3,5,8,13,21,34,55].

RecToSeq:=proc(C,n) local L,r,i,j:

L:=C[2]:
r:=nops(C[1]):

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

L:

end: 


#RecToPR(C,DEG,ORD,f,n): tries to find a (minimal) polynomial relation of degree DEG and order ORD,
#satisfied by the members of the C-finite sequence given by C. For example
#RecToPR([[-1,-1],[1,1]],4,1,f,n);
#should give
#(f(n)^2+f(n)*f(n+1)-f(n+1)^2-1) (f(n)^2+f(n)*f(n+1)-f(n+1)^2+1)=0.

RecToPR:=proc(C,DEG,ORD,f,n) local N,L,P,x:

N:=binomial(DEG+ORD+1,ORD+1):
L:=RecToSeq(C,N+20):
P:=GuessMPR(L,x,ORD):

if P=FAIL then
 RETURN(P):
fi:

subs({seq(x[i]=f(n+i),i=0..ORD)},P):

end:


#The Tribonacci numbers satisfy the following degree 3 and order 2 relation:
#-1+f(n+2)^3-2*f(n+2)^2*f(n+1)+2*f(n+1)^3-f(n)*f(n+2)^2-2*f(n)*f(n+2)*f(n+1)+2*f(n)*f(n+1)^2+f(n)^2*f(n+2)+2*f(n)^2*f(n+1)+f(n)^3

#The 4-bonacci numbers satisfy the following degree 4 and order 3 relation:
#(-1+4*f(n)^2*f(n+1)^2+f(n)*f(n+2)^3+3*f(n)^3*f(n+1)-f(n+3)^2*f(n)^2+f(n+3)*f(n)^3+f(n+3)^3*f(n)+3*f(n+3)^3*f(n+2)+5*f(n+1)*f(n+2)^3
#-2*f(n+3)*f(n+2)^3+7*f(n+1)^2*f(n+2)^2-2*f(n+3)^2*f(n+2)^2+3*f(n+1)^3*f(n+2)+2*f(n+1)*f(n+3)^3-3*f(n+1)^3*f(n+3)+2*f(n)*f(n+1)^3
#+2*f(n)^3*f(n+2)+7*f(n)^2*f(n+1)*f(n+2)-2*f(n+1)^2*f(n+3)*f(n+2)-f(n+1)*f(n+3)^2*f(n+2)-8*f(n+1)*f(n+3)*f(n+2)^2
#+8*f(n)*f(n+1)*f(n+2)^2-4*f(n+3)^2*f(n)*f(n+1)-5*f(n+3)*f(n)*f(n+2)^2-2*f(n+3)*f(n)^2*f(n+1)-2*f(n+3)*f(n)^2*f(n+2)
#-3*f(n+3)*f(n)*f(n+1)^2+12*f(n)*f(n+1)^2*f(n+2)-3*f(n+3)*f(n)*f(n+1)*f(n+2)-f(n+1)^4-f(n+3)^4+f(n)^4+3*f(n+2)^4
#+2*f(n+3)^2*f(n)*f(n+2))*(1+4*f(n)^2*f(n+1)^2+f(n)*f(n+2)^3+3*f(n)^3*f(n+1)-f(n+3)^2*f(n)^2+f(n+3)*f(n)^3+f(n+3)^3*f(n)
#+3*f(n+3)^3*f(n+2)+5*f(n+1)*f(n+2)^3-2*f(n+3)*f(n+2)^3+7*f(n+1)^2*f(n+2)^2-2*f(n+3)^2*f(n+2)^2+3*f(n+1)^3*f(n+2)+2*f(n+1)*f(n+3)^3
#-3*f(n+1)^3*f(n+3)+2*f(n)*f(n+1)^3+2*f(n)^3*f(n+2)+7*f(n)^2*f(n+1)*f(n+2)-2*f(n+1)^2*f(n+3)*f(n+2)-f(n+1)*f(n+3)^2*f(n+2)
#-8*f(n+1)*f(n+3)*f(n+2)^2+8*f(n)*f(n+1)*f(n+2)^2-4*f(n+3)^2*f(n)*f(n+1)-5*f(n+3)*f(n)*f(n+2)^2-2*f(n+3)*f(n)^2*f(n+1)
#-2*f(n+3)*f(n)^2*f(n+2)-3*f(n+3)*f(n)*f(n+1)^2+12*f(n)*f(n+1)^2*f(n+2)-3*f(n+3)*f(n)*f(n+1)*f(n+2)-f(n+1)^4-f(n+3)^4
#+f(n)^4+3*f(n+2)^4+2*f(n+3)^2*f(n)*f(n+2))

#The 5-bonacci numbers satisfy a degree 5 and order 4 relation.
#The 6-bonacci numbers: my program outputs FAIL for some reason.

#CONJECTURE: the k-bonacci numbers satisfy a degree k and order k-1 polynomial relation.