# restart; read `c:/Users/Asus/Google Drive/Aek/FiboConf/Cooper.txt`;

with(LinearAlgebra):


##############
# 4-by-4 Matrix case
##############

# Try: RunC(3,2,-1,5,6);
RunC := proc(a,b,c,d,N) option remember;
local n;

[seq(Cooper(a,b,c,d,n),n=0..N)];
end:

# Output: the coefficient of x^2
# Try: Cooper(3,2,-1,5,3);
# Cooper(a,b,c,d,2);
Cooper := proc(a,b,c,d,n) option remember;
local x,M,A,S;

M := Matrix([[a,b,c,d],[1,0,0,0],[0,1,0,0],[0,0,1,0]]);
A := Matrix([[x,0,0,0],[0,x,0,0],[0,0,x,0],[0,0,0,x]]);

S := expand(Determinant(A-M^n));

-coeff(S,x,2);
end:


##############
# 5-by-5 Matrix case
##############

# Try: RunC2(3,2,-1,5,7,6,0);
RunC2 := proc(a,b,c,d,e,N,deg) option remember;
local n;

[seq(Cooper2(a,b,c,d,e,n,deg),n=0..N)];
end:

# Try: Cooper2(3,2,-1,5,7,3,0);
# Cooper2(a,b,c,d,e,2,2);
Cooper2 := proc(a,b,c,d,e,n,deg) option remember;
local x,M,A,S;

M := Matrix([[a,b,c,d,e],[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0]
,[0,0,0,1,0]]);
A := Matrix([[x,0,0,0,0],[0,x,0,0,0],[0,0,x,0,0],[0,0,0,x,0]
,[0,0,0,0,x]]);

S := expand(Determinant(A-M^n));

coeff(S,x,deg);
end:





###########################################################

# Output: Solve the recurrence of S.
# Try: 4-by-4 case 
# S:=RunC(a,b,c,d,20);
# MySolve(S,6);
MySolve := proc(S,ord) option remember;
local A,eq,var;

eq := {seq(add(A[i]*S[i+j],i=0..ord-1)-S[ord+j]
,j=1..nops(S)-ord-1)};
var := {seq(A[i],i=0..ord-1)};

solve(eq,var);
end:

####################################################
# Check the recurrence with the actual coef. of
# x^2 from the 4-by-4 matrix.
# Input: N the number of terms of coefficeint
# Try: Check42(22);
Check42 := proc(N) option remember;
local j,a,b,c,d,S;
S:=RunC(a,b,c,d,N);


[seq(simplify(   
d^3*S[j]-b*d^2*S[j+1]+d*(d+a*c)*S[j+2]
+(c^2-2*b*d-a^2*d)*S[j+3]-(d+a*c)*S[j+4]
-b*S[j+5]-S[j+6]
),j=1..nops(S)-6)];

end:

# Check the recurrence with the actual coef. of
# x^1 from 5-by-5 matrix.
# Input: N the number of terms of coefficeint
# Try: Check51(22);
Check51 := proc(N) option remember;
local a,b,c,d,e,A,S;

S:=RunC2(a,b,c,d,e,N,1);
A := [0$6];

A[6] := -1;
A[5] := -d;
A[4] := -c*e;
A[3] := -b*e^2;
A[2] := -e^3*a;
A[1] := e^4;


[seq(simplify(add(A[i]*S[i+j],i=1..6)),j=0..nops(S)-6)];

end:

# Check the recurrence with the actual coef. of
# x^2 from 5-by-5 matrix.
# Input: N the number of terms of coefficeint
# Try: Check52(22);
Check52 := proc(N) option remember;
local a,b,c,d,e,A,S;

S:=RunC2(a,b,c,d,e,N,2);
A := [0$11];

A[11] := -1;
A[10] := c;
A[9] := a*e-b*d;
A[8] := e*(b^2+d)-a*(2*c*e-d^2);
A[7] := e^2*(a^2+b)+d^3-d*e*(a*b+3*c);
A[6] := e*(2*e^2+2*a*d*e+e*c*(a^2+2*b)-b*d^2);
A[5] := e^2*(d^2-c*e-a^3*e-3*a*b*e+a*c*d);
A[4] := -e^3*(d*a^2+a*e-c^2+2*b*d);
A[3] := -e^4*(a*c+d);
A[2] := -e^5*b;
A[1] := -e^6;


[seq(simplify(add(A[i]*S[i+j],i=1..11)),j=0..nops(S)-11)];

end:

# Check the recurrence with the actual coef. of
# x^3 from 5-by-5 matrix.
# Input: N the number of terms of coefficeint
# Try: Check53(22);
Check53 := proc(N) option remember;
local a,b,c,d,e,A,S;

S:=RunC2(a,b,c,d,e,N,3);
A := [0$11];

A[11]:=  -1  ;
A[10]:=  -b  ;
A[9]:=  -a*c-d  ;
A[8]:= -a*e+c^2-a^2*d-2*b*d   ;
A[7]:=  -a^3*e-3*a*b*e+a*c*d-c*e+d^2  ;
A[6]:= 2*e^2+e*(2*a*d+a^2*c+2*b*c)-b*d^2   ;
A[5]:= e^2*(a^2+b)+e*(-3*c*d-a*b*d)+d^3   ;
A[4]:= e*(b^2*e+a*d^2+d*e-2*a*c*e)   ;
A[3]:= e^2*(a*e-b*d)   ;
A[2]:= c*e^3   ;
A[1]:= -e^4   ;


[seq(simplify(add(A[i]*S[i+j],i=1..11)),j=0..nops(S)-11)];

end:




# Check the recurrence with the actual coef. of
# x^4 from 5-by-5 matrix.
# Input: N the number of terms of coefficeint
# Try: Check54(22);
Check54 := proc(N) option remember;
local a,b,c,d,e,A,S;

S:=RunC2(a,b,c,d,e,N,4);
A := [0$6];

A[6] := -1;
A[5] := a;
A[4] := b;
A[3] := c;
A[2] := d;
A[1] := e;

[seq(simplify(add(A[i]*S[i+j],i=1..6)),j=0..nops(S)-6)];
end: