/////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Invariant theory";
info="
LIBRARY: stratify.lib   Algorithmic Stratification for Unipotent Group-Actions
AUTHOR:  Anne Fruehbis-Krueger, anne@mathematik.uni-kl.de

PROCEDURES:
  prepMat(M,wr,ws,step);  list of submatrices corresp. to given filtration
  stratify(M,wr,ws,step); algorithmic stratifcation (main procedure)
";
////////////////////////////////////////////////////////////////////////////
// REQUIRED LIBRARIES
////////////////////////////////////////////////////////////////////////////

// first the ones written in Singular
LIB "general.lib";
LIB "primdec.lib";

// then the ones written in C/C++

////////////////////////////////////////////////////////////////////////////
// PROCEDURES
/////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////
// For the kernel of the Kodaira-Spencer map in the case of hypersurface
// singularities or CM codimension 2 singularities:
// * step=min{ord(x_i)}
// * wr corresponds to the weight vector of the d/dt_i (i.e. to -ord(t_i))
//   (since the entries should be non-negative it may be necessary to
//   multiply the whole vector by -1)
// * ws corresponds to the weight vector of the delta_i
// * M is the matrix delta_i(t_j)
/////////////////////////////////////////////////////////////////////////////

proc prepMat(matrix M, intvec wr, intvec ws, int step)
"USAGE:   prepMat(M,wr,ws,step);
         where M is a matrix, wr is an intvec of size ncols(M),
         ws an intvec of size nrows(M) and step is an integer
RETURN:  2 lists of submatrices corresponding to the filtrations
         specified by wr and ws:
         the first list corresponds to the list for the filtration
         of AdA, i.e. the ranks of these matrices will be the r_i,
         the second one to the list for the filtration of L, i.e.
         the ranks of these matrices will be the s_i
NOTE:    * the entries of the matrix M are M_ij=delta_i(x_j),
         * wr is used to determine what subset of the set of all dx_i is
           generating AdF^l(A):
           if (k-1)*step <= wr[i] < k*step, then dx_i is in the set of
           generators of AdF^l(A) for all l>=k and the i-th column
           of M appears in each submatrix starting from the k-th
         * ws is used to determine what subset of the set of all delta_i
           is generating Z_l(L):
           if (k-1)*step <= ws[i] < k*step, then delta_i is in the set
           of generators of Z_l(A) for l < k and the i-th row of M
           appears in each submatrix up to the (k-1)th
         * the entries of wr and ws as well as step should be positive
           integers
EXAMPLE: example prepMat; shows an example"
{
//----------------------------------------------------------------------
// Initialization and sanity checks
//----------------------------------------------------------------------
  int i,j;
  if ((size(wr)!=ncols(M)) || (size(ws)!=nrows(M)))
  {
    "size mismatch: wr should have " + string(ncols(M)) + "entries";
    "               ws should have " + string(nrows(M)) + "entries";
    return("ERROR");
  }
//----------------------------------------------------------------------
// Sorting the matrix to obtain nice structure
//----------------------------------------------------------------------
  list sortwr=sort(wr);
  list sortws=sort(ws);
  if(sortwr[1]!=wr)
  {
    matrix N[nrows(M)][ncols(M)];
    for(i=1;i<=size(wr);i++)
    {
      N[1..nrows(M),i]=M[1..nrows(M),sortwr[2][i]];
    }
    wr=sortwr[1];
    M=N;
    kill N;
  }
  if(sortws[1]!=ws)
  {
    matrix N[nrows(M)][ncols(M)];
    for(i=1;i<=size(ws);i++)
    {
      N[i,1..ncols(M)]=M[sortws[2][i],1..ncols(M)];
    }
    ws=sortws[1];
    M=N;
    kill N;
  }
//---------------------------------------------------------------------
// Forming the submatrices
//---------------------------------------------------------------------
  list R,S;
  i=1;
  j=0;
  while ((step*(i-1))<=wr[size(wr)])
  {
    while ((step*i)>wr[j+1])
    {
      j++;
      if(j==size(wr)) break;
    }
    if(j!=0)
    {
      matrix N[nrows(M)][j]=M[1..nrows(M),1..j];
    }
    else
    {
      matrix N=matrix(0);
    }
    R[i]=N;
    kill N;
    i++;
    if(j==size(wr)) break;
  }
  i=1;
  j=0;
  while ((step*i)<=ws[size(ws)])
  {
    while ((step*i)>ws[j+1])
    {
      j++;
      if(j==size(ws)) break;
    }
    if(j==size(ws)) break;
    if(j!=0)
    {
      matrix N[nrows(M)-j][ncols(M)]=M[j+1..nrows(M),1..ncols(M)];
      S[i]=N;
      kill N;
    }
    else
    {
      S[i]=M;
    }
    i++;
  }
  list ret=R,S;
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(t(1..3)),dp;
  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
  print(M);
  intvec wr=1,3,5;
  intvec ws=2,4;
  int step=2;
  prepMat(M,wr,ws,step);
}
/////////////////////////////////////////////////////////////////////////////
static
proc minorList (list matlist)
"USAGE:   minorList(l);
         where l is a list of matrices satisfying the condition that l[i]
         is a submatrix of l[i+1]
RETURN:  list of lists in which each entry of the returned list corresponds
         to one of the matrices of the list l and is itself the list of
         the minors (i.e. the 1st entry is the ideal generated by the
         1-minors of the matrix etc.)
EXAMPLE: example minorList(l); shows an example"
{
//---------------------------------------------------------------------------
// Initialization and sanity checks
//---------------------------------------------------------------------------
  int maxminor;
  int counter;
  if(size(matlist)==0)
  {
    return(matlist);
  }
  for(int i=1;i<=size(matlist);i++)
  {
    if(((typeof(matlist[i]))!="matrix") && ((typeof(matlist[i]))!="intmat"))
    {
      "The list should only contain matrices or intmats";
      return("ERROR");
    }
  }
  list ret,templist;
  int j;
  int k=0;
  ideal minid;
//---------------------------------------------------------------------------
// find the maximal size of the minors and compute all possible minors,
// and put a minimal system of generators into the list that will be returned
//---------------------------------------------------------------------------
  for(i=1;i<=size(matlist);i++)
  {
    if (nrows(matlist[i]) < ncols(matlist[i]))
    {
      maxminor=nrows(matlist[i]);
    }
    else
    {
      maxminor=ncols(matlist[i]);
    }
    if (maxminor < 1)
    {
      "The matrices should be of size at least 1 x 1";
      return("ERROR");
    }
    kill templist;
    list templist;
    for(j=1;j<=maxminor;j++)
    {
      minid=minor(matlist[i],j);
      if(size(minid)>0)
      {
        if (defined(watchdog_interrupt))
        {
          kill watchdog_interrupt;
        }
        string watchstring="radical(ideal(";
        for(counter=1;counter <size(minid);counter++)
        {
          watchstring=watchstring+"eval("+string(minid[counter])+"),";
        }
        watchstring=watchstring+"eval("+string(minid[size(minid)])+")))";
        def watchtempid=watchdog(180,watchstring);
        kill watchstring;
        if ((defined(watchdog_interrupt)) || (typeof(watchtempid)=="string"))
        {
          templist[j-k]=mstd(minid)[2];
        }
        else
        {
          templist[j-k]=mstd(watchtempid)[2];
        }
        kill watchtempid;
      }
      else
      {
        k++;
      }
    }
    k=0;
    ret[i]=templist;
  }
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(t(1..3)),dp;
  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
  intvec wr=1,3,5;
  intvec ws=2,4;
  int step=2;
  list l=prepMat(M,wr,ws,step);
  l[1];
  minorList(l[1]);
}
/////////////////////////////////////////////////////////////////////////////
static
proc strataList(list Minors, list d, ideal V, int r, int nl)
"USAGE:   strataList(Minors,d,V,r,nl);
         Minors: list of minors as returned by minorRadList
         d:      list of polynomials
                 the open set that we are dealing with is D(d[1])
                 d[2..size(d)]=list of the factors of d
         V:      ideal
                 the closed set we are dealing with is V(V)
         r:      offset of the rank
         nl:     nesting level of the recursion
RETURN:  list of lists, each entry of the big list corresponds to one
         locally closed set and has the following entries:
         1) intvec giving the corresponding r- resp. s-vector
         2) ideal determining the closed set (cf. 3rd parameter V)
         3) list of polynomials determining the open set (cf. 2nd
            parameter d)
NOTE:    * sensible default values are
           d[1]=1; (list of length 1)
           V=ideal(0);
           r=0;
           nl=0;
           these parameters are only important in the recursion
           (if you know what you are doing, you are free to set d, V
           and r, but setting nl to a value other than 0 may give
           unexpected results)
         * no sanity checks are performed, since the procedure is designed
           for internal use only
         * for use with the list of minors corresponding to the s-vectors,
           the list of minors has to be specified from back to front
EXAMPLE: example strataList; shows an example"
{
//---------------------------------------------------------------------------
// * No sanity checks, since the procedure is static
// * First reduce everything using the ideal V of which we know
//   that the desired stratum lies in its zero locus
// * Throw away zero ideals
//---------------------------------------------------------------------------
  poly D=d[1];
  int i,j,k,ll;
  int isZero,isEmpty;
  intvec rv=r;
  intvec delvec;
  list l=mstd(V);
  ideal sV=l[1];
  ideal mV=l[2];
  list Ml;
  for(i=1;i<=size(Minors);i++)
  {
    list templist;
    for(j=1;j<=size(Minors[i]);j++)
    {
      templist[j]=reduce(Minors[i][j],sV);
    }
    Ml[i]=templist;
    kill templist;
  }
  for(i=1;i<=size(Ml);i++)
  {
    list templist;
    isZero=1;
    for(j=size(Ml[i]);j>=1;j--)
    {
      if(size(Ml[i][j])!=0)
      {
        templist[j]=Ml[i][j];
        isZero=0;
      }
      else
      {
        if(isZero==0)
        {
          return("ERROR");
        }
      }
    }
    if(size(templist)!=0)
    {
      Ml[i]=templist;
    }
    else
    {
      rv=rv,r;
      delvec=delvec,i;
    }
    kill templist;
  }
  if(size(delvec)>=2)
  {
    intvec dummydel=delvec[2..size(delvec)];
    Ml=deleteSublist(dummydel,Ml);
    kill dummydel;
  }
//---------------------------------------------------------------------------
// We do not need to go on if Ml disappeared
//---------------------------------------------------------------------------
  if(size(Ml)==0)
  {
    list ret;
    list templist;
    templist[1]=rv;
    templist[2]=mV;
    templist[3]=d;
    ret[1]=templist;
    return(ret);
  }
//---------------------------------------------------------------------------
// Check for minors which cannot vanish at all
//---------------------------------------------------------------------------
  def rt=basering;
  ring ru=0,(U),dp;
  def rtu=rt+ru;
  setring rtu;
  def tempMl;
  def ML;
  def D;
  setring rt;
  int Mlrank=0;
  setring rtu;
  tempMl=imap(rt,Ml);
  ML=tempMl[1];
  D=imap(rt,D);
  while(Mlrank<size(ML))
  {
    if(reduce(1,std(ML[Mlrank+1]+poly((U*D)-1)))==0)
    {
      Mlrank++;
    }
    else
    {
      break;
    }
  }
  setring rt;
  if(Mlrank!=0)
  {
    kill delvec;
    intvec delvec;
    isEmpty=1;
    for(i=1;i<=size(Ml);i++)
    {
      if(Mlrank<size(Ml[i]))
      {
        list templi2=Ml[i];
        list templist=templi2[Mlrank+1..size(Ml[i])];
        kill templi2;
        Ml[i]=templist;
        isEmpty=0;
      }
      else
      {
        if(isEmpty==0)
        {
          return("ERROR");
        }
        rv=rv,(r+Mlrank);
        delvec=delvec,i;
      }
      if(defined(templist)>1)
      {
        kill templist;
      }
    }
    if(size(delvec)>=2)
    {
      intvec dummydel=delvec[2..size(delvec)];
      Ml=deleteSublist(dummydel,Ml);
      kill dummydel;
    }
  }
//---------------------------------------------------------------------------
// We do not need to go on if Ml disappeared
//---------------------------------------------------------------------------
  if(size(Ml)==0)
  {
    list ret;
    list templist;
    templist[1]=intvec(rv);
    templist[2]=mV;
    templist[3]=d;
    ret[1]=templist;
    return(ret);
  }
//---------------------------------------------------------------------------
// For each possible rank of Ml[1] and each element of Ml[1][rk]
// call this procedure again
//---------------------------------------------------------------------------
  ideal Did;
  ideal newV;
  ideal tempid;
  poly f;
  list newd;
  int newr;
  list templist,retlist;
  setring rtu;
  ideal newV;
  ideal Did;
  def d;
  poly f;
  setring rt;
  for(i=0;i<=size(Ml[1]);i++)
  {
// find the polynomials which are not allowed to vanish simulateously
    if((i<size(Ml[1]))&&(i!=0))
    {
      Did=mstd(reduce(Ml[1][i],std(Ml[1][i+1])))[2];
    }
    else
    {
      if(i==0)
      {
        Did=0;
      }
      else
      {
        Did=mstd(Ml[1][i])[2];
      }
    }
// initialize the rank
    newr=r + Mlrank + i;
// find the new ideal V
    for(j=0;j<=size(Did);j++)
    {
      if((i!=0)&&(j==0))
      {
        j++;
        continue;
      }
      if(i<size(Ml[1]))
      {
        newV=mV,Ml[1][i+1];
      }
      else
      {
        newV=mV;
      }
// check whether the intersection of V and the new D is empty
      setring rtu;
      newV=imap(rt,newV);
      Did=imap(rt,Did);
      D=imap(rt,D);
      d=imap(rt,d);
      if(j==0)
      {
        if(reduce(1,std(newV+poly(D*U-1)))==0)
        {
          j++;
          setring rt;
          continue;
        }
      }
      if(i!=0)
      {
        if(reduce(1,std(newV+poly(Did[j]*D*U-1)))==0)
        {
          j++;
          setring rt;
          continue;
        }
        f=Did[j];
        for(k=2;k<=size(d);k++)
        {
          while(((f/d[k])*d[k])==f)
          {
            f=f/d[k];
          }
          if(deg(f)==0) break;
        }
      }
      setring rt;
      f=imap(rtu,f);
// i==0 ==> f==0 ==> deg(f)<=0
// otherwise factorize f, if it does not take too long,
// and add its factors, resp. f itself, to the list d
      if(deg(f)>0)
      {
        f=cleardenom(f);
        if (defined(watchdog_interrupt))
        {
          kill watchdog_interrupt;
        }
        def watchtempid=watchdog(180,"factorize(eval(" + string(f) + "),1)");
        if (defined(watchdog_interrupt))
        {
          newd=d;
          newd[size(d)+1]=f;
          newd[1]=d[1]*f;
        }
        else
        {
          tempid=watchtempid;
          templist=tempid[1..size(tempid)];
          newd=d+templist;
          f=newd[1]*f;
          tempid=simplify(ideal(newd[2..size(newd)]),14);
          templist=tempid[1..size(tempid)];
          kill newd;
          list newd=f;
          newd=newd+templist;
        }
        kill watchtempid;
      }
      else
      {
        newd=d;
      }
// take the corresponding sublist of the list of minors
      list Mltemp=delete(Ml,1);
      for(k=1;k<=size(Mltemp);k++)
      {
        templist=Mltemp[k];
        if(i<size(Mltemp[k]))
        {
          Mltemp[k]=list(templist[(i+1)..size(Mltemp[k])]);
        }
        else
        {
          kill templist;
          list templist;
          Mltemp[k]=templist;
        }
      }
// recursion
      templist=strataList(Mltemp,newd,newV,newr,(nl+1));
      kill Mltemp;
// build up the result list
      if(size(templist)!=0)
      {
        k=1;
        ll=1;
        while(k<=size(templist))
        {
          if(size(templist[k])!=0)
          {
            retlist[size(retlist)+ll]=templist[k];
            ll++;
          }
          k++;
        }
      }
    }
  }
  kill delvec;
  intvec delvec;
// clean up of the result list
  for(i=1;i<=size(retlist);i++)
  {
    if(typeof(retlist[i])=="none")
    {
      delvec=delvec,i;
    }
  }
  if(size(delvec)>=2)
  {
    intvec dummydel=delvec[2..size(delvec)];
    retlist=deleteSublist(dummydel,retlist);
    kill dummydel;
  }
// set the intvec to the correct value
  for(i=1;i<=size(retlist);i++)
  {
    if(nl!=0)
    {
      intvec tempiv=rv,retlist[i][1];
      retlist[i][1]=tempiv;
      kill tempiv;
    }
    else
    {
      if(size(rv)>1)
      {
        intvec tempiv=rv[2..size(rv)];
        retlist[i][1]=tempiv;
        kill tempiv;
      }
    }
  }
  return(retlist);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(t(1..3)),dp;
  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
  intvec wr=1,3,5;
  intvec ws=2,4;
  int step=2;
  list l=prepMat(M,wr,ws,step);
  l[1];
  list l2=minorRadList(l[1]);
  list d=poly(1);
  strataList(l2,d,ideal(0),0,0);
}
/////////////////////////////////////////////////////////////////////////////
static
proc cleanup(list stratlist)
"USAGE:   cleanup(l);
         where l is a list of lists in the format which is e.g. returned
         by strataList
RETURN:  list in which entries to the same integer vector have been
         joined to one entry
         the changed entries may be identified by checking whether its
         3rd entry is an empty list, then all entries starting from the
         4th one give the different possibilities for the open set
NOTE:    use the procedure killdups first to kill entries which are
         contained in other entries to the same integer vector
         otherwise the result will be meaningless
EXAMPLE: example cleanup; shows an example"
{
  int i,j;
  list leer;
  intvec delvec;
  if(size(stratlist)==0)
  {
    return(stratlist);
  }
  list ivlist;
// sort the list using the intvec as criterion
  for(i=1;i<=size(stratlist);i++)
  {
    ivlist[i]=stratlist[i][1];
  }
  list sortlist=sort(ivlist);
  list retlist;
  for(i=1;i<=size(stratlist);i++)
  {
    retlist[i]=stratlist[sortlist[2][i]];
  }
  i=1;
// find duplicate intvecs in the list
  while(i<size(stratlist))
  {
    j=i+1;
    while(retlist[i][1]==retlist[j][1])
    {
      retlist[i][3+j-i]=retlist[j][3];
      delvec=delvec,j;
      j++;
      if(j>size(stratlist)) break;
    }
    if (j!=(i+1))
    {
      retlist[i][3+j-i]=retlist[i][3];
      retlist[i][3]=leer;
      i=j-1;
// retlist[..][3] is empty if and only if there was more than one entry to this intvec
    }
    if(j>size(stratlist)) break;
    i++;
  }
  if(size(delvec)>=2)
  {
    intvec dummydel=delvec[2..size(delvec)];
    retlist=deleteSublist(dummydel,retlist);
    kill dummydel;
  }
  return(retlist);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(t(1),t(2)),dp;
  intvec iv=1;
  list plist=t(1),t(1);
  list l1=iv,ideal(0),plist;
  plist=t(2),t(2);
  list l2=iv,ideal(0),plist;
  list l=l1,l2;
  cleanup(l);
}
/////////////////////////////////////////////////////////////////////////////
static
proc joinRS(list Rlist,list Slist)
"USAGE:   joinRS(Rlist,Slist);
         where Rlist and Slist are lists in the format returned by
         strataList
RETURN:  one list in the format returned by stratalist in which the
         integer vector is the concatenation of the corresponding vectors
         from Rlist and Slist
         (of course only non-empty locally closed sets are returned)
NOTE:    since Slist is a list returned by strataList corresponding to the
         s-vector, it corresponds to the list of minors read from back to
         front
EXAMPLE: no example available at the moment"
{
  int j,k;
  list retlist;
  list templist,templi2;
  intvec tempiv;
  ideal tempid;
  ideal dlist;
  poly D;
  def rt=basering;
  ring ru=0,(U),dp;
  def rtu=rt+ru;
  setring rtu;
  def Rlist=imap(rt,Rlist);
  def Slist=imap(rt,Slist);
  setring rt;
  for(int i=1;i<=size(Rlist);i++)
  {
    for(j=1;j<=size(Slist);j++)
    {
// skip empty sets
      if(Rlist[i][1][size(Rlist[i][1])]<Slist[j][1][size(Slist[j][1])])
      {
        j++;
        continue;
      }
      setring rtu;
      if(reduce(1,std(Slist[j][2]+poly(((Rlist[i][3][1])*U)-1)))==0)
      {
        j++;
        setring rt;
        continue;
      }
      if(reduce(1,std(Rlist[i][2]+poly(((Slist[j][3][1])*U)-1)))==0)
      {
        j++;
        setring rt;
        continue;
      }
      setring rt;
// join the intvecs and the ideals V
      tempiv=Rlist[i][1],Slist[j][1];
      kill templist;
      list templist;
      templist[1]=tempiv;
      if(size(Rlist[i][2]+Slist[j][2])>0)
      {
        templist[2]=mstd(Rlist[i][2]+Slist[j][2])[2];
      }
      else
      {
        templist[2]=ideal(0);
      }
// test again whether we are talking about the empty set
      setring rtu;
      def templist=imap(rt,templist);
      def tempid=templist[2];
      if(reduce(1,std(tempid+poly(((Slist[j][3][1])*(Rlist[i][3][1])*U)-1)))==0)
      {
        kill templist;
        kill tempid;
        j++;
        setring rt;
        continue;
      }
      else
      {
        kill templist;
        kill tempid;
        setring rt;
      }
// join the lists d
      if(size(Rlist[i][3])>1)
      {
        templi2=Rlist[i][3];
        dlist=templi2[2..size(templi2)];
      }
      else
      {
        kill dlist;
        ideal dlist;
      }
      if(size(Slist[j][3])>1)
      {
        templi2=Slist[j][3];
        tempid=templi2[2..size(templi2)];
      }
      else
      {
        kill tempid;
        ideal tempid;
      }
      dlist=dlist+tempid;
      dlist=simplify(dlist,14);
      D=1;
      for(k=1;k<=size(dlist);k++)
      {
        D=D*dlist[k];
      }
      if(size(dlist)>0)
      {
        templi2=D,dlist[1..size(dlist)];
      }
      else
      {
        templi2=list(1);
      }
      templist[3]=templi2;
      retlist[size(retlist)+1]=templist;
    }
  }
  return(retlist);
}
////////////////////////////////////////////////////////////////////////////

proc stratify(matrix M, intvec wr, intvec ws,int step)
"USAGE:   stratify(M,wr,ws,step);
         where M is a matrix, wr is an intvec of size ncols(M),
         ws an intvec of size nrows(M) and step is an integer
RETURN:  list of lists, each entry of the big list corresponds to one
         locally closed set and has the following entries:
         1) intvec giving the corresponding rs-vector
         2) ideal determining the closed set
         3) list d of polynomials determining the open set D(d[1])
            empty list if there is more than one open set
         4-n) lists of polynomials determining open sets which all lead
            to the same rs-vector
NOTE:    * ring ordering should be global, i.e. the ring should be a
           polynomial ring
         * the entries of the matrix M are M_ij=delta_i(x_j),
         * wr is used to determine what subset of the set of all dx_i is
           generating AdF^l(A):
           if (k-1)*step < wr[i] <= k*step, then dx_i is in the set of
           generators of AdF^l(A) for all l>=k
         * ws is used to determine what subset of the set of all delta_i
           is generating Z_l(L):
           if (k-1)*step <= ws[i] < k*step, then delta_i is in the set
           of generators of Z_l(A) for l < k
         * the entries of wr and ws as well as step should be positive
           integers
         * the filtrations have to be known, no sanity checks concerning
           the filtrations are performed !!!
EXAMPLE: example stratify; shows an example"
{
//---------------------------------------------------------------------------
// Initialization and sanity checks
//---------------------------------------------------------------------------
  int i,j;
  list submat=prepMat(M,wr,ws,step);
  if(defined(watchProgress))
  {
    "List of submatrices to consider:";
    submat;
  }
  if(ncols(submat[1][size(submat[1])])==nrows(submat[1][size(submat[1])]))
  {
    int symm=1;
    int nr=nrows(submat[1][size(submat[1])]);
    for(i=1;i<=nr;i++)
    {
      for(j=1;j<=nr-i;j++)
      {
        if(submat[1][size(submat[1])][i,j]!=submat[1][size(submat[1])][nr-j+1,nr-i+1])
        {
          symm=0;
          break;
        }
      }
      if(symm==0) break;
    }
  }
  if(defined(symm)>1)
  {
    if(symm==0)
    {
      kill symm;
    }
  }
  list Rminors=minorList(submat[1]);
  if(defined(watchProgress))
  {
    "minors corresponding to the r-vector:";
    Rminors;
  }
  if(defined(symm)<2)
  {
    list Sminors=minorList(submat[2]);
    if(defined(watchProgress))
    {
      "minors corresponding to the s-vector:";
      Sminors;
    }
  }
  if(size(Rminors[1])==0)
  {
    Rminors=delete(Rminors,1);
  }
//---------------------------------------------------------------------------
// Start the recursion and cleanup afterwards
//---------------------------------------------------------------------------
  list leer=poly(1);
  list Rlist=strataList(Rminors,leer,0,0,0);
  if(defined(watchProgress))
  {
    "list of strata corresponding to r-vectors:";
    Rlist;
  }
  Rlist=killdups(Rlist);
  if(defined(watchProgress))
  {
    "previous list after killing duplicate entries:";
    Rminors;
  }
  if(defined(symm)<2)
  {
// Sminors have the smallest entry as the last one
// In order to use the same routines as for the Rminors
// we handle the s-vector in inverse order
    list Stemp;
    for(i=1;i<=size(Sminors);i++)
    {
      Stemp[size(Sminors)-i+1]=Sminors[i];
    }
    list Slist=strataList(Stemp,leer,0,0,0);
    if(defined(watchProgress))
    {
      "list of strata corresponding to s-vectors:";
      Slist;
    }
//---------------------------------------------------------------------------
// Join the Rlist and the Slist to obtain the stratification
//---------------------------------------------------------------------------
    Slist=killdups(Slist);
    if(defined(watchProgress))
    {
      "previous list after killing duplicate entries:";
      Slist;
    }
    list ret=joinRS(Rlist,Slist);
    if(defined(watchProgress))
    {
      "list of strata corresponding to r- and s-vectors:";
      ret;
    }
    ret=killdups(ret);
    if(defined(watchProgress))
    {
      "previous list after killing duplicate entries:";
      ret;
    }
    ret=cleanup(ret);
  }
  else
  {
    list ret=cleanup(Rlist);
  }
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(t(1..3)),dp;
  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
  intvec wr=1,3,5;
  intvec ws=2,4;
  int step=2;
  stratify(M,wr,ws,step);
}
/////////////////////////////////////////////////////////////////////////////
static
proc killdups(list l)
"USAGE:   killdups(l);
         where l is a list in the form returned by strataList
RETURN:  list which is obtained from the previous list by leaving out
         entries which have the same intvec as another entry in which
         the locally closed set is contained
EXAMPLE: no example available at the moment"
{
  int i=1;
  int j,k,skip;
  while(i<size(l))
  {
    intvec delvec;
    for(j=i+1;j<=size(l);j++)
    {
// we do not need to check the V ideals, since the intvecs coincide
      if(l[i][1]==l[j][1])
      {
        if((l[i][3][1]/l[j][3][1])*l[j][3][1]==l[i][3][1])
        {

          delvec=delvec,i;
          break;
        }
        else
        {
          if((l[j][3][1]/l[i][3][1])*l[i][3][1]==l[j][3][1])
          {
            delvec=delvec,j;
            j++;
            continue;
          }
        }
      }
    }
    if(size(delvec)>=2)
    {
      delvec=sort(delvec)[1];
      intvec dummydel=delvec[2..size(delvec)];
      l=deleteSublist(dummydel,l);
      kill dummydel;
    }
    kill delvec;
    i++;
  }
  list ret=l;
  return(ret);
}
/////////////////////////////////////////////////////////////////////////////