version="$Id: control.lib,v 1.28 2005-04-29 14:53:50 levandov Exp $";
category="System and Control Theory";
info="
LIBRARY:  control.lib Procedures for System and Control Theory 
AUTHORS:  Oleksandr Iena       yena@mathematik.uni-kl.de
@*        Markus Becker        mbecker@mathematik.uni-kl.de
@*        Viktor Levandovskyy  levandov@mathematik.uni-kl.de

SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' (Project of Dr. E. Zerz
and V. Levandovskyy), Uni Kaiserslautern

NOTE: This library provides algebraic analysis tools for System and Control Theory

PROCEDURES:
control (R);            analysis of controllability-related properties of R (using Ext modules),
control2 (R);           analysis of controllability-related properties of R (using dimension),
autonom (R);            analysis of autonomy-related properties of R (using Ext modules),
autonom2 (R);           analysis of autonomy-related properties of R (using dimension),
LeftKernel (R);         a left kernel of R,
RightKernel (R);        a right kernel of R,
LeftInverse (R);        a left inverse of R,
RightInverse (R);       a right inverse of R,
smith (M);              a Smith form of a module M,
colrank (M);            a column rank of M as of matrix,
genericity (M);         analysis of the genericity of parameters,
canonize (L);           Groebnerification for modules in the output of control/autonomy procs,
iostruct (R);           computes an I/O-structure of behavior given by a module R
FindTorsion (R, I);     generators of the submodule of a module R, annihilated by the ideal I.


AUXILIARY PROCEDURES:
declare(NameOfRing,Variables [,Parameters, Ordering]);    define the ring easily
view();                      well-formatted output of lists, modules and matrices

NOTE (EXAMPLES): In order to use examples below, execute the commands
@* def A = exAntenna(); setring A; 
@* Thus A will become a basering from the example with the predefined module R (transposed), corresponding to the system. After that you can just type in 
@* control(R); // respectively autonom(R); 
and check the result.


EXAMPLES (AUTONOMY):
exCauchy();                  example of 1-dimensional Cauchy equation,
exCauchy2();                 example of 2-dimensional Cauchy equation,
exZerz();                    example from the lecture of Eva Zerz,

EXAMPLES (CONTROLLABILITY):
ex1();                       example of noncontrollable system,   
ex2();                       example of controllable system ,
exAntenna();                 Antenna,
exEinstein();                Einstein equations,
exFlexibleRod();             Flexible Rod,
exTwoPendula();              Two Pendula mounted on a cart,
exWindTunnel();              Wind Tunnel. 
";

// NOTE: static things should not be shown for end-user
// static Ext_Our(...)                  Copy of Ext_R from 'homolog.lib' in commutative case;
// static is_zero_Our(module R)         Copy of is_zero from 'poly.lib';
//  static space(int n)           Procedure used inside the procedure 'Print' to have a better formatted output
// static control_output();      Generating the output for the procedure 'control'
// static autonom_output();      Generating the output for the procedure 'autonom' and 'autonom2'
// static extgcd_Our(poly p, poly q)  Computes extgcd of p and q. for versions ealier than 2006 extgcd has a bug and is therefore not used
// static normalize_Our(matrix N, matrix Q) normalizes the columns of N and divides the columns of Q through the leading coefficients of the columns of N

LIB "homolog.lib"; 
LIB "poly.lib";
LIB "primdec.lib";
LIB "matrix.lib";

//---------------------------------------------------------------
static proc Opt_Our()
"USAGE: Opt_Our();
RETURN: intvec, where previous options are stored
PURPOSE: save previous options and set customized options
"
{
  intvec v;
  v=option(get);
  option(redSB);
  option(redTail);
  return (v);
}
//---------------------------------------------------------------
proc declare(string NameOfRing, string Variables, list #)
"USAGE: declare(NameOfRing, Variables,[Parameters, Ordering]);
@*          NameOfRing:  string with name of ring,
@*	    Variables:   string with names of variables separated by commas(e.g. \"x,y,z\"),
@*	    Parameters: string of parameters in the ring separated by commas(e.g. \"a,b,c\"),
@*	    Ordering:   string with name of ordering (by default, the ordering (dp,C) is used)
PURPOSE: define the ring easily
RETURN:  no return value
EXAMPLE:  example declare; shows an example 
"
{
  if(size(#)==0)
  {
    execute("ring "+NameOfRing+"=0,("+Variables+"),dp;");
  }
  else
  {
    if(size(#)==1)
    {
      execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),dp;" );
    }
    else
    { 
      if( (size(#[1])!=0)&&(#[1]!=" ") )
      {
        execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),("+#[2]+");" );
      }
      else
      {
        execute( "ring " + NameOfRing + "=0,("+Variables+"),("+#[2]+");" ); 
      };
    };
  };
  keepring(basering);
}
example
{"EXAMPLE:";echo = 2;
  string v="x,y,z";
  string p="q,p";
  string Ord ="c,lp";
  //----------------------------------
  declare("Ring_1",v);
  print(nameof(basering));
  print(basering);
  //----------------------------------
  declare("Ring_2",v,p);
  print(basering);
  print(nameof(basering));
  //----------------------------------
  declare("Ring_3",v,p,Ord);
  print(basering);
  print(nameof(basering));
  //----------------------------------
  declare("Ring_4",v,"",Ord);
  print(basering);
  print(nameof(basering));
  //----------------------------------
  declare("Ring_5",v," ",Ord);
  print(basering);
  print(nameof(basering));
};
//
// maybe reasonable to add this in declare
//
//  print("Please enter your representation matrix in the following form: 
//  module R=[1st row],[2nd row],...");
//  print("Type the command: R=transpose(R)");
//  print(" To compute controllability please enter: control(R)");
//  print(" To compute autonomy please enter: autonom(R)");
//-------------------------------------------------------------------------

static proc space(int n)
"USAGE:spase(n); 
         n: integer, number of needed spaces
RETURN:  string consisting of n spaces
NOTE:  the procedure is used in the procedure 'view' to have a better formatted output
"
{
  int i;
  string s="";
  for(i=1;i<=n;i++)
  {
    s=s+" ";
  };
  return(s);
};
//-----------------------------------------------------------------------------
proc view(M)
"USAGE:  view(M);
           M:  any type
RETURN:  no return value
PURPOSE:  procedure for ( well-) formatted output of modules, matrices, lists of modules, matrices;
          shows everything even if entries are long
NOTE:  in case of other types( not 'module', 'matrix', 'list') works just as standard 'print' procedure  
EXAMPLE:  example view; shows an example
"
{
  // to be replaced with something more feasible
  if ( (typeof(M)=="module")||(typeof(M)=="matrix") )
  {
  int @R=nrows(M);
  int @C=ncols(M);
  int i;
  int j;
  list MaxLength=list();
  int Size=0;
  int max;
  string s;
  
  for(i=1;i<=@C;i++)
  { 
    max=0; 
    
    for(j=1;j<=@R;j++)
    {
      Size=size( string( M[j,i] ) );
      if( Size>max )
      {
        max=Size;
      };
    };
    MaxLength[i] = max;
  };
  
  for(i=1;i<=@R;i++)
  {
    s="";
    for(j=1;j<@C;j++)
    {
      s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ) +",";
    };
    
    s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) );

    if (i!=@R)
    {
      s=s+",";
    };
    print(s);
  };

  return();    
  };
  
  if(typeof(M)=="list")
  {
    int sz=size(M);
    int i;
    for(i=1;i<=sz;i++)
    {
      view(M[i]);
      print("");
    }; 

    return();
  };
  print(M);
  return();
}
example
{"EXAMPLE:";echo = 2;
  ring r;
  matrix M[1][3] = x,y,z;
  print(M);
  view(M);
};
//--------------------------------------------------------------------------
proc RightKernel(matrix M)
"USAGE:  RightKernel(M);   M a matrix
PURPOSE: computes the right kernel of matrix M (a module of all elements v such that Mv=0)
RETURN:  module
EXAMPLE: example RightKernel; shows an example
"
{ 
  return(modulo(M,std(0)));
}
example
{"EXAMPLE:";echo = 2;
  ring r = 0,(x,y,z),dp;
  matrix M[1][3] = x,y,z;
  print(M);
  matrix R = RightKernel(M);
  print(R);
  // check:
  print(M*R);
};
//-------------------------------------------------------------------------
proc LeftKernel(matrix M)
"USAGE:  LeftKernel(M);   M a matrix
PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0)
RETURN:  module
EXAMPLE:  example LeftKernel; shows an example
"
{
  return( transpose( modulo( transpose(M),std(0) ) ) );
}
example
{"EXAMPLE:";echo = 2;
  ring r= 0,(x,y,z),dp;
  matrix M[3][1] = x,y,z;
  print(M);
  matrix L = LeftKernel(M);
  print(L);
  // check:
  print(L*M);
};
//------------------------------------------------------------------------
proc LeftInverse(module M)
"USAGE:  LeftInverse(M);  M a module
PURPOSE: computes such a matrix L, that LM == Id;
RETURN:  module
EXAMPLE: example LeftInverse; shows an example
NOTE: exists only in the case when Id \subseteq M!
"
{
  // it works also for the NC case;
  int NCols = ncols(M);
  module Id = freemodule(NCols);
  module N  = transpose(M);
  intvec old_opt=Opt_Our();
  Id = std(Id);
  matrix T;
  // check the correctness (Id \subseteq M)
  // via dimension: dim (M) = -1!
  int d = dim_Our(N);
  if (d != -1)
  {
    // No left inverse exists
    return(matrix(0));
  }
  matrix T2 = lift(N, Id);
  T2 =  transpose(T2);
  option(set,old_opt);  // set the options back
  return(T2);
}
example
{ // trivial example:
"EXAMPLE:";echo =2;
  ring r = 0,(x,z),dp;
  matrix M[2][1] = 1,x2z;
  print(M);
  print( LeftInverse(M) );  
  kill r;
  // derived from the exTwoPendula
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  matrix U[3][1];
  U[1,1]=(-L2)*Dt^4+(g)*Dt^2;
  U[2,1]=(-L1)*Dt^4+(g)*Dt^2;
  U[3,1]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2);
  module M = module(U);
  module L = LeftInverse(M);
  print(L);
};
//-----------------------------------------------------------------------
proc RightInverse(module R)
"USAGE:  RightInverse(M);  M a module
PURPOSE: computes such a matrix L, that ML == Id;
RETURN:  module
EXAMPLE: example RightInverse; shows an example
NOTE: exists only in the case when Id \subseteq M!
"
{
  return(transpose(LeftInverse(transpose(R))));
}
example
{ "EXAMPLE:";echo =2;
  // trivial example:
  ring r = 0,(x,z),dp;
  matrix M[1][2] = 1,x2+z;
  print(M);
  print( RightInverse(M) );  
  kill r;
  // derived from the exTwoPendula
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  matrix U[1][3];
  U[1,1]=(-L2)*Dt^4+(g)*Dt^2;
  U[1,2]=(-L1)*Dt^4+(g)*Dt^2;
  U[1,3]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2);
  module M = module(U);
  module L = RightInverse(M);
  print(L);
};
//-----------------------------------------------------------------------
static proc dim_Our(module R)
{
  int d;
  if (attrib(R,"isSB")<>1)
  {
    R = std(R);
  }
  d = dim(R);
  return(d);
}
//-----------------------------------------------------------------------
static proc Ann_Our(module R)
{
  return(Ann(R));
}
//-----------------------------------------------------------------------
static proc Ext_Our(int i, module R, list #)
{
  // mimicking 'Ext_R' from homolog.lib
  int ExtraArg = ( size(#)>0 );
  if (ExtraArg)
  {
    return( Ext_R(i,R,#[1]) );
  }
  else
  {
    return( Ext_R(i,R) );
  }
}
//------------------------------------------------------------------------
static proc is_zero_Our
{
  //just a copy of 'is_zero' from "poly.lib" patched with GKdim
  int d = dim_Our(std(#[1]));
  int a = ( d==-1 );
  if( size(#) >1 ) { list L=a,d; return(L); }
  return(a);
  //  return( is_zero(R) ) ;
};
//------------------------------------------------------------------------
static proc control_output(int i, int NVars, module R, module Ext_1, list Gen) 
//static proc control_output(int i, int NVars, module R, module Ext_1) 
"USAGE:  control_output(i, NVars, R, Ext_1), where
           i:  integer, number of first nonzero Ext or 
             just a number of variables in a base ring + 1 in case that all the Exts are zero,
	   NVars:  integer, number of variables in a base ring,
	   R:  module R (cokernel representation),
	   Ext_1:  module, the first Ext(its cokernel representation)      
RETURN:  list with all the contollability properties of the system which is to be returned in 'control' procedure
NOTE:  this procedure is used in 'control' procedure
"
{
  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
  // Is not clear what to do with gl.dim of qrings
  string DofS = "dimension of the system:";
  string Fn   = "number of first nonzero Ext:";
  string Gen_mes = "Parameter constellations which might lead to a non-controllable system:";

  module RK = RightKernel(R); 
  int d=dim_Our(std(transpose(R)));

  if (i==1)
  { 
    return( 
            list ( Fn,
	           i,
	          "not controllable , image representation for controllable part:",
         	   RK,	    
	          "kernel representation for controllable part:",
		   LeftKernel( RK ),
		  "obstruction to controllability",
		   Ext_1,
		  "annihilator of torsion module (of obstruction to controllability)",
		   Ann_Our(Ext_1),
		   DofS,
		   d
		 ) 
          );
  };
  
  if(i>NVars)
  { 
    return( list(  Fn,
                   -1, 
                  "strongly controllable(flat), image representation:",
		   RK,
		  "left inverse to image representation:",
		   LeftInverse(RK),
		   DofS,
		   d,
		   Gen_mes,
		   Gen) 
          );
  };
  
  //
  //now i<=NVars
  //
       
  if( (i==2) )
  {
    return( list( Fn,
                  i, 
                 "controllable, not reflexive, image representation:",
		  RK,
                  DofS,
                  d,
		  Gen_mes,
		  Gen) 
          ); 
  };
    
  if( (i>=3) )
  {
    return( list ( Fn,
                   i, 
                  "reflexive, not strongly controllable, image representation:",
		   RK,
		   DofS,
		   d,
   		   Gen_mes,
		   Gen) 
          );
  };  
};  
//-------------------------------------------------------------------------

proc control(module R)
"USAGE:  control(R);  R a module (R is the matrix of the system of equations to be investigated)
PURPOSE: compute the list of all the properties concerning controllability of the system (behavior), represented by the matrix R
RETURN:  list
EXAMPLE:  example control; shows an example
"
{
  int i;
  int NVars=nvars(basering); 
  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
  int ExtIsZero;
  intvec v=Opt_Our();
  module R_std=std(R);
  module Ext_1 = std(Ext_Our(1,R_std));
   
  ExtIsZero=is_zero_Our(Ext_1);
  i=1;
  while( (ExtIsZero) && (i<=NVars) )
  {
    i++;
    ExtIsZero = is_zero_Our( Ext_Our(i,R_std) );
  };
  matrix T=lift(R,R_std);
  list l=genericity(T);
  option(set,v);

  return( control_output( i, NVars, R, Ext_1, l ) );
}
example
{"EXAMPLE:";echo = 2; 
  //Wind Tunnel
  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
  module R;
  R = [D1+a, -k*a*delta, 0, 0],
      [0, D1, -1, 0],
      [0, omega^2, D1+2*zeta*omega, -omega^2];
  R=transpose(R);
  view(R);
  view(control(R));
};
//--------------------------------------------------------------------------
proc control2(module R)
"USAGE:  control2(R); R a module (R is the matrix of the system of equations to be investigated)
PURPOSE: computes list of all the properties concerning controllability of the system (behavior), represented by the  matrix R
RETURN: list
EXAMPLE:  example control2; shows an example
NOTE: same as control(R); but using dimensions, this procedure only works for full row rank matrices
"
{
  if( nrows(R) != colrank(transpose(R)) )
  {
    return ("control2 cannot be applied, since R does not have full row rank");
  }  
  intvec v=Opt_Our();
  module R_std=std(R);
  int d=dim_Our(R_std);
  int NVars=nvars(basering);
  int i=NVars-d;
  module Ext_1=std(Ext_Our(1,R_std));
  matrix T=lift(R,R_std);
  list l=genericity(T);
  option(set, v);
  return( control_output( i, NVars, R, Ext_1, l));
}
example
{"EXAMPLE:";echo = 2; 
  //Wind Tunnel
  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
  module R;
  R = [D1+a, -k*a*delta, 0, 0],
      [0, D1, -1, 0],
      [0, omega^2, D1+2*zeta*omega, -omega^2];
  R=transpose(R);
  view(R);
  view(control2(R));
};
//------------------------------------------------------------------------
proc colrank(module M)
"USAGE: proc colrank(M), M a matrix/module
PURPOSE: compute the column rank of M as of matrix
RETURN: int
NOTE: this procedure uses bareiss-algorithm which might not terminate in some cases
"
{
  module M_red  = bareiss(M)[1];
  int NCols_red = ncols(M_red);
  return (NCols_red);
}
example
{"EXAMPLE: ";echo = 2;
  // de Rham complex
  ring r=0,(D(1..3)),dp;
  module R;
  R=[0,-D(3),D(2)],
    [D(3),0,-D(1)],
    [-D(2),D(1),0];
  R=transpose(R);
  colrank(R); 
};
  
//------------------------------------------------------------------------
static proc autonom_output( int i, int NVars, module RC, int R_rank ) 
"USAGE:  proc autonom_output(i, NVars, RC, R_rank)
           i:  integer, number of first nonzero Ext or 
	       just number of variables in a base ring + 1 in case that all the Exts are zero
	   NVars:  integer, number of variables in a base ring  
	   RC: module, kernel-representation of controllable part of the system
	   R_rank: integer, column rank of the representation matrix
PURPOSE: compute all the autonomy properties of the system which is to be returned in 'autonom' procedure
RETURN:  list 
NOTE:  this procedure is used in 'autonom' procedure
"
{
  int d=NVars-i;//that is the dimension of the system
  string DofS="dimension of the system:";
  string Fn = "number of first nonzero Ext:";
  if(i==0)
  { 
    return( list(  Fn,
                   i,
                  "not autonomous",
		   "kernel representation for controllable part",
		   RC,
		   "column rank of the matrix",
		   R_rank,
		   DofS,
		   d ) 
	  );
  };
  
  if( i>NVars )
  {  
    return( list( Fn,
                  -1,
                  "trivial",
		  DofS,
		  d ) 
          );
  };
  
  //
  //now i<=NVars
  //
    
     
  if( i==1 )
  // in case that NVars==1 there is no sense to consider the notion
  // of strongly autonomous behavior, because it does not imply 
  // that system is overdetermined in this case
  {
    return( list ( Fn,
                   i, 
                  "autonomous, not overdetermined",
		   DofS,
		   d ) 
          ); 
  };
    
  if( i==NVars )
  {  
    return( list(  Fn,
                   i,
                  "strongly autonomous(fin. dimensional),in particular overdetermined",
		   DofS,
		   d) 
	  );
  };
  
  if( i<NVars )
  {
    return( list ( Fn,
                   i,
                  "overdetermined, not strongly autonomous",
		   DofS, 
		   d) 
          );
  };
};  
//--------------------------------------------------------------------------
proc autonom2(module R)
"USAGE:  autonom2(R);   R a module (R is a matrix of the system of equations which is to be investigated)
PURPOSE: computes the list of all the properties concerning autonomy of the system (behavior), represented by the matrix R
RETURN: list
NOTE:  this procedure is analogous to 'autonom' but uses dimension calculations
EXAMPLE:  example autonom2; shows an example
"
{ 
  int d;
  int NVars = nvars(basering);
  module RT = transpose(R);
  module RC;  //for computation of controllable part if if exists
  int R_rank = ncols(R);
  d     = dim_Our( std(RT) );  //this is the dimension of the system 
  int i = NVars-d;  //First non-zero Ext
  if( d==0 )
    {
      RC=LeftKernel(RightKernel(R));
      R_rank=colrank(R);
    }
  return( autonom_output(i,NVars,RC,R_rank) );
}
example
{"EXAMPLE:"; echo = 2;
  //Cauchy
  ring r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
	    [s3,-s4],
	    [s4, s3];	    
  R=transpose(R);
  view( R );
  view( autonom2(R) );
};  
//----------------------------------------------------------
proc autonom(module R)
"USAGE:  autonom(R);   R a module (R is a matrix of the system of equations which is to be investigated)
PURPOSE: find all the properties concerning autonomy of the system (behavior) represented by the  matrix R
RETURN:  list
EXAMPLE: example autonom; shows an example
"
{
  int NVars=nvars(basering);
  int ExtIsZero;
  module RT=transpose(R);
  module RC;
  int R_rank=ncols(R);
  ExtIsZero=is_zero_Our(Ext_Our(0,RT));     
  int i=0;
  while( (ExtIsZero)&&(i<=NVars) )
  {
    i++;
    ExtIsZero = is_zero_Our(Ext_Our(i,RT));
  };
  if (i==0)
  {
    RC=LeftKernel(RightKernel(R));
    R_rank=colrank(R);
  }
  return(autonom_output(i,NVars,RC,R_rank));
}
example
{"EXAMPLE:"; echo = 2;
  // Cauchy
  ring r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
	    [s3,-s4],
	    [s4, s3];	    
  R=transpose(R);
  view( R );
  view( autonom(R) );
};  

//----------------------------------------------------------
//
//Some example rings with defined systems
//----------------------------------------------------------
//autonomy:
//
//----------------------------------------------------------
proc exCauchy()
{
  ring @r=0,(s1,s2),dp;
  module R= [s1,-s2],
            [s2, s1];
  R=transpose(R);
  export R;
  return(@r);	    
};
//----------------------------------------------------------
proc exCauchy2()
{
  ring @r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
	    [s3,-s4],
	    [s4, s3];	    
  R=transpose(R);
  export R;
  return(@r);	    
}; 
//----------------------------------------------------------
proc exZerz()
{ 
  ring @r=0,(d1,d2),dp;
  module R=[d1^2-d2],
           [d2^2-1];
  R=transpose(R);
  export R;
  return(@r);	    
};  
//----------------------------------------------------------
//control
//----------------------------------------------------------
proc ex1()
{
  ring @r=0,(s1,s2,s3),dp;
  module R=[0,-s3,s2],
           [s3,0,-s1];
  R=transpose(R);	   
  export R;
  return(@r);
};
//----------------------------------------------------------
proc ex2()
{ 
  ring @r=0,(s1,s2,s3),dp;
  module R=[0,-s3,s2],
           [s3,0,-s1],
           [-s2,s1,0];
  R=transpose(R);	   
  export R;
  return(@r); 
};
//----------------------------------------------------------
proc exAntenna()
{
  ring @r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp;
  module R;
  R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
      [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
      [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
      [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
      [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
      [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];

  R=transpose(R);
  export R;
  return(@r);  
};

//----------------------------------------------------------

proc exEinstein()
{
  ring @r = 0,(D(1..4)),dp;
  module R =
  [D(2)^2+D(3)^2-D(4)^2, D(1)^2, D(1)^2, -D(1)^2, -2*D(1)*D(2), 0, 0, -2*D(1)*D(3), 0, 2*D(1)*D(4)],
  [D(2)^2, D(1)^2+D(3)^2-D(4)^2, D(2)^2, -D(2)^2, -2*D(1)*D(2), -2*D(2)*D(3), 0, 0, 2*D(2)*D(4), 0], 
  [D(3)^2, D(3)^2, D(1)^2+D(2)^2-D(4)^2, -D(3)^2, 0, -2*D(2)*D(3), 2*D(3)*D(4), -2*D(1)*D(3), 0, 0],
  [D(4)^2, D(4)^2, D(4)^2, D(1)^2+D(2)^2+D(3)^2, 0, 0, -2*D(3)*D(4), 0, -2*D(2)*D(4), -2*D(1)*D(4)],
  [0, 0, D(1)*D(2), -D(1)*D(2), D(3)^2-D(4)^2, -D(1)*D(3), 0, -D(2)*D(3), D(1)*D(4), D(2)*D(4)],
  [D(2)*D(3), 0, 0, -D(2)*D(3),-D(1)*D(3), D(1)^2-D(4)^2, D(2)*D(4), -D(1)*D(2), D(3)*D(4), 0],
  [D(3)*D(4), D(3)*D(4), 0, 0, 0, -D(2)*D(4), D(1)^2+D(2)^2, -D(1)*D(4), -D(2)*D(3), -D(1)*D(3)],
  [0, D(1)*D(3), 0, -D(1)*D(3), -D(2)*D(3), -D(1)*D(2), D(1)*D(4), D(2)^2-D(4)^2, 0, D(3)*D(4)],
  [D(2)*D(4), 0, D(2)*D(4), 0, -D(1)*D(4), -D(3)*D(4), -D(2)*D(3), 0, D(1)^2+D(3)^2, -D(1)*D(2)], 
  [0, D(1)*D(4), D(1)*D(4), 0, -D(2)*D(4), 0, -D(1)*D(3), -D(3)*D(4), -D(1)*D(2), D(2)^2+D(3)^2];

  R=transpose(R);
  export R;
  return(@r);  
};

//----------------------------------------------------------
proc exFlexibleRod()
{
  ring @r = 0,(D1, delta), dp;
  module R;
  R = [D1, -D1*delta, -1], [2*D1*delta, -D1-D1*delta^2, 0];
  
  R=transpose(R);
  export R;
  return(@r);  
};

//----------------------------------------------------------
proc exTwoPendula()
{  
  ring @r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  module R = [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2], 
             [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
             [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];

  R=transpose(R);
  export R;
  return(@r);  
};
//----------------------------------------------------------
proc exWindTunnel()
{ 
  ring @r = (0,a, omega, zeta, k),(D1, delta),dp;
  module R = [D1+a, -k*a*delta, 0, 0],
             [0, D1, -1, 0],
             [0, omega^2, D1+2*zeta*omega, -omega^2];

  R=transpose(R);
  export R;
  return(@r);  
};
//----------------------------------------------------------
proc genericity(matrix M)
"USAGE:  genericity(M), M is a matrix/module
PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis
RETURN:  list (of strings) 
NOTE: we strongly recommend to switch on the redSB and redTail options; 
@*    the procedure is effective with the lift procedure for modules with parameters
EXAMPLE:  example genericity; shows an example
"
{
  // returns "-", if there are no parameters!
  if (npars(basering)==0)
  {
    return("-");
  }
  list RT = evas_genericity(M); // list of strings
  if ((size(RT)==1) && (RT[1] == ""))
  {
    return("-");
  }
  return(RT);
}
example
{  // TwoPendula
  "EXAMPLE:"; echo = 2;
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  module RR = 
    [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2], 
    [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
    [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
  module R = transpose(RR);
  module SR = std(R);
  matrix T = lift(R,SR);
  genericity(T);
  //-- The result might be different when computing reduced bases:
  matrix T2;
  option(redSB);
  option(redTail);
  module SR2 = std(R);
  T2 =  lift(R,SR2);
  genericity(T2);
}
//---------------------------------------------------------------
static proc victors_genericity(matrix M)
{
  // returns "-", if there are no parameters!
  if (npars(basering)==0)
  {
    return("-");
  }
  int plevel = printlevel-voice+2;
  // M is a matrix over a ring with params and vars;
  ideal I = ideal(M); // a list of entries
  I = simplify(I,2); // delete 0's
  // decompose every coeff in every poly
  int i;
  int s = size(I);
  ideal NM;
  poly p;
  number num;
  int  cl=1;
  intvec ZeroVec; ZeroVec[nvars(basering)] = 0;
  intvec W;
  ideal Numero, Denomiro;
  int cNu=0; int cDe=0;
  for (i=1; i<=s; i++)
  {
    // remove contents and add them as polys
    p   = I[i];
    W   = leadexp(p);
    if (W == ZeroVec) // i.e. just a coef
    {
      num    = denominator(leadcoef(p)); // from poly.lib
      NM[cl] = numerator(leadcoef(p));
      dbprint(p,"numerator:");
      dbprint(p, string(NM[cl]));
      cNu++; Numero[cNu]= NM[cl];
      cl++;
      NM[cl] = num; // denominator
      dbprint(p,"denominator:");
      dbprint(p, string(NM[cl]));
      cDe++; Denomiro[cDe]= NM[cl];
      cl++;
      p = p - lead(p); // for the next cycle
    }      
    if ( p!= 0)
    {
      num = content(p);
      p   = p/num;
      NM[cl] = denominator(num);
      dbprint(p,"content denominator:");
      dbprint(p, string(NM[cl]));
      cNu++; Numero[cNu]= NM[cl];
      cl++;
      NM[cl] = numerator(num);
      dbprint(p,"content numerator:");
      dbprint(p, string(NM[cl]));
      cDe++; Denomiro[cDe]= NM[cl];
      cl++;
    }
    // it seems that the next elements will not have real influence
    while( p != 0)
    {
      NM[cl] = leadcoef(p); // should be all integer, i.e. non-rational
      dbprint(p,"coef:");
      dbprint(p, string(NM[cl]));
      cl++;
      p = p - lead(p);
    }
  }
  NM = simplify(NM,4); // delete identical
  string newvars = parstr(basering);
  def save = basering;
  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
  execute(NewRing);
  // get params as variables
  // create a list of non-monomials
  ideal @L;
  ideal F;
  ideal NM = imap(save,NM);
  NM = simplify(NM,8); //delete multiples
  poly p,q;
  cl = 1;
  int j, cf;
  for(i=1; i<=size(NM);i++)
  {
    p = NM[i] - lead(NM[i]);
    if (p!=0)
    {
      //      L[cl] = p;
      F  = factorize(NM[i],1); //non-constant factors only
      cf = 1;
      // factorize every polynomial
      // throw away every monomial from factorization (also constants from above ring)
      for (j=1; j<=size(F);j++)
      {
	q = F[j]-lead(F[j]);
	if (q!=0)
	{
	  @L[cl] = F[j];
	  cl++;
	}
      }
    }
  }
  // return the result [in string-format]
  @L = simplify(@L,2+4+8); // skip zeroes, doubled and entries, diff. by a constant
  list SL;
  for (j=1; j<=size(@L);j++)
  {
    SL[j] = string(@L[j]);
  }
  setring save;
  return(SL);
}
//---------------------------------------------------------------
static proc evas_genericity(matrix M)
{
  // called from the main genericity proc
  ideal I = ideal(M);
  I = simplify(I,2+4);
  int s = size(I);
  ideal Den;
  poly p;
  int i;
  for (i=1; i<=s; i++)
  {
    p = I[i];
    while (p !=0) 
    { 
      Den = Den, denominator(leadcoef(p)); 
      p   = p-lead(p); 
    }
  }
  Den = simplify(Den,2+4);  
  string newvars = parstr(basering);
  def save = basering;
  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
  execute(NewRing);
  ideal F;
  ideal Den = imap(save,Den);
  Den = simplify(Den,2);
  int s1 = size(Den);
  for (i=1; i<=s1; i++)
  { 
    if (Den[i] !=1)
    { 
      F= F, factorize(Den[i],1); 
    } 
  }
  F = simplify(F, 2+4+8); 
  ideal @L = F;
  list SL;
  int c,j;
  string Mono;
  c = 1;
  for (j=1; j<=size(@L);j++)
  {
    if (leadcoef(@L[j]) <0)
    {
      @L[j] = -1*@L[j];
    }
    if ( (@L[j] - lead(@L[j]))==0 ) //@L[j] is a monomial
    {
      Mono = Mono + string(@L[j])+ ","; // concatenation
    }
    else
    {
      c++;
      SL[c] = string(@L[j]);
    }
  }
  if (Mono!="")
  {
    Mono = Mono[1..size(Mono)-1]; // delete the last semicolon
  }
  SL[1] = Mono;
  setring save;
  return(SL);
}

//---------------------------------------------------------------
proc canonize(list L)
"USAGE:  canonize(L), L a list
PURPOSE: modules in the list are canonized by computing their reduced minimal (= unique up to constant factor w.r.t. the given ordering) Groebner bases
RETURN:  list 
ASSUME:  L is the output of control/autonomy procedures
EXAMPLE:  example canonize; shows an example
"
{
  list M = L;
  intvec v=Opt_Our();
  int s = size(L);
  int i;
  for (i=2; i<=s; i=i+2)
  {
    if (typeof(M[i])=="module")
    {
      M[i] = std(M[i]);
      //      M[i] = prune(M[i]); // mimimal embedding: no need yet
      //      M[i] = std(M[i]);
    }
  }
  option(set, v); //set old values back
  return(M);
}
example
{  // TwoPendula with L1=L2=L
  "EXAMPLE:"; echo = 2;
  ring r=(0,m1,m2,M,g,L),Dt,dp;
  module RR = 
    [m1*L*Dt^2, m2*L*Dt^2, -1, (M+m1+m2)*Dt^2], 
    [m1*L^2*Dt^2-m1*L*g, 0, 0, m1*L*Dt^2],
    [0, m2*L^2*Dt^2-m2*L*g, 0, m2*L*Dt^2];
  module R = transpose(RR);
  list C = control(R);
  list CC = canonize(C);
  view(CC);
}

//----------------------------------------------------------------

static proc elementof (int i, intvec v)
{
  int b=0;
  for(int j=1;j<=nrows(v);j++)
    {
      if(v[j]==i)
	{
	  b=1;
	  return (b);
	}
    }
  return (b);
}
//-----------------------------------------------------------------
proc iostruct(module R)
"USAGE: iostruct( R ); R a module 
RETURN:  list L with entries: string s, intvec v, module P and module Q
PURPOSE:  if R is the kernel-representation-matrix of some system, then we output a input-ouput representation Py=Qu of the system, the components that have been chosen as outputs(intvec v) and a comment s
NOTE:  the procedure uses Bareiss algorithm which might not terminate in some cases
EXAMPLE:  example iostruct; shows an example
"
{
  list L = bareiss(R);
  int R_rank = ncols(L[1]);
  int NCols=ncols(R);
  intvec v=L[2];
  int temp;
  int NRows=nrows(v);
  int i,j;
  int b=1;
  module P;
  module Q;
  int n=0;
  
  while(b==1)               //sort v through bubblesort
    {
      b=0;
      for(i=1;i<NRows;i++)
	{
	  if(v[i]>v[i+1])
	  {
	    temp=v[i];
	    v[i]=v[i+1];
	    v[i+1]=temp;
	    b=1;
	  }
	}
    }
  P=R[v];                     //generate P 
  for(i=1;i<=NCols;i++)       //generate Q
    {
      if(elementof(i,v)==1)
	{
	  i++;
	  continue;
	}
      Q=Q,R[i];
    }
  Q=simplify(Q,2);
  string s="The following components have been chosen as outputs: ";
  return (list(s,v,P,Q));
}
example
{"EXAMPLE:";echo = 2;
 //Example Antenna
 ring r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), (c,dp);
 
 module RR;
 RR = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
   [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
   [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
   [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
   [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
   [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];
 module R = transpose(RR);
 view(R);
 view(iostruct(R));
};	

//---------------------------------------------------------------
static proc smdeg(matrix N)  
// returns an intvec of length 2 with the index of an element of N with smallest degree
{
  int n = nrows(N);
  int m = ncols(N);
  int d,d_temp;
  intvec v;
  int i,j;            // counter

  if (N==0) 
  {
    v = 1,1;
    return(v);
  } 
  
  for (i=1; i<=n; i++) 
// hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet
  {
    for (j=1; j<=m; j++)
    {
      if( deg(N[i,j])!=-1 )
      {
	d=deg(N[i,j]);
	break;
      }
    }
    if (d != -1)
    {
      break;
    }  
  }
  for(i=1; i<=n; i++)
  {
    for(j=1; j<=m; j++)
    {
      d_temp = deg(N[i,j]);
      if ( (d_temp < d) && (N[i,j]!=0) )
      {
	d=d_temp;
      }
    }
  }
  for (i=1; i<=n; i++)
  {
    for (j=1; j<=m;j++)
    {
      if ( (deg(N[i,j]) == d) && (N[i,j]!=0) )
      {
	v = i,j;
	return(v);
      }
    }
  }
}
//---------------------------------------------------------------
static proc NoNon0Pol(vector v)
// returns 1, if there is only one non-zero element in v and 0 else
{
  int i,j;
  int n = nrows(v);
  for( j=1; j<=n; j++)
  {
    if (v[j] != 0)
    {
      i++;
    }
  }
  if ( i!=1 )
  {
    i=0;
  }
  return(i);
}
//---------------------------------------------------------------
static proc extgcd_Our(poly p, poly q)
{
  ideal J;   //for extgcd-computations
  matrix T; //----------"------------
  list L;
  // the extgcd-command has a bug in versions before 2-0-7
  if ( system("version")<=2006 )
  {
    J = p,q; // J = N[k-1,k-1],N[k,k]; //J is of type ideal
    L[1] = liftstd(J,T);  //T is of type matrix
    if(J[1]==p) //this is just for the case the SINGULAR swaps the 
    //      two elements due to ordering
    {  
      L[2] = T[1,1]; 
      L[3] = T[2,1]; 
    }
    else
    { 
      L[2] = T[2,1]; 
      L[3] = T[1,1]; 
    }
  }
  else
  {
    L=extgcd(p,q); 
    //    L=extgcd(N[k-1,k-1],N[k,k]);  
    //one can use this line if extgcd-bug is fixed
  }
  return(L);
}
static proc normalize_Our(matrix N, matrix Q)
"USAGE: normalize_Our(N,Q), N, Q are two matrices
PURPOSE: normalizes N and divides the columns of Q through the leading coefficients of the columns of N
RETURN: normalized matrix N and altered Q(according to the scheme mentioned in purpose). If number of columns of N and Q do not coincide, N and Q are returned unchanged
NOTE: number of columns of N and Q must coincide. 
"
{
  if(ncols(N) != ncols(Q))
    {
      return (N,Q);
    }	
  module M = module(N);
  module S = module(Q);
  int NCols = ncols(N);
  number n;
  for(int i=1;i<=NCols;i++)
    {
      n = leadcoef(M[i]);
      if( n != 0 )
        {
          M[i]=M[i]/n;
          S[i]=S[i]/n;
        }
     }
   N = matrix(M);
   Q = matrix(S);
   return (N,Q);
} 		
	
//---------------------------------------------------------------
proc smith( module M )
"USAGE: smith(M), M a module or a matrix,
PURPOSE: computes the Smith form of a matrix
RETURN: a list of length 4 with the following entries:
@*      [1]: The Smith-Form S of M, 
@*      [2]: the rank of M, 
@*      [3]: a unimodular matrix U,
@*      [4]: a unimodular matrix V, 
such that U*M*V=S. An warning is returned when no Smith Form exists. 
NOTE: The Smith form only exists over PIDs (principal ideal domains). Use global ordering for computations!
"
{
  if (nvars(basering)>1) //if more than one variable, return empty list
  {
    string s="The Smith-Form only exists for principal ideal domains";
    return (s);
  }
  matrix N = matrix(M);         //Typecasting
  int n = nrows(N);
  int m = ncols(N);
  matrix P = unitmat(n);       //left transformation matrix
  matrix Q = unitmat(m);       //right transformation matrix
  int k, i, j, deg_temp;
  poly tmp;
  vector v;
  list L;                      //for extgcd-computation
  intmat f[1][n];              //to save degrees
  matrix lambda[1][n];         //to save leadcoefficients
  intmat g[1][m];              //to save degrees
  matrix mu[1][m];             //to save leadcoefficients
  int ii;                       //counter
  
  while ((k!=n) && (k!=m) )    
  {
    k++;
    while ((k<=n) && (k<=m))  //outer while-loop for column-operations
    {
      while(k<=m )        //inner while-loop for row-operations
      {
	if( (n>m) && (k < n) && (k<m))
	{      
	  if( simplify((ideal(submat(N,k+1..n,k+1..m))),2)== 0)
	  {
	    return(N,k-1,P,Q);
	  }
	}
        i,j = smdeg(submat(N,k..n,k..m)); //choose smallest degree in the remaining submatrix
        i=i+(k-1);                        //indices adjusted to the whole matrix
	j=j+(k-1);
	if(i!=k)                    //take the element with smallest degree in the first position
        {
          N=permrow(N,i,k);
          P=permrow(P,i,k);
        }
        if(j!=k)
        {
          N=permcol(N,j,k);
          Q=permcol(Q,j,k);
        }
        if(NoNon0Pol(N[k])==1)
        {
          break;
        }
        tmp=leadcoef(N[k,k]);
        deg_temp=ord(N[k,k]);             //ord outputs the leading degree of N[k,k]
        for(ii=k+1;ii<=n;ii++)
        {
          lambda[1,ii]=leadcoef(N[ii,k])/tmp;     
          f[1,ii]=deg(N[ii,k])-deg_temp;
        }
        for(ii=k+1;ii<=n;ii++)
        {
          N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii);	         
          P = addrow(P,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
	  N,Q=normalize_Our(N,Q);
	}
      }
      if (k>n)
      {
	break;
      }
      if(NoNon0Pol(transpose(N)[k])==1)
      {
        break;
      }
      tmp=leadcoef(N[k,k]);
      deg_temp=ord(N[k,k]); //ord outputs the leading degree of N[k][k]
      
      for(ii=k+1;ii<=m;ii++)
      {
        mu[1,ii]=leadcoef(N[k,ii])/tmp;     
        g[1,ii]=deg(N[k,ii])-deg_temp;
      }
      for(ii=k+1;ii<=m;ii++)
      {
        N=addcol(N,k,-mu[1,ii]*var(1)^g[1,ii],ii);
        Q=addcol(Q,k,-mu[1,ii]*var(1)^g[1,ii],ii);
	N,Q=normalize_Our(N,Q);
      }
    }
    if( (k!=1) && (k<n) && (k<m) )
    {
      L = extgcd_Our(N[k-1,k-1],N[k,k]);      
      if ( N[k-1,k-1]!=L[1] )  //means that N[k-1,k-1] is not a divisor of N[k,k]
      {
	N=addrow(N,k-1,L[2],k);
	P=addrow(P,k-1,L[2],k);
        N,Q=normalize_Our(N,Q);	
	
        N=addcol(N,k,-L[3],k-1);
	Q=addcol(Q,k,-L[3],k-1);
	N,Q=normalize_Our(N,Q);
	k=k-2;
      }
    }
  }  
  if( (k<=n) && (k<=m) )
  {
    if( N[k,k]==0)
    {
      return(N,k-1,P,Q);
    }
  }
  return(N,k,P,Q);
}
example
{
  "EXAMPLE:";echo = 2;
  option(redSB);
  option(redTail);
  ring r=0,x,dp;
  // see what happens when the matrix is already in Smith-Form
  module M = [x,0,0],[0,x2,0],[0,0,x3];
  print(M);
  list L = smith(M);
  print(L[1]);
  matrix N=matrix(M);
  matrix B=L[3]*N*L[4];
  print(B);
  //------- and yet another example --------------
  module M2=[x2,x,3x3-4],[2x2-1,4x,5x2],[2x5,3x,4x];
  print(M2);
  list P=smith(M2);
  print(P[1]);
  matrix N2=matrix(M2);
  matrix B2=P[3]*N2*P[4];
  print(B2);
}
//---------------------------------------------------------------
proc list_tex(L, string name,link l,int nr_loop)
"USAGE: list_tex(L,name,l), where L is a list, name a string, l a link
	 writes the content of list L in a tex-file 'name'
RETURN: nothing
"
{
  if(typeof(L)!="list")  //in case L is not a list
  {
    texobj(name,L);
  }
  if(size(L)==0)  
  {
  }
  else
  {
    string t;
    for (int i=1;i<=size(L);i++)
    {
      while(1)
      {
	if(typeof(L[i])=="string")  //Fehler hier fuer normalen output->nur wenn string in liste dann verbatim
	{
	  t=L[i];
	  if(nr_loop==1)
	  {
	    write(l,"\\begin\{center\}");
	    write(l,"\\begin\{verbatim\}");
	  }
	  write(l,t);
	  if(nr_loop==0)
	  {
	    write(l,"\\par");
	  }
	  if(nr_loop==1)
	  {
	    write(l,"\\end\{verbatim\}");
	    write(l,"\\end\{center\}");
	  }
	  break;
	}
	if(typeof(L[i])=="module")
	{
	  texobj(name,matrix(L[i]));
	  break;
	}
	if(typeof(L[i])=="list")
	{
	  list_tex(L[i],name,l,1);
	  break;
	}
	write(l,"\\begin\{center\}");
	texobj(name,L[i]);
	write(l,"\\end\{center\}");
	write(l,"\\par");
	break;
      }
    }
  }
}
//---------------------------------------------------------------
proc verbatim_tex(string s, link l)
"USAGE: verbatim_tex(s,l), where s is a string and l a link
PURPOSE: writes the content of s in verbatim-environment in the file
	 specified by link
RETURN: nothing
"
{
  write(l,"\\begin{verbatim}");
  write(l,s);
  write(l,"\\end{verbatim}");
  write(l,"\\par");
}
//---------------------------------------------------------------
proc FindTorsion(module R, ideal TAnn)
"USAGE:  FindTorsion(R, I);   R an ideal/matrix/module, I an ideal
PURPOSE: computes the Groebner basis of the submodule of R, annihilated by I
ETURN:  module
NOTE: especially helpful, when I is the annihilator of the t(R) - the torsion submodule of R. In this case, the result is the explicit presentation of t(R) as
the submodule of R
EXAMPLE: example FindTorsion; shows an example
"
{
  // motivation: let R be a module,
  // TAnn is the annihilator of t(R)\subset R
  // compute the generators of t(R) explicitly
  ideal AS = TAnn;
  module S = R;
  if (attrib(S,"isSB")<>1)
  {
    S = std(S);
  }
  if (attrib(AS,"isSB")<>1)
  {
    AS = std(AS);
  }
  int nc  = ncols(S);
  module To = quotient(S,AS);
  To = std(NF(To,S));
  return(To);
}
example
{
  "EXAMPLE:";echo = 2;
  // Flexible Rod
  ring A = 0,(D1, D2), (c,dp);
  module R= [D1, -D1*D2, -1], [2*D1*D2, -D1-D1*D2^2, 0];
  module RR = transpose(R);
  list L = control(RR);
  // here, we have the annihilator:
  ideal LAnn = D1; // = L[10]
  module Tr  = FindTorsion(RR,LAnn);
  print(RR);  // the module itself
  print(Tr); // generators of the torsion submodule
}