##------------------------------------------------------------------------
##    Title: 	Mubasis
##    Created:	March 2010
##    Authors: 	Luu Ba Thang <lbathang@sophia.inria.fr>
##    Advisor:  Laurent Buse <Laurent.Buse@inria.fr>
##    
##    Description: Matrix of Koszul map, mu-basis of set univariate polynomial, mu-basis of set 
##                 homogeneous polynomial.
##----------------------------------------------------------------------                 
##    Version last modified 03/2010
## ----------------------------------------------------------------------
##	   Mu-basis: Basic functions to work with polynomials are developped. There is functions to
##         get basis of syzygies of set of univariates polynomials or set of homogeneous polynomials
##         of the same degree. This funtions allow our to compute representation matrices of parameterized 
##         plane curves and space curves.

## NB:The "warnings" obtained when loading representatiomatrices are due to the Maple 
## packages LinearAlgebra,MatrixPolynomialAlgebra,PolynomialTools and Groebner which are used here.

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


with(LinearAlgebra):
with(MatrixPolynomialAlgebra): 
with(PolynomialTools):
with(Groebner):

#----------------------------------------------------------------------
# addtohelp(SetPol):
#----------------------------------------------------------------------
## HELP SetPol
## CALLING SEQUENCE:
##     SetPol(f)
## DESCRIPTION:
## - f - set of polynomials with |f|=n;
## Return: Set of generators of module syzygies defining by f (Each generator is given under form matrix(1,n))
## EXAMPLES:
## > SetPol(<t, t*t+1, -1>);
## Matrix(1, 3, {(1, 1) = -t^2-1, (1, 2) = t, (1, 3) = 0}), 
## Matrix(1, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = t}), 
## Matrix(1, 3, {(1, 1) = 0, (1, 2) = 1, (1, 3) = t^2+1})

## SetPol(<t, t*t+1, 2*t, t+1>)
## Matrix(1, 4, {(1, 1) = -t^2-1, (1, 2) = t, (1, 3) = 0, (1, 4) = 0}), 
## Matrix(1, 4, {(1, 1) = -2*t, (1, 2) = 0, (1, 3) = t, (1, 4) = 0}), 
## Matrix(1, 4, {(1, 1) = -t-1, (1, 2) = 0, (1, 3) = 0, (1, 4) = t}), 
## Matrix(1, 4, {(1, 1) = 0, (1, 2) = -2*t, (1, 3) = t^2+1, (1, 4) = 0}), 
## Matrix(1, 4, {(1, 1) = 0, (1, 2) = -t-1, (1, 3) = 0, (1, 4) = t^2+1}), 
## Matrix(1, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = -t-1, (1, 4) = 2*t})   
                      
SetPol := proc (f) 
  local i, j, k, Z, C, n; 
  n := Dimension(f); 
  C := []; 
  Z := Matrix(1, n); 
  for i from 1 to n do 
      for j from i+1 to n do 
      Z := Matrix(1, n); 
      Z[1, i] := -f[j]; 
      Z[1, j] := f[i]; 
      C := [op(C), Z] end do; 
  end do; 
RETURN(C); 
end proc:
#----------------------------------------------------------------------
# addtohelp(LeadVector):
#----------------------------------------------------------------------
# DESCRIPTION:
# Return a leading coefficient matrix of univariate polynomial matrix of size (1,n).
LeadVector := proc (u, var)
  local m, n, Z, i, V; 
  n := ColumnDimension(u); 
  m := Degree(u, var); 
  Z := Matrix(1, n); 
  for i to n do 
    if degree(u[1, i], var) = m then V := CoefficientVector(expand(u[1, i]), var); 
                                Z[1, i] := V[m+1] 
     end if 
  end do; 
RETURN(Z);
end proc:

# Sort the matrix of leading vector and degree of each  polynomial matrix with arrangement
#----------------------------------------------------------------------
# addtohelp(Step12):
#----------------------------------------------------------------------
Step12 := proc (X, var) 
  local m, n, A, i, D, j, B;
 n := ColumnDimension(X[1]); 
 m := nops(X); A := Matrix(m, n); 
 D := Matrix(1, m); 
 for i to m do 
   for j to n do 
         B := LeadVector(X[i], var); 
         A[i, j] := B[1, j] 
   end do; 
   D[1, i] := Degree(X[i], var) 
  end do; RETURN(A, D); 
end proc:


# Arrange set of generator by increasing order
#----------------------------------------------------------------------
# addtohelp(Arrange):
#----------------------------------------------------------------------
Arrange := proc (u) 
   local i, j, s, n, Z; 
   Z := u; 
   n := ColumnDimension(Z); 
   for i to n do 
     for j from i+1 to n do 
      if Z[1, i] < Z[1, j] then
        s := Z[1, i]; 
        Z[1, i] := Z[1, j]; 
        Z[1, j] := s 
       end if 
     end do 
    end do; 
RETURN(Z); 
end proc:

#----------------------------------------------------------------------
# addtohelp(Step123):
#----------------------------------------------------------------------

Step123 := proc (T,var) 
local A, Z, Y, i, j, s, n, t, p, q, k, l,X; 
     X := T; 
  Y, A := Step12(T, var);
     p := RowDimension(Y); 
     q := ColumnDimension(Y); 
     Z := A; n := ColumnDimension(Z); 
    for i to n do 
    for j from i+1 to n do 
     if Z[1, i] < Z[1, j] then 
       s := Z[1, i]; 
       Z[1, i] := Z[1, j]; 
       Z[1, j] := s; 
     for k to q do 
     l := Y[i, k]; Y[i, k] := Y[j, k]; Y[j, k] := l 
     end do; 
     for k to q do 
     t := X[i][1, k]; 
    X[i][1, k] := X[j][1, k]; 
    X[j][1, k] := t 
    end do 
     end if 
     end do 
    end do; 
RETURN(X, Y, Z);
end proc:
#----------------------------------------------------------------------
# addtohelp(CompareMatrix):
#----------------------------------------------------------------------
CompareMatrix := proc (U, V) 
   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(Step1234):
#----------------------------------------------------------------------
Step1234 := proc (A) 
   local m, p, i, j, B, K, a, U, k, n, l, C; 
   m := ColumnDimension(A); 
   n := RowDimension(A); 
   C := Transpose(A); 
   p := n-Rank(C); 
   B := NullSpace(C); 
   K := Matrix(1, m); 
   for i to p do a := 0; 
   for j to n do U := Matrix(1, m); 
   for k to m do U[1, k] := A[j, k] end do; 
   if CompareMatrix(B[i][j]*U, K) = 0 then a := a+1 
   end if end do; 
   if 2 <= a then for l to n do if B[i][l] <> 0 then RETURN(B[i], l) 
   end if 
   end do 
   end if 
   end do 
end proc:
#----------------------------------------------------------------------
# addtohelp(Reduction):
#----------------------------------------------------------------------
Reduction := proc (T,var) 
  local A, B, C, H, X, E, G, m, s, j, p, k,Z; 
      X := T;
    (C,A,B) := Step123(X,var); 
       (E,m) := Step1234(A); 
          s := RowDimension(A); 
          k := ColumnDimension(A); 
          Z := Matrix(1, k); 
          G := Matrix(1, k); 
     for j from m to s do p := B[1, m]-B[1, j]; H := expand(E[j]*C[j]*var^p); G := G+H end do; 
     for j to k do G[1, j] := expand(G[1, j]) end do; 
     if CompareMatrix(G, Z) = 0 then X[m] := G else X := subsop(m = NULL, X) end if; 
    RETURN(X); end proc:

#----------------------------------------------------------------------
addtohelp(Mubase):
#----------------------------------------------------------------------

## HELP Mubase
## CALLING SEQUENCE:
##     Mubase(f,var,varimpl)
## DESCRIPTION:
## - f - set of univariate polynomials, f:=[f0,f1,...,fn] and GCD(f0,f1,...,fn)=1. 
## - var -: a variable
## - varimpl - Set of variable implicit.
## Return: A basis of module syzygies defining by f with decreasing degree of two form:
## List of polynomials and List of matrices..
## EXAMPLES:
## Mubase([t, t*t+1, -1], t,[x,y,z])
## Vector(2, {(1) = t*x-y-z, (2) = x+t*z})
##      Matrix(1, 3, {(1, 1) = t, (1, 2) = -1, (1, 3) = -1}), 
##      Matrix(1, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = t})

Mubase := proc (f, var,varimpl) 
  local l, k,X, U,V,i;
     l := nops(f); 
     X := SetPol(Vector(f)); 
     U := X; 
     k := nops(U); 
     while (l-1) < k do 
         U := Reduction(U,var); 
         k := nops(U) 
     end do; 
     V:=Vector(l-1);
     for i from 1 to (l-1) do
         V[i]:= Multiply(U[i], Vector(varimpl))[1];
         end do; 
RETURN(V,U); 
end proc:

#----------------------------------------------------------------------
addtohelp(MubaseHomogeneous):
#----------------------------------------------------------------------
## HELP MubaseHomogeneous
## CALLING SEQUENCE:
##     MubaseHomogeneous(G,var,varimpl)
## DESCRIPTION:
## - G - List of homogeneous polynomials of var1, var2, G:=[F0,F1,...,Fn] and GCD(F0,F1,...,Fn)=1.
## - varimpl - List of parameterized; 
## - var -: Two variables .
## RETURN: A basis of module syzygies defining by f under the form polynomial in K[var;varimpl] with decreasing degrees 
##         and a basis of module syzygies defining by f under the form matrix list.
## EXAMPLES: MubaseHomogeneous([s^2*t, (t*t)*s+(s*s)*s, -s^3, s*t^2], [s, t], [x, y, z, w])
##   Vector(3, {(1) = t*x-s*y-s*z, (2) = s*x+t*z, (3) = -y-z+w});
##   Matrix(1, 4, {(1, 1) = t, (1, 2) = -s, (1, 3) = -s, (1, 4) = 0}), 
##   Matrix(1, 4, {(1, 1) = s, (1, 2) = 0, (1, 3) = t, (1, 4) = 0}), 
##   Matrix(1, 4, {(1, 1) = 0, (1, 2) = -1, (1, 3) = -1, (1, 4) = 1})

MubaseHomogeneous := proc (G, var,varimpl) 
    local f, U1,U2, V, i, m, n,j,s,Z; 
     f := subs({var[1] = 1}, G); 
     (U1,U2) := Mubase(f, var[2],varimpl); 
     n := nops(G)-1; 
     V := Vector(n); 
     for i to n do V[i] := U1[i]; end do;
     for i to n do V[i] := Homogenize(V[i], var[1]); 
       end do; 
     for i from 1 to n do 
           for j from 1 to (n+1) do U2[i][1,j]:=coeff(V[i],varimpl[j]);
           end do;
      end do;
RETURN(V,U2); 
end proc:




[__];