##    Title: 	WDecomposition
##    Created:	April 2011
##    Authors: 	Luu Ba Thang <lbathang@sophia.inria.fr>
##    
## ----------------------------------------------------------------------
## Description: This part consists some functions to compute the W-Decomposition of polynomial
## matrix, the regular part of a bivariate pencil of matrices and the surface/surface intersection. 

##-------------------------------------------------------------------------
##   Version last modified 05/2011
##-------------------------------------------------------------------------

__ := NULL: 
addtohelp := proc() 
  global __; 
  __ := __, args: 
end:

with(linalg):
with(LinearAlgebra):
with(MatrixPolynomialAlgebra): 
#with(PolynomialTools):
#-----------------------------------------------------------------------------
addtohelp(Change):
addtohelp(ChangeColumn):
addtohelp(ChangeRow):
#------------------------------------------------------------------------------
Change:=proc(m,n,k) # a square matrix of size kxk obtained by permuting the columns number m with
# column number n.
local A,i;
A:=Matrix(k);
for i from 1 to k do
     A[i,i]:=1; end do;
A[m,m]:=0;
A[n,n]:=0;
A[m,n]:=1;
A[n,m]:=1;
RETURN(A);
end proc:

ChangeColumn:=proc(A,m,n) # permuting the column number m with column number n of a matrix A.
local k,B;
k:=ColumnDimension(A);
if m<=k and n<=k then
B:=Multiply(A,Change(m,n,k));
else RETURN("Do not change");
end if;
RETURN(B);
end proc:

ChangeRow:=proc(A,m,n) #permuting the row number m with row number n of a matrix A.
local k,B;
k:=RowDimension(A);
if m<=k and n<=k then
B:=Multiply(Change(m,n,k),A);
else RETURN("Do not change");
end if;
RETURN(B);
end proc:
#----------------------------------------------------------------------
addtohelp(CompareMatrix):
#----------------------------------------------------------------------
CompareMatrix := proc (U, V) # U<>V iff b=0.
   local a, n, b, i; 
     a := 0; n := ColumnDimension(V); 
    for i to n do 
      if V[1, i] <> U[1, i] then a := a+1
      end if end do; 
    if 1 <= a then b := 0 else b := 1 end if; 
RETURN(b);
end proc:
#----------------------------------------------------------------------
addtohelp(CompareColumn):
#----------------------------------------------------------------------
CompareColumn:=proc(M,N) # M<>N iff b:=1.
local a,n,b,i;
a:=0;
n:=Dimension(M);
for i from 1 to n do
if M[i]<>N[i] then a:=a+1; fi:
od:
if 1<=a then b:=1; else b:=0; fi:
RETURN(b);
end proc:
#----------------------------------------------------------------------
addtohelp(CompareRow):
#----------------------------------------------------------------------
CompareRow:=proc(M,N)  # M<>N iff b:=1.
local a,n,b,i;
a:=0;
n:=Dimension(M);
for i from 1 to n do
if M[i]<>N[i] then a:=a+1; fi:
od:
if 1<=a then b:=1; else b:=0; fi:
RETURN(b);
end proc:

#----------------------------------------------------------------------
addtohelp(PermuteMatrix):
#----------------------------------------------------------------------
## HELP PermuteMatrix 
## CALLING SEQUENCE:
##     PermuteMatrix(M)
## DESCRIPTION:
## - Given a matrix M. 
## RETURN: Matrix N, P and the number s where s is the number nonzero columns of M,
## N is a matrix that its s first columns are the nonzero columns of M and the rest columns are
## zero columns. M*P=N. 

PermuteMatrix:=proc(A)
local m,n,i,j,M,V,k,P,q;
M:=A;
n:=ColumnDimension(M);
m:=RowDimension(M); 
V:=Vector(m,0); 
P:=IdentityMatrix(n);
for i from 1 to n do
 for k from i to n do
   if CompareColumn(Column(M,k),V)=0 then 
      for j from 0 to n-1 do
       if CompareColumn(Column(M,n-j),V)=1 and k<=(n-j) then M:=ChangeColumn(M,k,n-j); P:=Multiply(P,Change(k,n-j,n)); break; fi: 
     od:
   fi: break;
od:
od:
q:=0;
for i from 1 to n do
   if CompareColumn(Column(M,i),V)=0 then q:=i-1;break;
else q:=n; fi: od:
RETURN(M,P,q);
end proc:

#----------------------------------------------------------------------
addtohelp(JoinDiagonalMatrix):
#----------------------------------------------------------------------
JoinDiagonalMatrix:=proc(A,B)
local M,N,K;
M:=Matrix(RowDimension(B), ColumnDimension(A),0);
N:=Matrix(RowDimension(A), ColumnDimension(B),0);
K:=<<A,M>|<N,B>>;
end proc:
#----------------------------------------------------------------------
addtohelp(DeleteCoeff):
#----------------------------------------------------------------------
# Delete the highest degree coefficients of a polynomial matrix.
DeleteCoeff:=proc(M,var)
local A,B,m,i;
A:=M;
m:=Degree(A,var);
B:=Coeff(A,var,0); 
if m=0 then RETURN([]); fi:
if m=1 then RETURN(B); fi:
if m>1 then 
for i from 1 to m-1 do
B:=B+Coeff(A,var,i)*var^i; 
od: fi:
RETURN(B);
end proc:

#......................................................................
addtohelp(MatrixN):
#......................................................................
MatrixN := proc(M,var) # Construct an auxiliary matrix N, permuted matrix and List index column.
local U,m,n,i,N,V,A,B,U1,B1,C,D,k,P,Lis,L;
A:=M;
V:=Vector(RowDimension(A),0);
n:=ColumnDimension(A);
L:=[];
A:=PermuteMatrix(A)[1];
for i from 1 to n do
   if CompareColumn(Column(A,i),V)=0 then L:=[op(L),i];
fi: od: 
if nops(L)>0 then A:=DeleteColumn(A,L[1]..L[nops(L)]); fi:
m:=Degree(A,var);
#print(A);
N:=Matrix(V); 
U:=PermuteMatrix(Coeff(A,var,m)); 
#print(U);
if U[3]=n then N:=U[1]; fi: 
if U[3]<n then  
     for i from 1 to U[3] do N:=<N|Column(U[1],i)>; 
od: N:=DeleteColumn(N,1);fi: 
#print(N);
k:=U[3];
P:=U[2];
Lis:=[k];
while ColumnDimension(N)<n do
U:=PermuteMatrix(Coeff(A,var,Degree(A,var)));
B1:=Multiply(A,U[2]);
C:=DeleteCoeff(B1,var); 
if C=[] then RETURN(N,P,Lis); fi:
B:=DeleteColumn(Coeff(C,var,Degree(C,var)),1..k); 
#print(B);
U1:=PermuteMatrix(B); 
if U1[3]=0 then Lis:=[op(Lis),0];fi:
if U1[3]>0 and ColumnDimension(N)+U1[3]<=n then Lis:=[op(Lis),U1[3]];
for i from 1 to U1[3] do N:=<N|Column(U1[1],i)>;od: 
else if U1[3]>0 and ColumnDimension(N)+U1[3]>n then Lis:=[op(Lis),n-ColumnDimension(N)];
for i from 1 to n-ColumnDimension(N) do
N:=<N|Column(U1[1],i)>;od: 
fi: 
fi:
D:=JoinDiagonalMatrix(IdentityMatrix(k),U1[2]);
A:=Multiply(C,D); 
#print(A);
P:=Multiply(P,D); 
k:=ColumnDimension(N); 
end do:
RETURN(N,P,Lis);
end proc:

