######################################################################
##QuasiDeterminants.txt: Save this file as  QuasiDeterminants.txt    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read QuasiDeterminants.txt<Enter>                  #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: June/July 2018
 
print(`Created:  June/July 2018`):
print(` This is QuasiDeterminants.txt `):
print(`It is a package dedicated to Vladimir Retakh on his 70th birthday `):
print(``):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Bdok, CheckCQD1, CheckHered1, inverse1, multiply1, RandMat `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  BU, CheckHered, EvalBM, EvalCom, qd `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Bdok then
print(`Bdok(n): checks that if all the entries commute, the |X|_pq = (-1)^(p+q) * det(X)/det(X^{pq})`):
print(`where X is the generic n by n matrix for all p,q, Try:`):
print(`Bdok(3);`):

elif nops([args])=1 and op(1,[args])=BU then
print(`BU(M,C): Inputs a square-matrix M of size 2 by 2, say, and a composition of n, C, outputs the correponding block matrix. Try:`):
print(`BU([seq([seq(a[i,j],j=1..4)],i=1..4)],[2,2]);`):

elif nops([args])=1 and op(1,[args])=CheckCQD1 then
print(`CheckCQD1(n,p,q): checks that if all the entries commute, the |X|_pq = (-1)^(p+q) * det(X)/det(X^{pq})`):
print(`where X is the generic n by n matrix`):
print(`Try: `):
print(`CheckCQD1(4,1,1); `):


elif nops([args])=1 and op(1,[args])=CheckHered then
print(`CheckHered(Mat,r): inputs a square matrix Mat, an integer r that is a divisor of nops(M)`):
print(`For each i,j between 1 and nops(Mat) computes the (i,j) quasi-determinant of M in two ways. First directly, and then`):
print(`via the quasi-determinant of the block matrix obtained by breaking it into r^2 parts. Checks that are all the same. Try:`):
print(`CheckHered([[a11,a12,a13,a14],[a21,a22,a23,a24],[a31,a32,a33,a34],[a41,a42,a43,a44]],2);`):
print(`CheckHered(RandMat(9,10^10),3);`):

elif nops([args])=1 and op(1,[args])=CheckHered1 then
print(`CheckHered1(Mat,r,i,j): inputs a square matrix Mat, an integer r that is a divisor of nops(M), and pos. integers i and j`):
print(`between 1 and nops(M), computes the (i,j) quasi-determinant of M in two ways. First directly, and then`):
print(`via the quasi-determinant of the block matrix obtained by breaking it into r^2 parts. Try:`):
print(`CheckHered1([[a11,a12,a13,a14],[a21,a22,a23,a24],[a31,a32,a33,a34],[a41,a42,a43,a44]],2,2,3);`):
print(`CheckHered1(RandMat(9,10^10),3,6,7);`):

elif nops([args])=1 and op(1,[args])=EvalBM then
print(`EvalBM(L,x,M,R,B): inputs an expression in non-coomutative variables of the form x[i,j], symbols M and R for`):
print(`multiplication and inverse respectively, and a block matrix B, given as a list of lists,`):
print(`evaluates it, where x[i,j] is replaced by B[i][j]`):
print(`Try:`):
print(`EvalBM(qd(x,2,M,R,1,1),x,M,R, BU([seq([seq(a[i,j],j=1..4)],i=1..4)],[2,2]));`):

elif nops([args])=1 and op(1,[args])=EvalCom then
print(`EvalCom(L,M,R): inputs an expression in non-coomutative variables returns the rational function, in normal form if all the variables commute.`):
print(`Try: `):
print(`EvalCom(qd(x,3,M,R,1,1),M,R);`):

elif nops([args])=1 and op(1,[args])=inverse1 then
print(`inverse1(A): the inveser of a square matrix given as a list of lists. Try: `):
print(`inverse1([[1,2],[4,5]]);`):

elif nops([args])=1 and op(1,[args])=multiply1 then
print(`multiply1(A,B): the matrix product of A and B given as list of lists. Try: `):
print(`multiply1([[1,2,3],[4,5,6]],[[1,2],[2,3],[4,5]]);`):

elif nops([args])=1 and op(1,[args])=qd then
print(`qd(x,n,M,R,i,j): inputs a symbol x, a positive integer n, symbols M and R, and positive integers i, j between 1 and n`):
print(`and outputs the (i,j)-quasi-determinant of the generic n by n matrix (x[a,b]) (1<=a,b<=n)  (with NON-COMMUTATIVE entries)`):
print(`It also inputs symbols`):
print(`M(a,b) means a TIMES b (in that order) and R(a) means a^(-1). Try:`):
print(` qd(x,3,M,R,1,1); `):

elif nops([args])=1 and op(1,[args])=RandMat then
print(`RandMat(n,K): a random n by n matrix with entries between -K and K. Try:`):
print(`RandMat(10,100000);`):

else
print(`There is no ezra for`,args):
fi:
 
end:

ez:=proc(): print(` , BU(M,C)`): end:


with(linalg):

#qd(x,n,M,R,i,j): inputs a symbol x, a positive integer n, symbols M and R, and positive integers i, j between 1 and n
#and outputs the (i,j)-quasi-determinant of the generic n by n matrix (x[a,b]) (1<=a,b<=n)  (with NON-COMMUTATIVE entries)
#It also inputs symbols
#M(a,b) means a TIMES b (in that order) and R(a) means a^(-1). Try:
#qd(x,n,1,1,M,R,1,1);
qd:=proc(x,n,M,R,i,j) local gu,lu,i1,j1,i1a,j1a,a,b,y:
option remember:

if not (type(n,integer) and type(i,integer) and type(j,integer) and 1<=i and i<=n  and 1<=j and j<=n ) then
 print(`bad input`):
 RETURN(FAIL):
fi:


if n=1  then
  RETURN(x[1,1]):
fi:

gu:=x[i,j]:

for i1 from 1 to n do
for j1 from 1 to n do

  if i1<>i and j1<>j then

  if i1<i then
    i1a:=i1:
   else
      i1a:=i1-1:
 fi:

  if j1<j then
    j1a:=j1:
   else
      j1a:=j1-1:
 fi:

 lu:=qd(x,n-1,M,R,i1a,j1a):

  for a from 1 to i-1 do
   
   for b from j to n-1 do
    lu:=subs(x[a,b]=y[a,b+1],lu):
   od:

  od:

  for a from i to n-1 do
   
   for b from 1 to j-1 do
    lu:=subs(x[a,b]=y[a+1,b],lu):
   od: 
 
   for b from j to n-1 do
    lu:=subs(x[a,b]=y[a+1,b+1],lu):
   od:

  od:

lu:=subs({seq(seq(y[a,b]=x[a,b],a=1..n),b=1..n)},lu):

gu:=gu -M(x[i,j1],R(lu),x[i1,j]):

fi:

od:
od:
 
gu:
    

        

end:


#EvalCom(L,M,R): inputs an expression in non-coomutative variables returns the rational function, in normal form if all the variables commute.
#Try:
#EvalCom(qd(x,3,1,1),M,R);
EvalCom:=proc(L,M,R) local i,mu,co:
option remember:

if type(L,indexed) then
 RETURN(L):
fi:


if type(L,`+`) then
 RETURN(normal(add(EvalCom(op(i,L),M,R),i=1..nops(L)))):
fi:


if type(L, `*`) then
 co:=op(1,L):
  if not type(co,numeric) then
    RETURN(FAIL):
  else
  mu:=op(2,L):

  RETURN(co*EvalCom(mu,M,R)):
fi:

fi:

if type(L,function) then

  if op(0,L)=R then
     RETURN(normal(1/ EvalCom(op(1,L),M,R))):
 elif op(0,L)=M then
   RETURN(normal(mul(EvalCom(op(i,L),M,R),i=1..nops(L)))):
 else
   RETURN(FAIL):
 fi:

fi:
 
FAIL:


end:

#CheckCQD1(n,p,q): checks that if all the entries commute, the |X|_pq = (-1)^(p+q) * det(X)/det(X^{pq})
#where X is the generic n by n matrix
CheckCQD1:=proc(n,p,q) local M,R,i,j,x,gu,mu,X,Xpq:

