######################################################################
## PeterL.txt Save this file as  PeterL.txt  to use it,              #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `Lacrombe.txt` <Enter>                       #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################

with(linalg):

print(`First Written: June 2022: tested for Maple 2018 `):
print():
print(`This is PeterL.txt, A Maple package`):
print(`generating Shalsoh B. Ekhad's  article: `):
print(` Extending the Larcombe Matrix for Cohort Sequences of Solutions of Linear Recurrences with Coefficients to Higher Order`):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/PeterL.txt .`):

ezra:=proc()
if args=NULL then
print(` : A Maple package for  discovering higher-order analogs of Larcombe relations`):

print(`The MAIN procedures are`):
print(`  Bmat, Paper, Paper11, Paper1a, Paper1 `):

elif nargs=1 and args[1]=Bmat then
print(`Bmat(P,A,V): Given vectors P,A,V of the same length,k, say, finds the k by k matrix such that`):
print(`BA^T=V^T and BH=HB. Try:`):
print(`Bmat([p,q,r],[c,b,a],[v2,v1,v0]);`):


elif nargs=1 and args[1]=Paper then
print(`Paper(P,A,V): Given vectors P,A,V of the same length,k, say, finds the k by k matrix such that`):
print(`BA^T=V^T and BH=HB. Verbose form of Bmat(P,A,V). Try:`):
print(`Paper([p,q,r],[c,b,a],[v2,v1,v0]);`):


elif nargs=1 and args[1]=Paper11 then
print(`Paper11(P,A,V): Given vectors P,A,V of the same length,k, say, finds the k by k matrix such that`):
print(`BA^T=V^T and BH=HB. Another verbose form of Bmat(P,A,V). Try:`):
print(`Paper11([p,q,r],[c,b,a],[v2,v1,v0]);`):


elif nargs=1 and args[1]=Paper1 then
print(`Paper1(P,A,V): combines Paper11 for orders 2,3,4,5. Try:`):
print(`Paper1();`):

else
 print(`There is no such thing as`, args):

fi:


end:

#Bmat(P,A,V): Given vectors P,A,V of the same length,k, say, finds the k by k matrix such that
#BA^T=V^T and BH=HB. Try:
#Bmat([p,q,r],[c,b,a],[v2,v1,v0]);
Bmat:=proc(P,A,V) local k,i,j,H,b,B,var,eq,k1:
k:=nops(P):
if not (nops(A)=k and nops(V)=k) then
 RETURN(FAIL):
fi:

H:=[P,seq([0$(i-1),1,0$(k-i)],i=1..k-1)]:

var:={seq(seq(b[i,j],j=1..k),i=1..k)}:

B:=[seq([seq(b[i,j],j=1..k)],i=1..k)]:

eq:={seq(add(B[i][j]*A[j],j=1..k)-V[i],i=1..k)}
union {seq(seq(add(B[i][k1]*H[k1][j],k1=1..k)-add(H[i][k1]*B[k1][j],k1=1..k),i=1..k),j=1..k)}:

var:=solve(eq,var):

normal(subs(var,B)):

end:


#Paper(P,A,V): verbose form of Bmat(P,A,V)
#Paper([p,q,r],[c,b,a],[v2,v1,v0]);
Paper:=proc(P,A,V) local k,w,i,n,H,B:

k:=nops(P):
if not (nops(A)=k and nops(V)=k) then
 RETURN(FAIL):
fi:

B:=Bmat(P,A,V):
print(`On the Linking matrix of Cohort Sequences Solving general linear recurrence equations of order`, k, `In the footsteps of Peter Larcombe`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`We are interested in solutions of the GENERAL recurrence of order`, k ):
print(``):
print(w[n+k]=add(P[i]*w[n+k-i],i=1..k)):
print(``):
print(`For n>=0`):
print(``):
print(`with generic initial conditions `):
print(``):
print(seq(w[i-1]=A[k+1-i],i=1..k)):
print(``):
print(`Let H be the square matrix of dimension`, k ):
print(``):
print(H=matrix([P,seq([0$(i-1),1,0$(k-i)],i=1..k-1)])):
print(``):
print(`Following Peter Larcombe, Palestinian J. Mathematics,  11(3) (2022), 659-663, we are interested in a square matrix of dimension`,k, `let's call it B such that `):
print(``):
print(`Property 1:  B times the transpose of`, A, ` equals the transpose of`, V):
print(``):
print(` Property 2: HB=BH `):
print(``):
print(`Solving we get the following square matrix of dimension`, k):
print(``):
print(matrix(B)):
print(``):
print(`and in Maple notation `):
print(``):
lprint(matrix(B)):
end:



#Paper11(P,A,V): verbose form of Bmat(P,A,V)
#Paper11([p,q,r],[c,b,a],[v2,v1,v0]);
Paper11:=proc(P,A,V) local k,w,i,j,n,H,B,alpha,b,beta,Alpha:

k:=nops(P):
if not (nops(A)=k and nops(V)=k) then
 RETURN(FAIL):
fi:

B:=Bmat(P,A,V):
print(`On the Linking matrix of Cohort Sequences Solving general linear recurrence equations of order`, k, `In the footsteps of Peter Larcombe`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`We are interested in solutions of the GENERAL recurrence of order`, k ):
print(``):
print(w[n+k]=add(P[i]*w[n+k-i],i=1..k)):
print(``):
print(`For n>=0`):
print(``):
print(`with generic initial conditions `):
print(``):
print(seq(w[i-1]=A[k+1-i],i=1..k)):
print(``):
print(`Let H be the square matrix of dimension`, k ):
print(``):
print(H=matrix([P,seq([0$(i-1),1,0$(k-i)],i=1..k-1)])):
print(``):
print(`Following Peter Larcombe, Palestinian J. Mathematics,  11(3) (2022), 659-663, we are interested in a square matrix of dimension`,k, `let's call it B such that `):
print(``):
print(`Property 1:  B times the transpose of`, A, ` equals the transpose of`, V):
print(``):
print(` Property 2: HB=BH `):
print(``):


alpha:={seq(seq(denom(B[i][j]),i=1..k),j=1..k)}:

if nops(alpha)<>1 then
 RETURN(FAIL):
else
 alpha:=alpha[1]:
fi:

print(` Let  `):
print(``):
lprint(Alpha=alpha):
print(``):

print(`Solving we get the following square matrix of dimension`, k):
print(``):
print(matrix([seq([seq(b[i,j],j=1..k)],i=1..k)])):
print(``):
print(` where `):
print(``):
print(b[i,j]=beta[i,j]/Alpha):
print(``):
print(`and the beta[i,j] are as follows `):
print(``):

for i from  1 to k do
 for j from 1 to k do
  lprint(beta[i,j]=numer(B[i][j])):
 od:
od:
print(``):
print(`----------------------------------------`):
print(``):

end:




#Paper1a(P,A,V): verbose form of Bmat(P,A,V)
#Paper1a([p,q,r],[c,b,a],[v2,v1,v0]);
Paper1a:=proc(P,A,V) local k,w,i,j,n,H,B,alpha,b,beta,Alpha:

k:=nops(P):
if not (nops(A)=k and nops(V)=k) then
 RETURN(FAIL):
fi:

B:=Bmat(P,A,V):
print(``):
print(`We are interested in solutions of the GENERAL recurrence of order`, k ):
print(``):
print(w[n+k]=add(P[i]*w[n+k-i],i=1..k)):
print(``):
print(`For n>=0`):
print(``):
print(`with generic initial conditions `):
print(``):
print(seq(w[i-1]=A[k+1-i],i=1..k)):
print(``):
print(`Let H be the square matrix of dimension`, k ):
print(``):
print(H=matrix([P,seq([0$(i-1),1,0$(k-i)],i=1..k-1)])):
print(``):
print(`Following Peter Larcombe, Palestinian J. Mathematics,  11(3) (2022), 659-663, we are interested in a square matrix of dimension`,k, `let's call it B such that `):
print(``):
print(`Property 1:  B times the transpose of`, A, ` equals the transpose of`, V):
print(``):
print(` Property 2: HB=BH `):
print(``):


alpha:={seq(seq(denom(B[i][j]),i=1..k),j=1..k)}:

if nops(alpha)<>1 then
 RETURN(FAIL):
else
 alpha:=alpha[1]:
fi:

print(` Let  `):
print(``):
lprint(Alpha=alpha):
print(``):

print(`Solving we get the following square matrix of dimension`, k):
print(``):
print(matrix([seq([seq(b[i,j],j=1..k)],i=1..k)])):
print(``):
print(` where `):
print(``):
print(b[i,j]=beta[i,j]/Alpha):
print(``):
print(`and the beta[i,j] are as follows `):
print(``):

for i from  1 to k do
 for j from 1 to k do
 print(``):
  lprint(beta[i,j]=numer(B[i][j])):
 print(``):
 od:
od:
print(``):
print(`----------------------------------------`):
print(``):

end:




#Paper1():
Paper1:=proc():



print(`On the Linking matrix of Cohort Sequences Solving general linear recurrence equations of order up to the fifth, in the footsteps of Peter Larcombe`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(``):
print(`-------------------------`):
print(``):
print(`For order two we have `):
print(``):

Paper1a([p,q],[b,a],[v1,v0]):



print(``):
print(`-------------------------`):
print(``):
print(`For order three we have `):
print(``):

Paper1a([p,q,r],[c,b,a],[v2,v1,v0]):




print(``):
print(`-------------------------`):
print(``):
print(`For order four we have `):
print(``):

Paper1a([p,q,r,s],[d,c,b,a],[v3,v2,v1,v0]):



print(``):
print(`-------------------------`):
print(``):
print(`For order five we have `):
print(``):

Paper1a([p,q,r,s,t],[e,d,c,b,a],[v4,v3,v2,v1,v0]):
print(``):
print(`----------------------------------`):
print(``):
print(`This ends this paper that took`, time(), `seconds to generate `):
end: