/////////////////////////////////////////////////////////////////////////////
//procedures examples comments
version="$Id$";
category="Noncommutative";
info="
LIBRARY: purityfiltration.lib     Algorithms for computing a purity filtration of a given module

AUTHORS: Christian Schilli,        christian.schilli@rwth-aachen.de
@*          Viktor Levandovskyy,   levandov@math.rwth-aachen.de


OVERVIEW:
Purity is a notion with several meanings. In our context it is equidimensionality
@* of a module (that is all M is pure iff any nonzero submodule of N has the same dimension as N).
@* Notably, one should define purity with respect to a given dimension function. In the context
@* of this library the corresponding function is the homological grade number j_A(M) of a module M over
@* an K-algebra A. j_A(M) is the minimal integer k, such that Ext^k_A(M,A) != 0.

REFERENCES: [AQ] Alban Quadrat: Grade filtration of linear functional systems, INRIA Report 7769 (2010), to appear in Acta Applicanda Mathematica.
@* [B93] Jan-Erik Bjoerk: Analytic D-modules and applications, Kluwer Acad. Publ., 1993.
@* [MB10] Mohamed Barakat: Purity Filtration and the Fine Structure of Autonomy. Proc. MTNS, 2010.

PROCEDURES:
projectiveDimension(matrix T,int i);        compute a shortest resolution of coker(T) and its projective dimension
purityFiltration(matrix R);                    compute the purity filtration of coker(R)
purityTriang(matrix R)                        compute a triangular blockmatrix T, such that coker(R) isomorphic to coker(T)
gradeNumber(matrix R);                              gives the grade number of the module coker(R)
showgrades(list T);                           gives all grade numbers of the modules represented by the elements of T
allExtOfLeft(matrix R);                              computes all right ext-modules ext^i(M,D) of a left module M=coker(R) over the ring D
allExtOfRight(matrix R);                      computes all left ext-modules ext^i(M,D) of a right module M=coker(R) over the ring D
doubleExt(matrix R, int i);                 computes the left module ext^i(ext^i(M,D),D) over the ring D, M=coker(R)
allDoubleExt(matrix R);                              computes all double ext modules ext^i(ext^j(M,D),D) of the left module coker(R) over the ring D
is_pure(matrix R);                              checks whether the module coker(R) is pure
purelist(list T);                              checks whether all the modules represented by the elements of T are pure

KEYWORDS: D-module; ext-module; filtration; projective dimension; resolution; purity
";

LIB "nctools.lib";
LIB "matrix.lib";
LIB "poly.lib";
LIB "general.lib";
LIB "control.lib";
LIB "nchomolog.lib";

//------------------- auxiliary procedures --------------------------

proc testPurityfiltrationLib()
{
example projectiveDimension;
example purityFiltration;
example purityTriang;
example gradeNumber;
example showgrades;
example allExtOfLeft;
example allExtOfRight;
example doubleExt;
example allDoubleExt;
example is_pure;
example purelist;
}

static proc iszero (matrix R)
"USAGE:  iszero(R); R a matrix
RETURN:  int, 1, if R is zero,
@*              or 0, if it's not
PURPOSE: checks, if the matrix R is zero or not
"
{
  ideal i=R;
  i=std(i);
  if (i==0)
  {
    return (1);
  }
  return (0);
}

proc lsyz (matrix R)
"USAGE:  lsyz(R), R a matrix
RETURN:         matrix, a left syzygy of R
PURPOSE: computes the left syzygy module of the module, generated by the rows of R, i.e.
@*         a matrix X with X*R=0
"
{
  matrix L=transpose(syz(transpose(R)));
  return(L);
}

proc rsyz (matrix R)
"USAGE:  rsyz(R), R a matrix
RETURN:         matrix, a rightsyzygy of R
PURPOSE: computes the right syzygy module of the module, generated by the rows of R, i.e.
@*         a matrix X with R*X=0
EXAMPLE: example rsyz; shows example
"
{
  def save = basering;            // with respect to non-commutative rings,
  def saveop = opposite(save);    // we have to switch to the oppose ring for a rightsyzygy
  setring saveop;
  matrix Rop = oppose(save,R);
  matrix Bop = syz(Rop);
  setring save;
  matrix B =oppose(saveop,Bop);
  kill saveop;
  return(B);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x,y,z),dp;
  matrix R[3][2]=x,0,0,x,y,-z;
  matrix X=rsyz(R);
  print(X);
  // check
  print(R*X);
}


static proc rinv (matrix R)
"USAGE:  rinv(R), R a matrix
RETURN:         matrix, a right inverse of R
PURPOSE: computes a right inverse matrix of R, if it exists
@*         if not, it returns the zero matrix
"
{
  return(rightInverse(R));
}

static proc linv (matrix R)
"USAGE:  linv(R), R a matrix
RETURN:         matrix, a left inverse of R
PURPOSE: computes a left inverse matrix of R, if it exists
@*         if not, it returns the zero matrix
"
{
  return (leftInverse(R));
}

proc rlift(matrix M, matrix N)
"USAGE:  rlift(M,N), M and N matrices, so that the module, generated by the columns of N
@*         is a submodule of the one, generated by the columns of M
RETURN:         matrix, a right lift of N in M
PURPOSE: computes a right lift matrix X of N in M,
@*         i.e. N=M*X
"
{
  def save = basering;           // with respect to non-commutative rings,
  def saveop = opposite(save);   // we have to change the ring for a rightlift
  setring saveop;
  matrix Mop = oppose(save,M);
  matrix Nop = oppose(save,N);
  matrix Bop = lift(Mop,Nop);
  setring save;
  matrix B =oppose(saveop,Bop);
  kill saveop;
  return(B);
}

proc llift(matrix M, matrix N)
"USAGE:  llift(M,N), M and N matrices, so that the module, generated by the rows of N
@*         is a submodule of the one, generated by the rows of M
RETURN:         matrix, a left lift of N in M
PURPOSE: computes a left lift matrix X of N in M,
@*         i.e. N=X*M
"
{
  matrix X=transpose(lift(transpose(M),transpose(N)));
  return(X);
}

static proc concatz(matrix M, matrix N)
"USAGE:  concatz(M,N), M and N matrices
RETURN:         matrix
PURPOSE: adds the rows of N under the rows of M, i.e. build the matrix (M^Tr,N^Tr)^Tr
"
{
  matrix X=transpose(concat(transpose(M),transpose(N)));
  return (X);
}

//------------------------- main procedures --------------------------

proc purityFiltration(matrix R)
"USAGE:  purityFiltration(S), S matrix with entries of an Auslander regular ring D
RETURN:         a list T of two lists, purity filtration of the module M=D^q/D^p(S^t)
PURPOSE: the first list T[1] gives a filtration {M_i} of M,
@*         where the i-th entry of T[1] gives the representation matrix of M_(i-1).
@*         the second list T[2] gives representations of the factor Modules,
@*         i.e. T[2][i] gives the repr. matrix for M_(i-1)/M_i
EXAMPLE: example purityFiltration; shows example
"
{
  int i,j;
  list re=projectiveDimension(R,0);
  list T=re[1];
  int di=re[2];
  list reres;             // Rji=reres[i][j+1], i=1,..,n+1; j=0,..,i
  for( i=1; i<=di+1; i++ )
  {
    list zw;
    zw[i+1]=T[i];
    for( j=i; j >= 1; j--)
    {
      zw[j]=rsyz(zw[j+1]);
    }
    reres[i]=zw;
    kill zw;
  }
  list F;    // Fij=F[j][i+1], j=2,..,n+1; i=0,..,j-1
  for(i=2;i<=di;i++)
  {
    list ehm;
    matrix I[nrows(T[i-1])][nrows(T[i-1])];
    I=I+1;
    ehm[i]=I;
    kill I;
    for (j=1; j<=i-1; j++)
    {
      ehm[i-j]=rlift(reres[i][i-j+1],ehm[i-j+1]*reres[i-1][i-j+1]);
    }
    F[i]=ehm;
    kill ehm;
  }
//  list M;       // Mi=M[i+1], i=0,...,n+1
//  M[1]=R1;
//  matrix Ti=lsyz(reres[1][1]);
//  matrix P[ncols(Ti)][ncols(Ti)];
//  P=P+1;
//  for (i=1;i<=di; i++)
//  {
//    M[i+1]=transpose(modulo(transpose(Ti*P),transpose(reres[i][2])));
//    P=F[i+1][1]*P;
//    Ti=lsyz(reres[i+1][1]);
//  }
//  M[di+2]=transpose(modulo(transpose(Ti*P),transpose(reres[di+1][2])));
//  list I;
//  for (i=1;i<=di+1;i++)
//  {
//    I[i]=transpose(modulo(transpose(M[i]),transpose(M[i+1])));
//  }
  list Rs,Rss;
  for(i=1; i<=di; i++)
  {
    list zw;
    zw[1]=lsyz(reres[i][1]);
    zw[2]=lsyz(zw[1]);
    Rss[i]=llift(zw[1],reres[i][2]);
    Rs[i]=zw;
    kill zw;
  }
  list Fs;
  for(i=2;i<=di;i++)
  {
    Fs[i]=llift(Rs[i-1][1],Rs[i][1]*F[i][1]);
  }
  list K,U;
  K[1]=transpose(R);
  U[1]=Rs[1][1];
  for(i=2;i<=di;i++)
  {
    K[i]=transpose(std(transpose(concatz(Rss[i-1], Rs[i-1][2]))));
    U[i]=transpose(std(transpose(concatz(concatz(Fs[i],Rss[i-1]),Rs[i-1][2]))));
  }
  K[di+1]=transpose(std(transpose(concatz(Rss[di], Rs[di][2]))));
  U[di+1]=K[di+1];
  list erg=(K,U);
  return (erg);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x1,x2,d1,d2),dp;
  def S=Weyl();
  setring S;
  int i;
  matrix R[3][3]=0,d2-d1,d2-d1,d2,-d1,-d1-d2,d1,-d1,-2*d1;
  print(R);
  list T=purityFiltration(transpose(R));
  // the purity filtration of coker(M)
  print(T[1][1]);
  print(T[1][2]);
  print(T[1][3]);
  // factor modules of the filtration
  print(T[2][1]);
  print(T[2][2]);
  print(T[2][3]);
}