if not ( type(n,integer) and n>=2 and p>=1 and p<=n and q>=1 and q<=n ) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=EvalCom(qd(x,n,M,R,p,q),M,R):

X:=[seq([seq(x[i,j],j=1..n)],i=1..n)]:

Xpq:=[
seq(
[seq(x[i,j],j=1..q-1),seq(x[i,j],j=q+1..n) ],i=1..p-1), 
 seq([seq(x[i,j],j=1..q-1),seq(x[i,j],j=q+1..n) ] ,i=p+1..n)
]:
mu:=(-1)^(p+q)*det(X)/det(Xpq):

evalb(normal(gu-mu)=0);

end:


#Bdok(n): checks that if all the entries commute, the |X|_pq = (-1)^(p+q) * det(X)/det(X^{pq})
#where X is the generic n by n matrix for all p,q
Bdok:=proc(n) local p,q:
evalb({seq(seq(CheckCQD1(n,p,q) ,p=1..n),q=1..n)}={true}):
end:

#BU(M,C): Inputs a square-matrix M of size 2 by 2, say, and a composition of n, C, outputs the correponding block matrix. Try:
#BU(seq([seq(a[i,j],j=1..4)],i=1..4)],[2,2]);
BU:=proc(M,C) local n,i1,j1,r,s,T,iTop,iBot,jLeft,jRight:

if not (type(M,list) and {seq(type(M[i1],list),i1=1..nops(M))}={true} and   {seq(nops(M[i1]),i1=1..nops(M))}={nops(M)}) then
print(`Bad input`):
RETURN(FAIL):
fi:

n:=nops(M):

if not (type(C,list) and convert(C,`+`)=n) then
print(`Bad input`):
RETURN(FAIL):
fi:

for r from 1 to nops(C) do
for s from 1 to nops(C) do
  iTop:=1+add(C[i1],i1=1..r-1):
  iBot:=add(C[i1],i1=1..r):
  jLeft:=1+add(C[i1],i1=1..s-1):
  jRight:=add(C[i1],i1=1..s):
 T[r,s]:= [seq([seq(M[i1][j1],j1=jLeft..jRight)],i1=iTop..iBot)]:
od:
od:

[seq([seq(T[r,s],s=1..nops(C))],r=1..nops(C))]:


end:


#multiply1(A,B): the matrix product of A and B given as list of lists. Try
#multiply1([[1,2,3],[4,5,6]],[[1,2],[2,3],[4,5]]);
multiply1:=proc(A,B) local i1,i,j,T:
if not (type(A,list) and type(B,list) and {seq(type(A[i1],list),i1=1..nops(A))}={true} and {seq(type(B[i1],list),i1=1..nops(B))}={true}
  and nops(A[1])=nops(B) and nops({seq(nops(A[i1]),i1=1..nops(A))})=1  and nops({seq(nops(B[i1]),i1=1..nops(B))})=1  )
then
 RETURN(FAIL):
fi:

T:=multiply(matrix(A),matrix(B)):

normal([seq([seq(T[i,j],j=1..nops(B[1]))],i=1..nops(A))]):
end:

#inverse1(A): the inveser of a square matrix given as a list of lists. Try
#inverse1([[1,2],[4,5]]);
inverse1:=proc(A) local gu,i1,i,j:
if normal(det(A))=0 then
  RETURN(FAIL):
fi:

if not (type(A,list)  and {seq(type(A[i1],list),i1=1..nops(A))}={true} 
  and nops({seq(nops(A[i1]),i1=1..nops(A))})=1  and nops(A)=nops(A[1]) ) then
 RETURN(FAIL):
fi:

  gu:=inverse(convert(A,matrix)):
   normal([seq([seq(gu[i,j],j=1..nops(A))],i=1..nops(A))]):


end:


#EvalBM(L,x,M,R,B): inputs an expression in non-coomutative variables of the form x[i,j], symbols M and R for
#multiplication and inverse respectively, and a block matrix B, given as a list of lists,
#evaluates it, where x[i,j] is replaced by B[i][j]
#Try:
#EvalBM(qd(x,2,1,1),x,M,R,BU(seq([seq(a[i,j],j=1..4)],i=1..4)],[2,2]));
EvalBM:=proc(L,x,M,R,B) local lu,i,mu,co,gu:
option remember:

if type(L,indexed) then
 lu:=[op(L)]:
 if not (type(lu[1],integer) and type(lu[2],integer) and lu[1]<=nops(B) and lu[2]<=nops(B[1]) ) then
  RETURN(FAIL):
 else
 RETURN(B[lu[1]][lu[2]]):
 fi:
fi:


if type(L,`+`) then
 RETURN(normal(add(EvalBM(op(i,L),x,M,R,B),i=1..nops(L)))):
fi:


if type(L, `*`) then
 co:=op(1,L):
  if not type(co,numeric) then
    RETURN(FAIL):
  else
  mu:=op(2,L):

  RETURN(co*EvalBM(mu,x,M,R,B)):
fi:

fi:

if type(L,function) then

  if op(0,L)=R then
  if nops(R)<>1 then
  RETURN(FAIL):
 else
 RETURN(inverse1(EvalBM(op(1,L),x,M,R,B))):
 fi:
 elif op(0,L)=M then

 gu:=EvalBM(op(1,L),x,M,R,B):


 for i from 2 to nops(L) do
  gu:=multiply1(gu,EvalBM(op(i,L),x,M,R,B)):
 od:
   RETURN(gu):
 else
   RETURN(FAIL):
 fi:

fi:
 
FAIL:


end:


#RandMat(n,K): a random n by n matrix with entries between -K and K. Try:
#RandMat(10,100000);
RandMat:=proc(n,K) local ra,i,j:
ra:=rand(-K..K):
[seq([seq(ra(),i=1..n)],j=1..n)]:
end:

#CheckHered1(Mat,r,i,j): inputs a square matrix M, an integer r that is a divisor of nops(M), and pos. integers i and j
#between 1 and nops(M), computes the (i,j) quasi-determinant of M in two ways. First directly, and then
#via the quasi-determinant of the block matrix obtained by breaking it into r^2 parts. Try:
#CheckHered1([[a11,a12,a13,a14],[a21,a22,a23,a24],[a31,a32,a33,a34],[a41,a42,a43,a44]],2,2,3);
CheckHered1:=proc(Mat,r,i,j) local M,R,x,gu,i1,j1,Mat1,q,mu,iNew,jNew,gu2,iNew0,jNew0,gu1,mu1,mu2:

if not (type(Mat,list)  and {seq(type(Mat[i1],list),i1=1..nops(Mat))}={true} 
  and nops({seq(nops(Mat[i1]),i1=1..nops(Mat))})=1  and nops(Mat)=nops(Mat[1]) ) then
 RETURN(FAIL):
fi:

if not type(nops(Mat)/r,integer) then
 RETURN(FAIL):
fi:

gu:=qd(x,nops(Mat),M,R,i,j):

mu:=EvalCom(gu,M,R) :

mu:=subs({seq(seq(x[i1,j1]=Mat[i1][j1],j1=1..nops(Mat)),i1=1..nops(Mat))},mu):


q:=nops(Mat)/r:
Mat1:=BU(Mat,[q$r]):

iNew:=trunc((i-1)/q)+1:
jNew:=trunc((j-1)/q)+1:

gu1:=qd(x,r,M,R,iNew,jNew):
mu1:=EvalBM(gu1,x,M,R,Mat1):


iNew0:=i mod q:

if iNew0=0 then
 iNew0:=q:
fi:

jNew0:=j mod q:

if jNew0=0 then
 jNew0:=q:
fi:

gu2:=EvalCom(qd(x,q,M,R,iNew0,jNew0),M,R):
mu2:=subs({seq(seq(x[i1,j1]=mu1[i1][j1],j1=1..q),i1=1..q)},gu2):

evalb(normal(mu-mu2)=0):

end:

#CheckHered(Mat,r): Does CheckHered1(Mat,r,i,j) for all i and j. Try:
#CheckHered1([[a11,a12,a13,a14],[a21,a22,a23,a24],[a31,a32,a33,a34],[a41,a42,a43,a44]],2);
CheckHered:=proc(Mat,r) local i,j:
evalb({seq(seq(CheckHered1(Mat,r,i,j),i=1..nops(Mat)),j=1..nops(Mat))}={true}):
end: