source: git/Singular/LIB/spcurve.lib @ 38c165

// $Id: spcurve.lib,v 1.9 2000-05-12 12:17:18 krueger Exp $
// (anne, last modified 31.5.99)
/////////////////////////////////////////////////////////////////////////////

version="$Id: spcurve.lib,v 1.9 2000-05-12 12:17:18 krueger Exp $";
info="
LIBRARY:  spcurve.lib    PROCEDURES FOR CM CODIMENSION 2 SINGULARITIES
AUTHOR:   Anne Fruehbis-Krueger, anne@mathematik.uni-kl.de
last modified: 31.5.99
SEE ALSO: sing_lib, deform_lib

PROCEDURES:
 isCMcod2(i);            presentation matrix of the ideal i, if i is CM
 CMtype(i);              Cohen-Macaulay type of the ideal i
 matrixT1(M,n);          1st order deformation T1 in matrix description
 semiCMcod2(M,T1);       semiuniversal deformation of maximal minors of M
 discr(sem,n);           discriminant of semiuniversal deformation
 qhmatrix(M);            weights if M is quasihomogeneous
 relweight(N,W,a);       relative matrix weight of N w.r.t. weights (W,a)
 posweight(M,T1,i);      deformation of coker(M) of non-negative weight
 KSpencerKernel(M);      kernel of the Kodaira-Spencer map
 mod2id(M,iv);           conversion of a module M to an ideal
 id2mod(i,iv);           conversion inverse to mod2id
 subrInterred(i1,i2,iv)  interred w.r.t. a subset of variables
";

LIB "elim.lib";
LIB "homolog.lib";
LIB "inout.lib";
LIB "poly.lib";
/////////////////////////////////////////////////////////////////////////////

proc isCMcod2(ideal kurve)
"USAGE:   @code{isCMcod2(i)};   @code{i} ideal
RETURN:  presentation matrix of @code{i}, if @code{i} is Cohen-Macaulay of codimension 2
         @* a zero matrix otherwise
EXAMPLE: @code{example isCMcod2}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------------
// Compute a minimal free resolution of the ideal and check if the
// resolution has the expected structure
//---------------------------------------------------------------------------
  list kurveres=mres(kurve,0);
  matrix M=kurveres[2];
  if ((size(kurveres)>3) &&
         ((size(kurveres[3])>1) ||
          ((size(kurveres[3])<=1) && (kurveres[3][1,1]!=0))))
  {
    dbprint(p,"//not Cohen-Macaulay, codim 2");
    matrix ret=0;
    return(ret);
  }
  return(M);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=32003,(x,y,z),ds;
  ideal i=xz,yz,x^3-y^4;
  print(isCMcod2(i));
}
/////////////////////////////////////////////////////////////////////////////

proc CMtype(ideal kurve)
"USAGE:   @code{CMtype(i)};  @code{i} ideal, CM of codimension 2
RETURN:  Cohen-Macaulay type of @code{i} (integer)
         @*(-1, if i is not Cohen-Macaulay of codimension 2)
EXAMPLE: @code{example CMtype}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
  int gt = -1;
//---------------------------------------------------------------------------
// Compute a minimal free resolution of the ideal and check if the
// resolution has the expected structure
//---------------------------------------------------------------------------
  list kurveres;
  kurveres=mres(kurve,0);
  if ((size(kurveres)>3) &&
         ((size(kurveres[3])>1) ||
          ((size(kurveres[3])<=1) && (kurveres[3][1,1]!=0))))
  {
    dbprint(p,"//not Cohen-Macaulay, codim 2");
    return(gt);
  }
//---------------------------------------------------------------------------
// Return the Cohen-Macaulay type of i
//---------------------------------------------------------------------------
  matrix M = matrix(kurveres[2]);
  gt = ncols(M);
  return(gt);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=32003,(x,y,z),ds;
  ideal i=xy,xz,yz;
  CMtype(i);
}
/////////////////////////////////////////////////////////////////////////////

proc matrixT1(matrix M ,int n)
"USAGE:   @code{matrixT1(M,n)};  @code{M} matrix, @code{n} integer
ASSUME:  @code{M} is a presentation matrix of an ideal i, CM of codimension 2;
         @* consider i as a family of ideals in a ring in the first @code{n}
         variables where the remaining variables are considered as
         parameters
RETURN:  list consisting of the k x (k+1) matrix @code{M} and a 
         module @math{K_M} such that @math{T1=Mat(k,k+1;R)/K_M} is the space 
         of first order deformations of i
EXAMPLE: @code{example matrixT1}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//--------------------------------------------------------------------------
// Initialization and sanity checks
//--------------------------------------------------------------------------
  int nr=nrows(M);
  int nc=ncols(M);
  if ( nr < nc )
  {
    M=transpose(M);
    temp=nc;
    nc=nr;
    nr=temp;
    tra=1;
  }
  if ( nr != (nc+1) )
  {
    dbprint(p,"//not a k x (k+1) matrix");
    return("ERROR");
  }
//---------------------------------------------------------------------------
// Construct the denominator - step by step
// step 1: initialization
//---------------------------------------------------------------------------
  int gt=nc;
  int i,j;
  ideal m = M;
  ideal dx;
  ideal rv;
  ideal lv;
  matrix R[gt][gt]=0;
  matrix L[gt+1][gt+1]=0;
  matrix T1[n+gt*gt+(gt+1)*(gt+1)][gt*(gt+1)] = 0;
//---------------------------------------------------------------------------
// step 2: the derivatives of the matrix are generators of the denominator
//---------------------------------------------------------------------------
  for( i=1; i<= n; i++ )
  {
    dx=diff(m,var(i));
    T1[i,1..gt*(gt+1)] = dx;
  }
//---------------------------------------------------------------------------
// step 3: M*R is a generator as well
//---------------------------------------------------------------------------
  for( i=1; i <= gt; i++ )
  {
    for ( j=1 ; j <= gt ; j++ )
    {
      R[i,j]=1;
      rv = M * R;
      T1[n+(i-1)*gt+j,1..gt*(gt+1)] = rv;
      R[i,j]=0;
    }
  }
//---------------------------------------------------------------------------
// step 4: so is L*M
//---------------------------------------------------------------------------
  for( i=1; i <= (gt+1); i++)
  {
    for( j=1 ; j <= (gt+1);j++ )
    {
      L[i,j]=1;
      lv = L * M;
      T1[n+gt*gt+(i-1)*(gt+1)+j,1..gt*(gt+1)] = lv;
      L[i,j]=0;
    }
  }
//---------------------------------------------------------------------------
// Compute the vectorspace basis of T1
//---------------------------------------------------------------------------
  module t1 = module(transpose(T1));
  list result=M,t1;
  return(result);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=32003,(x(1),x(2),x(3)),dp;
  ideal curve=x(1)*x(2),x(1)*x(3),x(2)*x(3);
  matrix M=isCMcod2(curve);
  matrixT1(M,3);
}
/////////////////////////////////////////////////////////////////////////////

proc semiCMcod2(matrix M, module t1)
"USAGE:   @code{semiCMcod2(M,t1)}; @code{M} matrix, @code{t1} module
ASSUME:  @code{M} is a presentation matrix of an ideal i, CM of codimension 2,
         and @code{t1} is a presentation of the space of first order 
         deformations of i (@code{(M,t1)} as returned by the procedure 
         @code{matrixT1})
CREATE:  new basering with name rneu
RETURN:  ideal in rneu describing the semiuniversal deformation of i
NOTE:    The current basering should not contain any variables named
         A(j) where j is some integer!
EXAMPLE: @code{example semiCMcod2}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------------
// Initialization
//---------------------------------------------------------------------------
  module t1erz=kbase(std(t1));
  int tau=vdim(t1);
  int gt=ncols(M);
  int i;
  def r=basering;
  if(size(M)!=gt*(gt+1))
  {
    gt=gt-1;
  }
  for(i=1; i<=size(t1erz); i++)
  {
    if(rvar(A(i)))
    {
      int jj=-1;
      break;
    }
  }
  if (defined(jj)>1)
  {
    if (jj==-1)
    {
      "Your ring contains a variable T(i)!";
      return("ERROR");
    }
  }
//---------------------------------------------------------------------------
// Definition of the new ring and the image of M and t1 in the new ring
//---------------------------------------------------------------------------
  ring rtemp=0,(A(1..tau)),dp;
  def rneu=r+rtemp;
  setring rneu;
  matrix M=imap(r,M);
  ideal m=M;
  module t1erz=imap(r,t1erz);
//---------------------------------------------------------------------------
// Construction of the presentation matrix of the versal deformation
//---------------------------------------------------------------------------
  matrix N=matrix(m);
  matrix Mtemp[gt*(gt+1)][1];
  for( i=1; i<=tau; i++)
  {
    Mtemp=t1erz[i];
    N=N+A(i)*transpose(Mtemp);
  }
  ideal n=N;
  matrix O[gt+1][gt]=n;
//---------------------------------------------------------------------------
// Construction of the return value
//---------------------------------------------------------------------------
  ideal result=minor(O,gt);
  export rneu;
  keepring rneu;
  return(result);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=32003,(x(1),x(2),x(3)),dp;
  ideal curve=x(1)*x(2),x(1)*x(3),x(2)*x(3);
  matrix M=isCMcod2(curve);
  list l=matrixT1(M,3);
  semiCMcod2(l[1],std(l[2]));
}
/////////////////////////////////////////////////////////////////////////////

proc discr(ideal kurve, int n)
"USAGE:   @code{discr(sem,n)};  @code{sem} ideal, @code{n} integer
ASSUME:  @code{sem} is the versal deformation of an ideal of codimension 2
         @*the first n variables of the ring are treated as variables
         all the others as parameters
RETURN:  ideal describing the discriminant
NOTE:    This is not a powerful algorithm!
EXAMPLE: @code{example discr}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------------
// some sanity checks and initialization
//---------------------------------------------------------------------------
  int i;
  ideal sem=std(kurve);
  ideal semdiff;
  ideal J2;
  int ncol=ncols(matrix(sem));
  matrix Jacob[n][ncol];
//---------------------------------------------------------------------------
// compute the Jacobian matrix
//---------------------------------------------------------------------------
  for (i=1; i<=n; i++)
  {
    semdiff=diff(sem,var(i));
    Jacob[i,1..ncol]=semdiff;
  }
//---------------------------------------------------------------------------
// eliminate the first n variables in the ideal generated by
// the versal deformation and the 2x2 minors of the Jacobian
//---------------------------------------------------------------------------
  semdiff=minor(Jacob,2);
  J2=sem,semdiff;
  J2=std(J2);
  poly eli=1;
  for(i=1; i<=n; i++)
  {
    eli=eli*var(i);
  }
  ideal dis=eliminate(J2,eli);
  return(dis);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=32003,(x(1),x(2),x(3)),dp;
  ideal curve=x(1)*x(2),x(1)*x(3),x(2)*x(3);
  matrix M=isCMcod2(curve);
  list l=matrixT1(M,3);
  def sem=semiCMcod2(l[1],std(l[2]));
  basering;
  discr(sem,3);
}
/////////////////////////////////////////////////////////////////////////////

proc qhmatrix(matrix M)
"USAGE:   @code{qhmatrix(M)};   @code{M} a k x (k+1) matrix
RETURN:  list, consisting of an integer vector containing the weights of
         the variables of the basering and an integer matrix giving the
         weights of the entries of @code{M}, if @code{M} is quasihomogeneous;
         @*zero integer vector and zero integer matrix, if @code{M} is not
         quasihomogeneous, i.e. does not allow row and column weights
EXAMPLE: @code{example qhmatrix}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------------
// Initialization and sanity checks
//---------------------------------------------------------------------------
  def r=basering;
  int i,j,temp;
  int tra=0;
  int nr=nrows(M);
  int nc=ncols(M);
  if ( nr > nc )
  {
    M=transpose(M);
    temp=nc;
    nc=nr;
    nr=temp;
    tra=1;
  }
  if ( nc != (nr+1) )
  {
    dbprint(p,"//not a k x (k+1) matrix");
    return("ERROR");
  }
  ideal m=minor(M,nr);
//---------------------------------------------------------------------------
// get the weight using the fact that the matrix is quasihomogeneous, if
// its maximal minors are, and check, whether M is really quasihomogeneous
//---------------------------------------------------------------------------
  intvec a=weight(m);
  string tempstr="ring rneu=" + charstr(r) + ",(" + varstr(r) + "),Ws(" + string(a) + ");";
  execute(tempstr);
  def M=imap(r,M);
  int difset=0;
  list l;
  int dif;
  int donttest=0;
  int comprow=0;
  intmat W[nr][nc];
//---------------------------------------------------------------------------
// find a row not containing a 0
//---------------------------------------------------------------------------
  for(i=1; i<=nr; i++)
  {
    if(comprow==0)
    {
      comprow=i;
      for(j=1; j<=nc; j++)
      {
        if(M[i,j]==0)
        {
          comprow=0;
          break;
        }
      }
    }
  }
//---------------------------------------------------------------------------
// get the weights of the comprow'th row or use emergency exit
//---------------------------------------------------------------------------
  if(comprow==0)
  {
    intvec v=0;
    intmat V=0
    list ret=v,V;
    return(ret);
  }
  else
  {
    for(j=1; j<=nc; j++)
    {
      l[j]=deg(lead(M[comprow,j]));
    }
  }
//---------------------------------------------------------------------------
// do the checks
//---------------------------------------------------------------------------
  for(i=1; i<=nr; i++)
  {
    if ( i==comprow )
    {
// this row should not be tested against itself
      donttest=1;
    }
    else
    {
// initialize the difference of the rows, but ignore 0-entries
      if (M[i,1]!=0)
      {
        dif=deg(lead(M[i,1]))-l[1];
        difset=1;
      }
      else
      {
        list memo;
        memo[1]=1;
      }
    }
// check column by column
    for(j=1; j<=nc; j++)
    {
      if(M[i,j]==0)
      {
        if(defined(memo)!=0)
        {
          memo[size(memo)+1]=j;
        }
        else
        {
          list memo;
          memo[1]=j;
        }
      }
      temp=deg(lead(M[i,j]));
      if((difset!=1) && (donttest!=1) && (M[i,j]!=0))
      {
// initialize the difference of the rows, if necessary - still ignore 0s
        dif=deg(lead(M[i,j]))-l[j];
        difset=1;
      }
// is M[i,j] quasihomogeneous - else emergency exit
      if(M[i,j]!=jet(M[i,j],temp,a)-jet(M[i,j],temp-1,a))
      {
        intvec v=0;
        intmat V=0;
        list ret=v,V;
        return(ret);
      }
      if(donttest!=1)
      {
// check row and column weights - else emergency exit
        if(((temp-l[j])!=dif) && (M[i,j]!=0) && (difset==1))
        {
          intvec v=0;
          intmat V=0;
          list ret=v,V;
          return(ret);
        }
      }
// set the weight matrix entry
      W[i,j]=temp;
    }
// clean up the 0's we left out
    if((difset==1) && (defined(memo)!=0))
    {
      for(j=1; j<=size(memo); j++)
      {
        W[i,memo[j]]=dif+l[memo[j]];
      }
      kill memo;
    }
    donttest=0;
  }
//---------------------------------------------------------------------------
// transpose, if M was transposed during initialization, and return the list
//---------------------------------------------------------------------------
  if ( tra==1 )
  {
    W=transpose(W);
  }
  setring r;
  list ret=a,W;
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),ds;
  matrix M[3][2]=z,0,y,x,x^3,y;
  qhmatrix(M);
  pmat(M);
}
/////////////////////////////////////////////////////////////////////////////

proc relweight(matrix N, intmat W, intvec a)
"USAGE:   @code{relweight(N,W,a)}; @code{N} matrix, @code{W} intmat, @code{a} intvec
ASSUME:  @code{N} is a non-zero matrix
         @* @code{W} is an integer matrix of the same size as @code{N}
         @* @code{a} is an integer vector giving the weights of the variables
RETURN:  integer, @math{max(a-weighted order(N_ij) - W_ij | all entries ij)}
         @* string \"ERROR\" if sizes do not match
EXAMPLE: @code{example relweight}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------------
// Initialization and sanity checks
//---------------------------------------------------------------------------
  if ((size(N)!=size(W)) || (ncols(N)!=ncols(W)))
  {
    dbprint(p,"//matrix size does not match");
    return("ERROR");
  }
  if (size(a)!=nvars(basering))
  {
    dbprint(p,"//length of weight vector != number of variables");
    return("ERROR");
  }
  int i,j,temp;
  def r=basering;
//---------------------------------------------------------------------------
// Comparision entry by entry
//---------------------------------------------------------------------------
  for(i=1; i<=nrows(N); i++)
  {
    for(j=1; j<=ncols(N); j++)
    {
      if (N[i,j]!=0)
      {
        temp=mindeg1(N[i,j],a)-W[i,j];
        if (defined(ret))
        {
          if(temp > ret)
          {
            ret=temp;
          }
        }
        else
        {
          int ret=temp;
        }
      }
    }
  }
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=32003,(x,y,z),ds;
  matrix N[2][3]=z,0,y,x,x^3,y;
  intmat W[2][3]=1,1,1,1,1,1;
  intvec a=1,1,1;
  relweight(N,W,a);
}
/////////////////////////////////////////////////////////////////////////////

proc posweight(matrix M, module t1, int choose, list #)
"USAGE:   @code{posweight(M,t1,n[,s])}; @code{M} matrix, @code{t1} module, @code{n} int, @code{s} string
         @*@code{n}=0 : all deformations of non-negative weight
         @*@code{n}=1 : only non-constant deformations of non-negative weight
         @*@code{n}=2 : all deformations of positive weight
         @*As an optional parameter the name of a new ring may be
         specified.
ASSUME:  @code{M} is a presentation matrix of a Cohen-Macaulay codimension 2
         ideal and @code{t1} is a presentation of its T1 space in matrix 
         notation
CREATE:  new basering (default name: rneu); a different name for this ring
         may be given as a 4th parameter
RETURN:  list, consisting of a presentation matrix describing the deformation
         given by the generators of T1 of non-negative/positive weight
         and the weight vector for the new variables
NOTE:   The current basering should not contain any variables named
         T(i) where i is some integer!
EXAMPLE: @code{example posweight}; shows an example
"
{
//---------------------------------------------------------------------------
// Initialization and sanity checks
//---------------------------------------------------------------------------
  if (size(#)>0)
  {
    if (typeof(#[1])=="string")
    {
      string newname=#[1];
    }
  }
  if (attrib(t1,"isSB"))
  {
    module t1erz=kbase(t1);
    int tau=vdim(t1);
  }
  else
  { module t1erz=kbase(std(t1));
    int tau=vdim(std(t1));
  }
  for(int i=1; i<=size(t1erz); i++)
  {
    if(rvar(T(i)))
    {
      int jj=-1;
      break;
    }
  }
  kill i;
  if (defined(jj))
  {
    if (jj==-1)
    {
      "Your ring contains a variable T(i)!";
      return("ERROR");
    }
  }
  int pw=0;
  int i;
  def r=basering;
  list l=qhmatrix(M);
  int gt=ncols(M);
  if(size(M)!=gt*(gt+1))
  {
    gt=gt-1;
  }
  matrix erzmat[gt+1][gt];
  list erz;
  if ((size(l[1])==1) && (l[1][1]==0) && (size(l[2])==1) && (l[2][1,1]==0))
  {
    return("ERROR");
  }
//---------------------------------------------------------------------------
// Find the generators of T1 of non-negative weight
//---------------------------------------------------------------------------
  int relw;
  list rlw;
  for(i=1; i<=tau; i++)
  {
    erzmat=t1erz[i];
    kill relw;
    def relw=relweight(erzmat,l[2],l[1]);
    if(typeof(relw)=="int")
    {
      if (((choose==0) && (relw>=0))
           || ((choose==1) && (relw>=0) && (CMtype(minor(M+erzmat,gt))==gt))
           || ((choose==2) && (relw > 0)))
      {
        pw++;
        rlw[pw]=relw;
        erz[pw]=erzmat;
      }
    }
    else
    {
      return("ERROR");
    }
  }
//---------------------------------------------------------------------------
// Definition of the new ring and the image of M and erz in the new ring
//---------------------------------------------------------------------------
  if(size(rlw)==0)
  {
    return("ERROR");
  }
  intvec iv=rlw[1..size(rlw)];
  ring rtemp=0,(T(1..pw)),dp;
  def rneu=r+rtemp;
  setring rneu;
  matrix M=imap(r,M);
  ideal m=M;
// we cannot imap erz, if its size=0
  if(pw==0)
  {
    list erz1;
  }
  else
  {
    list erz1=imap(r,erz);
  }
//---------------------------------------------------------------------------
// Construction of the presentation matrix of the deformation
//---------------------------------------------------------------------------
  matrix N=matrix(m);
  ideal mtemp;
  matrix Mtemp[gt*(gt+1)][1];
  for( i=1; i<=pw; i++)
  {
    mtemp=erz1[i];
    Mtemp=mtemp;
    N=N+T(i)*transpose(Mtemp);
  }
  ideal n=N;
  matrix O[gt+1][gt]=n;
//---------------------------------------------------------------------------
// Keep the ring and return the matrix
//---------------------------------------------------------------------------
  if (defined(newname)>1)
  {
    def `newname`=rneu;
    setring `newname`;
    export `newname`;
    keepring `newname`;
  }
  else
  {
    export rneu;
    keepring rneu;
  }
  list ret=O,iv;
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=32003,(x(1),x(2),x(3)),dp;
  ideal curve=(x(3)-x(1)^2)*x(3),(x(3)-x(1)^2)*x(2),x(2)^2-x(1)^7*x(3);
  matrix M=isCMcod2(curve);
  list l=matrixT1(M,3);
  list li=posweight(l[1],std(l[2]),0);
  pmat(li[1]);
  li[2];
}
/////////////////////////////////////////////////////////////////////////////

proc KSpencerKernel(matrix M,list #)
"USAGE:     @code{KSpencerKernel(M[,s][,v])};  @code{M} matrix, @code{s} string, @code{v} intvec
           @*optional parameters (please specify in this order, if both are
           present):
           @* *  @code{s} = first of the names of the new rings
           @*   e.g. \"R\" leads to ring names R and R1
           @* *  @code{v} of size n(n+1) leads to the following module ordering
           @*   @code{gen(v[1]) > gen(v[2]) > ... > gen(v[n(n+1)])}
           @*   where the matrix entry ij corresponds to @code{gen((i-1)*n+j)}
ASSUME:    @code{M} is a quasihomogeneous n x (n+1) matrix where the n minors 
           define an isolated space curve singularity
CREATE:    2 new rings (default names: rneu and reneu)
           different ring names may be specified as a 2nd parameter
RETURN:    coefficient matrix representing the kernel of the Kodaira-
           Spencer map of the family of non-negative deformations
           having the given singularity as special fibre
NOTE:      * the initial basering should not contain variables with name
             e(i) or T(i), since those variable names will be internally
             used by the script
           * setting an intvec with 5 entries and name watchProgress
             shows the progress of the computations:
           @*  watchProgress[1]>0  => option(prot) in groebner commands
           @*  watchProgress[2]>0  => trace output for highcorner
           @*  watchProgress[3]>0  => output of deformed matrix
           @*  watchProgress[4]>0  => result of elimination step
           @*  watchProgress[4]>1  => trace output of multiplications with xyz
           @*                         and subsequent reductions
           @*  watchProgress[5]>0  => matrix representing the kernel using print
EXAMPLE:   @code{example KSpencerKernel}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------------
// Initialization and sanity checks
//---------------------------------------------------------------------------
  intvec optvec=option(get);
  if (size(#)>0)
  {
    if (typeof(#[1])=="string")
    {
      string newname=#[1];
    }
    if (typeof(#[1])=="intvec")
    {
      intvec desiredorder=#[1];
    }
    if (size(#)>1)
    {
      if (typeof(#[2])=="intvec")
      {
        intvec desiredorder=#[2];
      }
    }
  }
  if (defined(watchProgress))
  {
    if ((typeof(watchProgress)!="intvec") || (size(watchProgress)<5))
    {
      "watchProgress should be an intvec with at least 5 entries";
      "ignoring watchProgress";
      def kksave=watchProgress;
      kill watchProgress;
    }
  }
  option(redTail);
  if (nvars(basering) != 3 )
  {
    "It should be a curve in 3 space";
    return("ERROR");
  }
//---------------------------------------------------------------------------
// change to a basering with the correct weihted order
//---------------------------------------------------------------------------
  def rt=basering;
  list wl=qhmatrix(M);
  if ((size(wl)!=2) || ((wl[1]==0) && (wl[2]==0)))
  {
    "The matrix was not n x (n+1) or not quasihomogenous";
    return("ERROR");
  }
  string ringre=" ring r=" + charstr(rt) + ",(x,y,z), Ws(" + string(wl[1]) + ");";
  execute(ringre);
  matrix M=imap(rt,M);
  int ne=size(M);
  if (defined(desiredorder)>1)
  {
    intvec iv;
    for(int i=1;i<=size(desiredorder);i++)
    {
      iv[desiredorder[i]]=i;
    }
  }
  else
  {
    intvec iv=1..ne;
  }
  list l=matrixT1(M,3);
  if (dim(std(l[2])) != 0)
  {
    "The matrix does not define an isolated space curve singularity";
    return("ERROR");
  }
  module t1qh=l[2];
//--------------------------------------------------------------------------
// Passing to a new ring with extra variables e(i) corresponding to
// the module generators gen(i) for weighted standard basis computation
// accepting weights for the gen(i)
//--------------------------------------------------------------------------
  int jj=0;
  for(int i=1; i<=ne; i++)
  {
    if(rvar(e(i)))
    {
      jj=-1;
    }
  }
  if (jj==-1)
  {
    "Your ring contains a variable e(" +string(i)+")!";
    return("ERROR");
  }
  if(defined(desiredorder)>1)
  {
    ringre="ring re=" + charstr(r) +",(e(1.." + string(ne) + "),"+
                        varstr(basering) + "),Ws(";
    intvec tempiv=intvec(wl[2]);
    for(i=1;i<=ne;i++)
    {
      ringre=ringre + string((-1)*tempiv[desiredorder[i]]) + ",";
    }
    ringre= ringre  + string(wl[1]) + ");";
  }
  else
  {
    ringre="ring re=" + charstr(r) +",(e(1.." + string(ne) + "),"+ varstr(basering)
                        + "),Ws(" + string((-1)*intvec(wl[2])) + ","
                        + string(wl[1]) + ");";
  }
  execute(ringre);
  module temp=imap(r,t1qh);
  ideal t1qh=mod2id(temp,iv);
  if (defined(watchProgress))
  {
    if (watchProgress[1]!=0)
    {
      option(prot);
      "Protocol output of the groebner computation (quasihomogenous case)";
    }
  }
  ideal t1qhs=std(t1qh);
  if (defined(watchProgress))
  {
    if (watchProgress[1]!=0)
    {
      "groebner computation finished";
      option(noprot);
    }
  }
  ideal t1qhsl=lead(t1qhs);
  module mo=id2mod(t1qhsl,iv);
//--------------------------------------------------------------------------
// Return to the initial ring to compute the kbase and noether there
// (in the new ring t1qh is of course not of dimension 0 but of dimension 3
// so we have to go back)
//--------------------------------------------------------------------------
  setring r;
  module mo=imap(re,mo);
  attrib(mo,"isSB",1);                // mo is monomial ==> SB
  attrib(mo,"isHomog",intvec(wl[2])); // highcorner has to respect the weights
  vector noe=highcorner(mo);
  if (defined(watchProgress))
  {
    if (watchProgress[2]!=0)
    {
      "weights corresponding to the entries of the matrix:";
      wl;
      "leading term of the groebner basis (quasihomogeneous case)";
      mo;
      "noether";
      noe;
    }
  }
//--------------------------------------------------------------------------
// Define the family of curves with the same quasihomogeneous initial
// matrix M, compute T1 and pass again to the ring with the variables e(i)
//--------------------------------------------------------------------------
  if (defined(newname)>1)
  {
    list li=posweight(M,mo,2,newname+"1");
    def rneu=basering;
  }
  else
  {
    list li=posweight(M,mo,2);
  }
  if (size(li)<=1)
  {
    return("ERROR");
  }
  if (defined(watchProgress))
  {
    if(watchProgress[3]!=0)
    {
      "perturbed matrix and weights of the perturbations:";
      li;
    }
  }
  list li2=matrixT1(li[1],3);
  module Mpert=transpose(matrix(ideal(li2[1])));
  module t1pert=li2[2];
  int nv=nvars(rneu)-nvars(r);
  ring rtemp=0,(T(1..nv)),wp(li[2]);
  def reneu=re+rtemp;
  setring reneu;
  module noe=matrix(imap(r,noe));
  ideal noet=mod2id(noe,iv);
  module temp=imap(rneu,t1pert);
  ideal t1pert=mod2id(temp,iv);
//--------------------------------------------------------------------------
// Compute the standard basis and select those generators with leading term
// divisible by some T(i)
//--------------------------------------------------------------------------
  noether=noet[size(noet)];
  if (defined(watchProgress))
  {
    if (watchProgress[1]!=0)
    {
      "protocol output of the groebner command (perturbed case)";
      option(prot);
    }
  }
  ideal t1perts=std(t1pert);
  noether=noet[size(noet)];
  t1perts=interred(t1perts);
  if (defined(Debug))
  {
    if (watchProgress[1]!=0)
    {
      "groebner computation finished (perturbed case)";
      option(noprot);
    }
  }
  ideal templ=lead(t1perts);
  for(int j=1;j<=nv;j++)
  {
    templ=subst(templ,T(j),0);
  }
  ideal mx;
  ideal mt;
  for(j=1;j<=size(t1perts);j++)
  {
    if(templ[j]!=0)
    {
      mx=mx,t1perts[j];
    }
    else
    {
      mt=mt,t1perts[j];
    }
  }
//--------------------------------------------------------------------------
// multiply by the initial ring variables to shift the generators with
// leading term divisible by some T(i) and reduce afterwards
//--------------------------------------------------------------------------
                       // This is obviously no SB, but we have to reduce by
  attrib(mx,"isSB",1); // it and setting isSB suppresses error messages
  noether=noet[size(noet)];
  ideal ker_gen=reduce(mt,mx);
  ideal ovar=var(ne+1),var(ne+2),var(ne+3);
  j=1;
  noether=noet[size(noet)];
  if (defined(watchProgress))
  {
    if (watchProgress[4]!=0)
    {
      "generators of the kernel as a C[T]{x} module:";
      mt;
      "noether:";
      noether;
    }
  }
  int zeros;
  templ=ker_gen;
  while(zeros==0)
  {
    zeros=1;
    templ=templ*ovar;
    templ=reduce(templ,mx);
    if(defined(watchProgress))
    {
      if(watchProgress[4]>1)
      {
        templ;
      }
    }
    if (size(templ)!= 0)
    {
      zeros=0;
      ker_gen=ker_gen,templ;
    }
  }
//-------------------------------------------------------------------------
// kill zero entries, keep only one of identical entries
//-------------------------------------------------------------------------
  ovar=var(1);
  for(i=2;i<=ne;i++)
  {
    ovar=ovar,var(i);
  }
  ker_gen=ker_gen,ovar^2;
  noether=noet[size(noet)];
  ker_gen=simplify(ker_gen,10);
//-------------------------------------------------------------------------
// interreduce ker_gen as a k[T]-module
//-------------------------------------------------------------------------
  intvec mgen=1..(ne+3);
  ideal Mpert=mod2id(imap(rneu,Mpert),iv);
  templ=0;
  for(i=1;i<=nv;i++)
  {
    templ[i]=diff(Mpert[size(Mpert)],T(i));
  }
  templ=templ,ovar^2;
  list retl=subrInterred(templ,ker_gen,mgen);
// Build up the matrix representing L
  module retlm=transpose(retl[2]);
  for(i=1;i<=size(retl[1]);i++)
  {
    if(reduce(retl[1][1,i],std(ovar^2))==0)
    {
      retlm[i]=0;
    }
  }
  retlm=simplify(transpose(simplify(transpose(retlm),10)),10);
  if(defined(watchProgress))
  {
    if(watchProgress[5]>0)
    {
      print(retlm);
    }
  }
  ker_gen=retl[3];
// we define ret=i(L),(delta_j(t_k))_jk
  list ret=id2mod(ker_gen,iv),matrix(retlm);
// cleanups - define what we previously killed
  if(defined(kksave)>1)
  {
    def watchProgress=kksave;
    export watch Progress;
  }
  option(set,optvec);
// make sure that the exported ring has the right name
  if (defined(newname)>1)
  {
    def `newname`=reneu;
    setring `newname`;
    export `newname`;
    keepring `newname`;
  }
  else
  {
    export reneu;
    keepring reneu;
  }
  return(ret[2]);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),ds;
  matrix M[3][2]=z-x^7,0,y^2,z,x^9,y;
  KSpencerKernel(M,"ar");
  nameof(basering);
  basering;
}
///////////////////////////////////////////////////////////////////////////

proc mod2id(matrix M,intvec vpos)
"USAGE:     @code{mod2id(M,vpos)}; @code{M} matrix, @code{vpos} intvec
ASSUME:    @code{vpos} is an integer vector such that @code{gen(i)} corresponds
           to @code{var(vpos[i])}
           @*the basering contains variables @code{var(vpos[i])} which do not 
           occur in @code{M}
NOTE:      this procedure should be used in the following situation:
           one wants to pass to a ring with new variables, say e(1),..,e(s),
           which correspond to the components gen(1),..,gen(s) of the
           module @code{M} such that e(i)*e(j)=0 for all i,j
           the new ring should already exist and be the current ring
RETURN:    ideal i in which each gen(i) from the module is replaced by 
           @code{var(vpos[i])} and all monomials @code{var(vpos[i])*var(vpos[j])} 
           have been added to the generating set of i
EXAMPLE:   @code{example mod2id}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//----------------------------------------------------------------------
// define the ideal generated by the var(vpos[i]) and set up the matrix
// for the multiplication
//----------------------------------------------------------------------
  ideal vars=var(vpos[1]);
  for(int i=2;i<=size(vpos);i++)
  {
    vars=vars,var(vpos[i]);
  }
  matrix varm[1][size(vpos)]=vars;
  if (size(vpos) > nrows(M))
  {
    matrix Mt[size(vpos)][ncols(M)];
    Mt[1..nrows(M),1..ncols(M)]=M;
    kill M;
    matrix M=Mt;
  }
//----------------------------------------------------------------------
// define the desired ideal
//----------------------------------------------------------------------
  ideal ret=vars^2,varm*M;
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3));
  module mo=x*gen(1)+y*gen(2);
  intvec iv=2,1;
  mod2id(mo,iv);
}
////////////////////////////////////////////////////////////////////////

proc id2mod(ideal i,intvec vpos)
"USAGE:     @code{id2mod(I,vpos)}; @code{I} ideal, @code{vpos} intvec
NOTE:      * use this procedure only makes sense if the ideal contains
             all @code{var(vpos[i])*var(vpos[j])} as monomial generators and
             all other generators are linear combinations of the
             @code{var(vpos[i])} over the ring in the other variables
           @* * this is the inverse procedure to @code{mod2id} and should be 
             applied only to ideals created by @code{mod2id} using the same 
             intvec @code{vpos} (possibly after a standard basis computation)
RETURN:    module corresponding to the ideal by replacing @code{var(vpos[i])} by
           gen(i) and omitting all generators @code{var(vpos[i])*var(vpos[j])}
EXAMPLE:   @code{example id2mod}; shows an example"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------
// Initialization
//---------------------------------------------------------------------
  int n=size(i);
  int v=size(vpos);
  matrix tempmat;
  matrix mm[v][n];
//---------------------------------------------------------------------
// Conversion
//---------------------------------------------------------------------
  for(int j=1;j<=v;j++)
  {
    tempmat=coeffs(i,var(vpos[j]));
    mm[j,1..n]=tempmat[2,1..n];
  }
  for(j=1;j<=v;j++)
  {
    mm=subst(mm,var(vpos[j]),0);
  }
  module ret=simplify(mm,10);
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3));
  ideal i=e(2)^2,e(1)*e(2),e(1)^2,e(1)*y+e(2)*x;
  intvec iv=2,1;
  id2mod(i,iv);
}
///////////////////////////////////////////////////////////////////////

proc subrInterred(ideal mon, ideal sm, intvec iv)
"USAGE:    @code{subrInterred(mon,sm,iv)};
          @*@code{sm}:   ideal in a ring r with n + s variables,
                e.g. x_1,..,x_n and t_1,..,t_s
          @*@code{mon}:  ideal with monomial generators (not divisible by
                one of the t_i) such that sm is contained in the module
                @code{k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)]}
          @*@code{iv}:   intvec listing the variables which are supposed to be 
                used as x_i
RETURN:   interreduced system of generators of sm seen as a submodule
          of @code{k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)]}
EXAMPLE:  @code{example subrInterred}; shows an example
"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//-----------------------------------------------------------------------
// check that mon is really generated by monomials
// and sort its generators with respect to the monomial ordering
//-----------------------------------------------------------------------
  int err;
  for(int i=1;i<=ncols(mon);i++)
  {
    if ( size(mon[i]) > 1 )
    {
      err=1;
    }
  }
  if (err==1)
  {
    "mon has to be generated by monomials";
    return("ERROR");
  }
  intvec sv=sortvec(mon);
  ideal mons;
  for(i=1;i<=size(sv);i++)
  {
    mons[i]=mon[sv[i]];
  }
  ideal itemp=maxideal(1);
  for(i=1;i<=size(iv);i++)
  {
    itemp=subst(itemp,var(iv[i]),0);
  }
  itemp=simplify(itemp,10);
  def r=basering;
  string tempstr="ring rtemp=" + charstr(basering) + ",(" + string(itemp)
                               + "),(C,dp);";
//-----------------------------------------------------------------------
// express m in terms of the generators of mon and check whether m
// can be considered as a submodule of k[t_1,..,t_n]*mon
//-----------------------------------------------------------------------
  module motemp;
  motemp[ncols(sm)]=0;
  poly ptemp;
  int j;
  for(i=1;i<=ncols(mons);i++)
  {
    for(j=1;j<=ncols(sm);j++)
    {
      ptemp=sm[j]/mons[i];
      motemp[j]=motemp[j]+ptemp*gen(i);
    }
  }
  for(i=1;i<=size(iv);i++)
  {
    motemp=subst(motemp,var(iv[i]),0);
  }
  matrix monmat[1][ncols(mons)]=mons;
  ideal dummy=monmat*motemp;
  for(i=1;i<=size(sm);i++)
  {
    if(sm[i]-dummy[i]!=0)
    {
      "the second argument is not a submodule of the assumed structure";
      return("ERROR");
    }
  }
//----------------------------------------------------------------------
// do the interreduction and convert back
//----------------------------------------------------------------------
  execute(tempstr);
  def motemp=imap(r,motemp);
  motemp=interred(motemp);
  setring r;
  kill motemp;
  def motemp=imap(rtemp,motemp);
  list ret=monmat,motemp,monmat*motemp;
  for(i=1;i<=ncols(ret[2]);i++)
  {
    ret[2][i]=cleardenom(ret[2][i]);
  }
  for(i=1;i<=ncols(ret[3]);i++)
  {
    ret[3][i]=cleardenom(ret[3][i]);
  }
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),dp;
  ideal i=x^2+x*y^2,x*y+x^2*y,z;
  ideal j=x^2+x*y^2,x*y,z;
  ideal mon=x^2,z,x*y;
  intvec iv=1,3;
  subrInterred(mon,i,iv);
  subrInterred(mon,j,iv);
}
////////////////////////////////////////////////////////////////////////