#----------------------------------------------------------------------
addtohelp(AddColumn):
#----------------------------------------------------------------------
AddColumn:=proc(M,k1,k2,u1,v1) #` Replace the column  k1 of M with  u1*column k1 + v1* column k2`  
local k;  
for k from 1 to RowDimension(M) do  
M[k,k1]:=u1*M[k,k1]+v1*M[k,k2]; 
od:  
RETURN(M);
end proc:
#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
GaussMatrix:=proc(M) #
local N,i,b;
N:=M;
b:=0;
for i from 1 to ColumnDimension(N) do
if N[1,i]<>0 then b:=i; break; fi: od:
if b>0 then 
N:=ChangeColumn(N,1,b); 
for i from 2 to ColumnDimension(N) do N:=AddColumn(N,i,1,N[1,1],-N[1,i]) od:
fi:
RETURN(N,b);
end proc:

#-------------------------------------------------------------------------
#-------------------------------------------------------------------------

ColumnZero:=proc(M)
local A,V,b,i;
A:=M;
V:=Vector(RowDimension(A),0);
b:=0;
for i from 1 to ColumnDimension(A) do
if CompareColumn(Column(A,i),V)=0 then b:=b+1; fi: od:
RETURN(b);
end proc:

#-------------------------------------------------------------------------
#-------------------------------------------------------------------------

AddIdenColumn:=proc(M,i)
local A,m,V;
A:=M;
m:=RowDimension(A);
V:=Vector[row](ColumnDimension(A),0); if CompareRow(Row(A,i),V)=0 then 
if i=1 then A:=<Matrix(Column(IdentityMatrix(m),i))|SubMatrix(A,1..m,i..ColumnDimension(A))>; fi:
if i>1 then A:=<SubMatrix(A,1..m,1..(i-1))|Matrix(Column(IdentityMatrix(m),i))|SubMatrix(A,1..m,i..ColumnDimension(A))>; 
fi:
fi:
RETURN(A);
end proc:
#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
WReduction := proc (M) 
local N, A, B, C, B1, C1, D, D1, b, i; 
A := M; 
N := GaussMatrix(A); 
if N[2] <> 0 then A := N[1] 
else A := AddIdenColumn(A, 1) end if; 
b := 2; 
while b <= RowDimension(A) do 
B := SubMatrix(A, b .. RowDimension(A), b .. ColumnDimension(A)); 
C := GaussMatrix(B); 
if C[2] <> 0 then B := C[1] else B := AddIdenColumn(B, 1) end if; 
B := `<,>`(Matrix(b-1, ColumnDimension(B)), B); 
A := `<|>`(SubMatrix(A, 1 .. RowDimension(A), 1 .. b-1), B); 
b := b+1 
end do; 
for i to ColumnDimension(A) do 
A := AddColumn(A, i, 1, 1/A[i, i], 0) end do; 
RETURN(A); end proc:
#-------------------------------------------------------------------------
addtohelp(ExpandMatrix):
#-------------------------------------------------------------------------
ExpandMatrix := proc (M) 
local U, i, j; U := M; 
for i to RowDimension(U) do 
for j to ColumnDimension(U) do 
U[i, j] := expand(U[i, j]);
 end do; 
end do; 
RETURN(U) end proc:
#-------------------------------------------------------------------------
addtohelp(MatrixNullSpace):
#-------------------------------------------------------------------------
MatrixNullSpace := proc (M)  ## Output: NullSpace of a matrix under the matrix form
local U, V, i; U := NullSpace(M); 
if nops(U) = 0 then RETURN(U) end if; 
if 0 < nops(U) then V := Matrix(U[1]) end if; 
for i from 2 to nops(U) do V := `<|>`(V, Matrix(U[i])) end do; 
RETURN(V) end proc:

