## Kristen Lew, Homework 8
## ExpMath 2012
## Post if you wish.  

Help:=proc():
print(` New procedures: GuessMPR(L,x,ORDER), RecToPR(C,DEG,ORD,f,n), MinTrib()`):
print(` Previous procedures:  GuessPR(L,x,ORDER,DEGREE), GP(X,d,a),`):
print(` GenHOMPol(X,d,a,co), RecToSeq(C,n), GuessRec(L), GuessRec1(L,r)` ):
end:

### Problem 1
## GuessMPR(L,x,ORDER) finds the (non-zero) polynomial relation of order ORDER 
# of smallest degree, compatible with the amount of data in L, or returns FAIL.  
GuessMPR:=proc(L,x,ORDER) local m,i,G:
	m:=nops(L):
	for i from 1 to m do
		G:=GuessPR(L,x,ORDER,i):
		if G<> 0 and G<>FAIL then
			RETURN(G):
		fi:
	od:
	FAIL:
end:

		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( {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:

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

### Problem 2
## 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,2,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 C1,G,i:
	C1:=RecToSeq(C,10*DEG*ORD):
	G:=GuessMPR(C1,f,ORD,DEG):
	for i from 0 to ORD do
		G:=subs(f[i]=f[n+i],G):
	od:
	if G=FAIL then
		RETURN(FAIL):
	fi:
	G=0:
end:
		RecToSeq:=proc(C,n) local L, L1, L2, i, j, t:
			L1:=C[2]: L2:=C[1]: t:=nops(L1)+1: 
			for i from t to n do
				L1:=[op(L1), add(-L2[j]*L1[i-j], j=1..t-1)]:
			od:
			L1:
		end:


### Problem 3
## What is the minimal order and degree, and what is the actual relation, satisfied
# by the Tribonacci numbers.  

# Using GuessRec, from C6.txt, we find that the C-finite sequence of the Tribonacci 
# numbers is [[-1, -1, -1], [1, 1, 2]].  
# Using the procedure MinBonacci(L) 
# with L=[0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504]
# returns a[4]*(-f[n+3]+f[n+2]+f[n+1]+f[n]) = 0, 1, 3.
# The recursive relation is (-f[n+3]+f[n+2]+f[n+1]+f[n]) = 0
# The minimal order is 3 and minimal degree is 1.  

MinBonacci:=proc(L) local L1, i, j, Guess:
	L1:=GuessRec(L):
	print(L1):
	for i from 1 to 6 do
		for j from 4 to 6 do
			Guess:=RecToPR(L1,i,j,f,n):
			if Guess <> FAIL then
				RETURN(Guess, i, j):
			fi:
		od:
	od:

end:
		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:
			FAIL:
		end:
		GuessRec1:=proc(L,r) local c,i,var,eq:
			if nops(L)<=2*r+6 then
 				RETURN(FAIL):
			fi:
			var:={seq(c[i],i=1..r)}:
			eq:={seq( L[n]+add(c[i]*L[n-i],i=1..r)=0 , n=r+1..nops(L))}:
			var:=solve(eq,var):
			if var=NULL then
  				RETURN(FAIL):
			fi:
			[[seq(subs(var, c[i]),i=1..r)], L[1..r]   ] ;
		end:

### Problem 4
## What is the minimal order and degree and what is the actual relation, 
# satisfied by the K-bonacci numbers (the sequence f(n) defined by
# f(0)=…=f(K-2)=0, f(K-1)=1, f(n)=f(n-1)+..+f(n-k) for n>=K for K=4,5,6.  

# K4:=[0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773, 1490, 2872, 5536]:
# GuessRec(K4);
# [[-1, -1, -1, -1], [0, 0, 0, 1]]
# RecToPR(LK4,1,4,f,n);
# a[5] (-f[n + 4] + f[n + 3] + f[n + 2] + f[n + 1] + f[n]) = 0
### Minimal order = 4, minimal degree = 1
### Recursive relation: (-f[n + 4] + f[n + 3] + f[n + 2] + f[n + 1] + f[n]) = 0

# K5:=[0, 0, 0, 0, 1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525]:
# GuessRec(K5);
# [[-1, -1, -1, -1, -1], [0, 0, 0, 0, 1]]
# RecToPR(LK5,1,5,f,n);
# a[6] (-f[n + 5] + f[n + 4] + f[n + 3] + f[n + 2] + f[n + 1] + f[n]) = 0
### Minimal order = 5, minimal degree = 1
### Recursive relation: 
### (-f[n + 5] + f[n + 4] + f[n + 3] + f[n + 2] + f[n + 1] + f[n]) = 0

# K6:=[0, 0, 0, 0, 0, 1, 1, 2, 4, 8, 16, 32, 63, 125, 248, 492, 976, 1936, 3840, 7617]:
# GuessRec(K6);
# [[-1, -1, -1, -1, -1, -1], [0, 0, 0, 0, 0, 1]]
# RecToPR(LK6,1,6,f,n);
# a[7] (-f[n + 6] + f[n + 5] + f[n + 4] + f[n + 3] + f[n + 2] + f[n + 1] + f[n]) = 0
### Minimal order = 5, minimal degree = 1
### Recursive relation: 
### (-f[n + 6] + f[n + 5] + f[n + 4] + f[n + 3] + f[n + 2] + f[n + 1] + f[n]) = 0