proc purityTriang(matrix R)
"USAGE:  purityTriang(S), S matrix with entries of an Auslander regular ring D
RETURN:         a matrix T
PURPOSE: compute a triangular block matrix T, such that M=D^p/D^q(S^t) is isomorphic to M'=D^p'/D^q(T^t)
EXAMPLE: example purityTriang; shows example
"
{
  int i,j;
  list re=projectiveDimension(R,0);
  list T=re[1];
  int di=re[2];
  list reres;             // Rji=reres[i][j+1], i=1,..,n+1; j=0,..,i
  for( i=1; i<=di+1; i++ )
  {
    list zw;
    zw[i+1]=T[i];
    for( j=i; j >= 1; j--)
    {
      zw[j]=rsyz(zw[j+1]);
    }
    reres[i]=zw;
    kill zw;
  }
  list F;    // Fij=F[j][i+1], j=2,..,n+1; i=0,..,j-1
  for(i=2;i<=di;i++)
  {
    list ehm;
    matrix I[nrows(T[i-1])][nrows(T[i-1])];
    I=I+1;
    ehm[i]=I;
    kill I;
    for (j=1; j<=i-1; j++)
    {
      ehm[i-j]=rlift(reres[i][i-j+1],ehm[i-j+1]*reres[i-1][i-j+1]);
    }
    F[i]=ehm;
    kill ehm;
  }

  list Rs,Rss;
  for(i=1; i<=di; i++)
  {
    list zw;
    zw[1]=lsyz(reres[i][1]);
    zw[2]=lsyz(zw[1]);
    Rss[i]=llift(zw[1],reres[i][2]);
    Rs[i]=zw;
    kill zw;
  }
  list Fs;
  for(i=2;i<=di;i++)
  {
    Fs[i]=llift(Rs[i-1][1],Rs[i][1]*F[i][1]);
  }


  int sp; list spnr;
  spnr[1]=ncols(Rs[1][1]);
  for (i=2;i<=di;i++)
    {
      spnr[i]=ncols(Fs[i]);
    }
  spnr[di+1]=ncols(Rss[di]);
  sp=sum(spnr);

  matrix E[nrows(Rs[1][1])][nrows(Rs[1][1])]; E=E-1;
  list Z; int sumh;
  Z[1]=concat(Rs[1][1],E);
  sumh=ncols(Rs[1][1]);
  kill E;

  for(i=2;i<=di;i++)
  {
    matrix A;
    matrix B[1][sumh];
    matrix E[nrows(Fs[i])][nrows(Fs[i])]; E=E-1;

    A=Fs[i];
    if (i>2)
    {
      if (iszero(Rss[i-1])==0)
      {
        A=concatz(A,Rss[i-1]);
      }
    }
    if (iszero(Rs[i-1][2])==0)
    {
      A=concatz(A,Rs[i-1][2]);
    }
    A=concat(B,A,E);
    Z[i]=A;
    sumh=sumh+spnr[i];
    kill A,B,E;
  }


  matrix hi,his;
  matrix N[1][sumh];

  if (iszero(Rss[di])==0)
  {
    hi=concat(N,Rss[di]);
  }

  if (iszero(Rs[di][2])==0)
  {
    his=concat(N,Rs[di][2]);
    if (iszero(hi)==1)
    {
      hi=his;
    }
    if (iszero(hi)==0)
    {
    hi=concatz(hi,his);
    }
  }

  kill his;

  matrix ges=Z[1];
  for (i=2;i<=di;i++)
  {
    ges = concatz(ges,Z[i]);
  }

  if (iszero(hi)==0)
  {
    ges=concatz(ges,hi);
  }
return (ges);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x1,x2,d1,d2),dp;
  def S=Weyl();
  setring S;
  int i;
  matrix R[3][3]=0,d2-d1,d2-d1,d2,-d1,-d1-d2,d1,-d1,-2*d1;
  print(R);
  matrix T=purityTriang(transpose(R));
  // a triangular blockmatrix representing the module coker(R)
  print(T);
}


proc gradeNumber(matrix R)
"USAGE:  gradeNumber(R), R matrix, representing M=D^p/D^q(R^t) over a ring D
RETURN:         int, grade number of M
PURPOSE: computes the grade number of M, i.e. the first i, with ext^i(M,D) !=0
@*         returns -1 if M=0
EXAMPLE: example gradeNumber; shows examples
"
{
  matrix M=transpose(R);
  if (is_zero(transpose(M))==1)
  {
    return (-1);
  }
  list ext = allExtOfLeft(transpose(M));
  int i=1;
  matrix L=ext[i];
  while (is_zero(transpose(L))==1)
  {
    i=i+1;
    L=ext[i];
  }
  return (i-1);
}
example
{"EXAMPLE:";echo = 2;
  // trivial example
  ring D=0,(x,y,z),dp;
  matrix R[2][1]=1,x;
  gradeNumber(R);
  // R has left inverse, so M=D/D^2R=0
  gradeNumber(transpose(R));
  print(ncExt_R(0,R));
  // so, ext^0(coker(R),D) =! 0)
  //
  // a little bit more complex
  matrix R1[3][1]=x,-y,z;
  gradeNumber(transpose(R1));
  print(ncExt_R(0,transpose(R1)));
  print(ncExt_R(1,transpose(R1)));
  print(ncExt_R(2,transpose(R1)));
  // ext^i are zero for i=0,1,2
  matrix ext3=ncExt_R(3,transpose(R1));
  print(ext3);
  // not zero
  is_zero(ext3);
}

proc allExtOfLeft(matrix Ps)
"USAGE:  allExtOfLeft(M),
RETURN:         list, entries are ext-modules
ASSUME: M presents a left module of finite left projective dimension n
PURPOSE: For a left module presented by M over the basering D,
@*           compute a list T, whose entry T[i+1] is a matrix, presenting the right module Ext^i_D(M,D) for i=0..n
EXAMPLE: example allExtOfLeft; shows example
"
{
  // old doc: ... T[i] gives the repr. matrix of ext^(i-1)(M,D), i=1,.., n+1
  list ext, Phi;
  ext[1]=ncHom_R(Ps);
  Phi = mres(Ps,0);
  int di = size(Phi);
  Phi[di+1]= transpose(lsyz(transpose(Phi[di])));
  int i;
  def Rbase = basering;
  for(i=1;i<=di;i++)
  {
    module f   = transpose(matrix(Phi[i+1]));
    module Im2 = transpose(matrix(Phi[i]));
    def Rop   = opposite(Rbase);
    setring Rop;
    module fop    = oppose(Rbase,f);
    module Im2op  = oppose(Rbase,Im2);
    module ker_op = modulo(fop,std(0));
    module ext_op = modulo(ker_op,Im2op);
    setring Rbase;
    ext[i+1] = oppose(Rop,ext_op); // a right module!
    kill f, Im2, Rop;
  }
  return(ext);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x,y,z),dp;
  matrix R[6][4]=
  0,-2*x,z-2*y-x,-1,
  0,z-2*x,2*y-3*x,1,
  z,-6*x,-2*y-5*x,-1,
  0,y-x,y-x,0,
  y,-x,-y-x,0,
  x,-x,-2*x,0;
  // coker(R) consider the left module M=D^6/D^4R
  list T=allExtOfLeft(transpose(R));
  print(T[1]);
  print(T[2]);
  print(T[3]);
  print(T[4]);
  // right modules coker(T[i].)!!
}

proc allExtOfRight(matrix Ps)
"USAGE:  allExtOfRight(R), R matrix representing the right Module M=D^q/RD^p over a ring D
@*           M module with finite right projective dimension n
RETURN:         list, entries are ext-modules
PURPOSE: computes a list T, which entries are representations of the left modules ext^i(M,D)
@*         T[i] gives the repr. matrix of ext^(i-1)(M,D), i=1,..,n+1
EXAMPLE: example allExtOfRight; shows example
"
{
  // matrix Ps=transpose(Y);
  list ext, Phi;
  def Rbase = basering;
  def Rop   = opposite(Rbase);
  setring Rop;
  matrix Psop=oppose(Rbase,Ps);
  matrix ext1_op = ncHom_R(Psop);
  setring Rbase;
  ext[1]=oppose(Rop,ext1_op);
  kill Rop;
  list zw = rightreso(transpose(Ps)); // right resolution
  int di = size(zw);
  zw[di+1]=lsyz(zw[di]);
  Phi = zw;
  kill zw;
  int i;
  for(i=1;i<=di;i++)
  {
    module f   = Phi[i+1];
    module Im2 = Phi[i];
    module ker = modulo(f,std(0));
    ext[i+1] = modulo(ker,Im2);  // a left module!
    kill f, Im2, ker;
  }
  return(ext);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x,y,z),dp;
  matrix R[6][4]=
  0,-2*x,z-2*y-x,-1,
  0,z-2*x,2*y-3*x,1,
  z,-6*x,-2*y-5*x,-1,
  0,y-x,y-x,0,
  y,-x,-y-x,0,
  x,-x,-2*x,0;
  // coker(R) considered as right module
  projectiveDimension(R,1)[2];
  list T=allExtOfRight(R);
  print(T[1]);
  print(T[2]);
  // left modules coker(.T[i])!!
}

static proc rightreso(matrix T)
"USAGE:  rightreso(T), T matrix representing the right module M=D*/TD*
RETURN:         list L, a right resolution of M
PURPOSE: computes a right resolution of M, using mres
@*         the i-th entry of L gives the (i-1)th right syzygy module of M
"
{
  int j;
  matrix M=transpose(T);
  list res;
  def save = basering;            // with respect to non-commutative rings,
  def saveop = opposite(save);    // we have to change the ring for a rightresolution
  setring saveop;
  matrix Mop=oppose(save,M);
  list aufl=mres(Mop,0);
  list resop=aufl;
  kill aufl;
  for (j=1; j<=size(resop); j++)
  {
    matrix zw=resop[j];
    setring save;
    res[j]=transpose(oppose(saveop,zw));
    setring saveop;
    kill zw;
  }
  setring save;
  kill saveop;
  return(res);
}

proc showgrades(list T)
"USAGE:  showgrades(T), T list, which includes representation matrices of modules
RETURN:         list, gradenumbers of the entries in T
PURPOSE: computes a list L with L[i]=gradenumber(M), M=D^p/D^qT[i]
EXAMPLE: example showgrades; shows example
"
{
  list grades;
  int gr=size(T);
  int i;
  for (i=1;i<=gr;i++)
  {
    grades[i]=gradeNumber(transpose(T[i]));
  }
  return (grades);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x,y,z),dp;
  matrix R[6][4]=
  0,-2*x,z-2*y-x,-1,
  0,z-2*x,2*y-3*x,1,
  z,-6*x,-2*y-5*x,-1,
  0,y-x,y-x,0,
  y,-x,-y-x,0,
  x,-x,-2*x,0;
  list T=purityFiltration(transpose(R))[2];
  showgrades(T);
  // T[i] are i-1 pure (i=1,3,4) or zero (i=2)
}

proc doubleExt(matrix R, int i)
"USAGE:  doubleExt(R,i), R matrix representing the left Module M=D^p/D^q(R^t) over a ring D
@*         int i, less or equal the left projective dimension of M
RETURN:         matrix P, representing the double ext module
PURPOSE: computes a matrix P, which represents the left module ext^i(ext^i(M,D))
EXAMPLE: example doubleExt; shows example
"
{
  return (allExtOfRight(   allExtOfLeft(R)[i+1]   )[i+1]);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x,y,z),dp;
  matrix R[7][3]=
  0 ,0,1,
  1 ,-4*x+z,-z,
  -1,8*x-2*z,z,
  1 ,0  ,0,
  0 ,x-y,0,
  0 ,x-y,y,
  0 ,0  ,x;
  // coker(R) is 2-pure, so all doubleExt are zero
  print(doubleExt(transpose(R),0));
  print(doubleExt(transpose(R),1));
  print(doubleExt(transpose(R),3));
  // except of the second
  print(doubleExt(transpose(R),2));
}

proc allDoubleExt(matrix R)
"USAGE:  allDoubleExt(R), R matrix representing the left Module M=D^p/D^q(R^t) over a ring D
RETURN:         list T, double indexed, which include all double-ext modules
PURPOSE: computes all double ext-modules
@*         T[i][j] gives a representation matrix of ext^(j-1)(ext(i-1)(M,D))
EXAMPLE: example allDoubleExt; shows example
"
{
list ext=allExtOfLeft(transpose(R));
list extext;
int i;
for(i=1;i<=size(ext);i++)
  {
    extext[i]=allExtOfRight(ext[i]);
  }
kill ext;
return (extext);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x1,x2,x3,d1,d2,d3),dp;
  def S=Weyl();
  setring S;
  matrix R[6][4]=
  0,-2*d1,d3-2*d2-d1,-1,
  0,d3-2*d1,2*d2-3*d1,1,
  d3,-6*d1,-2*d2-5*d1,-1,
  0,d2-d1,d2-d1,0,
  d2,-d1,-d2-d1,0,
  d1,-d1,-2*d1,0;
  list T=allDoubleExt(transpose(R));
  // left projective dimension of M=coker(R) is 3
  // ext^i(ext^0(M,D)), i=0,1,2,3
  print(T[1][1]);
  print(T[1][2]);
  print(T[1][3]);
  print(T[1][4]);
  // ext^i(ext^1(M,D)), i=0,1,2,3
  print(T[2][1]);
  print(T[2][2]);
  print(T[2][3]);
  print(T[2][4]);
  // ext^i(ext^2(M,D)), i=0,1,2,3  (all zero)
  print(T[3][1]);
  print(T[3][2]);
  print(T[3][3]);
  print(T[3][4]);
  // ext^i(ext^3(M,D)), i=0,1,2,3  (all zero)
  print(T[4][1]);
  print(T[4][2]);
  print(T[4][3]);
  print(T[4][4]);
}

proc is_pure(matrix R)
"USAGE:  is_pure(R), R representing the module M=D^p/D^q(R^t)
RETURN:         int, 0 or 1
PURPOSE: checks pureness of M.
@*         returns 1, if M is pure, or 0, if it's not
@*         remark: if M is zero, is_pure returns 1
EXAMPLE: example is_pure; shows example
"
{
  matrix M=transpose(R);
  int gr=gradeNumber(transpose(M));
  int di=projectiveDimension(transpose(M),0)[2];
  int i=0;
  while(i<=di)
  {
    if (i!=gr)
    {
      if (  is_zero( doubleExt(transpose(M),i) ) == 0 )
      {
        return (0);
      }
    }
    i=i+1;
  }
  return (1);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x,y,z),dp;
  matrix R[3][2]=y,-z,x,0,0,x;
  list T=purityFiltration(transpose(R));
  print(transpose(std(transpose(T[2][2]))));
  // so the purity filtration of coker(R) is trivial,
  // i.e. coker(R) is already pure
  is_pure(transpose(R));
  // we can also have non-pure modules:
  matrix R2[6][4]=
  0,-2*x,z-2*y-x,-1,
  0,z-2*x,2*y-3*x,1,
  z,-6*x,-2*y-5*x,-1,
  0,y-x,y-x,0,
  y,-x,-y-x,0,
  x,-x,-2*x,0;
  is_pure(transpose(R2));
}

proc purelist(list T)
"USAGE:  purelist(T), T list, in which the i-th entry R=T[i] represents M=D^p/D^q(R^t)
RETURN:         list M, entries of M are 0 or 1
PURPOSE: if T[i] is pure, M[i] is 1, else M[i] is 0
EXAMPLE: example purelist; shows example
"
{
  int i;
  list erg;
  for(i=1;i<=size(T);i++)
  {
    erg[i]=is_pure(transpose(T[i]));
  }
  return (erg);
}
example
{"EXAMPLE:";echo = 2;
  ring D = 0,(x,y,z),dp;
  matrix R[6][4]=
  0,-2*x,z-2*y-x,-1,
  0,z-2*x,2*y-3*x,1,
  z,-6*x,-2*y-5*x,-1,
  0,y-x,y-x,0,
  y,-x,-y-x,0,
  x,-x,-2*x,0;
  is_pure(transpose(R));
  // R is not pure, so we do the purity filtration
  list T=purityFiltration(transpose(R));
  // all Elements of T[2] are either zero or pure
  purelist(T[2]);
}


proc projectiveDimension(matrix T, list #)
"USAGE:  projectiveDimension(R,i,j), R matrix representing the Modul M=coker(R)
@*       int i, with i=0 or i=1, j a natural number
RETURN:  list T, a projective resolution of M and its projective dimension
PURPOSE: if i=0 (and by default), T[1] gives a shortest left resolution of M=D^p/D^q(R^t) and T[2] the left projective dimension of M
@*         if i=1, T[1] gives a shortest right resolution of M=D^p/RD^q and T[2] the right projective dimension of M
@*          in both cases T[1][j] is the (j-1)-th syzygy module of M
NOTE: The algorithm is due to A. Quadrat, D. Robertz, Computation of bases of free modules over the Weyl algebras, J.Symb.Comp. 42, 2007.
EXAMPLE: example projectiveDimension; shows examples
"
{
  int i = 0; // default
  if (size(#) >0)
  {
    i = int(#[1]);
    if ( (i!=0) and (i!=1) )
    {
      printf("Unaccepted second argument. Use 0 to get a left resolution, 1 for a right one.");
    }
  }
  if (i==0)
  {
    return(prodim(T));
  }
    int j;
    matrix M=T;
    list res;
    def save = basering;            // with respect to non-commutative rings,
    def saveop = opposite(save);    // we have to change the ring for a rightresolution
    setring saveop;
    matrix Mop=oppose(save,M);
    list aufl=prodim(Mop);
    int k=aufl[2];
    list resop=aufl[1];
    kill aufl;
    for (j=1; j<=size(resop); j++)
    {
      matrix zw=resop[j];
      setring save;
      res[j]=transpose(oppose(saveop,zw));
      setring saveop;
      kill zw;
    }
    setring save;
    list Y;
    Y[1]=res;
    Y[2]=k;
    kill saveop;
    kill res;
    return(Y);

}
example
{"EXAMPLE:";echo = 2;
  // commutative example
  ring D = 0,(x,y,z),dp;
  matrix R[6][4]=
  0,-2*x,z-2*y-x,-1,
  0,z-2*x,2*y-3*x,1,
  z,-6*x,-2*y-5*x,-1,
  0,y-x,y-x,0,
  y,-x,-y-x,0,
  x,-x,-2*x,0;
  // compute a left resolution of M=D^4/D^6*R
  list T=projectiveDimension(transpose(R),0);
  // so we have the left projective dimension
  T[2];
  //we could also compute a right resolution of M=D^6/RD^4
  list T1=projectiveDimension(R,1);
  // and we have right projective dimension
  T1[2];
  // check, that a syzygy matrix of R has left inverse:
  print(leftInverse(syz(R)));
  // so lpd(M) must be 1.
  // Non-commutative example
  ring D1 = 0,(x1,x2,x3,d1,d2,d3),dp;
  def S=Weyl();  setring S;
  matrix R[3][3]=
  1/2*x2*d1, x2*d2+1, x2*d3+1/2*d1,
  -1/2*x2*d2-3/2,0,1/2*d2,
  -d1-1/2*x2*d3,-d2,-1/2*d3;
  list T=projectiveDimension(R,0);
  // left projective dimension of coker(R) is
  T[2];
  list T1=projectiveDimension(R,1);
  // both modules have the same projective dimension, but different resolutions, because D is non-commutative
  print(T[1][1]);
  // not the same as
  print(transpose(T1[1][1]));
}

static proc prodim(matrix M)
"USAGE:  prodim(R), R matrix representing the Modul M=coker(R)
RETURN:  list T, a left projective resolution of M and its left projective dimension
PURPOSE: T[1] gives a shortest left resolution of M and T[2] the left projective dimension of M
@*         it is T[1][j] the (j-1)-th syzygy module of M
"
{
  matrix T=transpose(M);
  list R,zw;
  R[1]=T;
  if (rinv(R[1])==0)
  {
    R[2]=transpose(std(transpose(lsyz(R[1]))));
  }
  else
  {
    matrix S[1][ncols(T)];
    R[1]=S;
    zw[1]=R;
    zw[2]=0;
    return (zw);
  }
  if (iszero(R[2])==1)
  {
    zw[1]=R;
    zw[2]=1;
    return (zw);
  }
  int i=1;
  matrix N;
  while (iszero(R[i+1])==0)
  {
    i=i+1;
    N=rinv(R[i]);
    if (iszero(N)==0)
    {
      if (i==2)
      {
        R[i-1]=concat(R[i-1],N);
        matrix K[1][nrows(R[1])];
        R[2]=K;
        zw[1]=R;
        zw[2]=i-1;
        return (zw);
      }
      if (i>2)
      {
        R[i-1]=concat(R[i-1],N);
        matrix K[ncols(N)][1];
        R[i-2]=concatz(R[i-2],K);
        R[i]=0;
        zw[1]=R;
        zw[2]=i-1;
        return(zw);
      }
    }
    R[i+1]=transpose(std(transpose(lsyz(R[i]))));
  }
  zw[1]=R;
  zw[2]=i;
  return (zw);
}