#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
BlockMatrix := proc (M, var, L, k) 
local A, B, m, j; 
m := nops(L); 
if L[k]=0 then
for j from k+1 to m do 
if L[j]<>0 then k:=j; break; fi: od: fi:
A := SubMatrix(M, 1 .. L[k], 1 .. ColumnDimension(M)); 
for j from k+1 to m do 
B := ExpandMatrix(var^(j-k)*SubMatrix(M, L[j-1]+1 .. L[j], 1 .. ColumnDimension(M))); 
A := `<,>`(A, B);
end do;
RETURN(A) end proc:

#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
 
WMatrix := proc (M, var, L) # M is a square matrix.
local List, k, a, i, A, B; 
k := nops(L); 
List := []; 
a := 0; 
for i to k do a := a+L[i]; List := [op(List), a] end do; 
A := BlockMatrix(SubMatrix(M, 1 .. RowDimension(M), 1 .. List[1]), var, List, 1); 
for i to k-1 do 
B := BlockMatrix(SubMatrix(M, 1 .. RowDimension(M), List[i]+1 .. List[i+1]), var, List, i+1); 
A := `<|>`(A, B); 
end do; 
RETURN(A);
end proc:

#-------------------------------------------------------------------------
addtohelp(WDecompositionStep):
#-------------------------------------------------------------------------

WDecompositionStep := proc (M, var) 
local b, A, B, B1, B2, B3, C, W; 
A := M; 
if Rank(A) = ColumnDimension(A) then 
 RETURN(A,IdentityMatrix(ColumnDimension(A)),ColumnDimension(A) ) end if; 
b := ColumnZero(A); 
W := IdentityMatrix(ColumnDimension(A)); 
while b = 0 do 
B := MatrixN(A, var); 
B1 := MatrixNullSpace(B[1]); 
B2 := WReduction(B1); 
B3 := WMatrix(B2, var, B[3]); 
A := ExpandMatrix(Multiply(Multiply(A, B[2]), B3)); 
W := ExpandMatrix(Multiply(W, Multiply(B[2], B3))); 
b := ColumnZero(A) 
end do; 
C := PermuteMatrix(A); 
W := ExpandMatrix(Multiply(W, C[2])); 
RETURN(C[1], W, C[3]);
end proc:

#-------------------------------------------------------------------------
addtohelp(WDecompositionMatrix):
#-------------------------------------------------------------------------

## HELP WDecompositionMatrix
## CALLING SEQUENCE:
##     WDecompositionMatrix(M,var);
## DESCRIPTION:
## - var: variable of polynomial.
## - M - Univariable polynomial matrix. 
## - var: variable of polynomial.
## RETURN: A polynomial matrix A under the form [B,0] where B is of full column rank and unimodular polynomial matrix W such that M*W=A.
 
WDecompositionMatrix := proc (M, var) 
local A, m, n, WD1, WD, W1, B, W;
A := M; 
m := RowDimension(A); 
n := ColumnDimension(A); 
if Rank(A) = n then RETURN(A, IdentityMatrix(n), n) end if; 
B := A; 
WD := WDecompositionStep(A, var); 
A:=WD[1];
W:=WD[2];
B := SubMatrix(WD[1], 1 .. m, 1 .. WD[3]); 
while Rank(B) < ColumnDimension(B) do 
WD1 := WDecompositionStep(B, var);
W1 := JoinDiagonalMatrix(WD1[2], IdentityMatrix(n-WD[3])); 
W := ExpandMatrix(Multiply(W, W1)); 
A := `<|>`(WD1[1], Matrix(m, n-WD[3], 0)); 
WD := WDecompositionStep(WD1[1], var); 
B := SubMatrix(WD[1], 1 .. m, 1 .. WD[3]);
end do; 
RETURN(A, W); end proc:

##################################
##################################
RegularPencilVarStep:= proc(M,N,var)
local A1,A,B,m,n,k1,k2,B0,B1,B2,B3,B4,A3,A4,B5 ;
A:=M;
B:=N;
#if Rank(A)=ColumnDimension(A) then A:=Transpose(A); B:=Transpose(B); fi:
while Rank(A)< ColumnDimension(A) do
n:=RowDimension(A);
m:=ColumnDimension(A);
A1:= WDecompositionMatrix(A,var); 
k1:=ColumnZero(A1); 
B0:= ExpandMatrix(Multiply(B,A1[2])); 
B1:= SubMatrix(B0,1..n,m-k1+1..m); 
B5:=Transpose(B1);
B2:=WDecompositionMatrix(B5,var); 
B3:= ExpandMatrix(Multiply(Transpose(B2[2]),B0)); 
k2:=ColumnZero(B2[1]); 
A3:= ExpandMatrix(Multiply(Transpose(B2[2]),A1[1])); 
A4:=SubMatrix(A3,n-k2+1..n,1..m-k1);
B4:=SubMatrix(B3,n-k2+1..n,1..m-k1);
A:=A4;
B:=B4; 
end do:
RETURN(A,B);
end proc:



#...........................................................................
addtohelp(RegularPencilVar):
#...........................................................................
## HELP RegularPencil
## CALLING SEQUENCE:
##     RegularPencilVar(M,N,var);
## DESCRIPTION:
## - M,N - Univariable polynomial matrix of the same size.  
## var: variable of polynomial.
## - var - variable of M, N.
## RETURN: A pair square polynomial matrix A,B which are regular part of M,N.

RegularPencilVar:= proc(M,N,var)
local A1,A,B,m,n,k1,k2,B0,B1,B2,B3,B4,A3,A4,B5 ;
A:=M;
B:=N;
if Rank(A)=ColumnDimension(A) then A:=Transpose(A); B:=Transpose(B); fi:
while Rank(A)< ColumnDimension(A) do
n:=RowDimension(A);
m:=ColumnDimension(A);
A1:= WDecompositionMatrix(A,var); 
k1:=ColumnZero(A1); 
B0:= ExpandMatrix(Multiply(B,A1[2])); 
B1:= SubMatrix(B0,1..n,m-k1+1..m); 
if k1=ColumnDimension(A) then A,B:=RegularPencilVarStep(B1,A1[1],var); RETURN(B,A); end if;
B5:=Transpose(B1); 
B2:=WDecompositionMatrix(B5,var); 
B3:= ExpandMatrix(Multiply(Transpose(B2[2]),B0)); 
k2:=ColumnZero(B2[1]);
A3:= ExpandMatrix(Multiply(Transpose(B2[2]),A1[1])); 
A4:=SubMatrix(A3,n-k2+1..n,1..m-k1);
B4:=SubMatrix(B3,n-k2+1..n,1..m-k1);
A:=A4;
B:=B4; 
print(A3,B3);
end do:
RETURN(A,B);
end proc:

#-----------------------------------------------------------------------------
addtohelp(IntersectionSurface):
#-----------------------------------------------------------------------------
## HELP IntersectionSurface
## CALLING SEQUENCE:
##     IntersectionSurf(Mat,varmat,Surface,varsurf)
## DESCRIPTION:
## - Mat - Matrix representation of a  parametric surface S_1;  
## - varmat - set of implicit variables of Mat.
## - Surface - a parametric surface S_2.
## - varsurf - set of parametric of S_2.
## RETURN: Intersection curve in the parameter space coresponds to the intersection S_1 and S_2.

     
IntersectionSurface:= proc(M,varmat,S,varsurf) 
local U,V,A,B,A1,B1,Curve;
U:=S;
V:=M;
U:=subs({varsurf[3]=1},U); 
V:=subs({varmat[1]=U[1],varmat[2]=U[2],varmat[3]=U[3],varmat[4]=U[4]},V);
A:= FirstCompanionMatrix(V,varsurf[2]); 
B:= SecondCompanionMatrix(V,varsurf[2]);
A1,B1:= RegularPencilVar(A,B,varsurf[1]);
Curve:= Determinant(A1-varsurf[2]*B1);
RETURN(Curve);
end proc:
[__];