# restart; read `c:/Users/Asus/Google Drive
# /Aek/Fibo/FiboMatrix.txt`;

#with(Quaternions):
with(LinearAlgebra,Transpose):
with(LinearAlgebra,Determinant):

# Generalize Fibonacci into the complex plane.
# First Version: June 19, 2014
# Second Version: August 17, 2015
# Third Version: April 10, 2016
# Fourth Version: December 25, 2018

ez:=proc(): 
print(` Section 1: Sequence`): 
print(``):
print(` Section 2: One Matrix Multiply`):
print(``):
print(` Section 3: Two Matrix Multiply`):
print(``):
print(` Section 4: Multiplication`): 
print(` of Three (or more) Matrices`):
print(``):
print(` Section 5: Numbers Related`): 
print(` to 3-by-3 Matrix`):
print(``):
print(` Section 6: Complex Fibonacci`):
print(` numbers by Harman`):
end:



######################
# Section 1: Sequence
######################

# C-finite of order nops(C):
# Input: constants, C,A.
# Try: [seq(CF(n,[1,1],[0,1]),n=-10..10)];
CF := proc(n,C,A) option remember;
local i,tt,temp;

if C[nops(C)]=0 then
	ERROR(NotDefined);
fi:

tt := nops(C):
if n >= 0 and n <= nops(A)-1 then
	return(A[n+1]);
fi:

if n >= nops(A) then
	add( C[i]*CF(n-i,C,A), i=1..tt);
else
	temp := add( -C[i]*CF(n+(tt-i),C,A), i=1..tt-1);
	(CF(n+tt,C,A)+temp)/C[tt];
fi:
end:




################################
# Section 2: One Matrix Multiply
################################

# Fibonacci alone.

# Try: [seq(Fib(n),n=0..10)];
Fib := proc(n) option remember;
CF(n,[1,1],[0,1]);
end:


# Fibonacci Matrix
# Try: [seq(FibMat(n),n=1..10)];
FibMat := proc(n) option remember;
Matrix([[Fib(n+1), Fib(n)]
,[Fib(n), Fib(n-1)]]);
end:


# Lucas Matrix
# Try: [seq(LucasMat(n),n=1..10)];
LucasMat := proc(n) option remember;
Matrix([[Lucas(n+1), Lucas(n)]
,[Lucas(n), Lucas(n-1)]]);
end:


################################
# Section 2.1.1: 
# Identities of Fibonacci numbers  
################################
# Try: [seq(Id3(n),n=0..30)];
Id3 := proc(n) option remember;
local j,k,f;

f := k->CF(k,[1,1],[0,1]);

[add(f(j),j=1..n),f(n+2)-1];
end:



# Intermediate calculation for
# identity 4
# Try: [seq(TestId4(n),n=1..10)];
TestId4 := proc(n) option remember;
local k,L,R;

L := FibMat(n);
R := add(FibMat(n-1-2*k),k=0..(n-1)/2);

if n mod 2 = 0 then
	R := R+1;
else
	R := R+Matrix([[0,1],[1,-1]]);
fi:
[L,simplify(R-L)];
end:


# Try: [seq(Id16(2,3,k),k=1..10)];
Id16 := proc(n,m,k) option remember;
local j,L,R;

L := add(Fib(n+j)*Fib(m+j),j=1..2*k);
R := Fib(n+2*k+1)*Fib(m+2*k)
-Fib(n+1)*Fib(m);
[L,R]:
end:


# Try: seq(Id18(2,3,k),k=1..10);
Id18 := proc(n,m,k) option remember;
local j,t,P,L,R;
t := add(Fib(n+j)*Fib(m+j),j=1..2*k);
P := Matrix([[1,1],[1,0]]):
L := simplify(add(FibMat(n+m+4*j-1),j=1..k));
R := P^(n+m+2).(P^(4*k)-1).(1+P^2)^(-1);

[t,L,simplify(R)];
end:



################################
# Section 2.1.2: Sum of Products
################################
# Theorem 2.1.1, 2.1.2, 2.1.3, 2.1.4
# Sum of k-product of Fib 
# page 370 of Koshy
# [seq(HFib(5,j),j=1..10)];
# [seq(ProdF(m,5),m=1..10)];
ProdF := proc(m,k) option remember;
local i,a;

add( mul(Fib(a[i]),i=1..nops(a)),
a in AddT(m,k));
end:

# add to m with k tuples
# Try: AddT(6,3);
AddT := proc(m,k) option remember;
local i,a,S;

if k =1 then
	return({[m]});
fi:

S := {};
for i from 1 to m-k+1 do
	S := S union 
	{seq([i,op(a)] ,a in AddT(m-i,k-1))};
od:
S;
end:

# Try: [seq(HFib(1,j),j=1..10)];
HFib := proc(i,j) option remember;
if i = 0 and j >=0 then
	0;
elif i=j and j >=1 then
	1;
elif i > j then
	0;
elif i >=1 and j>=2 then
	HFib(i,j-2)+HFib(i,j-1)+HFib(i-1,j-1);
else
	ERROR("NoCase");
fi:
end:

# Try: [seq(HFibMat(n),n=2..10)];
HFibMat := proc(n) option remember;
Matrix([[HFib(2,n+1), HFib(2,n)]
,[HFib(2,n), HFib(2,n-1)]]);
end:

# Try: [seq(TestHFib(n),n=2..10)];
TestHFib := proc(n) option remember;
local k,P,R,T1,T2,T3,T,S;

P := Matrix([[1,1],[1,0]]);
R := Matrix([[1,0],[0,0]]);
T1 := HFibMat(n).P+P^(n).R;
T2 := HFibMat(n-1).P^2
+P^(n-1).R.P+P^(n).R;
T3 := HFibMat(n-2).P^3
+P^(n-2).R.P^2+P^(n-1).R.P+P^(n).R;

T := simplify(add( P^(n-k).R.P^k,k=0..n));

[T1,T2,T3,T,T-HFibMat(n+1)];
end:




################################
# Section 2.2: Guassian Fibonacci
# by Jordan
################################

# Try: seq(Jor(n),n=-5..10);
Jor := proc(n) option remember;

if n=0 then
	I;
elif n=1 then
	1
elif n >0 then;
	Jor(n-1)+Jor(n-2); 
else
	Jor(n+2)-Jor(n+1);
fi:
end:

# Try: seq(JorMat(n),n=-1..5);
JorMat := proc(n) option remember;
Matrix([[Jor(n+1),Jor(n)],[Jor(n),Jor(n-1)]]):
end:

CJorMat := proc(n) option remember;
local P,B,C;
P := Matrix([[1,1],[1,0]]):
B := Matrix([[1,I],[I,1-I]]);
C := Matrix([[1,-I],[-I,1+I]]):
P^n.C;
end:


################################
# Section 2.2.1: Identities of
# Gaussian Fibonacci
################################

# Identity 5
# Try: seq(Jor5(n),n=1..10);
Jor5 := proc(n) option remember;
local L,R;

L := Jor(n+1)^2-Jor(n-1)^2;
R := Fib(2*n-1)*(1+2*I);

[L,R];
end:


# Matrix of identity 5
# Try: seq(JorMat5(n),n=1..5);
JorMat5 := proc(n) option remember;
local P,R,S,A,B,C;

P := Matrix([[1,1],[1,0]]):
R := Matrix([[1,I],[I,1-I]]):
S := Matrix([[1,0],[0,-1]]):

A := P^(n).R.S.P^(n+1).R
+P^(n-1).R.S.P^(n).R;
B := (1+2*I)*P^(2*n-1);
C := JorMat(n)^2;

[A,B,C];
end:


# Identity 7
# Try: [seq(Jor7(n),n=1..16)];
# [seq(Jor7(2*n)/(1+2*I),n=1..6)];
Jor7 := proc(n) option remember;
local j,L,R;

L := add(Jor(j)^2,j=1..n);
R := (1+2*I)*Fib(n)^2+(-1)^(n)*I-I;

[L,R];
end:



# Matrix for identity 7
# Try: seq(JorMat7(2*n),n=1..5);
JorMat7 := proc(n) option remember;
local j,P,R,A,B;
if n mod 2 <> 0 then
	ERROR("n must be even");
fi:
P := Matrix([[1,1],[1,0]]):
R:=Matrix([[1,I],[I,1-I]]);

P.(P^(2*n)-1).(P^2+1)^(-1);

A := Matrix([[Fib(2*n+1)-1,Fib(2*n)]
,[Fib(2*n),Fib(2*n-1)-1]]);
B := Matrix([[1,2],[2,-1]]);

1/5*A.B;
end:



# Matrix form for identity 8
# Try: seq(Jor8(n),n=1..10);
Jor8 := proc(n) option remember;
local j,P,R,A,B;

P := Matrix([[1,1],[1,0]]):
R := Matrix([[1,I],[I,1-I]]):

A := P^(-n).R;
B := (P^(n+1).R)^(-1).P.R^2;

#print(P.R^2);
[A,B];
end:

# Fib(2*n+1)= Fib(n+1)^2+Fib(n)^2
# Try: seq(TestF(n),n=0..6);
# seq(Fib(n+1)^2+Fib(n)^2,n=1..10);
# seq(Fib(2*n+1),n=1..10);
TestF := proc(n) option remember;
local j,P,A,B;

A:= Matrix([[1,1],[0,0]]):
P:= Matrix([[1,1],[1,0]]):
B:= Matrix([[0,1],[0,1]]):

A.P^(2*n).B;
end:



################################
# Section 2.3: 
# Fibonacci Quaternions by Iyer
################################

# Try: [seq(FQ(n),n=1..10)];
FQ := proc(n) option remember;
Fib(n)+Fib(n+1)*I+Fib(n+2)*J
+Fib(n+3)*K;
end:

# Try: [seq(FQMat(n),n=1..5)];
FQMat := proc(n) option remember;
local Q,R;

Q := Matrix([[Fib(n+1),Fib(n)]
,[Fib(n),Fib(n-1)]]):
R := Matrix([[1+I+2*J+3*K,I+J+2*K]
,[I+J+2*K,1+J+K]]):

MulQ(Q,R);
end:

# inverse of noncommuative entry 2x2 matrix
# Try: InverseR();
InverseR := proc() option remember;
local A,B,C,D;

A := 1+I+2*J+3*K;
B := I+J+2*K;
C := I+J+2*K;
D := 1+J+K;

Matrix([[(A-B*D^(-1)*C)^(-1),
(C-D*B^(-1)*A)^(-1)],
[-D^(-1)*C*(A-B*D^(-1)*C)^(-1),
-B^(-1)*A*(C-D*B^(-1)*A)^(-1)]]);

end:


#########################
# Section 2.3.1: 
# Identities by Iyer
#########################

# Have to turn on Quaternions package
# Identity 2
# Try: seq(Iyer17(n),n=1..10);
Iyer17 := proc(n) option remember;
(FQ(n)-I*FQ(n+1)-J*FQ(n+2)-K*FQ(n+3))
/Lucas(n+3);
end:

# Have to turn on Quaternions package
# Identity 4
# Try: seq(Iyer19(n),n=0..5);
Iyer19 := proc(n) option remember;
local L,R;
L := FQ(n)^2+FQ(n-1)^2;
R := 2*FQ(2*n-1)-3*Lucas(2*n+2);

[L,R]
end:



# Matrix form of Identity 4
# Try: IyerMat19();
IyerMat19 := proc() option remember;
local L,R;
L := MulQ(FQMat(-1),FQMat(0));
R := AddQ(CQ(2,FQMat(-1)),
CQ(-3,LucasMat(2)));

[L,R]
end:


# Identity 5
# Try: seq(Iyer20(n),n=0..5);
Iyer20 := proc(n) option remember;
local L,R;
L := FQ(n+1)^2-FQ(n-1)^2;
R := 2*FQ(2*n)-3*Lucas(2*n+3);

[L,L-R]
end:



# Matrix form of Identity 5
# Try: IyerMat20();
IyerMat20 := proc() option remember;
local L,R;
L := MulQ(FQMat(0),FQMat(0));
R := AddQ(CQ(2,FQMat(0)),
CQ(-3,LucasMat(3)));

[L,R]
end:


# Identity 6
# Try: seq(Iyer21(n),n=0..5);
Iyer21 := proc(n) option remember;
local L,R;
L := FQ(n)*FQ(n+1)+FQ(n-1)*FQ(n-2);
R := 6*Fib(n)*FQ(n-1)-9*Fib(2*n+2)
+2*(-1)^(n+1)*(1-I-K);

[L,L-R]
end:


# Matrix form of Identity 6
# Try: seq(IyerMat21(n),n=0..4);
IyerMat21 := proc(n) option remember;
local P,A,B;
P := Matrix([[1,1],[1,0]]):
A := Matrix([[Fib(1),Fib(-2)],[Fib(0),Fib(-1)]]):
B := Matrix([[Fib(1),Fib(0)],[Fib(-2),Fib(-3)]]):

P^n.B.P^n
end:



# Identity 7
# Try: seq(Iyer22(n),n=-3..7);
Iyer22 := proc(n) option remember;
local L,R;
L := FQ(n-1)*FQ(n+3)-FQ(n+1)^2;
R := (-1)^n*(2+4*I+6*J+K);

[L,L-R]
end:


# Identity 8
# Try: seq(Iyer23(n),n=-3..7);
Iyer23 := proc(n) option remember;
local L,R;
L := FQ(n+2)*FQ(n-2)-FQ(n+1)*FQ(n-1);
R := 2*(-1)^(n)*(-2+I-J-6*K);

[L,L-R];
end:


# Identity 9
# Try: seq(Iyer24(n),n=-3..7);
Iyer24 := proc(n) option remember;
local L,R;
L := FQ(n)*FQ(n+1)+FQ(n-2)*FQ(n-3);
R := 4*FQ(2*n-2)-6*Lucas(2*n+1)
+2*(-1)^(n+1)*(I+J-K);

[L,L-R];
end:


# Identity 10
# Try: seq(Iyer25(n),n=-3..7);
Iyer25 := proc(n) option remember;
local L,R;
L := FQ(n+1)^2+FQ(n-1)^2;
R := 6*Fib(n+1)*FQ(n-1)-9*Fib(2*n+3)
+2*(-1)^(n+1)*(1-I-K);

[L,L-R];
end:


# Matrix form of Identity 10
# Try: IyerMat25();
IyerMat25 := proc() option remember;
local A,B;
A := Matrix([[FQ(1),FQ(-1)],[FQ(0),FQ(-2)]]):
B := Matrix([[FQ(1),FQ(0)],[FQ(-1),FQ(-2)]]):
MulQ(A,B);
end:

# Identity 11
# Try: seq(Iyer26(n,3),n=-3..7);
Iyer26 := proc(n,r) option remember;
local L,R;
L := FQ(n+r)+(-1)^r*FQ(n-r);
R := Lucas(r)*FQ(n);

[L,L-R];
end:



################################
# Section 2.3.2: 
# Fibonacci Quaternions, Halici
################################

# Try: seq(Halicib(n),n=0..7);
Halicib := proc(n) option remember;
local i,L,R;
L := 0;
for i from 0 to n do
	L := L+FQ(2*i);
od:
R := FQ(2*n+1)-(1+J+K);

[L,simplify(L-R)];
end:


# Try: seq(Halici13(n),n=0..7);
Halici13 := proc(n) option remember;
local i,L,R;

L := add(binomial(n,i)*FQ(i),i=0..n);
R := FQ(2*n);

[L,simplify(L-R)];
end:



#########################################
# Multiply two matrix
# Try: 
# R := Matrix([[1+I+2*J+3*K,I+J+2*K]
# ,[I+J+2*K,1+J+K]]):
# MulQ(R,R);
MulQ := proc(A,B) option remember;
local a,b,c,d;

a := A[1,1]*B[1,1]+A[1,2]*B[2,1];
b := A[1,1]*B[1,2]+A[1,2]*B[2,2];
c := A[2,1]*B[1,1]+A[2,2]*B[2,1];
d := A[2,1]*B[1,2]+A[2,2]*B[2,2];

Matrix([[a,b],[c,d]]):
end:

# Output: A^n
# Try: R := Matrix([[1,I],[I,1-I]]):
# PowQ(R,3);
PowQ := proc(A,n) option remember;
if n = 0 then
	Matrix([[1,0],[0,1]]);
else
	MulQ(A,PowQ(A,n-1));
fi:
end:

# Try: AddQ(R,R);
AddQ := proc(A,B) option remember;
local a,b,c,d;

a := A[1,1]+B[1,1];
b := A[1,2]+B[1,2];
c := A[2,1]+B[2,1];
d := A[2,2]+B[2,2];

Matrix([[a,b],[c,d]]):
end:

CQ := proc(c,A) option remember;
Matrix([[c*A[1,1],c*A[1,2]]
,[c*A[2,1],c*A[2,2]]]):
end:





################################
# Section 3: Two Matrix Multiply
################################

#############################
# Section 3.1:
# Gaussian Fibonacci Numbers
# by Berzsenyi
#############################

# Output: Value of F_{n+m*I}
# Try: [seq(Berz(n,1),n=-2..7)];
Berz := proc(n,m) option remember;
local k;

add(binomial(m,k)*I^k*CF(n-k,[1,1],[0,1])
,k=0..m);
end:


# Matrix stuffs

# Generate GFib with 2 
# constant Matrices
# Try: PMat(3);
PMat := proc(n) option remember;
Matrix([[1,1],[1,0]])^n;
end:

# Generate GFib with 2 
# constant Matrices
# Try: RMat(3);
RMat := proc(m) option remember;
Matrix([[1,I],[I,1-I]])^m;
end:

# Try: BerMat(2,3);
BerMat := proc(n,m) option remember;
Matrix([[Berz(n+1,m),Berz(n,m)]
,[Berz(n,m),Berz(n-1,m)]]);

end:

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

# Corollary 3.1.2
# Try: Berz312(n,s,m,t);
Berz312 := proc(n,s,m,t) option remember;
local F;
[F(n+s,m+t),
F(n+1,m)*F(s,t)+F(n,m)*F(s-1,t)];
end:


# Try: seq(Am(m),m=0..20);
Am := proc(m) option remember;

if m = 0 then
	return(1);
fi:

Am(m-1)-2*Bm(m-1);
end:
	
# Try: seq(Bm(m),m=0..20);
Bm := proc(m) option remember;

if m = 0 then
	return(0);
fi:

Bm(m-1)+2*Am(m-1);
end:


# Identity 2
# Try: Ber9a(5,F,n);
Ber9a := proc(m,F,n) option remember;
add( binomial(2*m,2*k)*(-1)^k*F(n-2*k)
,k=0..m)=Am(m)*F(n-m);
end:


# Identity 2
# Try: Ber9b(5,F,n);
Ber9b := proc(m,F,n) option remember;
add( binomial(2*m+2,2*k+1)*(-1)^k*F(n-2*k)
,k=0..m)=Bm(m+1)*F(n-m);
end:


# Corollary 3.1.3, 
# Generalize Cassini's identity
# Try: seq(GCas(5,m),m=0..6);
GCas := proc(n,m) option remember;

[Berz(n+1,m)*Berz(n-1,m)-Berz(n,m)^2,
(-1)^n*(2-I)^m];
end:


# Corollary 3.1.4, 
# Generalize Vajda's identity
# Try: seq(seq(GVaj(5,-4,5,m,k,5),m=0..6)
# ,k=0..4);
GVaj := proc(n,i,j,m,k,l) option remember;
local A,B;

A :=factor(Berz(n+i,m+k)*Berz(n+j,m+l)
-Berz(n,m)*Berz(n+i+j,m+k+l));
B :=(-1)^n*(2-I)^m*Berz(i,k)*Berz(j,l);

[A,B/A];
end:




######################################
# Section 3.2 Generalize Iyer,
# (I call it Berzenyi-Lucas, BL)
# Case a=1, b=2
######################################

# S_{n,m} as define in the definition
# Try: BL(3,2);
BL := proc(n,m) option remember;
local a,b,N,M,R;

if m = 0 then
	Fib(n);
elif m < 0 then
	ERROR("BAD INPUT");
else
	2*BL(n+1,m-1)-BL(n,m-1);
fi
end:


# Corollary 3.2.2
# Try: [seq(Cor322BL(n),n=0..13)];
# [seq(BL(n,1),n=0..13)];
Cor322BL := proc(n) option remember;
Fib(n)+2*Fib(n-1);
end:


# Theorem 3.2.5
# Try: [seq(Thm325BL(n,2),n=0..13)];
# [seq(BL(n,2),n=0..13)];
Thm325BL := proc(n,m) option remember;
local j;
add( binomial(m,j)*2^j*Fib(n-j),j=0..m);
end:


# Corollary 3.2.6 (First one)
# Generate identity of F_n
# Try: Cor3261BL(n,5);
Cor3261BL := proc(n,m) option remember;
local F;
5^m*F(n)=
add( binomial(2*m,k)*2^k*F(n-k),k=0..2*m);

end:

# Corollary 3.2.6 (Second one)
# Generate identity of L_n and F_n
# Try: Cor3262BL(n,5);
Cor3262BL := proc(n,m) option remember;
local L,F;
5^m*L(n)=
add( binomial(2*m+1,k)*2^k*F(n-k),k=0..2*m+1);
end:


##################################
# Section 3.4 Fibonacci as Sum of
# Binomial 
##################################

# Table
# Search for b and c for given
# exponents A,B
# Try: FindID(1,2);
# FindID(3,2);
FindID := proc(A,B) option remember;
local a,b,Q,R,C,S;
Q := Matrix([[1,1],[1,0]]);
R := Matrix([[1,b],[b,1-b]]);

C := simplify(Q^A.R^B);
S := {};

for a in {solve(C[1,2]=0,b)} do
	S := S union {[a
	,factor(simplify(subs(b=a,C)))]};
od:
S;
end:



# Big Theorem 3.4.1 
# Try: 
# seq(seq(Thm341(n,m,0,2,2,5),n=3..5),m=6..9);
# seq(seq(Thm341(n,m,1,2,I,1+2*I)
# ,n=3..5),m=6..9);
# seq(seq(Thm341(n,m,2,1,-1,1),n=3..5),m=6..9);
# seq(seq(Thm341(n,m,2,2,-1/2,5/4)
# ,n=3..5),m=6..9);
# seq(seq(Thm341(n,m,3,1,-2,-1)
# ,n=3..5),m=6..9);
Thm341 := proc(n,m,A,B,b,c) option remember;
local j,L,R;

L:=c^m*Fib(n);
R:=add( binomial(B*m,j)*b^j*Fib(n+A*m-j)
,j=0..B*m);

[L,L-R];
end:


# Theorem 3.4.1 write out in n
# Try: 
# seq(print(Thm341W(m,0,2,2,5)),m=1..4);
# seq(print(Thm341W(m,1,2,I,1+2*I)),m=1..4);
# seq(print(Thm341W(m,2,1,-1,1)),m=1..4);
# seq(print(Thm341W(m,2,2,-1/2,5/4)),m=1..4);
# seq(print(Thm341W(m,3,1,-2,-1)),m=1..4);
Thm341W := proc(m,A,B,b,c) option remember;
local k,n,F;

c^m*F(n-A*m)
-add( binomial(B*m,k)*b^k*F(n-k)
,k=0..B*m)=0;

end:



# Theorem 3.4.3
# Try: [seq([seq(Thm343(m,n,-3,0,2,2,5)[3]
# ,m=1..7)],n=0..5)];
# [seq([seq(Thm343(m,n,-3,2,1,-1,1)[3]
# ,m=1..7)],n=0..5)];
Thm343 := proc(m,n,s,A,B,b,c) 
option remember;
local L,R;

L := c^(n*s)*MyT(B*m*n,0,c);
R := add(binomial(B*n,j)*MyT(m,s,b)^j
*MyT(m-1,s,b)^(B*n-j)*MyT(A*n*s+j,0,c)
,j=0..B*n);

[L,R,L-R];
end:


# Try: MyT(4,1,2);
MyT := proc(n,m,b) option remember;
local P,R;

P := Matrix([[1,1],[1,0]]);
R := Matrix([[1,b],[b,1-b]]);

((P^n).(R^m))[1,2];
end:



# Proposition 3.4.7
# where T0,T1,T2,T3 are for 
# s=0,1,2,3 respectively
# Try:
# seq(Prop347(n,0,4),n=1..10);
# seq(Prop347(n,2,b),n=1..10);
Prop347 := proc(n,m,b) option remember;
local t,j,G,C1,C2,T0,T1,T2,T3;

G := Matrix([[1,b],[b,1-b]])^m;

C1 := Fib(n+1)*G[1,2]+Fib(n)*G[2,2];
C2 := Fib(n)*G[1,1]+Fib(n-1)*G[1,2];
T0 := add(binomial(m,j)*(1-b)^(m-j)*
b^j*Fib(n+j),j=0..m);
T1 := add(binomial(m,j)*b^(m-j)*Fib(n-m+j)
,j=0..m);
T2 := add(binomial(m,j)*(b+1)^j*Fib(n-2*m+j)
,j=0..m);
T3 := add(binomial(m,j)*(b+1)^(m-j)*
(b+2)^j*Fib(n-3*m+j)
,j=0..m);

[simplify(C1),simplify(C2),simplify(T0)
,simplify(T1),simplify(T2),simplify(T3)];
end:


# Proposition 3.4.7
# New Identities that arise from
# combining the two binomials on
# the right hand sides
# This one we compare for s=2,3
# Try:
# seq(NewId347(m,n,b,F),m=1..10);
NewId347 := proc(m,n,b,F) option remember;
local t,T2,T3,T;

#T1 := add(binomial(m,j)*b^j*F(n-j),j=0..m);

T2 := add(binomial(m,j)*(b+1)^j*F(n-2*m+j)
,j=0..m);
T3 := add(binomial(m,j)*(b+1)^(m-j)*
(b+2)^j*F(n-3*m+j),j=0..m);

T := simplify((T2-T3)/(b+1));

add(
factor(simplify(coeff(T,F(n-t))))*F(n-t)
,t=0..3*m)=0;
end:



##############################
# Section 4: Multiplication
# of Three (or more) Matrices
##############################

##############################
# Section 4.1: Multiplication
# of Three Matrices
##############################

# Definition 4.1
# Generalize Berzsenyi to three matrix
# Try: B3(1,3,0,1,I,3,1);
B3 := proc(n,m,p,a,b,A,B) option remember;
local P,Q,R;
P := Matrix([[1, 1],[1, 0]]);
Q := Matrix([[a, b],[b, a-b]]);
R := Matrix([[A, B],[B, A-B]]);

P^n.Q^m.R^p;
end:


# The nth entry of the matrix
# Try: NumB3(1,2,1,1,I,3,1);
NumB3 := proc(n,m,p,a,b,A,B) option remember;
B3(n,m,p,a,b,A,B)[1,2];
end:



# Theorem 4.1.4
# Try: VajdaB3(1,2,1,1,I,3,1);
VajdaB3 := proc(n,m,p,a,b,A,B) option remember;
local L,R,i1,i2,i3,j1,j2,j3;

i1:=2:i2:=3:i3:=1:j1:=1:j2:=2:j3:=12:
L := NumB3(n+i1,m+i2,p+i3,a,b,A,B)*
NumB3(n+j1,m+j2,p+j3,a,b,A,B)
-NumB3(n,m,p,a,b,A,B)
*NumB3(n+i1+j1,m+i2+j2,p+i3+j3,a,b,A,B);

R := (-1)^n*(a^2-a*b-b^2)^m*(A^2-A*B-B^2)^p
*NumB3(i1,i2,i3,a,b,A,B)
*NumB3(j1,j2,j3,a,b,A,B);

[L,R];
end:


# Test the matrix of theorem 4.1.4
# Try: VajdaMatTest(1,2,1,1,I,3,1);
VajdaMatTest := proc(n,m,p,a,b,A,B) 
option remember;
local L,R1,R2,i1,i2,i3,j1,j2,j3;

i1:=2:i2:=3:i3:=1:j1:=1:j2:=2:j3:=12:
L := Matrix([[
NumB3(n+i1,m+i2,p+i3,a,b,A,B),
NumB3(n,m,p,a,b,A,B)]
,[ NumB3(n+i1+j1,m+i2+j2,p+i3+j3,a,b,A,B),
NumB3(n+j1,m+j2,p+j3,a,b,A,B)]]);

R1 := Matrix([[
NumB3(n+1,0,0,a,b,A,B),
NumB3(n,0,0,a,b,A,B)]
,[ NumB3(n+j1+1,j2,j3,a,b,A,B),
NumB3(n+j1,j2,j3,a,b,A,B)]]);

R2 := Matrix([[
NumB3(i1,m+i2,p+i3,a,b,A,B),
NumB3(0,m,p,a,b,A,B)]
,[ NumB3(i1-1,m+i2,p+i3,a,b,A,B),
NumB3(-1,m,p,a,b,A,B)]]);
[L,R1.R2];
end:


# Test some Matrix multiplications
# Try: TestHom(2,1);
TestHom := proc(A,B) option remember;
local P,X,Y,C,a,b;

P := Matrix([[1,1],[1,0]]);
X := Matrix([[1,a],[a,1-a]]);
Y := Matrix([[1,b],[b,1-b]]);

C := simplify((X^A).(Y^B));
solve(C[1,2]=0);
[simplify(P.X.Y),
simplify(P+a*X+b*Y)
,simplify(X.Y)];
end:





#############################
# Section 4.2: Quaternions
#############################

# Generalize Berzsenyi to Quaternions

# Definition 4.2
# Try: BQ(1,3,0,0);
BQ := proc(n,m,p,q) option remember;
local N,M,P,Q;
N := Matrix([[1, 1],[1, 0]]);
M := Matrix([[1, I],[I, 1-I]]);
P := Matrix([[1, J],[J, 1-J]]);
Q := Matrix([[1, K],[K, 1-K]]);

MulQ(
MulQ(PowQ(N,n),PowQ(M,m)),
MulQ(PowQ(P,p),PowQ(Q,q))
);
end:


# Number Q_{n,m,p,q} from the matrix
# Try: NumBQ(1,2,1,0);
NumBQ := proc(n,m,p,q) option remember;
if n =-1 then
	BQ(0,m,p,q)[2,2];
elif m =-1 then
	BQ(n,0,p,q)[2,2];
elif p =-1 then
	BQ(n,m,0,q)[2,2];
elif q =-1 then
	BQ(n,m,p,0)[2,2];
else
	BQ(n,m,p,q)[1,2];
fi:
end:


# Have to read Quaternions package
# Corollary 4.2.1 line 2
# Try: Thm11BQ(13,3,5,7);
Thm11BQ := proc(n,m,p,q) option remember;

NumBQ(n,m,p,q)-I*NumBQ(n+1,m-1,p,q)
-(1-I)*NumBQ(n,m-1,p,q);
end:

# Similar to Corollary 4.2.1 line 2
# Try: Thm12BQ(13,3,5,7);
Thm12BQ := proc(n,m,p,q) option remember;

NumBQ(n,m,p,q+1)-NumBQ(n+1,m,p,q)*K
-NumBQ(n,m,p,q)*(1-K);
end:


# Similar to Corollary 4.2.1 line 2
# Try: Thm13BQ(13,3,5,7);
Thm13BQ := proc(n,m,p,q) option remember;
local L,R;

L := I*(NumBQ(n,m,p,q+1)-NumBQ(n,m,p,q));
R := (NumBQ(n,m+1,p,q)-NumBQ(n,m,p,q))*K;
[L,L-R];
end:



# Corollary 4.2.1 line 4,5
# Try: Cor421(13,3,5,7);
Cor421 := proc(n,m,p,q) option remember;
local L,R1,R2;

L := NumBQ(n,m,p,q);
R1 := NumBQ(n+1,0,0,0)*NumBQ(0,m,p,q)+
NumBQ(n,0,0,0)*NumBQ(-1,m,p,q);
R2 := NumBQ(n+1,m,0,0)*NumBQ(0,0,p,q)+
NumBQ(n,m,0,0)*NumBQ(-1,0,p,q);
[L,R1,R2];
end:



# Theorem 4.2.3, identity 2
# Try: 
# seq(seq(Thm2BQ(13,m,5,q),m=0..5),q=4..9);
Thm2BQ := proc(n,m,p,q) option remember;
local s1,s2,s3,L,R;

L := NumBQ(n,m,p,q);
R := 0;

for s1 from 0 to m do
for s2 from 0 to p do
for s3 from 0 to q do
	R:=R+binomial(m,s1)*binomial(p,s2)*
	binomial(q,s3)*I^s1*J^s2*K^s3*
	Fib(n-s1-s2-s3);
od:od:od:

[L,L-R];
end:



# Theorem 4.2.3, identity 5
# Try: 
# seq(seq(Thm5BQ(13,m,5,q),m=0..5),q=4..9);
Thm5BQ := proc(n,m,p,q) option remember;
local L,R;

L:= NumBQ(n,2*m,p,q);
R:= (1+2*I)^m*NumBQ(n-m,0,p,q);
[L,L-R];
end:



# Theorem 4.2.3, identity 6
# Try: 
# seq(seq(Thm6BQ(13,m,5,q),m=0..5),q=4..9);
Thm6BQ := proc(n,m,p,q) option remember;
local L,R;

L:=NumBQ(n,2*m+1,p,q);
R:= (1+2*I)^m*
(NumBQ(n-m,0,p,q)+I*NumBQ(n-m-1,0,p,q));
[L,L-R];
end:



# Corollary 4.2.4
# It could be that p and q have to be 0
# Try: 
# seq(seq(CassiniBQ(13,m,5,q),m=0..5),q=4..9);
CassiniBQ := proc(n,m,p,q) option remember;
local L,R;

L := NumBQ(n+1,m,p,q)*NumBQ(n-1,m,p,q)
-NumBQ(n,m,p,q)^2;
R := (-1)^n*(2-I)^m*(2-J)^p*(2-K)^q;  

[L,R];
end:





################################
# Section 5: Numbers Related
# to 3-by-3 Matrix
################################

##################################
# Section 5.1: Tribonacci Numbers
##################################

# Tribonacci numbers
# Try: [seq(Tribo(n),n=-7..7)];
Tribo := proc(n) option remember;
CF(n,[1,1,1],[0,0,1]);
end:


# Tribonacci with flexible
# initial conditions
# Try: [seq(TriGen(n,3,4,5),n=-7..7)];
TriGen := proc(n,a,b,c) option remember;
CF(n,[1,1,1],[a,b,c]);
end:


# Matrix of the form P^n.R
# Try: [seq(TriMat(n,3,4,5),n=-7..3)];
# Trucas
# [seq(TriMat(n,3,1,3),n=-7..3)];
TriMat := proc(n,a,b,c) option remember;
local N,M;
N := Matrix([[1,1,1],[1,0,0],[0,1,0]]);
M := Matrix([[c,b+a,b],[b,c-b,a],[a,b-a,c-b-a]]);
(N^n).M;
end:



# Some calculation I did awhile ago
# Try: Id52();
Id52 := proc() option remember;
local N,M;
N := Matrix([[1,1,1],[1,0,0],[0,1,0]]);
M := Matrix([[3,4,1],[1,2,3],[3,-2,-1]]);
print(Determinant(M)):
[N+2+3*N^(-1),M]:
end:


##############################
# Section 5.2: Trucas Numbers
##############################

# Generalize Iyer, Lucas
# to 3-by-3 matrix
# Iyer, Lucas: L^3=4(L+1):
# Try: Iyer3(2,3);
Iyer3 := proc(A,B) option remember;
local M,L;

M := Mat3(0,1,1);
L := Matrix([[1,2,1],[1,0,1],[1,0,-1]]);

M^A.L^B;
end:


# Trucas matrix
# Trucas: T^3=4(T^2+11):
# Try: Trucas3(2,3);
Trucas3 := proc(A,B) option remember;
local M,T;

M := Mat3(0,1,1);
T := Mat3(1,1/3,1);

M^A.T^B;
end:



#############################
# Section 5.3 Multiplication
# of Two Matrices
#############################

# The general form of 3-by-3 
# generator matrix
# Try: Mat3(0,I,1);
Mat3 := proc(a,b,c) option remember;
Matrix([[c,b+a,b],[b,c-b,a],[a,b-a,c-b-a]]);
end:


# Main Table
# Search for a,b and c for given
# exponents A,B of 3x3 matrix 
# Try: seq(print(FindID3(A,1)),A=1..12);
FindID3 := proc(A,B) option remember;
local a,b,c,M,R,C,S,T;

M := Mat3(0,1,1);
R := Mat3(a,b,1);

C := simplify(M^A.R^B);
S := solve({C[1,2]=0,C[1,3]=0});

if nops([S]) = 0 then
	return();
elif nops([S])= 1 then 
	T := S union {c=subs(S,C[1,1])};
else
	T := S[1] union {c=subs(S,C[1,1])};
fi:

#print("Exponent",[A,B]);
subs(T,simplify([a,b,c]));
end:


# Double check M^A.R^B = cI
# Try: CheckID3(7,1,-4/3,-1/3,-1/3); 
CheckID3 := proc(A,B,a,b,c) option remember;
local M,R,T;

M := Mat3(0,1,1);
R := Mat3(a,b,1);
T := Matrix([[1,0,0],[0,1,0],[0,0,1]]);

[simplify(M^A.R^B),c*T];
end:


# Theorem 5.3.2
# Try: Thm532(15,19,6,1);
# [seq(Thm532(n,19,6,1),n=1..15)];
Thm532 := proc(n,m,A,B) option remember;
local i,j,L,R,S,a,b,c;

S := FindID3(A,B);
a := S[1];
b := S[2];
c := S[3];

L := c^m*Tribo(n);
R := add(add(
(B*m)!/i!/j!/(B*m-i-j)!*(a+b)^i*b^j*
Tribo(n+A*m-i-2*j)
,j=0..B*m-i),i=0..B*m);
[L,R];
end:




#######################
# Section 6: Complex 
# Fibonacci numbers by
# Harman
#######################

# Lucas numbers
# Try: seq(Lucas(n),n=0..10);
Lucas := proc(n) option remember;
CF(n,[1,1],[2,1]);
end:


# Complex Fibonacci numbers
# as defined by Harman
# Try: Har(5,4);
Har := proc(n,m) option remember;
if m = 0 then 
	Fib(n);
elif m= 1 then 
	Fib(n)+I*Fib(n+1);
elif m >= 2 then
	Har(n,m-1)+Har(n,m-2);
elif m < 0 then
	Har(n,m+2)-Har(n,m+1);
else
	ERROR(NoCase);
fi:

end:


# Harman Matrix
# Try: HarMat(2,3);
HarMat := proc(n,m) option remember;
Matrix([[Har(n+1,m+1), Har(n,m+1)]
,[Har(n+1,m), Har(n,m)]]);

end:


# Generate HarMat with 3
# contant Matrices
# Try: HarPG(14,23); HarMat(14,23);
HarPG := proc(n,m) option remember;
local P,R;

P := Matrix([[1,1],[1,0]]);
R := Matrix([[1+I, I],[1,0]]);

P^m.R.P^n;
end:


##########################
# Section 6.1: Identities
##########################

# Identity 6, Harman(2.7)
# Try: Har27(4,3); Har(4,3);
Har27 := proc(n,m)
Fib(n)*Fib(m+1)+I*Fib(n+1)*Fib(m);
end:

# Identity 8
# Try: [seq(Har8(5,m),m=1..10)];
Har8 := proc(n,m) option remember;
[Har(n,m)+Har(n-1,m-1),(1+I)*Fib(n+m)];
end:


# Intermediate calculation 
# for Identity 8.
# Try: LeftHar(2,3);
LeftHar := proc(n,m) option remember;
Matrix([[Har(n+1,m+1)+Har(n,m)
, Har(n,m+1)+Har(n-1,m)]
,[Har(n+1,m)+Har(n,m-1)
, Har(n,m)+Har(n-1,m-1)]]);
end:


# Try: RightHar(2,3);
RightHar := proc(n,m) option remember;
local id1,P,S;

id1 := Matrix([[1, 0],[0, 1]]);
P := Matrix([[1, 1],[1, 0]]);
S := Matrix([[1+I, I],[I, 1]]);

P^(m-1+n).(id1+P.S);
end:



# Identity 9, Harman(3.2)
# Try: [seq(Har32(5,m),m=1..10)];
Har32 := proc(n,m) option remember;
[Har(n,m)-Har(m-1,n-1),
(1+I)*Fib(n)*Fib(m)];
end:



##########################
# Section 6.2: Constant
# Matrix Representation
##########################

# Matrix of Harman number as defined
# Try: HarG(2);
HarG := proc(n) option remember;
Matrix([[Har(n,1), Har(n,0)]
,[Har(n,0), Har(n,-1)]]);

end:

# Golden ratio Matrix
# Try: FibS(3);
FibS := proc(n) option remember;
Matrix([[(1+sqrt(5))/2,0]
,[0,(1+sqrt(5))/2]])^n;
end:

# The beta Matrix part
# Try: simplify(FibS(3)-FibT(3));
FibT := proc(n) option remember;
Matrix([[(1-sqrt(5))/2,0]
,[0,(1-sqrt(5))/2]])^n;
end:


# Another Diagonal matrix
# Try: seq(FibDig(n),n=0..10);
FibDig := proc(n) option remember;
simplify(
Matrix([[1/sqrt(5),0]
,[0,1/sqrt(5)]]).(FibS(n)-FibT(n))
);
end:


# Harman as a sum of three matrices
# Try: seq([SplitHar(n),HarG(n)],n=0..10);
SplitHar := proc(n) option remember;
local a,b,na,nb;

a := 1/10*(-(5-sqrt(5))+(5+sqrt(5))*I);
na := Matrix([[a,0],[0,a]]);
b := 1/10*(-(5+sqrt(5))+(5-sqrt(5))*I);
nb := Matrix([[b,0],[0,b]]);

simplify(
PMat(n)+na.FibS(n)+nb.FibT(n)
);

end:


############################
# Section 6.3: Third Order 
# Recurrence
############################

# Matrix definition of 
# this third order numbers
# Try: ThirdHarMat(0,4);
ThirdHarMat := proc(n,m) option remember;
local M,N,A;

M := Matrix([[1,1,1],[1,0,0],[0,1,0]]);
N := Matrix([[1,1,0],[1,0,1],[1,0,0]]);
A := matrix([[2+2*I,2+I,I],[1+2*I,1+I,I],[1,1,0]]); 

M^m.A.N^n;
end:


# Hardamard third order numbers
# Try: ThirdHar(10,12);
ThirdHar := proc(n,m) option remember;
ThirdHarMat(n,m)[3,3];
end:


# Identities 1 and 2
# Try: ThirdHar12(10,12); 
ThirdHar12 := proc(n,m) option remember;
local A;

if n in {0,1,2} and m in {0,1,2} then
	A := matrix([[2+2*I,2+I,I],[1+2*I,1+I,I],[1,1,0]]); 
	return(A[3-m,3-n]);
elif n <0 or m<0 then
	ERROR(NoCaseYet);
	
fi:

if n >= 3 then 
	ThirdHar(n-1,m)+ThirdHar(n-2,m)+ThirdHar(n-3,m);
elif m >= 3 then 
	ThirdHar(n,m-1)+ThirdHar(n,m-2)+ThirdHar(n,m-3);
else
	ERROR(NoCaseYet);
fi:	
end:


# Identity 5
# Try:   ThirdHar5(13,24);
ThirdHar5 := proc(n,m)   option remember;
[ThirdHar(n,m), 
Tribo(m+1)*ThirdHar(n,1)+
(Tribo(m)+Tribo(m-1))*ThirdHar(n,0),
Tribo(n+1)*ThirdHar(1,m)+
(Tribo(n)+Tribo(n-1))*ThirdHar(0,m)];
end:


# Identity 6
# Try: ThirdHar6(13,24);
ThirdHar6 := proc(n,m)   option remember;
[ThirdHar(n,m), 
Tribo(n+1)*Tribo(m+2)+I*Tribo(n+2)*Tribo(m+1)];
end:



############################
# Section 6.4: Higher 
# Dimensions
############################

# Harman numbers for 3 dimensions
# Try: Har3d(3,5,4);
Har3d := proc(l,m,n) option remember;
if {l,m,n} subset {0,1} then 
	[l,m,n];

elif l >= 2 then
	Har3d(l-1,m,n)+Har3d(l-2,m,n);
elif l < 0 then
	Har3d(l+2,m,n)-Har3d(l+1,m,n);

elif m >= 2 then
	Har3d(l,m-1,n)+Har3d(l,m-2,n);
elif m < 0 then
	Har3d(l,m+2,n)-Har3d(l,m+1,n);

elif n >= 2 then
	Har3d(l,m,n-1)+Har3d(l,m,n-2);
elif n < 0 then
	Har3d(l,m,n+2)-Har3d(l,m,n+1);
else
	ERROR(NoCase);
fi:

end:


# Theorem 6.4.1
# Try: ForHar3d(3,5,4);
ForHar3d := proc(l,m,n) option remember;

[Fib(l)*Fib(m+1)*Fib(n+1),
Fib(l+1)*Fib(m)*Fib(n+1),
Fib(l+1)*Fib(m+1)*Fib(n)];

end: