# restart; read `d:/Aek/Fibo/DetFibo.txt`;

with(LinearAlgebra):

# First Version: June 2, 2014.
# Second Version: April 16, 2015
# Third Version: December 7, 2015
# Note: generalize F to order 3  
# seems hopeless.

ez:=proc(): 
print(` Section 1: Catalan's Identities:`): 
print(` CF, GCatalan,`):
print(` Eq23, Eq24, Eq25, Eq26`):
print(` `):
print(` Section 2: Main Result with Fibo.`): 
print(` PolMat, LPolMat,`):
print(` GenCarlitz, Carlitz`):
print(` `):
print(` Section 3: Product Form`):
print(` Krat3, VersionKrat3,`):
print(` Force, Alfred, PMat `):
print(` `):
print(` Section 4: Other functions`):
print(` AramMat`):
print(` `):
end:

##########################
# Section 1: Catalan's Identities
# Functions are
# CF, GCatalan,
# Eq23, Eq24, Eq25, Eq26
##########################

# C-finite of order nops(C):
# Input: constants, C,A.
# Try: [seq(CF(n,[1,1],[0,1]),n=0..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:

# Proposition 2.1
# Try: 
# [seq([seq([seq(GCatalan(s,i,j,-2,5,4,-3,5,2)
# ,i=-3..-1)],j=1..3)],s=-3..3)];
GCatalan := proc(s,i,j,c1,c2,a0,a1,b0,b1) 
option remember;
local W,Y,U,n,L,R;

W := n -> CF(n,[c1,c2],[a0,a1]):
Y := n -> CF(n,[c1,c2],[b0,b1]):
U := n -> CF(n,[c1,c2],[0,1]):

L := W(s+i)*Y(s+j)-W(s)*Y(s+i+j);
R := (-c2)^(s)*(W(1)*Y(j)-W(0)*Y(j+1))
*U(i); 

L-R;

#Note: 
#[W(0)*Y(j+1)-W(1)*Y(j),
#(b0*a1*c1+c2*a0*b0-a1*b1)*U(j)
#+(a0*b1-b0*a1)*U(j+1)]:
end:

# Equation 2.3 in the paper
# Try: 
# [seq([seq([seq(Eq23(s,i,j,-2,5)
# ,j=-3..-1)],i=1..3)],s=-3..3)];
Eq23 := proc(s,i,j,c1,c2) 
option remember;
local U,n,L,R;

U := n -> CF(n,[c1,c2],[0,1]):

L := U(s+i)*U(s+j)-U(s)*U(s+i+j);
R := (-c2)^(s)*U(j)*U(i); 

L-R;
end:

# Try: 
# [seq([seq([seq(Eq24(s,i,j,-2,5,4,-3)
# ,i=-3..-1)],j=1..3)],s=-3..3)];
Eq24 := proc(s,i,j,c1,c2,a0,a1) 
option remember;
local W,U,n,L,R;

W := n -> CF(n,[c1,c2],[a0,a1]):
U := n -> CF(n,[c1,c2],[0,1]):

L := U(s+i)*W(s+j)-U(s)*W(s+i+j);
R := (-c2)^(s)*W(j)*U(i); 

L-R;
end:

# Try: 
# [seq([seq([seq(Eq25(s,i,j,-2,5,4,-3)
# ,i=-3..-1)],j=1..3)],s=-3..3)];
Eq25 := proc(s,i,j,c1,c2,a0,a1) 
option remember;
local W,U,n,L,R;

W := n -> CF(n,[c1,c2],[a0,a1]):
U := n -> CF(n,[c1,c2],[0,1]):

L := W(s+i)*W(s+j)-W(s)*W(s+i+j);
R := (-c2)^(s)*(W(1)^2-W(0)*W(2))*U(j)*U(i); 

L-R;
end:


# Equaiton 2.6, Restate Catalan
# Try: 
# [seq([seq(Eq26(s,k,2,-4,12,3,-5,2,3)
# ,k=1..5)],s=-3..3)];
Eq26 := proc(s,k,n,t,r,c1,c2,a0,a1) 
option remember;
local del,N,U,W,R,L;

del := a1^2-a1*a0*c1-c2*a0^2;

U := N -> CF(N,[c1,c2],[0,1]):
W := N -> CF(N,[c1,c2],[a0,a1]):
L := W(s+k*(n+t));
R := W(k*(t-r))/W(-k*r)*W(s+k*n)
-(-c2)^(-k*r)*del*U(k*t)/W(-k*r)*U(s+k*(n+r));

L-R;
end:


##########################
# Section 2: Main result with Fibo.
# Functions are
# PolMat, LPolMat,
# GenCar,Carlitz
##########################

# Proposition 3.1
# Try: PolMat(3);
PolMat := proc(r) option remember;
local c,x,M;

M := Matrix([seq([seq( 
(c[j]*x[i]+1)^r
,j=0..r)],i=0..r)]);
print(M);
factor(Determinant(M));
end:


# Verstion of proposition 3.1
# Try: LPolMat(2);
LPolMat := proc(r) option remember;
local M,A,B,x,y;

M := Matrix([seq([seq( 
(A[j]*x[i]+B[j]*y[i])^r
,j=0..r)],i=0..r)]);
print(M);
factor(Determinant(M));
end:





# Theorem 3.2: Main result
# Main Result: Determinant Power 
# on (r+1)x(r+1) matrix
# Try: 
# [seq([seq([seq(GenCarlitz(r,s,k,-3,3,-5,2,3)
# ,k=1..5)],s=-3..3)],r=1..4)];
GenCarlitz := proc(r,s,k,n,c1,c2,a0,a1) 
option remember;
local i,j,N,U,W,M,del,L,R;

U := N -> CF(N,[c1,c2],[0,1]):
W := N -> CF(N,[c1,c2],[a0,a1]):
M := Matrix([seq([seq( 
W(s+k*(n+i+j))^r,i=0..r)],j=0..r)]);
L:= Determinant(M);

del := a1^2-a1*a0*c1-c2*a0^2;
R := (-1)^((s+k*n+1)*r*(r+1)/2)
*c2^((s+k*n)*r*(r+1)/2+k*(r-1)*r*(r+1)/3)
*del^(r*(r+1)/2)
*mul(U(i*k)^(r+1-i),i=1..r)^2
*mul(binomial(r,i),i=0..r);

L-R:
end:



# Remark 3.3: Main result with Fibo.
# Input: matrix of size r+1 with entries
# F_(s+k*(n+i+j))^r
# Try: 
# {seq(seq(seq(seq(Carlitz(r,s,k,n)
# ,s=-5..5),k=-5..5),n=-3..3),r=1..7)};
Carlitz := proc(r,s,k,n) option remember;
local i,F,N,L,R;

F := N -> CF(N,[1,1],[0,1]);
L := AramMat(s,k,n,r,r+1,[1,1],[0,1]):

R := (-1)^((s+k*n+1)*r*(r+1)/2)
*mul(F((i+1)*k)^(r-i),i=0..r-1)^2
*mul(binomial(r,i),i=0..r);

L-R;
end:





########################
# Section 3: Product Form
# Krat3, VersionKrat3,
# Force, Alfred, PMat
########################

# Lemma 4.1
# Try: seq(Krat3(r),r=0..3);
Krat3 := proc(r) option remember;
local i,j,l,m,X,c,d,M,a,b;

M := Matrix([seq([seq( 
mul(X[i]+c[l],l=j+1..r)
*mul(X[i]+d[m],m=1..j)
,j=0..r)],i=0..r)]);

print(M);
a:=factor(Determinant(M));
b:=mul(mul(X[j]-X[i],i=0..j-1),j=0..r)
*mul(mul(c[j]-d[i],i=1..j),j=1..r);

simplify(a/b)
end:


# Alternative way of Lemma 4.1
# Try: seq(VersionKrat3(r),r=0..3);
VersionKrat3 := proc(r) option remember;
local i,j,l,m,X,Y,A,B,c,d,M,a,b;

M := Matrix([seq([seq( 
mul(A[c[l]]*X[i]+B[c[l]]*Y[i],l=j+1..r)
*mul(A[d[m]]*X[i]+B[d[m]]*Y[i],m=1..j)
,j=0..r)],i=0..r)]);

#print(M);
a:=factor(Determinant(M));
#b:=mul(mul(X[i]*Y[j]-X[j]*Y[i]
#,i=0..j-1),j=0..r)
#*mul(mul(B[d[i]]*A[c[j]]-A[d[i]]*B[c[j]]
#,i=1..j),j=1..r);

#simplify(a/b);
end:




# Theorem 4.2: main result
# Try:
# r:=3: Force(r,3,2,0,1,1,0,1,[0$r],[1$r]);
# seq(Force(r,2,3,0,1,1,0,1,[0$r],[1$r]),r=0..10);
# Force(3,2,3,0,1,1,0,1,[2,3,4],[0,1,2]);
# Force(3,2,3,0,5,3,0,1,[2,3,4],[0,1,2]);
Force := proc(r,s,k,n,c1,c2,a0,a1,C,D) 
option remember;
local i,j,l,m,N,U,W,M,del,L,R,A,B,t;

if nops(C) <> r or nops(D) <> r then
	ERROR("Wrong dimension");
fi:

U := N -> CF(N,[c1,c2],[0,1]):
W := N -> CF(N,[c1,c2],[a0,a1]):
M := Matrix([seq([seq( 
mul(W(s+k*(n+i+C[l])),l=j+1..r)
*mul(W(s+k*(n+i+D[m])),m=1..j)
,j=0..r)],i=0..r)]);

L:= Determinant(M);

del := a1^2-a1*a0*c1-c2*a0^2;
#A := t-> W(k*(t-r))/W(-k*r);
#B := t-> -(-c2)^(-k*r)*del*W(k*t)/W(-k*r);

R := (-del)^binomial(r+1,2)
*(-c2)^((s+k*n)*binomial(r+1,2)
+k*binomial(r+1,3))
*mul(U(i*k)^(r+1-i),i=1..r)
*mul(mul((-c2)^(k*C[j])*
U(k*(D[i]-C[j])),i=1..j),j=1..r);

L/R:
end:


# Corollary 4.3: Alfred determinant
# Try: Alfred(5,0,3);
# {seq(seq(seq(Alfred(r,s,k)
# ,k=-5..5),s=-3..3),r=0..8)};
Alfred := proc(r,s,k) option remember;
local i,j,F,N,M,L,R;

F := N -> CF(N,[1,1],[0,1]):
M := Matrix([seq([seq( 
F(s+k*i)^(r-j)*F(s+k*(i+1))^j
,j=0..r)],i=0..r)]);

L:= Determinant(M);

R := (-1)^((s+1)*r*(r+1)/2+k*binomial(r+1,3))
*F(k)^(r*(r+1)/2)
*mul(F(i*k)^(r+1-i),i=1..r);
L-R:
end:

# Corollary 4.4
# Try:
#seq(seq(PMat(3,s,3,n,5,3,0,1,2),n=0..5),s=-3..4);
PMat := proc(r,s,k,n,c1,c2,a0,a1,p) 
option remember;
local i,j,l,m,N,U,W,M,del,L,R,C,D;

C := [seq(i+p,i=0..r-1)];
D := [seq(i,i=0..r-1)];

U := N -> CF(N,[c1,c2],[0,1]):
W := N -> CF(N,[c1,c2],[a0,a1]):
M := Matrix([seq([seq( 
mul(W(s+k*(n+i+C[l])),l=j+1..r)
*mul(W(s+k*(n+i+D[m])),m=1..j)
,j=0..r)],i=0..r)]);

L:= Determinant(M);

del := a1^2-a1*a0*c1-c2*a0^2;

R := del^binomial(r+1,2)
*(-c2)^((s+k*n)*binomial(r+1,2)
+2*k*binomial(r+1,3))
*mul(U(i*k)^(r+1-i),i=1..r)
*mul(U(k*(p+i))^(r-i),i=0..r-1);

L/R:
end:




########################
# Section 4: Other functions
# AramMat
########################

# General Form of Matrix
# Try: AramMat(0,5,3,1,2,[1,1],[0,1]);
# [seq(AramMat(0,1,n,2,3,[1,1],[0,1]),n=1..10)];
# [seq(2*(-1)^(n+1)*CF(1,[1,1],[0,1])^4
# *CF(2,[1,1],[0,1])^2 ,n=1..10)];
AramMat := proc(s,k,n,pow,t,C,A) 
option remember;
local i,j,M;

M := Matrix([seq([seq( 
CF(s+k*(n+i+j),C,A)^pow
,j=0..t-1)],i=0..t-1)]);
Determinant(M);
end: