source: git/Singular/LIB/linalg.lib @ 0b59f5

//
//      last modified: 
//////////////////////////////////////////////////////////////////////////////

version="$Id: linalg.lib,v 1.4 1999-12-13 15:33:48 obachman Exp $";
info="
LIBRARY:  linalg.lib    PROCEDURES FOR ALGORITHMIC LINEAR ALGEBRA
AUTHOR:   Ivor Saynisch (ivs@math.tu-cottbus.de)

PROCEDURES:
 inverse(A);        matrix, the inverse of A (for constant matrizes)
 inverse_B(A);      list(matrix Inv,poly p),Inv*A=p*En ( using busadj(A) )
 inverse_L(A);      list(matrix Inv,poly p),Inv*A=p*En ( using lift )
 sym_gauss(A);      symmetric gauss algorithm
 orthogonalize(A);  Gram-Schmidt orthogonalization
 diag_test(A);      test whether A can be diagnolized
 busadj(A);         coefficients of Adj(E*t-A) and coefficients of det(E*t-A)
 charpoly(A,v);     characteristic polynomial of A. ( using busadj(A) ) 
 adj(A);            ( using busadj(A) )
 det_B(A);          determinant of A. ( uses busadj(A) )
 gaussred(A);       gaussian reduction: P*A=U*S, S is a row reduced form of A
 gaussred_pivot(A); gaussian reduction: P*A=U*S, S is a row reduced form of A
 gauss_nf(A);            gauss normal form of A
 mat_rk(A);             rank of A
 U_D_O(A);            P*A=U*D*O  
 pos_def(A,i);      test of positive symmetric matrizes.        
";

LIB "matrix.lib";
LIB "ring.lib";

//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
// help functions
//////////////////////////////////////////////////////////////////////////////

proc abs(poly c)
"RETURN: absolut value of c, c must be a constant"
{
  if(c>=0){ return(c);}
  else{ return(-c);}
}
 
proc const_mat(matrix A)
"RETURN:   1 (0) if A is (is not) a constant matrix"
{
  int i;
  int n=ncols(A);
  def BR=basering;
  changeord("@R","dp,c",BR);
  matrix A=fetch(BR,A);
  for(i=1;i<=n;i=i+1){
    if(deg(lead(A)[i])>=1){
      //"input is not a constant matrix";
      kill @R;
      setring BR;
      return(0);
    }
  }
  kill @R;
  setring BR;
  return(1);  
}

proc red(matrix A,int i,int j)
"USAGE:    red(A,i,j);  A = constant matrix
          reduces column j with respect to A[i,i] and column i
          reduces row j with respect to A[i,i] and row i
RETURN:   matrix
"
{       
  module m=module(A);
  
  if(A[i,i]==0){
    m[i]=m[i]+m[j];
    m=module(transpose(matrix(m)));
    m[i]=m[i]+m[j];
    m=module(transpose(matrix(m)));
  }
  
  A=matrix(m);
  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
  m=module(transpose(matrix(m)));
  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
  m=module(transpose(matrix(m)));
  
  return(matrix(m));
}

proc sp(vector v1,vector v2)
"RETURN:   inner product <v1,v2> "
{
  int k;
  if (nrows(v2)>nrows(v1)) { k=nrows(v2); } else { k=nrows(v1); } 
  return ((transpose(matrix(v1,k,1))*matrix(v2,k,1))[1,1]);
}


/////////////////////////////////////////////////////////////////////////////
// user functions
/////////////////////////////////////////////////////////////////////////////

proc inverse(matrix A)
"USAGE:    inverse(A);  A = constant, square matrix
RETURN:   matrix, the inverse matrix of A
EXAMPLE:  example inverse; shows an example"
{
  //erlaubt Parameterringe und minpoly

  int i;
  int n=nrows(A);
  string mp=string(minpoly);

  if (ncols(A)!=n){
    "input is not a square matrix";
    matrix invmat;
    return(invmat);
  }
  if(!const_mat(A)){
    "input is not a constant matrix";
    matrix invmat;
    return(invmat);
  }
          
  def BR=basering;
  changeord("@R","c,dp",BR);
  execute("minpoly="+mp+";");
  matrix A=fetch(BR,A);
  
  matrix B=transpose(concat(A,unitmat(n)));
  matrix D=transpose(std(B));
  matrix D1=submat(D,1..n,1..n);
  matrix D2=submat(D,1..n,(n+1)..2*n);
  B=transpose(concat(D2,D1));

  changeord("@@R","C,dp",@R);
  execute("minpoly="+mp+";");
  matrix B=imap(@R,B);
 
  for(i=1;i<=n/2;i=i+1){
    B=permcol(B,i,n-i+1);
  }
  B=std(B);
  module C=module(B);
  
  for (i=1;i<=n;i=i+1){
    if (deg(B[n+i,i])<0){
      "matrix is not invertible";
      setring BR;
      kill @R;kill @@R;
      matrix invmat;
      return(invmat);
    }
    else{ C[i] = C[i]/B[n+i,i];}
  }
  
  B=transpose(matrix(C));       
  matrix invmat[n][n]=B[1..n,1..n];
  setring BR;
  matrix invmat=fetch(@@R,invmat);
  kill @R;kill @@R;     
  return(invmat);
}
example{
  "EXAMPLE";
  ring r=0,(x),lp;
  matrix A[4][4]=5,7,3,8,4,6,9,2,3,7,5,7,2,5,6,12;
  "A=";print(A);
  "inverse(A)=";print(inverse(A));
  "A*inverse(A)=";print(inverse(A)*A)
}

//////////////////////////////////////////////////////////////////////////////
proc sym_gauss(matrix A)
"USAGE:    sym_gauss(A);  A = symmetric matrix
          diagonalisation with symmetric gauss algorithm
RETURN:   matrix
EXAMPLE:  example sym_gauss; shows an example"
{
  int i;
  int j;
  int n=nrows(A);

  if (ncols(A)!=n){
    "input is not a square matrix";;
    return(A);
  } 
  
  if(!const_mat(A)){
    "input is not a constant matrix";
    return(A);
  }
  
  if(deg(std(A-transpose(A))[1])!=-1){ 
    "input is not a symmetric matrix";
    return(A);
  }

  for(i=1;i<n;i++){
    for(j=i+1;j<=n;j++){
      if(A[i,j]!=0){ A=red(A,i,j); }
    }
  }
  
  return(A);
}
example{
"EXAMPLE:";
ring r=0,(x),lp;
matrix A[2][2]=1,4,4,15;
"A="; print(A);
"sym_gauss(A)="; print(sym_gauss(A));
}

//////////////////////////////////////////////////////////////////////////////
proc orthogonalize(matrix A) 
"USAGE:    orthogonalize(A); A = constant matrix
RETURN:    matrix, orthogonal basis of the colum space of A
EXAMPLE:   example orthogonalize; shows an example "
{    
  int i;
  int j;
  int n=ncols(A);
  poly k;

  if(!const_mat(A)){
    "input is not a constant matrix";
    matrix B;
    return(B);
  }  

  module B=module(interred(A));
  
  for(i=1;i<=n;i=i+1) {
    for(j=1;j<i;j=j+1) {
      k=sp(B[j],B[j]);
      if (k==0) { "Error: vector with zero length"; return(matrix(B)); }
      B[i]=B[i]-(sp(B[i],B[j])/k)*B[j];
    }
  }
  
  return(matrix(B));
}
example{
"EXAMPLE";
ring r=0,(x),lp;
matrix A[4][4]=5,6,12,4,7,3,2,6,12,1,1,2,6,4,2,10;
"A=";print(A);
"orthogonalize(A)=";print(orthogonalize(A));
}

////////////////////////////////////////////////////////////////////////////
proc diag_test(matrix A)
"USAGE:   diag_test(A); A = const square matrix
NOTE:     Der Test ist nur fuer zerfallende Matrizen aussagefaehig.
          Parameterringe werden nicht unterstuetzt (benutzt factorize,gcd).
RETURNS:  int,  1 falls A diagonalisierbar, 0 nicht diagonalisierbar
               -1 keine Aussage moeglich, da A nicht zerfallend.
EXAMPLE:  example diag_test; shows an example"
{
  int i;
  int j;
  int n=nrows(A);       
  string mp=string(minpoly);
  string cs=charstr(basering);
  int np=0;

  if(ncols(A) != n) { 
    "input is not a square matrix";
    return(-1);
  }             

   if(!const_mat(A)){
    "input is not a constant matrix";
    return(-1);
  }  

  //Parameterring wegen factorize nicht erlaubt
  for(i=1;i<size(cs);i=i+1){
    if(cs[i]==","){np=np+1;} //Anzahl der Parameter
  }
  if(np>0){
        "parametric rings not allowed";
        return(-1);
  }

  //speichern des aktuellen Rings  
  def BR=basering;
  //setze R[t]
  execute("ring rt=("+charstr(basering)+"),(@t,"+varstr(basering)+"),lp;");
  execute("minpoly="+mp+";");
  matrix A=imap(BR,A);  

  intvec z;
  intvec s;
  poly X;         //characteristisches Polynom 
  poly dXdt;      //Ableitung von X nach t
  ideal g;        //ggT(X,dXdt)
  poly b;         //Komponente der Busadjunkten-Matrix
  matrix E[n][n]; //Einheitsmatrix
        
  E=E+1;
  A=E*@t-A;
  X=det(A);

  matrix Xfactors=matrix(factorize(X,1));       //zerfaellt die Matrtix ?
  int nf=ncols(Xfactors);

  for(i=1;i<=nf;i++){
    if(lead(Xfactors[1,i])>=@t^2){
      //"Die Matrix ist nicht zerfallend";
      setring BR;
      return(-1);
    }
  }     

  dXdt=diff(X,@t);
  g=std(ideal(gcd(X,dXdt)));

  //Busadjunkte
  z=2..n;       
  for(i=1;i<=n;i++){
    s=2..n;  
    for(j=1;j<=n;j++){
      b=det(submat(A,z,s));
    
      if(0!=reduce(b,g)){
        //"Die Matrix ist nicht diagonalisierbar!";
        setring BR;
        return(0);
      }
      
      s[j]=j;
    }
    z[i]=i;
  }
  
  //"Die Matrix ist diagonalisierbar";          
  setring BR;
  return(1);
}
example{
"EXAMPLE:";
ring r=0,(x),dp;
matrix A[4][4]=6,0,0,0,0,0,6,0,0,6,0,0,0,0,0,6;
"A="; print(A);
"diag_test(A);"; diag_test(A);
}

//////////////////////////////////////////////////////////////////////////////
proc busadj(matrix A)
"USAGE:   busadj(A);  A =  square matrix
RETURN:  list L; 
         L[1] contains the (n+1)coefficients of the characteristic polynomial of A
              X=L[1][1]+..+L[1][k]*t^(k-1)+..+(L[1][n+1]=1)*t^n
         L[2] contains the n matrices Hk of type nxn with
              bA=(Hn-1)*t^(n-1)+...+H0   ,( Hk=L[2][k+1] ) 
              this matrix polynomial forms the busadjoint bA=adj(E*t-A)
EXAMPLE: example busadj; shows an example"
{
  int k;
  int n=nrows(A);
  matrix E=unitmat(n);
  matrix H[n][n];
  matrix B[n][n];
  list bA;
  list X;
  list L;
  poly a;

  if(ncols(A) != n) { 
    "input is not a square matrix";
    return(L);
  }     

  bA=E;
  X[1]=1;
  for(k=1;k<n;k=k+1){
    B=A*bA[1]; //bA[1] is the last H
    a=-trace(B)/k;
    H=B+a*E;
    bA=insert(bA,H);
    X=insert(X,a);
  }
  B=A*bA[1];
  a=-trace(B)/n;
  X=insert(X,a);

  L=insert(L,bA);
  L=insert(L,X);
  return(L);
}
example{
"EXAMPLE";
ring r=0,(t,x),lp;
matrix A[2][2]=1,x2,x,x2+3x;
"A=";print(A);
"list L=busadj(A);";list L=busadj(A);
"poly X=L[1][1]+L[1][2]*t+L[1][3]*t2;";
poly X=L[1][1]+L[1][2]*t+L[1][3]*t2; X;
"matrix bA[2][2]=L[2][1]+L[2][2]*t;";
matrix bA[2][2]=L[2][1]+L[2][2]*t; print(bA);
"(t*unitmat(2)-A)*bA=X*E";
print(bA*(t*unitmat(2)-A));
}

//////////////////////////////////////////////////////////////////////////////
proc charpoly(matrix A,string v)
"USAGE:   charpoly(A,v); A = square matrix, v = name of a ringvariable
NOTE:    A must be constant in the variable v. The computation uses busadj(A).
RETURN:  poly, characteristic polynomial det(E*v-A)
EXAMPLE: example charpoly; shows an example"
{
  int n=nrows(A);
  int i;
  int j;
  poly X=0;
  def BR=basering;
  string s;
  list L;
  string mp=string(minpoly);

  if(ncols(A) != n){    
    "input is not a square matrix";
    return(X);
  }     

  //test for corect variable
  j=-1;
  for(i=1;i<=nvars(BR);i=i+1){
    if(varstr(i)==v){j=i;}
  }
  if(j==-1){
    v+" is not a variable in your basering";
    return(X);
  }
  
  //var can not be in A
  s="Wp(";
  for(i=1;i<=nvars(BR);i=i+1){
    if(i!=j){ s=s+"0";}
    else{ s=s+"1";}
    if(i<nvars(BR)){s=s+",";}
  }
  s=s+")";

  changeord("@R",s);
  execute("minpoly="+mp+";");
  matrix A=imap(BR,A);
  for(i=1;i<=n;i=i+1){
    if(deg(lead(A)[i])>=1){
      "matrix cannot contain the variable "+v;
      kill @R;
      setring BR;
      return(X);
    }
  }

  //get coefficients and build the char. poly
  kill @R;
  setring BR;
  L=busadj(A);
  for(i=1;i<=n+1;i=i+1){
    execute("X=X+L[1][i]*"+v+"^"+string(i-1)+";");
  }
  
  return(X); 
}
example{
"EXAMPLE";
ring r=0,(x,t),dp;
matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
"A="; print(A);
"charpoly(A,\"t\")="; charpoly(A,"t");
}

//////////////////////////////////////////////////////////////////////////////
proc adj(matrix A)
"USAGE:    adj(A);  A = square matrix
NOTE:     computation uses busadj(A)
RETURN:   adjoint 
EXAMPLE:  example adj; shows an example"
{
  int n=nrows(A);
  matrix Adj[n][n];
  list L;

  if(ncols(A) != n) {   
    "input is not a square matrix";
    return(Adj);
  }     
  
  L=busadj(A);
  Adj=(-1)^(n-1)*L[2][1];
  return(Adj);
  
}
example{
"EXAMPLE";
ring r=0,(t,x),lp;
matrix A[2][2]=1,x2,x,x2+3x;
"A="; print(A);
"matrix Adj[2][2]=adj(A);"; matrix Adj[2][2]=adj(A);
print(Adj);
"Adj*A=det(A)*E";
print(Adj*A);
}


//////////////////////////////////////////////////////////////////////////////
proc inverse_B(matrix A)
"USAGE:     inverse_B(A);  A = square matrix
RETURN:    list Inv with Inv[1]=matrix I and Inv[2]=poly p 
           and I*A = unitmat(n)*p; 
NOTE:      p=1 if 1/det(A) is computable and  p=det(A) if not
EXAMPLE:   example inverse_B; shows an example"
{
  int i;
  int n=nrows(A);
  matrix I[n][n];
  poly factor;
  list L;
  list Inv;
  
  if(ncols(A) != n) {   
    "input is not a square matrix";
    return(I);
  }     
  
  L=busadj(A);
  I=module(-L[2][1]);        //+-Adj(A)   

  if(reduce(1,std(L[1][1]))==0){ 
    I=I*lift(L[1][1],1)[1][1];      
    factor=1;
  }
  else{ factor=L[1][1];}     //=+-det(A) or 1
  Inv=insert(Inv,factor);
  Inv=insert(Inv,matrix(I));

  return(Inv);
}
example{
"EXAMPLE";
ring r=0,(x,y),lp;
matrix A[3][3]=x,y,1,1,x2,y,x+y,6,7;
"A="; print(A);
"list Inv=inverse_B(A);";list Inv=inverse_B(A);
"Inv[1]=";print(Inv[1]);
"Inv[2]=";print(Inv[2]);
"Inv[1]*A=";print(Inv[1]*A);
}


//////////////////////////////////////////////////////////////////////////////
proc det_B(matrix A)
"USAGE:     det_B(A);  A any matrix  
RETURN:    returns the determinant of A
NOTE:      the computation uses the busadj algorithm
EXAMPLE:   example det_B; shows an example"
{
  int n=nrows(A); 
  list L;

  if(ncols(A) != n){ return(0);}

  L=busadj(A);
  return((-1)^n*L[1][1]);
  
  
}
example{
"EXAMPLE";
ring r=0,(x),dp;
matrix A[10][10]=random(2,10,10)+unitmat(10)*x;
"A=";print(A);
"det_B(A)=";det_B(A); 
}


//////////////////////////////////////////////////////////////////////////////
proc inverse_L(matrix A)
"USAGE:     inverse_L(A);  A = square matrix
RETURN:    list Inv with Inv[1]=matrix I and Inv[2]=poly p 
           and I*A = unitmat(n)*p; 
NOTE:      p=1 if 1/det(A) is computable and  p=det(A) if not
EXAMPLE:   example inverse_L; shows an example"
{
  int n=nrows(A);
  matrix I;
  matrix E[n][n]=unitmat(n);
  poly factor;
  poly d=1;
  list Inv;

  if (ncols(A)!=n){
    "input is not a square matrix";
    return(I);
  }

  d=det(A);
  if(d==0){
    "matrix is not invertible";
    return(Inv);
  }

  // test if 1/det(A) exists
  if(reduce(1,std(d))!=0){ E=E*d;}
  
  I=lift(A,E);
  if(I==unitmat(n)-unitmat(n)){ //catch error in lift 
    "matrix is not invertible";
    return(Inv);
  }

  factor=d;      //=det(A) or 1
  Inv=insert(Inv,factor);
  Inv=insert(Inv,I);
  
  return(Inv);
}
example{
"EXAMPLE";
ring r=0,(x,y),lp;
matrix A[3][3]=x,y,1,1,x2,y,x+y,6,7;
"A="; print(A);
"list Inv=inverse_L(A);";list Inv=inverse_L(A);
"Inv[1]=";print(Inv[1]);
"Inv[2]=";print(Inv[2]);
"Inv[1]*A=";print(Inv[1]*A);
}

//////////////////////////////////////////////////////////////////////////////
proc gaussred(matrix A)
"USAGE:     gaussred(A);   any constant matrix A
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
           gives a row reduced matrix S, a permutation matrix P and a
           normalized lower triangular matrix U , with P*A=U*S            
EXAMPLE:   example gaussred; shows an example"
{
  int i;
  int j;
  int l;
  int k=0;
  int jp;
  int n=nrows(A);
  int m=ncols(A);
  int mr=n; //max. rang
  matrix P[n][n]=unitmat(n);
  matrix U[n][n]=P;
  poly c;
  poly pivo;
  int rang;
  list Z;

  if(!const_mat(A)){
    "input is not a constant matrix";
    return(Z);
  }

  if(n>m){mr=m;} //max. rang

  for(i=1;i<=mr;i=i+1){
    if((i+k)>m){break}; 
    
    //Test: Diagonalelement=0
    if(A[i,i+k]==0){
      jp=i;pivo=0;
      for(j=i+1;j<=n;j=j+1){
        c=abs(A[j,i+k]);
        if(pivo<c){ pivo=c;jp=j;}
      }
      if(jp != i){ //Zeilentausch
        for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
          c=A[i,j];         
          A[i,j]=A[jp,j];
          A[jp,j]=c;
        }
        for(j=1;j<=n;j=j+1){ //Zeilentausch in P
          c=P[i,j];        
          P[i,j]=P[jp,j];
          P[jp,j]=c;  
        }
      }
      if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
    }                          //i sollte im naechsten Lauf nicht erhoeht sein
    
    //Eliminationsschritt   
    for(j=i+1;j<=n;j=j+1){ 
      c=A[j,i+k]/A[i,i+k];
      for(l=i+k+1;l<=m;l=l+1){
        A[j,l]=A[j,l]-A[i,l]*c;
      }
      A[j,i+k]=0;  // nur wichtig falls k>0 ist
      A[j,i]=c;    // bildet U
    }
  rang=i;    
  }
 
  for(i=1;i<=mr;i=i+1){
    for(j=i+1;j<=n;j=j+1){
      U[j,i]=A[j,i];
      A[j,i]=0;
    }
  }
  
  Z=insert(Z,rang);
  Z=insert(Z,A);
  Z=insert(Z,U);
  Z=insert(Z,P);
  
  return(Z);
}
example{
"EXAMPLE";echo=2;
ring r=0,(x),dp;
matrix A[5][4]=1,3,-1,4,2,5,-1,3,1,3,-1,4,0,4,-3,1,-3,1,-5,-2;
print(A);
list Z=gaussred(A);
"P=";print(Z[1]);
"U=";print(Z[2]);
"S=";print(Z[3]);
"rank=";print(Z[4]);
"P*A=";print(Z[1]*A);
"U*S=";print(Z[2]*Z[3]);
}

//////////////////////////////////////////////////////////////////////////////
proc gaussred_pivot(matrix A)
"USAGE:     gaussred(A);   any constant matrix A
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
           gives n row reduced matrix S, a permutation matrix P and a
           normalized lower triangular matrix U , with
                    P*A=U*S
NOTE:      with row pivoting            
EXAMPLE:   example gaussred; shows an example"
{
  int i;
  int j;
  int l;
  int k=0;
  int jp;
  int n=nrows(A);
  int m=ncols(A);
  int mr=n; //max. rang
  matrix P[n][n]=unitmat(n);
  matrix U[n][n]=P;
  poly c;
  poly pivo;
  int rang;
  list Z;

  if(!const_mat(A)){
    "input is not a constant matrix";
    return(Z);
  }

  if(n>m){mr=m;} //max. rang

  for(i=1;i<=mr;i=i+1){
    if((i+k)>m){break}; 
    
    //Pivotisierung
    pivo=abs(A[i,i+k]);jp=i;
    for(j=i+1;j<=n;j=j+1){
      c=abs(A[j,i+k]);
      if(pivo<c){ pivo=c;jp=j;}
    }
    if(jp != i){ //Zeilentausch
      for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
        c=A[i,j];         
        A[i,j]=A[jp,j];
        A[jp,j]=c;
      }
      for(j=1;j<=n;j=j+1){ //Zeilentausch in P
        c=P[i,j];        
        P[i,j]=P[jp,j];
        P[jp,j]=c;  
      }
    }
    if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
                               //i sollte im naechsten Lauf nicht erhoeht sein
    //Eliminationsschritt   
    for(j=i+1;j<=n;j=j+1){ 
      c=A[j,i+k]/A[i,i+k];
      for(l=i+k+1;l<=m;l=l+1){
        A[j,l]=A[j,l]-A[i,l]*c;
      }
      A[j,i+k]=0;  // nur wichtig falls k>0 ist
      A[j,i]=c;    // bildet U
    }
  rang=i;    
  }
 
  for(i=1;i<=mr;i=i+1){
    for(j=i+1;j<=n;j=j+1){
      U[j,i]=A[j,i];
      A[j,i]=0;
    }
  }
  
  Z=insert(Z,rang);
  Z=insert(Z,A);
  Z=insert(Z,U);
  Z=insert(Z,P);
  
  return(Z);
}
example{
"EXAMPLE";echo=2;
ring r=0,(x),dp;
matrix A[5][4]=1,3,-1,4,2,5,-1,3,1,3,-1,4,0,4,-3,1,-3,1,-5,-2;
print(A);
list Z=gaussred(A);
"P=";print(Z[1]);
"U=";print(Z[2]);
"S=";print(Z[3]);
"rank=";print(Z[4]);
"P*A=";print(Z[1]*A);
"U*S=";print(Z[2]*Z[3]);
}


//////////////////////////////////////////////////////////////////////////////
proc gauss_nf(matrix A)
"USAGE:    gauss_nf(A); any constant matrix A
RETURN:   matrix; gauss normal form of A
EXAMPLE:  example gauss_nf; shows an example"
{
  list Z;
  if(!const_mat(A)){
    "input is not a constant matrix";
    return(A);
  }
  Z=gaussred(A);
  return(Z[3]);
}
example{
"EXAMPLE";echo=2;
ring r=7,(x),dp;
matrix A[4][4]= 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
print(gauss_nf(A));
}

//////////////////////////////////////////////////////////////////////////////
proc mat_rk(matrix A)
"USAGE:    mat_rk(A); any constant matrix A
RETURN:   int; rank of A 
EXAMPLE:  example mat_rk; shows an example"
{
  list Z;
  if(!const_mat(A)){
    "input is not a constant matrix";
    return(-1);
  }
  Z=gaussred(A);
  return(Z[4]);
}
example{
"EXAMPLE";echo=2;
ring r=7,(x),dp;
matrix A[4][4]= 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
mat_rk(A);
}

//////////////////////////////////////////////////////////////////////////////
proc U_D_O(matrix A)
"USAGE:     gaussred(A);   constant invertible matrix A
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=D , Z[4]=O
           gives a permutation matrix P, a
           normalized lower triangular matrix U ,
           a diagonal matrix D, and a normalized upper triangular matrix O
           withP*A=U*D*O
NOTE:      Z[1]=-1 means A is not regular
EXAMPLE:   example gaussred; shows an example"
{
  int i;
  int j;
  int n=nrows(A);
  list Z;
  list L;
  matrix O[n][n]=unitmat(n);
  matrix D[n][n];

  if (ncols(A)!=n){
    "input is not a square matrix";
    return(Z);
  }
  if(!const_mat(A)){
    "input is not a constant matrix";
    return(Z);
  }
  
  L=gaussred(A);

  if(L[4]!=n){
    "input is not a invertible matrix";
    Z=insert(Z,-1);//hint for calling procedures 
    return(Z);
  }

  D=L[3];

  for(i=1;i<=n;i=i+1){
    for(j=i+1;j<=n;j=j+1){
      O[i,j]=D[i,j]/D[i,i];
      D[i,j]=0;
    }
  }

  Z=insert(Z,O);
  Z=insert(Z,D);
  Z=insert(Z,L[2]);
  Z=insert(Z,L[1]);
  return(Z);
}
example{
"EXAMPLE";echo=2;
ring r=0,(x),dp;
matrix A[5][5]=10,4,0,-9,8,-3,6,-6,-4,9,0,3,-1,-9,-8,-4,-2,-6,
               -10,10,-9,5,-1,-6,5; 
print(A);
list Z=U_D_O(A);
"P=";print(Z[1]);
"U=";print(Z[2]);
"D=";print(Z[3]);
"O=";print(Z[4]);
"P*A=";print(Z[1]*A);
"U*D*O=";print(Z[2]*Z[3]*Z[4]);
}


//////////////////////////////////////////////////////////////////////////////
proc pos_def(matrix A)
"USAGE:     pos_def(A); A = const, symmetric ;       
RETURN:    int; 1=positiv definit , 0=not positiv definit , -1= I dont know 
EXAMPLE:   example pos_def; shows an example"
{
  int n=nrows(A);
  int j;
  list Z;
  matrix H[n][n];

  if (ncols(A)!=n){
    //"input is not a square matrix";
    return(0);
  }
  if(!const_mat(A)){
    "input is not a constant matrix";
    return(-1);
  }     
  if(deg(std(A-transpose(A))[1])!=-1){ 
    "input is not a hermitesche (symetric) matrix";
    return(-1);
  }
  
  Z=U_D_O(A);

  if(Z[1]==-1){return(0);} //A not regular , therefor not pos.definit

  H=Z[1];
  //es fand Zeilentausch statt: also nicht positiv definit 
  if(deg(std(H-unitmat(n))[1])!=-1){return(0);}
  
  H=Z[3];
  for(j=1;j<=n;j=j+1){
    if(H[j,j]<=0){return(0);}//eigenvalue<=0 ,not pos.definit
  }

  return(1); //positiv definit;
}
example{
"EXAMPLE";
}


//////////////////////////////////////////////////////////////////////////////
proc solve(matrix A,matrix b)
"USAGE:
RETURN:
EXAMPLE:"
{
  int i;
  int j;
  int k;
  int rc;
  int r;
  int n=nrows(A);
  int m=ncols(A);
  int n_b=nrows(b);
  matrix Ab[n][m+1];
  matrix X[m][1];
  poly c;

  list Z;
        
  if(ncols(b)!=1){
    "right side b is not a nx1 matrix";
    return(X);
  }

  if(!const_mat(A)){
    "input is not a constant matrix";
    return(X);
  }     
  
  if(n_b>n){
    for(i=n;i<=n_b;i=i+1){
      if(b[i,1]!=0){
        "right side b not in Image(A)";
        return X;
      }
    }    
  }
  
  if(n_b<n){ 
    matrix copy[n_b][1]=b;
    matrix b[n][1]=0;
    for(i=1;i<=n_b;i=i+1){
      b[i,1]=copy[i,1];
    }
  }
  
  r=mat_rk(A);
        
  //1. b constant vector
  if(const_mat(b)){     
    //extend A with b
    for(i=1;i<=n;i=i+1){
      for(j=1;j<=m;j=j+1){
        Ab[i,j]=A[i,j];
      }
      Ab[i,m+1]=b[i,1];
    }
      
    //Gaussreduction
    Z=gaussred(Ab);
    Ab=Z[3]; //normal form
    rc=Z[4];  //rank(Ab)  
 print(Ab);

    if(r<rc){
        "no solution";
        return(X); 
    }
    k=m;        
    for(i=r;i>=1;i=i-1){
      
      j=1;      
      while(Ab[i,j]==0){j=j+1;}// suche Ecke 
   
      for(;k>j;k=k-1){ X[k]=0;}//springe zur Ecke
     

      c=Ab[i,m+1]; //i-te Komponene von b
      for(j=m;j>k;j=j-1){
        c=c-X[j,1]*Ab[i,j];
      }
      if(Ab[i,k]==0){
        X[k,1]=1; //willkuerlich
      }
      else{ 
        X[k,1]=c/Ab[i,k];
      }
      k=k-1;
      if(k==0){break;}
    }
    
    
  }//endif (const b)
  else{  //b not constant
    "!not implementet!";
 
  }

  return(X);
}
example{
"EXAMPLE";echo=2;
ring r=0,(x),dp;
int n=random(1,10);
int m=random(1,10);
matrix A[n][m]=random(10,n,m);
matrix b[n][1]=random(10,n,1);
"A=";print(A);
"b=";print(transpose(b));
matrix X[n][1]=solve(A,b);
"X=";print(transpose(X));
print(A*X);
}