######################################################################
## HANKEL: Save this file as HANKEL.txt to use it,                   #
## open Maple and type 						     #
##	read `Dir/HANKEL.txt`;				             #
## where Dir is the directory where you saved the HANKEL.txt,        #
## and then hit <enter>.                                             #
##                                                                   #
## Written by Emilie Hogan, Rutgers University ,                     #
##  eahogan at math dot rutgers dot edu.                             # 
######################################################################


print(`This is the Maple package HANKEL created by`);
print(`Emilie Hogan`);
print(`to accompany Chapter 3 of her PhD thesis.`);
print(``);
print(`For general help type Help()`);
print(`For help on a specific procedure type Help(<procedure name>)`);
print(`Help on supplementary procedures can be found using Help1`);
print(``);
print(`This package has been tested in Maple 13`);
print(``);

with(LinearAlgebra):

#=============================================================
Hn := proc(n,K) option remember;

if n <= 2*K+1 then
	RETURN(1);
else
	RETURN( (Hn(n-1,K)*Hn(n-2*K,K)+Hn(n-K,K)+Hn(n-(K+1),K) )/Hn(n-(2*K+1),K) );
fi;
end:

#=============================================================
An := proc(n,K,i) option remember;

if n <= (2*K+1)*i then
	RETURN(1);
else
	RETURN( (An(n-i,K,i)*An(n-2*K*i,K,i)+An(n-K*i,K,i)+An(n-(K+1)*i,K,i) )/An(n-(2*K+1)*i,K,i) );
fi;
end:

#=============================================================
En := proc(n,K) option remember;

if n <= 2*K then
	RETURN(1);
else
	RETURN( (En(n-1,K)*En(n-2*K+1,K)+2*En(n-K,K) )/En(n-(2*K),K) );
fi;
end:

#=============================================================
Bn := proc(n,K,i) option remember;

if n <= (2*K)*i then
	RETURN(1);
else
	RETURN( (Bn(n-i,K,i)*Bn(n-(2*K-1)*i,K,i)+2*Bn(n-K*i,K,i) )/Bn(n-(2*K)*i,K,i) );
fi;
end:

#=============================================================
Xn := proc(n,i,j,k) option remember;

if n <= k then
	RETURN(1);
else
	RETURN( (Xn(n-i,i,j,k)*Xn(n-k+i,i,j,k)+Xn(n-j,i,j,k)+Xn(n-k+j,i,j,k) )/Xn(n-k,i,j,k) );
fi;
end:

#=============================================================
SeqToHankel := proc(S,n,m)

RETURN([seq([seq( S[i+j] ,i=n..n+m)],j=0..m)]);

end:

#=============================================================
SeqToHankelDet := proc(S,n,m)

Determinant(Matrix(SeqToHankel(S,n,m)));

end:

#=============================================================
ConjLinOrder := proc(S,m)
local i,n,DH;

n := nops(S);

for i from 1 to min(n-2*m,m) do
	DH := SeqToHankelDet(S,1,i); #print(i,DH);
	if DH=0 then
		if {seq(SeqToHankelDet(S,j,i),j=1..n-2*i)}={0} then
			RETURN(i);
		fi;
	fi;
od;

RETURN(FAIL);

end:

#=============================================================
EigenVec0 := proc(M)
local EV, Evals, cols;

EV := [Eigenvectors(M)];
Evals := convert(EV[1],list); #print(Evals);

cols := select(i->Evals[i]=0,[seq(i,i=1..nops(Evals))]); #print(cols);

RETURN([seq([seq(EV[2][i][c],i=1..nops(Evals))], c in cols)]);

end:

#=============================================================
ConjLinRecur := proc(S,m)
local ord, M, Vs;

ord := ConjLinOrder(S,m);

if ord=FAIL then
	RETURN(FAIL);
fi;

M := Matrix(SeqToHankel(S,1,ord));
Vs := EigenVec0(M);

RETURN(Vs);

end:

#=============================================================
ConjLinRecurVerbose := proc(S,m)
local ord, M, Vs, i, j, v;

ord := ConjLinOrder(S,m);
if ord=FAIL then
	print(`There is no linear recurrence of order <= `,m);
	RETURN(FAIL);
fi;

print(`The conjectured order of the linear recurrence that annihilates S is`, ord);

M := Matrix(SeqToHankel(S,1,ord));
Vs := EigenVec0(M);

print(`The conjectured linear recurrence is given by the following vector(s):`);
for i from 1 to nops(Vs) do
	print(Vs[i]);
od;
print(`The recurrence is then:`);

for j from 1 to nops(Vs) do
	v := Vs[j];
	print(add(v[k]*s[n+k-1],k=1..nops(v))=0);
od;

RETURN(Vs);

end:

Help := proc()
if nargs=0 then
	print(`The main procedures in this package are:`);
	print(`Xn(n,i,j,k), ConjLinOrder(S,m), ConjLinRecur(S,m), ConjLinRecurVerbose(S,m)`);
	print(`For help on a specific procedure type Help(<procedure name>)`);
elif nargs=1 then
if args[1]=Xn then
	lprint(`Xn(n,i,j,k):`);
	lprint(`Inputs:`);
	lprint(` n,i,j,k - positive integers`);
	lprint(`Outputs:`);
	lprint(` The nth term in the sequence produced by the recurrence:`);
	lprint(` x(n)*x(n-k)=x(n-i)*x(n-k+i)+x(n-j)+x(n-k+j) `);
	lprint(` with initial conditions x(m)=1 for m <= k.`);
	lprint(` This recurrence was studied in Emilie Hogan's PhD thesis. Special cases were`);
	lprint(` also studied (h,e,a,b recurrences). Help on these are found in the Help1 procedure`);
	lprint(`For example, try:`);
	lprint(` [seq(Xn(n,1,2,5),n=1..100)];`);
elif args[1]=ConjLinOrder then
	lprint(`ConjLinOrder(S,m):`);
	lprint(`Inputs:`);
	lprint(` S - List `);
	lprint(` m - positive integer`);
	lprint(`Outputs:`);
	lprint(` The conjectured order of a linear recurrence that annihilates the sequence S.`);
	lprint(`For example, try:`);
	lprint(` ConjLinOrder([seq(Xn(n,1,2,5),n=1..100)], 13); `);
elif args[1]=ConjLinRecur then
	lprint(`ConjLinRecur(S,m):`);
	lprint(`Inputs:`);
	lprint(` S - List `);
	lprint(` m - positive integer`);
	lprint(`Outputs:`);
	lprint(` The conjectured linear recurrence that annihilates the sequence S.`);
	lprint(`For example, try:`);
	lprint(` ConjLinRecur([seq(Xn(n,1,1,3),n=1..100)], 13); `);
elif args[1]=ConjLinRecurVerbose then
	lprint(`ConjLinRecurVerbose(S,m):`);
	lprint(`Inputs:`);
	lprint(` S - List `);
	lprint(` m - positive integer`);
	lprint(`Outputs:`);
	lprint(` The conjectured linear recurrence that annihilates the sequence S.`);
	lprint(` Also prints out statements to help interpret the output.`);
	lprint(`For example, try:`);
	lprint(` ConjLinRecurVerbose([seq(Xn(n,1,1,3),n=1..100)], 13); `);
fi;
fi;
end:

Help1 := proc()
if nargs=0 then
	print(`All procedures in this package are:`);
	print(`Hn(n,K), An(n,K,i), En(n,K), Bn(n,K,i), Xn(n,i,j,k)`);
	print(`SeqToHankel(S,n,m), SeqToHankelDet(S,n,m)`);
	print(`ConjLinOrder(S,m), ConjLinRecur(S,m), ConjLinRecurVerbose(S,m)`);
	print(`EigenVec0(M)`);
	print(`For help on a specific procedure type Help1(<procedure name>)`);
elif nargs=1 then
if args[1]=Xn then
	Help(Xn);
elif args[1]=ConjLinRecur then
	Help(ConjLinRecur);
elif args[1]=ConjLinRecurVerbose then
	Help(ConjLinRecurVerbose);
elif args[1]=ConjLinOrder then
	Help(ConjLinOrder);
elif args[1]=Hn then
	lprint(`Hn(n,K):`);
	lprint(`Inputs:`);
	lprint(` n,K - positive integers`);
	lprint(`Outputs:`);
	lprint(` The nth term in the sequence produced by the recurrence:`);
	lprint(` h(n)*h(n-(2*K+1))=h(n-1)*h(n-2*K)+h(n-K)+h(n-(K+1)) `);
	lprint(` with initial conditions h(m)=1 for m <= 2*K+1.`);
	lprint(` This recurrence was studied in Emilie Hogan's PhD thesis.`);
	lprint(`For example, try:`);
	lprint(` [seq(Hn(n,2),n=1..100)];`);
elif args[1]=En then
	lprint(`En(n,K):`);
	lprint(`Inputs:`);
	lprint(` n,K - positive integers`);
	lprint(`Outputs:`);
	lprint(` The nth term in the sequence produced by the recurrence:`);
	lprint(` e(n)*e(n-2*K)=e(n-1)*e(n-(2*K-1))+2*e(n-K) `);
	lprint(` with initial conditions e(m)=1 for m <= 2*K.`);
	lprint(` This recurrence was studied in Emilie Hogan's PhD thesis.`);
	lprint(`For example, try:`);
	lprint(` [seq(En(n,2),n=1..100)];`);
elif args[1]=An then
	lprint(`An(n,K,i):`);
	lprint(`Inputs:`);
	lprint(` n,K,i - positive integers`);
	lprint(`Outputs:`);
	lprint(` The nth term in the sequence produced by the recurrence:`);
	lprint(` a(n)*a(n-(2*K+1)*i)=a(n-i)*a(n-2*K*i)+a(n-K*i)+a(n-(K+1)*i) `);
	lprint(` with initial conditions a(m)=1 for m <= (2*K+1)*i.`);
	lprint(` This recurrence was studied in Emilie Hogan's PhD thesis.`);
	lprint(`For example, try:`);
	lprint(` [seq(An(n,2,2),n=1..100)];`);
elif args[1]=Bn then
	lprint(`Bn(n,K,i):`);
	lprint(`Inputs:`);
	lprint(` n,K,i - positive integers`);
	lprint(`Outputs:`);
	lprint(` The nth term in the sequence produced by the recurrence:`);
	lprint(` b(n)*b(n-2*K*i)=b(n-i)*b(n-(2*K-1)*i)+2*b(n-K*i) `);
	lprint(` with initial conditions b(m)=1 for m <= 2*K*i.`);
	lprint(` This recurrence was studied in Emilie Hogan's PhD thesis.`);
	lprint(`For example, try:`);
	lprint(` [seq(Bn(n,2,2),n=1..100)];`);
elif args[1]=SeqToHankel then
	lprint(`SeqToHankel(S,n,m):`);
	lprint(`Inputs:`);
	lprint(` S - List`);
	lprint(` n,m - positive integers`);
	lprint(`Outputs:`);
	lprint(` The Hankel matrix (list of lists) of size (m+1)x(m+1) starting at`);
	lprint(` position n in S, i.e., the matrix whos (i,j) entry is S[n+i+j-2].`);
	lprint(`For example, try:`);
	lprint(` SeqToHankel([seq(Hn(i,2),i=1..7)],3,2); `);
elif args[1]=SeqToHankelDet then
	lprint(`SeqToHankelDet(S,n,m):`);
	lprint(`Inputs:`);
	lprint(` S - List`);
	lprint(` n,m - positive integers`);
	lprint(`Outputs:`);
	lprint(` The determinant of the matrix SeqToHankel(S,n,m).`);
	lprint(`For example, try:`);
	lprint(` SeqToHankelDet([seq(Hn(i,2),i=1..7)],3,2); `);
elif args[1]=EigenVec0 then
	lprint(`EigenVec0(M):`);
	lprint(`Inputs:`);
	lprint(` M - a matrix (formatted for the LinearAlgebra package)`);
	lprint(`      [use Matrix(L) where L is a list of lists to put]`);
	lprint(`      [the output of SeqToHankel in the correct format]`);
	lprint(`Outputs:`);
	lprint(` The eigenvector corresponding to the 0 eigenvalue `);
	lprint(` of the matrix M (if one exists).`);
	lprint(`For example, try:`);
	lprint(` EigenVec0(Matrix([[1,1,1],[2,2,3],[0,0,1]])); `);
fi;
fi;

end: