#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 8
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

# Stuff copied from class 8:

Help := proc(): print(`Help(), GenPol(X,d,a,co), GenHOMPol(X,d,a,co), GP(X,d,a), GenHOMPol(X,d,a,co), GuessPR(L,x,ORDER,DEGREE), GuessMPR(L, x, ORDER), RecToSeq(Seq, n), RecToPR(C, deg, ord, f, n)`): 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,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:

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

GuessMPR := proc(L, x, ORDER) local d, rv:
    for d from 1 to nops(L)/2 do:
        rv := GuessPR(L, x, ORDER, d):
        if rv <> FAIL and rv <> 0 then:
            return rv:
        fi:
    od:
    FAIL:
end:

# From HW 6:

RecToSeq := proc(Seq, n) local ini, co, out, i, j:
    ini := Seq[2]:
    co := Seq[1]:
    out := ini:
    for i from nops(ini)+1 to n do:
        out := [op(out), add(-co[j]*out[i-j], j=1..nops(co))]:
    od:
    return out:
end:

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

RecToPR := proc(C, deg, ord, f, n) local L, PR, x, i:
    L := RecToSeq(C, (deg+1)^(ord+1)+6):
    PR := GuessPR(L, x, ord, deg):
    return subs({seq(x[i] = f[n+i], i=0..ord)}, PR):
end:

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

# The minimal order is 2. The minimal degree is 3. The relation is:

# f(n+2)^3 - 2f(n+2)^2f(n+1) + 2f(n+1)^3 - f(n)f(n+2)^2  -
# 2f(n+2)f(n+1)f(n) + 2 2f(n+1)^2f(n) + f(n)^2f(n+2) + 2f(n)^2f(n+1) +
# f(n)^3

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

# For the sequence of the Kbonacci numbers, the minimal order appears
# to be k-1 (this seems likely), and the minimal degree appears to be
# k if k is odd, or 2k if k is even. This is because when k is even,
# you get a relation of the form
#         P(f(n), f(n+1), ..., f(n+k-1)) = (-1)^n,
# where P is a polynomial of degree k.

# The first couple beyond 3 (with no attempt made to format them nicely):

# k = 4:
#(-27-3*f[n+3]*f[n]*f[n+1]*f[n+2]-2*f[n+3]*f[n]^2*f[n+2]+f[n]^4-f[n+1]^4+3*f[n+2]^4-f[n+3]^4-2*f[n+3]*f[n]^2*f[n+1]-5*f[n+3]*f[n]*f[n+2]^2-3*f[n+
#3]*f[n]*f[n+1]^2+2*f[n+3]^2*f[n]*f[n+2]-4*f[n+3]^2*f[n]*f[n+1]+8*f[n]*f[n+1]*f[n+2]^2+7*f[n]^2*f[n+1]*f[n+2]+12*f[n]*f[n+1]^2*f[n+2]-f[n+1]*f[n+3]^2*f[n+2]-\
#8*f[n+1]*f[n+3]*f[n+2]^2-2*f[n+1]^2*f[n+3]*f[n+2]+f[n+3]*f[n]^3-f[n+3]^2*f[n]^2+2*f[n]*f[n+1]^3+4*f[n]^2*f[n+1]^2+2*f[n]^3*f[n+2]+3*f[n]^3*f[n+1]+5*f[n+1]*f
#[n+2]^3-2*f[n+3]*f[n+2]^3+3*f[n+1]^3*f[n+2]+7*f[n+1]^2*f[n+2]^2+3*f[n+3]^3*f[n+2]-2*f[n+3]^2*f[n+2]^2-3*f[n+1]^3*f[n+3]+2*f[n+1]*f[n+3]^3+f[n+3]^3*f[n]+f[n]
#*f[n+2]^3)*(27-3*f[n+3]*f[n]*f[n+1]*f[n+2]-2*f[n+3]*f[n]^2*f[n+2]+f[n]^4-f[n+1]^4+3*f[n+2]^4-f[n+3]^4-2*f[n+3]*f[n]^2*f[n+1]-5*f[n+3]*f[n]*f[n+2]^2-3*f[n+3]
#*f[n]*f[n+1]^2+2*f[n+3]^2*f[n]*f[n+2]-4*f[n+3]^2*f[n]*f[n+1]+8*f[n]*f[n+1]*f[n+2]^2+7*f[n]^2*f[n+1]*f[n+2]+12*f[n]*f[n+1]^2*f[n+2]-f[n+1]*f[n+3]^2*f[n+2]-8*
#f[n+1]*f[n+3]*f[n+2]^2-2*f[n+1]^2*f[n+3]*f[n+2]+f[n+3]*f[n]^3-f[n+3]^2*f[n]^2+2*f[n]*f[n+1]^3+4*f[n]^2*f[n+1]^2+2*f[n]^3*f[n+2]+3*f[n]^3*f[n+1]+5*f[n+1]*f[n
#+2]^3-2*f[n+3]*f[n+2]^3+3*f[n+1]^3*f[n+2]+7*f[n+1]^2*f[n+2]^2+3*f[n+3]^3*f[n+2]-2*f[n+3]^2*f[n+2]^2-3*f[n+1]^3*f[n+3]+2*f[n+1]*f[n+3]^3+f[n+3]^3*f[n]+f[n]*f
#[n+2]^3)


# k = 5:
# (-256+14*f[n+1]^3*f[n+3]^2+10*f[n]^2*f[n+2]^3+f[n]^2*f[n+4]^3-2*f[n]^2*f[n+3]^3-f[n]^3*f[n+4]^2+6*f[n]^3*f[n+2]^2+8*f[n]*f[n+2]^4-f[n]*f[n+3]^4-
# f[n]*f[n+4]^4+3*f[n]*f[n+1]^4+7*f[n]^3*f[n+1]^2-2*f[n+1]^3*f[n+2]^2-f[n+1]^4*f[n+4]+8*f[n+1]^4*f[n+3]+f[n+1]^4*f[n+2]-4*f[n+4]^4*f[n+3]+5*f[n+4]^3*f[n+3]^2-
# 5*f[n+4]*f[n+3]^4-3*f[n+2]*f[n+4]^4+9*f[n+2]*f[n+3]^4+2*f[n+2]^2*f[n+4]^3+16*f[n+2]^2*f[n+3]^3+10*f[n+2]^3*f[n+3]^2+4*f[n+2]^4*f[n+3]-2*f[n+1]*f[n+4]^4+2*f[
# n+1]*f[n+3]^4+4*f[n+1]*f[n+2]^4+4*f[n]^4*f[n+1]+3*f[n]^4*f[n+2]+2*f[n]^4*f[n+3]+4*f[n+1]^2*f[n+4]^3+8*f[n+1]^2*f[n+3]^3+f[n]^4*f[n+4]+6*f[n]^2*f[n+1]^3+4*f[
# n+2]^5+2*f[n+1]^5+f[n]^5+4*f[n+3]^5+f[n+4]^5-26*f[n+1]^2*f[n+2]*f[n+4]*f[n+3]-32*f[n+1]*f[n+2]^2*f[n+4]*f[n+3]-22*f[n+1]*f[n+2]*f[n+4]*f[n+3]^2+32*f[n]*f[n+
# 1]^2*f[n+2]*f[n+3]+2*f[n]*f[n+1]*f[n+4]^2*f[n+3]-24*f[n]*f[n+2]^2*f[n+4]*f[n+3]-14*f[n]*f[n+2]*f[n+4]*f[n+3]^2-22*f[n]*f[n+1]^2*f[n+2]*f[n+4]-10*f[n]*f[n+1]
# ^2*f[n+4]*f[n+3]-12*f[n]*f[n+1]*f[n+2]^2*f[n+4]+2*f[n+1]*f[n+2]*f[n+4]^2*f[n+3]-6*f[n]^2*f[n+1]*f[n+4]*f[n+3]+26*f[n]^2*f[n+1]*f[n+2]*f[n+3]-10*f[n]^2*f[n+1
# ]*f[n+2]*f[n+4]+32*f[n]*f[n+1]*f[n+2]^2*f[n+3]+32*f[n]*f[n+1]*f[n+2]*f[n+3]^2+2*f[n]*f[n+1]*f[n+2]*f[n+4]^2-16*f[n]^2*f[n+2]*f[n+4]*f[n+3]-10*f[n]*f[n+1]*f[
# n+4]*f[n+3]^2+2*f[n]^3*f[n+2]*f[n+3]-14*f[n+1]*f[n+4]*f[n+3]^3+10*f[n+1]*f[n+4]^2*f[n+3]^2-16*f[n+2]^3*f[n+4]*f[n+3]-2*f[n]*f[n+1]^2*f[n+4]^2+2*f[n+1]*f[n+4
# ]^3*f[n+3]+6*f[n+2]*f[n+4]^3*f[n+3]+6*f[n+2]^2*f[n+4]^2*f[n+3]-20*f[n+2]*f[n+4]*f[n+3]^3+5*f[n+2]*f[n+4]^2*f[n+3]^2-16*f[n+2]^2*f[n+4]*f[n+3]^2+2*f[n+1]*f[n
# +2]*f[n+4]^3+10*f[n+1]*f[n+2]^2*f[n+4]^2+22*f[n+1]*f[n+2]*f[n+3]^3-4*f[n+1]*f[n+2]^3*f[n+4]-2*f[n]^3*f[n+1]*f[n+4]+26*f[n+1]*f[n+2]^2*f[n+3]^2-4*f[n+1]^2*f[
# n+4]^2*f[n+3]-11*f[n+1]^2*f[n+4]*f[n+3]^2+6*f[n+1]^2*f[n+2]*f[n+4]^2+35*f[n+1]^2*f[n+2]*f[n+3]^2-6*f[n+1]^2*f[n+2]^2*f[n+4]-10*f[n+1]^3*f[n+4]*f[n+3]+10*f[n
# +1]^2*f[n+2]^2*f[n+3]-10*f[n+1]^3*f[n+2]*f[n+4]+14*f[n+1]^3*f[n+2]*f[n+3]+21*f[n]*f[n+1]^2*f[n+3]^2-8*f[n]*f[n+1]^3*f[n+4]+10*f[n]*f[n+1]^2*f[n+2]^2+4*f[n]*
# f[n+1]*f[n+4]^3+2*f[n]*f[n+1]^3*f[n+2]+18*f[n]*f[n+1]^3*f[n+3]-2*f[n]*f[n+1]*f[n+3]^3+10*f[n]*f[n+2]*f[n+3]^3+6*f[n]*f[n+2]*f[n+4]^3-6*f[n]*f[n+4]*f[n+3]^3+
# 16*f[n]*f[n+1]*f[n+2]^3-2*f[n]*f[n+2]^2*f[n+4]^2+16*f[n]*f[n+2]^3*f[n+3]-8*f[n]*f[n+2]^3*f[n+4]+26*f[n]*f[n+2]^2*f[n+3]^2+2*f[n]^2*f[n+4]^2*f[n+3]+9*f[n]*f[
# n+4]^2*f[n+3]^2-2*f[n]*f[n+4]^3*f[n+3]+f[n]^2*f[n+2]*f[n+4]^2+4*f[n]^2*f[n+2]*f[n+3]^2-4*f[n]^2*f[n+2]^2*f[n+4]+12*f[n]^2*f[n+2]^2*f[n+3]-6*f[n]^2*f[n+1]*f[
# n+4]^2+10*f[n]^2*f[n+1]*f[n+3]^2+22*f[n]^2*f[n+1]*f[n+2]^2-7*f[n]^2*f[n+1]^2*f[n+4]+22*f[n]^2*f[n+1]^2*f[n+3]+15*f[n]^2*f[n+1]^2*f[n+2]-2*f[n]^3*f[n+4]*f[n+
# 3]-2*f[n]^3*f[n+2]*f[n+4]+12*f[n]^3*f[n+1]*f[n+2]+10*f[n]^3*f[n+1]*f[n+3]-36*f[n]*f[n+1]*f[n+2]*f[n+4]*f[n+3])