Changeset 6f2edc for Singular/LIB – Singular

Singular/LIB/deform.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: deform.lib,v 1.1.1.1 1997-04-25 15:13:25 obachman Exp $
 //(BM+GMG)
+// $Id: deform.lib,v 1.2 1997-04-28 19:27:15 obachman Exp $
+//(BM/GMG, last modified 22.06.96)
 ///////////////////////////////////////////////////////////////////////////////
 LIBRARY:  deform.lib    PROCEDURES FOR COMPUTING MINIVERSAL DEFORMATION
  miniversal(id[,deg]);  miniversal deformation of an isolated singularity id
   SUB-PROCEDURES        used by main procedure:
   apply_col(A,B);       put A into col-nf and apply same col-operations to B
 …
 LIB "inout.lib";
 LIB "general.lib";
 LIB "sing.lib";
+LIB "sing.lib";
 LIB "matrix.lib";
 ///////////////////////////////////////////////////////////////////////////////
 proc miniversal (ideal id,list #)
 USAGE:   miniversal(id[,d,na,va,o,iv]); id=ideal, d=integer,
+USAGE:   miniversal(id[,d,na,va,o,iv]); id=ideal, d=integer,
          na,va,o=strings, iv=intvec of positive integers
 COMUPTE: miniversal deformation of id up to degree d (default d=100)
 …
          the basering by new variables given by va (deformation parameters).
          -- The new vars come before the old vars
          -- The characteristic of `na`is the characteristic of the basering.
          -- The new vars are derived from va. If va is a single letter, say
             va="T", and if n<=26 then T and the following n-1 letters from
             T..Z..T (resp. T(1..n) if n>26) are taken as additional variables.
             If va is a single letter followed by (, say va="x(", the new
+         -- The characteristic of `na`   is the characteristic of the basering.
+         -- The new vars are derived from va. If va is a single letter, say
+            va="T", and if n<=26 then T and the following n-1 letters from
+            T..Z..T (resp. T(1..n) if n>26) are taken as additional variables.
+            If va is a single letter followed by (, say va="x(", the new
             variables are x(1),...,x(n) (default va="A").
          -- The ordering is the product ordering between the ordering of r and
             an ordering derived from `o`, which has to be local!! (default:
             o="ds") [and iv (a weight vector)].
+            an ordering derived from `o`, which has to be local!! (default:
+            o="ds") [and iv (a weight vector)].
             Type 'help extendring' for a more detailed explanation of the
             ordering
          -- Even if na,va,o are given, d and/or iv may be ommited. Then the
+            ordering
+         -- Even if na,va,o are given, d and/or iv may be ommited. Then the
             default values d=100, iv=0 (i.e. all weights = 1) are used.
          The procedure creates also two ideals:
             ideal jetJ - defining the miniversal base space (in `na`)
             ideal jetF - defining miniversal total space (in `na`)
+NOTE:    int printlevel=2;  shows what is going on
+         int printlevel=3;  shows also memory usage
+         This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
+NOTE:    printlevel >=0: display dimT1,T2 and explain created objects (default)
+         printlevel >=1: show partial + final result during computation
+         printlevel >=2: show also memory and time usage
+         printlevel >=3: test and show obstructions
+         printlevel >=4: create a file 'minbaseout' and (over) write part of
+                         ideal of miniversal base up to current degree into it
+         This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
          start with @ (see the file HelpForProc)
 EXAMPLE: example miniversal; shows an example
+EXAMPLE: example miniversal; shows an example
+{
 //------- initialisation ------------------------------------------------------
+   int @d,@deg,@t1,@t2,@colR,@noObstr;
+   int @d,@deg,@t1,@t2,@colR,@noObstr,@j;
+   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
    intvec @iv,@jv;
    string @na,@va,@o;
 …
    if( size(#)==4 ) { @deg=#[1]; @na=#[2];  @va=#[3]; @o=#[4];}
    if( size(#)==5 ) { @deg=#[1]; @na=#[2];  @va=#[3]; @o=#[4]; @iv=#[5]; }
    if( find(@o,"s")==0 )
+   if( find(@o,"s")==0 )
    { "// ordering must be an s-ordering, please change!"; return();}
   def @Pn = basering;
    string @ords = ordstr(@Pn);
+  def @Pn = basering;
+   string @ords = ordstr(@Pn);
    id = simplify(id,10);
    int @rowR = size(id);
    if( @rowR<=1 )
+   {
+   {
       "// hypersurface, use proc deform from sing.lib";
       return();
+   }
+   }
 //------- change ordering if not correct --------------------------------------
    @t1=1;
    for( @d=1;@d<=nvars(@Pn);@d++ ) { @t1=@t1*(lead(1+var(@d))==var(@d)); }
+   for( @d=1;@d<=nvars(@Pn);@d=@d+1 ) { @t1=@t1*(lead(1+var(@d))==var(@d)); }
    if( @t1==0 )
+   {
       if( @ords[size(@ords)]!="c" and @ords[size(@ords)]!="C" )
+      {
+      if( @ords[size(@ords)]!="c" and @ords[size(@ords)]!="C" )
+      {
          if( @ords[1]=="c" ) { @ords=@ords[3,size(@ords)-2]+",c"; @t1=1;}
          if( @ords[1]=="C" ) { @ords=@ords[3,size(@ords)-2]+",C"; @t1=1;}
+      }
       if( @t1==1 )
+      {
+      if( @t1==1 )
+      {
          changeord("@On",@ords,@Pn);
          ideal id  = imap(@Pn,id);
 …
    if( defined(@On)==0 ) { def @On=@Pn; setring @On; }
 //-------  reproduce T12 -------------------------------------------------------
    list   Ls   = T12(id,1);
    matrix Ro   = Ls[4];                         //syz(i)
    matrix InfD = Ls[3];                         //matrix of inf. deformations
    matrix PreO = Ls[5];                         //present. mat of Syz/Kos^*
    module PreT = Ls[2];                         //present. module of modT2
    @t1 = Ls[8];                                 //vdim of T1
    @t2 = Ls[9];                                 //vdim of T2
+   list   Ls   = T12(id,1);
+   matrix Ro   = Ls[6];                         //syz(i)
+   matrix InfD = Ls[5];                         //matrix of inf. deformations
+   matrix PreO = Ls[7];                         //present. mat of Syz/Kos^*
+   module PreT = Ls[9];                         //present. module of modT2
+   @t1 = Ls[3];                                 //vdim of T1
+   @t2 = Ls[4];                                 //vdim of T2
    kill Ls;
    dbpri(2,"","// ___ matrix of infinitesimal deformations:",InfD);
    @colR = ncols(Ro);
+   dbprint(p-1,"","// ___ matrix of infinitesimal deformations:",InfD);
+   @colR = ncols(Ro);
    ideal i0 = std(id);
   qring @Ox = i0;                               //ring of singularity to deform
    matrix Cup,lCup; module PreT;
+   matrix Cup,lCup;
    ideal testid;
    matrix Ro   = fetch(@On,Ro);
    matrix PreO = fetch(@On,PreO);
+   module PreT = fetch(@On,PreT);
 //---- create new ring with @t1=dim T1 additional variables and initialize ----
 …
    export jetF;
    matrix Fo = matrix(jetF);                    //initial equations
+   matrix Rs = imap(@On,Ro);                    //deformed syzygies
+   matrix Ro = imap(@On,Ro);
+   matrix Rs = imap(@On,Ro);                    //deformed syzygies
    ideal  jetJ;                                 //(jet)ideal of minversal defor
    export jetJ;
 …
    jetF= Fs;
    F_R = Fs*Rs;
    if (@t2<=0) { @d=0; }                         //finished, if "T2=0"
+   if (@t2<=0) { @d=0; }                         //finished, if "T2=0"
 //------- start the loop ------------------------------------------------------
    for (@d=1;@d<=@deg;@d++)
+   {
       dbpri(2,"","// ___ start computation in degree "+string(@d)+":");
       dbpri(3,"memory="+string(kmemory())+"k");
+   for (@d=1;@d<=@deg;@d=@d+1)
+   {
+      dbprint(p-1,"","// ___ start computation in degree "+string(@d)+":");
+      dbprint(p-2,"// memory = "+string(kmemory())+"k");
 //------- lift relation to next degree ----------------------------------------
       F_r = reduce_s(F_R,Jo,@d+1);
+      F_r = reduce_s(F_R,Jo,@d+1);
       Cup = matrix(jet(F_r,@d,@jv),1,@colR);
       Rn  = (-1)*lift(Fo,Cup);
       Rs  = Rs + Rn;
       F_R = F_R + Fs*Rn;
+      F_R = F_R + Fs*Rn;
 //------- test: already finished? ---------------------------------------------
       testid = simplify(reduce(ideal(F_R),Jo),10);
       if (testid[1]==0)
+      {
+         "// computation finished in degree "+string(@d);
+         if( @d==@deg )
+         {"// degree bound reached, result may not yet be complete!";}
+      {  dbprint(p,"// ___ computation finished in degree "+string(@d));
+         if( @d==@deg )
+         { dbprint(p,"// ___ degree bound reached, result may not yet be complete!");}
          break;
+      }
 //------- compute obstruction-matrix  -----------------------------------------
       F_r = reduce_s(F_R,Jo,@d+1);
       Cup = matrix(jet(F_r,@d+1,@jv),1,@colR);
+      F_r = reduce_s(F_R,Jo,@d+1);
+      Cup = matrix(jet(F_r,@d+1,@jv),1,@colR);
       Test= Cup;
       dbpri(2,"","// ___ obstruction vector:",ideal(Cup));
+      dbprint(p-3,"","// ___ obstruction vector:",ideal(Cup));
       Cup,Mon = coef_ideal(Cup,@t1);
 //------- express obstructions in kbase of T2  --------------------------------
 …
       Cup  = imap(`@na`,Cup);
       lCup = lift(PreO,Cup);
-      PreT = fetch(@On,PreT);
       lCup = lift_kbase(lCup,PreT);
       @t2   = nrows(lCup);
       dbpri(2,"","// ___ obstructions in kbase of T2:",lCup);
       testid = simplify(ideal(lCup),10);           // test no obstructions
+      dbprint(p-3,"","// ___ obstructions in kbase of T2:",lCup);
+      testid = simplify(ideal(lCup),10);               // test no obstructions
       if (testid[1]==0)
+      { @noObstr=1; } else { @noObstr=0; }
+//------- compute ideal of minversal(its k-jet) -------------------------------
+      { @noObstr=1;dbprint(p-3,"// ___ no obstruction"); } else { @noObstr=0; }
+      @j=size(module(gauss_col(lCup)));                // test:full obstruction
+      if (@j==ncols(lCup))
+      {  dbprint(p,"","// nothing to lift!",
+         "// ___ miniversal base, defined by jetJ, is a fat point!");
+         break;
+      }
+//------- compute ideal of minversal base (its k-jet) -------------------------
    setring `@na`;
       if (@noObstr==0)                              //case of non-zero obstr.
+      if (@noObstr==0)                              //case of non-zero obstr.
+      {
          lCup = imap(@Ox,lCup);
          Jo   = lCup*transpose(Mon);
+         jetJ = matrix(jetJ,1,@t2)+matrix(Jo,1,@t2);
+         dbpri(2,"","// ___ degree-"+string(@d+1)+"-part of ideal of miniversal base"+":",Jo);
+         jetJ = matrix(jetJ,1,@t2)+matrix(Jo,1,@t2);
+         dbprint(p-1,"","// ___ degree-"+string(@d+1)+"-part of ideal of miniversal base"+":",Jo);
+         if( p-1>=4 )
+         { write (">minbaseout","// part of ideal of miniversal base up to degree <= "+string(@d+1)+":",jetJ); }
          Jo = std(jetJ);
+//------- choose a defining system --------------------------------------------
+         @iv,Cup = defining_system(lCup,Cup);
+         dbpri(2,"","// ___ number of cols of defining system:",@iv);
+//------- lift the equations --------------------------------------------------
+         if (sum(@iv)==0)
+         {
+            "// nothing to lift!";
+            "// miniversal base, defined by jetJ, is a fat point!";break;
+         }
+   setring @Ox;
+         Cup = imap(`@na`,Cup);
+         Cup = submat(Cup,1..nrows(Cup),@iv);
+         dbpri(2,"","// ___ matrix of defining system:",Cup);
+      }
+      else                                         // case of zero obstructions
+      {
+   setring @Ox;
+         Cup = imap(`@na`,Cup);
+      }
+      Cup = lift(transpose(Ro),module(Cup));
+   setring `@na`;
+      Cup = imap(@Ox,Cup);
+      if (@noObstr==0)
+      {  Mon = submat(Mon,1..nrows(Mon),@iv); }
+      Fn  = (-1)*transpose(Cup*transpose(Mon));
+      }
+      F_r = reduce_s(F_R,Jo,@d+1);
+      Cup = matrix(jet(F_r,@d+1,@jv),1,@colR);
+//---------------- repeat test: jetJ ok in deg d+1? --------------------------
+      if( (p-1>=3) && (@noObstr==0) )
+      {
+         lCup,Mon = coef_ideal(Cup,@t1);
+      setring @Ox;
+         Cup  = imap(`@na`,Cup);
+         lCup = lift(PreO,Cup);
+         lCup = lift_kbase(lCup,PreT);
+         dbprint(p-3,"","// ____ test: jetJ ok iff all entries are 0",lCup);
+      setring `@na`;
+      }
+//---------------- lift equations F -----------------------------------------
+      if (defined(Qrg)) {kill Qrg;}
+  qring Qrg = std(ideal(Fo));
+      def Ro=fetch(`@na`,Ro);
+      def Cup=fetch(`@na`,Cup);
+      def Fn = lift(transpose(Ro),transpose(Cup));
+      Fn=(-1)*transpose(Fn);
+  setring `@na`;
+      Fn  = fetch(Qrg,Fn);
       Fs  = Fs+Fn;
       F_R = F_R+Fn*Rs;
+      F_R = F_R+Fn*Rs;
       jetF = matrix(Fs);
       dbpri(2,"","// ___ degree-"+string(@d+1)+"-part of deformed equations:",Fn);
+   }
+      dbprint(p-1,"","// ___ degree-"+string(@d+1)+"-part of deformed equations:",Fn);
+   }
 //---------  end loop and final output ---------------------------------------
+  "";
   "// ___ Equations of miniversal base space ___";jetJ; "";
+  "// ___ Equations of miniversal total space ___";jetF; "";
   "// Result belongs to ring",@na,"(total space of miniversal deformation).";
   "// Make",@na,"the basering and list objects defined in",@na,"by typing:";
   "   setring",@na,"; show("+@na+");";
   "   listvar(ideal);";
   kill @On;
+dbprint(p-1,"","// ___ Equations of miniversal base space ___",jetJ,
+            "","// ___ Equations of miniversal total space ___",jetF);
+dbprint(p,"","// Result belongs to ring "+@na+".",
+       "// Equations of total space of miniversal deformation are ",
+       "// given by jetF, equations of miniversal base space by jetJ.",
+       "// Make "+@na+" the basering and list objects defined in "+@na+" by typing:",
+       "   setring "+@na+"; show("+@na+");","   listvar(ideal);");
+  kill @On;
   return();
+}
 example
+example
 { "EXAMPLE:"; echo = 2;
+   ring r1=0,(x,y,z,u,v),ds;
+   matrix m[2][4]=x,y,z,u,y,z,u,v;
+   ideal i=minor(m,2);          //cone over rational normal curve of degree 4
+   int p          = printlevel;
+   ring r1        = 0,(x,y,z,u,v),ds;
+   matrix m[2][4] = x,y,z,u,y,z,u,v;
+   ideal i        = minor(m,2);          //cone over rational normal curve of degree 4
    miniversal(i,"R","T(");
+   // hit return-key to continue;
+   // pause;
+   ring  r = 0,(x,y,z),ds;
+   ideal i = x2,xy,yz,zx;
+   printlevel = 2;
+   miniversal(i);"";
+   kill printlevel;
+// NOTE: rings R and Ont are still alive!
+}
+///////////////////////////////////////////////////////////////////////////////
+proc apply_col (matrix A, matrix B)
+   setring R;"";
+   // ___ Equations of miniversal base space ___:
+   jetJ;"";
+   // ___ Equations of miniversal total space ___:
+   jetF;"";
+   ring  r        = 0,(x,y,z),ds;
+   ideal i        = x2,xy,yz,zx;
+   printlevel     = 3;
+   miniversal(i);"";
+   printlevel     = p;
+   // NOTE: rings R and Ont are still alive!
+}
+///////////////////////////////////////////////////////////////////////////////
+proc apply_col (matrix A, matrix B)
 USAGE:   apply_col(A,B);  A,B=matrices
 ASSUME:  A = constant matrix in row-reduced (upper triangular) normal form,
+ASSUME:  A = constant matrix in row-reduced (upper triangular) normal form,
          B = matrix of same size
 COMUPTE: apply to B those col-operations which reduce A into col-reduced nf
 RETURN:  two transformed matrices: col-reduced A, transformed B
 EXAMPLE: example apply_col; shows an example
+EXAMPLE: example apply_col; shows an example
+{
    int i,j,k;
 …
    matrix C  = concat(transpose(A),transpose(B));
    module mC = transpose(C);
    for( k=1;k<=r;k++ )
+   for( k=1;k<=r;k=k+1 )
+   {
       j=1;
       while( C[j,k]==0 && j<c ) { j++; }
       for( i=j+1;i<=c;i++ )
+      while( C[j,k]==0 && j<c ) { j=j+1; }
+      for( i=j+1;i<=c;i=i+1 )
+      {
          m = C[i,k];
          mC[i] = mC[i]-m*mC[j];
+         mC[i] = mC[i]-m*mC[j];
+      }
+   }
 …
    return(transpose(a),transpose(b));
+}
 example
+example
 { "EXAMPLE:"; echo = 2;
    ring R=0,(x,y,z),dp;
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc defining_system (matrix A,matrix B)
+proc defining_system (matrix A,matrix B)
 USAGE:   defining_system(A,B);  A,B=matrices
 ASSUME:  A a constant matrix
 …
 RETURN:  two objects: intvec iv, matrix M (the transformed matrix B)
          The columns of M with index from iv are a defining sytem
 EXAMPLE: example defining_system; shows an example
+EXAMPLE: example defining_system; shows an example
+{
    int    k,l;
 …
    int rg = ncols(A);
    A,B    = apply_col(A,B);                // special columne-reduction
    for( k=1;k<=rg;k++ )                  // collect zero-cols of B
+   {
       if( A[k]==0) {l++;iv[l]=k;}        // test if kth column is 0
+   for( k=1;k<=rg;k=k+1 )                  // collect zero-cols of B
+   {
+      if( A[k]==0) {l=l+1;iv[l]=k;}        // test if kth column is 0
    }                                       // collect indices of 0-columns in iv
    return(iv,B);
+}
 example
+example
 { "EXAMPLE:"; echo = 2;
    ring R=0,(x,y,z),dp;
 …
   int d,k;
   ideal j0 = std(j);
   for (k=1;k<=m;k++)
+  for (k=1;k<=m;k=k+1)
+  {
     if (deg(i[k])>=0)
 …
   return(i);
+}
 example
+example
 { "EXAMPLE:"; echo = 2;
    ring r = 0,(x,y),ds;
    poly f = x7+y7+(x-y)^2*x^2*y^2;
    ideal j = jacob(f);
    reduce_s(f,j,10);
+   reduce_s(f,j,10);
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 proc lift_kbase (N, M)
 USAGE:   lift_kbase(N,M); N,M=poly/ideal/vector/module
 RETURN:  matrix A, coefficient matrix expressing N as linear combination of
          k-basis of M. Let the k-basis have k elements and A c columns.
+RETURN:  matrix A, coefficient matrix expressing N as linear combination of
+         k-basis of M. Let the k-basis have k elements and size(N)=c columns.
          Then A satisfies:
              matrix(reduce(N,std(M)),k,c) = matrix(kbase(std(M)))*A
 ASSUME:  dim(M)=0 and the monomial ordering is a well ordering or the last
+ASSUME:  dim(M)=0 and the monomial ordering is a well ordering or the last
          block of the ordering is c or C
 EXAMPLE: example lift_kbase; shows an example
 …
 //------- check wether ordering is correct ------------------------------------
    k=1;
    for( l=1;l<=nvars(basering);l++ ) { k=k*(lead(1+var(l))==var(l)); }
+   for( l=1;l<=nvars(basering);l=l+1 ) { k=k*(lead(1+var(l))==var(l)); }
    if( k==0 )
+   {
       if( ords[size(ords)]!="c" and ords[size(ords)]!="C" )
+      {
+      if( ords[size(ords)]!="c" and ords[size(ords)]!="C" )
+      {
          "// change ordering!";
          "// ordering "+ordstr(basering)+" not implemented for this proc";
          return();
+         "// ordering "+ordstr(basering)+" not implemented for this proc";
+         return();
+      }
+   }
 …
    if( d<1 )
    { "// second argument in `lift_kbase` has vdim",d; return(); }
 //----------  compute kbase and reduce(N,M) -----------------------------------
+//----------  compute kbase and reduce(N,M) -----------------------------------
    kb = kbase(M);
    col = ncols(N);
    N = reduce(N,M);
+   N = matrix(N,nrows(N),col);
 //----------  collecting coefficients of reduce(N,M) --------------------------
    matrix result[d][col];
 …
+      {
          for( k=1;k<=d;k=k+1 )
+         {
+         {
             p = kb[k];
             q = lead(v);
+            q = lead(v);
             if( size(p-q)<2 )
+            {
 …
                else { k=0; }
+            }
+         }
+         }
+      }
+   }
 //---------  final test -------------------------------------------------------
    testm = matrix(N,nrows(kb),ncols(result))- matrix(kb)*result;
    if( size(module(testm))!=0 )
+   {
+   if( size(module(testm))!=0 )
+   {
       "// proc `lift_kbase` did'nt work correctly!";
       "// Please inform tthe authors";
 …
+}
 example
+{
+  "EXAMPLE:";     echo=2;
+{"EXAMPLE:";     echo=2;
   ring R=0,(x,y),ds;
   module M=[x2,xy],[y2,xy],[0,xx],[0,yy];
 …
 ASSUME:  M=matrix with only one row and without any constant term
 COMPUTE: coef_matrices with respect to first s variables
 RETURN:  2 matrices:
+RETURN:  2 matrices:
          matrix of coefficients (each column is formed by the coefficients
            of M with respect to some monomial)
          row-matrix of corresponding monomials
+           of M with respect to some monomial)
+         row-matrix of corresponding monomials
 EXAMPLE: example coef_ideal; shows an example
+{
   int   k,l,n,z;
   int   cM = ncols(M);
   ideal flatM = M;
+  ideal flatM = M;
   ideal monId,flat;
   poly  mon = product(maxideal(1),1..s);
 //--------- collect all monomials (!=1) ---------------------------------------
   for (k=1;k<=cM;k++)
+  for (k=1;k<=cM;k=k+1)
+  {
     matrix mci(k) = coef(flatM[k],mon);
 …
     if (flat[1]!=1)
     { monId = monId,flat;}
+  }
+  }
   monId = monId+ideal(0);
   k=size(monId);
   matrix BIG[cM][k];
 //---------  create coef_matrices  --------------------------------------------
   for (n=1;n<=k;n++)
+  for (n=1;n<=k;n=n+1)
+  {
     for (l=1;l<=cM;l++)
+    for (l=1;l<=cM;l=l+1)
+    {
       for (z=1;z<=ncols(mci(l));z++)
+      for (z=1;z<=ncols(mci(l));z=z+1)
+      {
         if(mci(l)[1,z]==monId[n])
         { BIG[l,n] = mci(l)[2,z];}
+      }
+      }
+    }
+  }
+  }
   return(BIG,matrix(monId));
+}
 example
+}
+example
 { "EXAMPLE:"; echo = 2;
   ring Z = 0,(A,B,C,x,y,z),ds;
   int  s = 3;
+  int  s = 3;
   matrix M[1][4]=A+yB,2C,3AA,4BB+5CC;
   print(M);
 …
+}
 ///////////////////////////////////////////////////////////////////////////////
-----------
-"example in r1: i=cone rational normal curve d=4";
-int d=4;
-ring r1=0,(x,y,z,u,v),ds;
-matrix m[2][4]=x,y,z,u,y,z,u,v;
-ideal i=minor(m,2);
-i=minbase(i);
-i;pause;
-int t=timer;miniversal(i);timer-t;
-----------
-"example: in r4 i=cone rational normal curve d=5";
-int d=5;
-ring s=0,(x,y,z,u,v,w),ds;
-matrix m[2][5]=x,y,z,u,v,y,z,u,v,w;
-ideal i=minor(m,2);
-i=minbase(i);
-i;pause;
-----------
-"Example: in r5 i=L_n^n,   n=4;";
-ring r5=0,(x,y,z,u),ds;
-ideal i;
-i=xy,xz,xu,yz,yu,zu;
-i;pause;
-----------
-"Example 1 : cyclic quotient  in ws
-      (type setring r1;sud(i);)";
-ring r1=0,(x,y,z,u,v),ws(4,3,2,3,4);
-ideal i=xz-y2,yz2-xu,xv-yzu,yu-z3,z2u-yv,zv-u2;
-i;
-"Example 2: same in wp
-      (ringr r2)";
-ring r2=0,(x,y,z,u,v),wp(4,3,2,3,4);
-ideal i=xz-y2,yz2-xu,xv-yzu,yu-z3,z2u-yv,zv-u2;
-i;
-"Example 3: same in ls";
-ring r3=0,(x,y,z,u,v),ls;
-ideal i=xz-y2,yz2-xu,xv-yzu,yu-z3,z2u-yv,zv-u2;
-i;
-"Example 4: by chance for testing";
-ring r4=0,(x,y,z),ds;
-ideal i=xy,yz,xz+y3,x2+y2+z3;
-i;

Singular/LIB/elim.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: elim.lib,v 1.1.1.1 1997-04-25 15:13:25 obachman Exp $
 //system("random",787422842);
 //(GMG)
+// $Id: elim.lib,v 1.2 1997-04-28 19:27:16 obachman Exp $
+// system("random",787422842);
+// (GMG, last modified 22.06.96)
 ///////////////////////////////////////////////////////////////////////////////
 LIBRARY:  elim.lib      PROCEDURES FOR ELIMINATIOM, SATURATION AND BLOWING UP
  blowup0(j[,s1,s2]);    create presentation of blownup ring of ideal j
+ blowup0(j[,s1,s2]);    create presentation of blownup ring of ideal j
  elim(id,n,m);          variable n..m eliminated from id (ideal/module)
  elim1(id,p);           p=product of vars to be eliminated from id
+ elim1(id,p);           p=product of vars to be eliminated from id
  nselect(id,n[,m]);     select generators not containing nth [..mth] variable
  sat(id,j);             saturated quotient of ideal/module id by ideal j
+ sat(id,j);             saturated quotient of ideal/module id by ideal j
  select(id,n[,m]);      select generators containing nth [..mth] variable
            (parameters in square brackets [] are optional)
 …
 proc blowup0 (ideal j,list #)
 USAGE:   blowup0(j[,s1,s2]); j ideal, s1,s2 nonempty strings
 CREATE:  Create a presentation of the blowup ring of j
+USAGE:   blowup0(j[,s1,s2]); j ideal, s1,s2 nonempty strings
+CREATE:  Create a presentation of the blowup ring of j
 RETURN:  no return value
 NOTE:    s1 and s2 are used to give names to the blownup ring and the blownup
          ideal (default: s1="j", s2="A")
          Assume R = char,x(1..n),ord is the basering of j, and s1="j", s2="A"
          then the procedure creates a new basering with name Bl_jR
+         then the procedure creates a new ring with name Bl_jR
          (equal to R[A,B,...])
                Bl_jR = char,(A,B,...,x(1..n)),(dp(k),ord)
          with k=ncols(j) new variables A,B,... and ordering wp(d1..dk) if j is
+               Bl_jR = char,(A,B,...,x(1..n)),(dp(k),ord)
+         with k=ncols(j) new variables A,B,... and ordering wp(d1..dk) if j is
          homogeneous with deg(j[i])=di resp. dp otherwise for these vars.
          If k>26 or size(s2)>1, say s2="A()", the new vars are A(1),...,A(k).
+         If k>26 or size(s2)>1, say s2="A()", the new vars are A(1),...,A(k).
          Let j_ be the kernel of the ring map Bl_jR -> R defined by A(i)->j[i],
          x(i)->x(i), then the quotient ring Bl_jR/j_ is the blowup ring of j
          in R (being isomorphic to R+j+j^2+...). Moreover the procedure creates
+         a std basis of j_ with name j_ and Bl_jR as basering.
+EXAMPLE: example blowup0; shows an example
+         a std basis of j_ with name j_ in Bl_jR.
+         This proc uses 'execute' or calls a procedure using 'execute'.
+DISPLAY: printlevel >=0: explain created objects (default)
+EXAMPLE: example blowup0; shows examples
+{
    string bsr = nameof(basering);
    def br = basering;
    string cr, vr, o = charstr(br), varstr(br), ordstr(br);
    int n, k = nvars(br), ncols(j);
    int i;
+   string cr,vr,o = charstr(br),varstr(br),ordstr(br);
+   int n,k,i = nvars(br),ncols(j),0;
+   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
 //---------------- create coordinate ring of blown up space -------------------
    if( size(#)==0 ) { #[1] = "j"; #[2] = "A"; }
 …
    if( k<=26 and size(#[2])==1 ) { string nv = A_Z(#[2],k)+","; }
    else { string nv = (#[2])[1]+"(1.."+string(k)+"),"; }
    if( ishomog(j) )
+   {
+   if( is_homog(j) )
+   {
       intvec v=1;
       for( i=1; i<=k; i++) { v[i+1]=deg(j[i]); }
+      for( i=1; i<=k; i=i+1) { v[i+1]=deg(j[i]); }
       string nor = "),(wp(v),";
+   }
 …
 //---------- map to new ring, eliminate and create blown up ideal -------------
    ideal j=imap(br,j);
    for( i=1; i<=k; i++) { j[i]=var(1+i)-t*j[i]; }
+   for( i=1; i<=k; i=i+1) { j[i]=var(1+i)-t*j[i]; }
    j=eliminate(j,t);
    v=v[2..size(v)];
 …
    export basering;
    export `#[1]+"_"`;
+   keepring basering;
+   //keepring basering;
+   setring br;
 //------------------- some comments about usage and names  --------------------
+   if( voice ==2 )
+   {
+"// NOTE:";
+"// basering is now Bl_"+#[1]+bsr+" (equal to "+bsr+"["+nv[1,size(nv)-1]+"])";
+"// it contains the ideal "+#[1]+"_ , such that";
+"//             Bl_"+#[1]+bsr+"/"+#[1]+"_ is the blowup ring";
+"// For blowing-up another ideal in "+bsr+" type first:";
+"//             setring "+bsr+";";
+   }
+dbprint(p,"",
+"// The proc created the ring Bl_"+#[1]+bsr+" (equal to "+bsr+"["+nv[1,size(nv)-1]+"])",
+"// it contains the ideal "+#[1]+"_ , such that",
+"//             Bl_"+#[1]+bsr+"/"+#[1]+"_ is the blowup ring",
+"// show(Bl_"+#[1]+bsr+"); shows this ring.",
+"// Make Bl_"+#[1]+bsr+" the basering and see "+#[1]+"_ by typing:",
+"   setring Bl_"+#[1]+bsr+";","   "+#[1]+"_;");
    return();
+}
 example
+example
 { "EXAMPLE:"; echo = 2;
    ring R=0,(x,y),dp;
    poly f=y2+x3; ideal j=jacob(f);
+   poly f=y2+x3; ideal j=jacob(f);
    blowup0(j);
+   type basering; "";
+// NOTE:
+// basering is now Bl_jR (equal to R[A,B])
+// it contains the ideal j_ , such that
+//             Bl_jR/j_ is the blowup ring
+// For blowing-up another ideal in R type first:
+//             setring R;
+   type j_; "";
+   show(Bl_jR);
+   setring Bl_jR;
+   j_;"";
    ring r=32003,(x,y,z),ds;
    blowup0(maxideal(1),"m","T()");
+   type basering; "";
+// NOTE:
+// basering is now Bl_mr (equal to r[T(1..3)])
+// it contains the ideal m_ , such that
+//             Bl_mr/m_ is the blowup ring
+// For blowing-up another ideal in r type first:
+//             setring r;
+   show(Bl_mr);
+   setring Bl_mr;
    m_;
    kill Bl_jR, Bl_mr;
+   kill Bl_jR, Bl_mr;
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc elim (id, int n, int m)
+USAGE:   elim(id,n,m);  id ideal/module, n,m integers
+RETURNS: ideal/module obtained from id by eliminating variables n..m
+NOTE:    no special monomial ordering is required, result is a SB with
+         respect to ordering dp (resp. ls) if the first var not to be
+         eliminated belongs to a -p (resp. -s) blockordering
+EXAMPLE: example elim; shows an example
+USAGE:   elim(id,n,m);  id ideal/module, n,m integers
+RETURNS: ideal/module obtained from id by eliminating variables n..m
+NOTE:    no special monomial ordering is required, result is a SB with
+         respect to ordering dp (resp. ls) if the first var not to be
+         eliminated belongs to a -p (resp. -s) blockordering
+         This proc uses 'execute' or calls a procedure using 'execute'.
+EXAMPLE: example elim; shows examples
+{
 //---- get variables to be eliminated and create string for new ordering ------
    int ii; poly vars=1;
    for( ii=n; ii<=m; ii=ii+1 ) { vars=vars*var(ii); }
+   for( ii=n; ii<=m; ii=ii+1 ) { vars=vars*var(ii); }
    if( n>1 ) { poly p = 1+var(1); }
    else { poly p = 1+var(m+1); }
 …
 //-------------- create new ring and map objects to new ring ------------------
    def br = basering;
    string str = "ring newr = ("+charstr(br)+"),("+varstr(br)+ordering;
+   string str = "ring @newr = ("+charstr(br)+"),("+varstr(br)+ordering;
    execute(str);
    def i = imap(br,id);
 …
    i = eliminate(i,vars);
    setring br;
    return(imap(newr,i));
+}
 example
+   return(imap(@newr,i));
+}
+example
 { "EXAMPLE:"; echo = 2;
    ring r=0,(x,y,u,v,w),dp;
    ideal i=x-u,y-u2,w-u3,v-x+y3;
    elim(i,3,4);
+   elim(i,3,4);
    module m=i*gen(1)+i*gen(2);
    m=elim(m,3,4);show(m);
+}
 ///////////////////////////////////////////////////////////////////////////////
+   m=elim(m,3,4);show(m);
+}
+///////////////////////////////////////////////////////////////////////////////
 proc elim1 (id, poly vars)
 …
          respect to ordering dp (resp. ls) if the first var not to be
          eliminated belongs to a -p (resp. -s) blockordering
+EXAMPLE: example elim1; shows an example
+         This proc uses 'execute' or calls a procedure using 'execute'.
+EXAMPLE: example elim1; shows examples
+{
 //---- get variables to be eliminated and create string for new ordering ------
+   int ii;
+   for( ii=1; ii<=nvars(basering); ii++ )
+   {
+      if( vars/var(ii)==0 ) { poly p = 1+var(ii); }
+      break;
+   int ii;
+   for( ii=1; ii<=nvars(basering); ii=ii+1 )
+   {
+      if( vars/var(ii)==0 ) { poly p = 1+var(ii); break;}
+   }
    if( ord(p)==0 ) { string ordering = "),ls;"; }
 …
 //-------------- create new ring and map objects to new ring ------------------
    def br = basering;
    string str = "ring newr = "+charstr(br)+",("+varstr(br)+ordering;
+   string str = "ring @newr = ("+charstr(br)+"),("+varstr(br)+ordering;
    execute(str);
    def id = fetch(br,id);
 …
    id = eliminate(id,vars);
    setring br;
    return(imap(newr,id));
+}
 example
+   return(imap(@newr,id));
+}
+example
 { "EXAMPLE:"; echo = 2;
    ring r=0,(x,y,t,s,z),dp;
    ideal i=x-t,y-t2,z-t3,s-x+y3;
    elim1(i,ts);
+   elim1(i,ts);
    module m=i*gen(1)+i*gen(2);
    m=elim1(m,st); show(m);
+}
 ///////////////////////////////////////////////////////////////////////////////
+   m=elim1(m,st); show(m);
+}
+///////////////////////////////////////////////////////////////////////////////
 proc nselect (id, int n, list#)
 USAGE:   nselect(id,n[,m]); id a module or ideal, n, m integers
 RETURN:  generators of id not containing the variable n [up to m]
+EXAMPLE: example nselect; shows an example
+{
+EXAMPLE: example nselect; shows examples
+{
+   int j,k;
    if( size(#)==0 ) { #[1]=n; }
+   int j,k;
+   for( k=1; k<=ncols(id); k++ )
+   {
+      for( j=n; j<=#[1]; j++ )
+      {
+         if( size(id[k]/var(j))!=0) { id[k]=0; break; }
+   for( k=1; k<=ncols(id); k=k+1 )
+   {  for( j=n; j<=#[1]; j=j+1 )
+      {  if( size(id[k]/var(j))!=0) { id[k]=0; break; }
+      }
+   }
    return(simplify(id,2));
+}
 example
+example
 { "EXAMPLE:"; echo = 2;
    ring r=0,(x,y,t,s,z),(c,dp);
    ideal i=x-y,y-z2,z-t3,s-x+y3;
    nselect(i,3);
    module m=i*(gen(1)+gen(2));
+   module m=i*(gen(1)+gen(2));
    show(m);
    show(nselect(m,3,4));
 …
 proc sat (id, ideal j)
+USAGE:   sat(id,j);  id ideal or module, j ideal
+RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
+NOTE:    result is a std basis in the basering
+USAGE:   sat(id,j);  id=ideal/module, j=ideal
+RETURN:  list of an ideal/module [1] and an integer [2]:
+         [1] = saturation of id with respect to j (= union_(k=1...) of id:j^k)
+         [2] = saturation exponent (= min( k | id:j^k = id:j^(k+1) ))
+NOTE:    [1] is a standard basis in the basering
+DISPLAY: saturation exponent during computation if printlevel >=1
 EXAMPLE: example sat; shows an example
+{
    int ii,kk;
+   def i=id; id=std(id);
+   def i=id;
+   id=std(id);
+   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
    while( ii<=size(i) )
+   {
+      if( voice==2 )
+      {"// start quotient:",kk+1;}
+   {
+      dbprint(p-1,"// compute quotient "+string(kk+1));
       i=quotient(id,j);
       for( ii=1; ii<=size(i); ii++ )
+      for( ii=1; ii<=size(i); ii=ii+1 )
+      {
          if( reduce(i[ii],id)!=0 ) break;
 …
       id=std(i); kk=kk+1;
+   }
+   if( voice==2 )
+   {"//  saturation becomes stable after",kk-1,"iteration(s)";"";}
+   return (id);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=2,(x,y,z),dp;
+   poly F=x5+y5+(x-y)^2*xyz;
+   ideal j= jacob(F);
+   sat(j,maxideal(1));
+   dbprint(p-1,"// saturation becomes stable after "+string(kk-1)+" iteration(s)","");
+   list L = id,kk-1;
+   return (L);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   int p      = printlevel;
+   ring r     = 2,(x,y,z),dp;
+   poly F     = x5+y5+(x-y)^2*xyz;
+   ideal j    = jacob(F);
+   sat(j,maxideal(1));
+   printlevel = 2;
+   sat(j,maxideal(2));
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 USAGE:   select(id,n[,m]); id ideal/module, n, m integers
 RETURN:  generators of id containing the variable n [up to m]
 EXAMPLE: example select; shows an example
+EXAMPLE: example select; shows examples
+{
    if( size(#)==0 ) { #[1]=n; }
    int j,k;
+   for( k=1; k<=ncols(id); k++ )
+   {
+      for( j=n; j<=#[1]; j++ )
+      {
+         if( size(id[k]/var(j))==0) { id[k]=0; break; }
+   for( k=1; k<=ncols(id); k=k+1 )
+   {  for( j=n; j<=#[1]; j=j+1 )
+      {   if( size(id[k]/var(j))==0) { id[k]=0; break; }
+      }
+   }
    return(id+id);
+}
 example
+   return(simplify(id,2));
+}
+example
 { "EXAMPLE:"; echo = 2;
    ring r=0,(x,y,t,s,z),(c,dp);
    ideal i=x-y,y-z2,z-t3,s-x+y3;
    ideal j=select(i,1);
+   module m=i*(gen(1)+gen(2)); show(m);
+   show(select(m,1,2));
+}
+///////////////////////////////////////////////////////////////////////////////
+   j;
+   module m=i*(gen(1)+gen(2));
+   m;
+   select(m,1,2);
+}
+///////////////////////////////////////////////////////////////////////////////

Singular/LIB/factor.lib

r6d09c56	r6f2edc
1		// $Id: factor.lib,v 1.~~1.1.1 1997-04-25 15:13:26~~ obachman Exp $
	1	// $Id: factor.lib,v 1.2 1997-04-28 19:27:17 obachman Exp $
2	2	//(RS)
3	3	///////////////////////////////////////////////////////////////////////////////

Singular/LIB/general.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: general.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
+// $Id: general.lib,v 1.2 1997-04-28 19:27:17 obachman Exp $
 //system("random",787422842);
 //(GMG, last modified 22.06.96)
 …
  sum(vector/id/..[,v]); add components of vector/ideal/...[with indices v]
            (parameters in square brackets [] are optional)
+ wich(command);         searches for command and returns absolute
+                        path, if found
 LIB "inout.lib";
 ///////////////////////////////////////////////////////////////////////////////
 …
 proc killall
 USAGE:   killall(); (no parameter)
+         killall("proc");
+         killall("type_name");
+         killall("not", "type_name");
 COMPUTE: killall(); kills all user-defined variables but not loaded procedures
+         killall("proc"); kills all loaded procedures
+         killall("type_name"); kills all user-defined variables, of type "type_name"
+         killall("not", "type_name"); kills all user-defined
+         variables, except those of type "type_name" and except loaded procedures
 RETURN:  no return value
 NOTE:    killall should never be used inside a procedure
 …
          if( L[joni]!="LIB" and typeof(`L[joni]`)!="proc" ) { kill `L[joni]`; }
+      }
+      return();
+   }
+   if( #[1] == "proc" )
+   {
+      for ( joni=size(L); joni>0; joni-- )
+      {
+         if( L[joni]=="LIB" or typeof(`L[joni]`)=="proc" ) { kill `L[joni]`; }
+      }
+   }
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring rtest; ideal i=x,y,z; number n=37; string str="hi";
+   export rtest,i,n,str;     //this makes the local variables global
+   }
+   else
+   {
+     if( size(#)==1 )
+     {
+       if( #[1] == "proc" )
+       {
+          for ( joni=size(L); joni>0; joni-- )
+          {
+             if( L[joni]=="LIB" or typeof(`L[joni]`)=="proc" )
+               { kill `L[joni]`; }
+          }
+       }
+       else
+       {
+          for ( ; joni>2; joni-- )
+          {
+            if(typeof(`L[joni]`)==#[1] and L[joni]!="LIB" and typeof(`L[joni]`)!="proc") { kill `L[joni]`; }
+          }
+        }
+     }
+     else
+     {
+        for ( ; joni>2; joni-- )
+        {
+          if(typeof(`L[joni]`)!=#[2] and L[joni]!="LIB" and typeof(`L[joni]`)!="proc") { kill `L[joni]`; }
+        }
+     }
+  }
+  return();
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring rtest; ideal i=x,y,z; number n=37; string str="hi"; int j = 3;
+   export rtest,i,n,str,j;     //this makes the local variables global
    listvar(all);
+   killall();
+   killall("string"); // kills all string variables
+   listvar(all);
+   killall("not", "int"); // kills all variables except int's (and procs)
+   listvar(all);
+   killall(); // kills all vars except loaded procs
    listvar(all);
+}
 …
+}
 ///////////////////////////////////////////////////////////////////////////////
+proc which (command)
+USAGE:    which(command); command = string expression
+RETURN:   Absolute pathname of command, if found in search path.
+          Empty string, otherwise.
+NOTE:     Based on the Unix command 'which'.
+EXAMPLE:  example which; shows an example
+{
+   int rs;
+   int i;
+   string fn = "/tmp/which_" + string(system("pid"));
+   string pn;
+   if( typeof(command) != "string")
+   {
+        return pn;
+   }
+   i = system("sh", "which " + command + " > " + fn);
+   pn = read(fn);
+   pn[size(pn)] = "";
+   i = 1;
+   while ((pn[i] != " ") and (pn[i] != ""))
+   {
+     i = i+1;
+   }
+   if (pn[i] == " ") {pn[i] = "";}
+   rs = system("sh", "ls " + pn + " > " + fn + " 2>&1 ");
+   i = system("sh", "rm " + fn);
+   if (rs == 0) {return (pn);}
+   else
+   {
+     print (command + " not found ");
+     return ("");
+   }
+}
+example
+{  "EXAMPLE:"; echo = 2;
+    which("Singular");
+}
+///////////////////////////////////////////////////////////////////////////////

Singular/LIB/homog.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: homog.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
+// $Id: homog.lib,v 1.2 1997-04-28 19:27:18 obachman Exp $
 //(BM 11/95)
 ///////////////////////////////////////////////////////////////////////////////
 LIBRARY:  homog.lib           homological algebra porcedures
+LIBRARY:  homog.lib           homological algebra porcedures
 ext(<int>,<ideal>);             Ext^k(R/I,R),    I - ideal in basering R
 Ext(<int>,<mod>,<mod>);         Ext^k(M,N),      M,N modules
 …
 example
+{
   "EXAMPLE";   echo=2;  int printlevel=2;
+  "EXAMPLE";   echo=2;  printlevel=2;
   ring  r = 0,(x,y),(dp,C);
   ideal i = x2-y3;

Singular/LIB/homolog.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: homolog.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
 //(BM 11/95)
+// $Id: homolog.lib,v 1.2 1997-04-28 19:27:19 obachman Exp $
+//(BM/GMG, last modified 22.06.96)
 ///////////////////////////////////////////////////////////////////////////////
+LIBRARY:  homolog.lib   PROCEDURES FOR HOMOLOGICAL ALGEBRA
+Ext_R(k,id);            Ext^k(R/id,R), id  = ideal
+Ext(k,M,N);             Ext^k(M,N),    M,N = modules
+Hom(M,N);               Hom(M,N),      M,N = modules
+kohom(A,k);             Hom(R^k,A), A = matrix
+kontrahom(A,k);         Hom(A,R^k), A = matrix
+LIBRARY:  homolog.lib   PROCEDURES FOR HOMOLOGICAL ALGEBRA
+ cup(M);                cup: Ext^1(M',M') x Ext^1() --> Ext^2()
+ cupproduct(M,N,P,p,q); cup: Ext^p(M',N') x Ext^q(N',P') --> Ext^p+q(M',P')
+ Ext_R(k,M);            Ext^k(M',R),    M module, R basering, M'=coker(M)
+ Ext(k,M,N);            Ext^k(M',N'),   M,N modules, M'=coker(M), N'=coker(N)
+ Hom(M,N);              Hom(M',N'),     M,N modules, M'=coker(M), N'=coker(N)
+ homology(A,B,M,N)      ker(B)/im(A),   homology of complex R^k--A->M'--B->N'
+ kernel(A,M,N);         ker(M'--A->N')  M,N modules, A matrix
+ kohom(A,k);            Hom(R^k,A),     A matrix over basering R
+ kontrahom(A,k);        Hom(A,R^k),     A matrix over basering R
+LIB "general.lib";
+LIB "deform.lib";
 LIB "matrix.lib";
+LIB "poly.lib";
 ///////////////////////////////////////////////////////////////////////////////
+proc Ext_R (int k, ideal I, list #)
+USAGE:   Ext_R(k,id[,any]);  k=integer, id=ideal
+DISPLAY: degree of Ext^k(P/id,P)
+RETURN:  Ext^k(P/id,P), presentation as a quotient of a free module.
+         If a third argument is given (of any type) return a list of two
+         modules:
+           [1] = module Ext^k(P/id,P)
+           [2] = SB of  Ext^k(P/id,P)
+EXAMPLE: example Ext_R; shows an example
+proc cup (module M,list #)
+USAGE:   cup(M,[,any,any]);  M=module
+COMPUTE: cup-product Ext^1(M',M') x Ext^1(M',M') ---> Ext^2(M',M'),
+         where M':=R^m/M, if M in R^m, R basering (i.e. M':=coker(matrix(M)))
+         in case of a second argument: symmetrized cup-product
+ASSUME:  all Ext's  are finite dimensional
+RETURN:  matrix of the associated linear map,
+         i.e. the columns of <matrix> present the coordinates of b_i & b_j
+         (resp. (1/2)(b_i & b_j + b_j & b_i) in the symmetric version)
+         with respect to a kbase of Ext^2,
+         (where b_1,b_2,... is a kbase of Ext^1 and & denotes cup product)
+         in case of  a third argument return a list:
+               L[1] = matrix see above (and symmetric case)
+               L[2] = matrix of kbase of Ext^1
+               L[3] = matrix of kbase of Ext^2
+NOTE:    printlevel >=1;  shows what is going on
+         printlevel >=2;  shows result in another representation
+         For computing cupproduct of M itself, apply proc to syz(M)!
+EXAMPLE: example cup; shows examples
+{
+   module m1,m2,ret,ret0;
+//------------------ compute resolution of P/I --------------------------------
+   mres(I,k+1,resI);
+//-----------------  apply Hom(_,P) at k-th place -----------------------------
+   m2 = transpose(resI(k+1));
+   m1 = transpose(resI(k));
+//-----------------  ker(m2)/im(m1) -------------------------------------------
+   m2 = simplify(m2,10);
+   if (m2[1]==0) { ret = m1; }
+   else          { ret = modulo(syz(m2),m1); }
+   ret0 = std(ret);
+   "// degree of Ext^"+string(k)+"(P/I,P):";degree(ret0);
+   if (size(#)>0) { return(ret,ret0); }
+   return(ret);
+//---------- initialization ---------------------------------------------------
+   int    i,j,k,f0,f1,f2,f3,e1,e2;
+   module M1,M2,A,B,C,ker,ima,ext1,ext2,ext10,ext20;
+   matrix cup[1][0];
+   matrix kb1,lift1,kb2,mA,mB,mC;
+   ideal  tes1,tes2,null;
+   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
+//-----------------------------------------------------------------------------
+//take a resolution of M<--F(0)<--- ...  <---F(3)
+//apply Hom(-,M) and compute the Ext's
+//-----------------------------------------------------------------------------
+   list resM = res(M,3);
+   M1 = resM[2];
+   M2 = resM[3];
+   f0 = nrows(M);
+   f1 = ncols(M);
+   f2 = ncols(M1);
+   f3 = ncols(M2);
+   tes1 = simplify(ideal(M),10);
+   tes2=simplify(ideal(M1),10);
+   if ((tes1[1]*tes2[1]==0) or (tes1[1]==1) or (tes2[1]==1))
+   {
+      dbprint(p,"// Ext == 0 , hence 'cup' is the zero-map");
+      return(@cup);
+   }
+//------ compute Ext^1 --------------------------------------------------------
+   B     = kohom(M,f2);
+   A     = kontrahom(M1,f0);
+   C     = intersect(A,B);
+   C     = reduce(C,std(null));C = simplify(C,10);
+   ker   = lift(A,C)+syz(A);
+   ima   = kohom(M,f1);
+   ima   = ima + kontrahom(M,f0);
+   ext1  = modulo(ker,ima);
+   ext10 = std(ext1);
+   e1    = vdim(ext10);
+   dbprint(p-1,"// vdim (Ext^1) = "+string(e1));
+   if (e1 < 0)
+   {
+     "// Ext^1 not of finite dimension";
+     return(cup);
+   }
+   kb1 = kbase(ext10);
+   kb1 = matrix(ker)*kb1;
+   dbprint(p-1,"// kbase of Ext^1(M,M)",
+               "//  - the columns present the kbase elements in Hom(F(1),F(0))",
+               "//  - F(*) a free resolution of M",kb1);
+//------ compute the liftings of Ext^1 ----------------------------------------
+   mC = matrix(A)*kb1;
+   lift1 =lift(B,mC);
+   dbprint(p-1,"// lift kbase of Ext^1:",
+    "//  - the columns present liftings of kbase elements into Hom(F(2),F(1))",
+    "//  - F(*) a free resolution of M ",lift1);
+//------ compute Ext^2  -------------------------------------------------------
+   B    = kohom(M,f3);
+   A    = kontrahom(M2,f0);
+   C    = intersect(A,B);
+   C    = reduce(C,std(null));C = simplify(C,10);
+   ker  = lift(A,C)+syz(A);
+   ima  = kohom(M,f2);
+   ima  = ima + kontrahom(M1,f0);
+   ext2 = modulo(ker,ima);
+   ext20= std(ext2);
+   e2   = vdim(ext20);
+   if (e2<0)
+   {
+     "// Ext^2 not of finite dimension";
+     return(cup);
+   }
+   dbprint(p-1,"// vdim (Ext^2) = "+string(e2));
+   kb2 = kbase(ext20);
+   kb2 = matrix(ker)*kb2;
+   dbprint(p-1,"// kbase of Ext^2(M,M)",
+               "//  - the columns present the kbase elements in Hom(F(2),F(0))",
+               "//  - F(*) is a a free resolution of M ",kb2);
+//-------  compute: cup-products of base-elements -----------------------------
+   for (i=1;i<=e1;i=i+1)
+   {
+     for (j=1;j<=e1;j=j+1)
+     {
+       mA = matrix(ideal(lift1[j]),f1,f2);
+       mB = matrix(ideal(kb1[i]),f0,f1);
+       mC = mB*mA;
+       if (size(#)==0)
+       {                                          //non symmestric
+         mC = matrix(ideal(mC),f0*f2,1);
+         cup= concat(cup,mC);
+       }
+       else                                       //symmetric version
+       {
+         if (j>=i)
+         {
+           if (j>i)
+           {
+             mA = matrix(ideal(lift1[i]),f1,f2);
+             mB = matrix(ideal(kb1[j]),f0,f1);
+             mC = mC+mB*mA;mC=(1/2)*mC;
+           }
+           mC = matrix(ideal(mC),f0*f2,1);
+           cup= concat(cup,mC);
+         }
+       }
+     }
+   }
+   dbprint(p-1,"// matrix of cup-products (in Ext^2)",cup,"////// end level 2 //////");
+//------- comptute: presentation of base-elements -----------------------------
+   cup = lift(ker,cup);
+   cup = lift_kbase(cup,ext20);
+   if( p>2 )
+   {
+     "// the associated matrices of the bilinear mapping 'cup' ";
+     "// corresponding to the kbase elements of Ext^2(M,M) are shown,";
+     "//  i.e. the rows of the final matrix are written as matrix of";
+     "//  a bilinear form on Ext^1 x Ext^1";
+     matrix BL[e1][e1];
+     for (k=1;k<=e2;k=k+1)
+     {
+       "//_____"+string(k)+". component:";
+       for (i=1;i<=e1;i=i+1)
+       {
+         for (j=1;j<=e1;j=j+1)
+         {
+           if (size(#)==0) { BL[i,j]=cup[k,j+e1*(i-1)]; }
+           else
+           {
+             if (i<=j)
+             {
+               BL[i,j]=cup[k,j+e1*(i-1)-binomial(i,2)];
+               BL[j,i]=BL[i,j];
+             }
+           }
+         }
+       }
+       print(BL);
+     }
+     "////// end level 3 //////";
+   }
+   if (size(#)>2) { return(cup,kb1,kb2);}
+   else {return(cup);}
+}
 example
+{ "EXAMPLE:";     echo=2;
+   ring r=0,(x,y,z),dp;
+   ideal  i = x2y,y2z,z3x;
+   module E = Ext_R(1,i);"";
+   E = Ext_R(2,i);
+   show(E);"";
+   list LE = Ext_R(3,i,"");
+   kbase(LE[2]);"";
+   E = Ext_R(4,i);
+{"EXAMPLE";  echo=2;
+  int p      = printlevel;
+  ring  rr   = 32003,(x,y,z),(dp,C);
+  ideal  I   = x4+y3+z2;
+  qring  o   = std(I);
+  module M   = [x,y,0,z],[y2,-x3,z,0],[z,0,-y,-x3],[0,z,x,-y2];
+  print(cup(M));
+  print(cup(M,1));
+  // 2nd EXAMPLE  (shows what is going on)
+  printlevel = 3;
+  ring   r   = 0,(x,y),(dp,C);
+  ideal  i   = x2-y3;
+  qring  q   = std(i);
+  module M   = [-x,y],[-y2,x];
+  print(cup(M));
+  printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
+proc Ext (int k, module M, module N, list #)
+USAGE:   Ext(k,M,N[,any]);  k=integer, M,N=modules
+DISPLAY: degree of Ext^k(M,N)
+RETURN:  Ext^k(M,N), presentation as a quotient of a free module.
+         If a third argument is given (of any type) return a list of two
+         modules:
+           [1] = module Ext^k(M,N)
+           [2] = SB of  Ext^k(M,N)
+EXAMPLE: example Ext; shows an example
+proc cupproduct (module M,N,P,int p,q,list #)
+USAGE:   cupproduct(M,N,P,p,q[,any]);  M,N,P modules, p,q integers
+COMPUTE: cup-product Ext^p(M',N') x Ext^q(N',P') ---> Ext^p+q(M',P')
+         where M':=R^m/M, if M in R^m, R basering (i.e. M':=coker(matrix(M)))
+ASSUME:  all Ext's  are of finite dimension
+RETURN:  matrix of the associated linear map Ext^p(tensor)Ext^q-->Ext^p+q
+         i.e. the columnes of <matrix> present the coordinates of
+         the cup products (b_i & c_j) with respect to a kbase of Ext^p+q
+         (b_i resp. c_j are choosen bases of Ext^p resp. Ext^q)
+         in case of a 6th argument:
+            return a list
+            L[1] = matrix (see above)
+            L[2] = matrix of kbase of Ext^p(M',N')
+            L[3] = matrix of kbase of Ext^q(N',P')
+            L[4] = matrix of kbase of Ext^p+q(N',P')
+NOTE:    printlevel >=1;  shows what is going on
+         printlevel >=2;  shows result in another representation
+         For computing cupproduct of M,N itself, apply proc to syz(M),syz(N)!
+EXAMPLE: example cupproduct; shows examples
+{
+   module A,B,C,M1,M2,ker,imag,extMN,extMN0;
+   ideal  test1,test2;
+//----------  resolution of M  and N ------------------------------------------
+  if (k>0)
+  {
+     mres(M,k+1,resM);
+     M1 = resM(k);
+     M2 = resM(k+1);
+     test1 = simplify(ideal(M1),10);
+     test2 = simplify(ideal(N),10);
+     if ((test1[1]==0) or (test2[1]==0))
+        { "//Ext(M,N)=0";return(extMN); }
+     else
+     {
+       test1 = simplify(ideal(M2),10);
+       if (test1[1]==0)                             //Ker(Hom(m2,N))
+          { ker = freemodule(ncols(M1)*nrows(N));}
+       else
+       {
+          A   = kontrahom(M2,nrows(N));
+          B   = kohom(N,ncols(M2));
+          C   = intersect(A,B);
+          C   = reduce(C,std(ideal(0)));C=simplify(C,10);
+          ker = lift(A,C)+syz(A);
+       }
+       imag  = kohom(N,ncols(M1));
+       A     = kontrahom(M1,nrows(N));
+       imag  = imag+A;                              //im(Hom(m1,M))
+       extMN = modulo(ker,imag);
+       extMN0= std(extMN);
+       "// degree of Ext^"+string(k)+"(M,N):";degree(extMN0);
+//---------- initialization ---------------------------------------------------
+   int    e1,e2,e3,i,j,k,f0,f1,f2;
+   module M1,M2,N1,N2,P1,P2,A,B,C,ker,ima,extMN,extMN0,extMP,
+          extMP0,extNP,extNP0;
+   matrix cup[1][0];
+   matrix kbMN,kbMP,kbNP,lift1,mA,mB,mC;
+   ideal  test1,test2,null;
+   int pp = printlevel-voice+3;  // pp=printlevel+1 (default: p=1)
+//-----------------------------------------------------------------------------
+//compute resolutions of M and N
+//                     M<--F(0)<--- ...  <---F(p+q+1)
+//                     N<--G(0)<--- ...  <---G(q+1)
+//-----------------------------------------------------------------------------
+   list resM = res(M,p+q+1);
+   M1 = resM[p];
+   M2 = resM[p+1];
+   list resN = res(N,q+1);
+   N1 = resN[q];
+   N2 = resN[q+1];
+   P1 = resM[p+q];
+   P2 = resM[p+q+1];
+//-------test: Ext==0?---------------------------------------------------------
+   test1 = simplify(ideal(M1),10);
+   test2 = simplify(ideal(N),10);
+   if (test1[1]==0) { dbprint(pp,"//Ext(M,N)=0");return(cup); }
+   test1 = simplify(ideal(N1),10);
+   test2 = simplify(ideal(P),10);
+   if (test1[1]==0) { dbprint(pp,"//Ext(N,P)=0");return(cup); }
+   test1 = simplify(ideal(P1),10);
+   if (test1[1]==0) { dbprint(pp,"//Ext(M,P)=0");return(cup); }
+ //------ compute kbases of Ext's ---------------------------------------------
+ //------ Ext(M,N)
+   test1 = simplify(ideal(M2),10);
+   if (test1[1]==0) { ker = freemodule(ncols(M1)*nrows(N));}
+   else
+   {
+     A   = kontrahom(M2,nrows(N));
+     B   = kohom(N,ncols(M2));
+     C   = intersect(A,B);
+     C   = reduce(C,std(ideal(0)));C=simplify(C,10);
+     ker = lift(A,C)+syz(A);
+   }
+   ima   = kohom(N,ncols(M1));
+   A     = kontrahom(M1,nrows(N));
+   ima   = ima+A;
+   extMN = modulo(ker,ima);
+   extMN0= std(extMN);
+   e1    = vdim(extMN0);
+   dbprint(pp-1,"// vdim Ext(M,N) = "+string(e1));
+   if (e1 < 0)
+   {
+     "// Ext(M,N) not of finite dimension";
+     return(cup);
+   }
+   kbMN  = kbase(extMN0);
+   kbMN = matrix(ker)*kbMN;
+   dbprint(pp-1,"// kbase of Ext^p(M,N)",
+          "//  - the columns present the kbase elements in Hom(F(p),G(0))",
+          "//  - F(*),G(*) are free resolutions of M and N",kbMN);
+//------- Ext(N,P)
+   test1 = simplify(ideal(N2),10);
+   if (test1[1]==0) {  ker = freemodule(ncols(N1)*nrows(P)); }
+   else
+   {
+     A = kontrahom(N2,nrows(P));
+     B = kohom(P,ncols(N2));
+     C = intersect(A,B);
+     C = reduce(C,std(ideal(0)));C=simplify(C,10);
+     ker = lift(A,C)+syz(A);
+   }
+   ima   = kohom(P,ncols(N1));
+   A     = kontrahom(N1,nrows(P));
+   ima   = ima+A;
+   extNP = modulo(ker,ima);
+   extNP0= std(extNP);
+   e2    = vdim(extNP0);
+   dbprint(pp-1,"// vdim Ext(N,P) = "+string(e2));
+   if (e2 < 0)
+   {
+     "// Ext(N,P) not of finite dimension";
+     return(cup);
+   }
+   kbNP  = kbase(extNP0);
+   kbNP  = matrix(ker)*kbNP;
+   dbprint(pp-1,"// kbase of Ext(N,P):",kbNP,
+          "// kbase of Ext^q(N,P)",
+          "//  - the columns present the kbase elements in Hom(G(q),H(0))",
+          "//  - G(*),H(*) are free resolutions of N and P",kbNP);
+//------ Ext(M,P)
+   test1 = simplify(ideal(P2),10);
+   if (test1[1]==0) { ker = freemodule(ncols(P1)*nrows(P)); }
+   else
+   {
+     A = kontrahom(P2,nrows(P));
+     B = kohom(P,ncols(P2));
+     C = intersect(A,B);
+     C = reduce(C,std(ideal(0)));C=simplify(C,10);
+     ker = lift(A,C)+syz(A);
+   }
+   ima   = kohom(P,ncols(P1));
+   A     = kontrahom(P1,nrows(P));
+   ima   = ima+A;
+   extMP = modulo(ker,ima);
+   extMP0= std(extMP);
+   e3    = vdim(extMP0);
+   dbprint(pp-1,"// vdim Ext(M,P) = "+string(e3));
+   if (e3 < 0)
+   {
+     "// Ext(M,P) not of finite dimension";
+     return(cup);
+   }
+   kbMP  = kbase(extMP0);
+   kbMP  = matrix(ker)*kbMP;
+   dbprint(pp-1,"// kbase of Ext^p+q(M,P)",
+          "//  - the columns present the kbase elements in Hom(F(p+q),H(0))",
+          "//  - F(*),H(*) are free resolutions of M and P",kbMP);
+//----- lift kbase of Ext(M,N) ------------------------------------------------
+   lift1 = kbMN;
+   for (i=1;i<=q;i=i+1)
+   {
+     mA = kontrahom(resM[p+i],nrows(resN[i]));
+     mB = kohom(resN[i],ncols(resM[p+i]));
+     lift1 = lift(mB,mA*lift1);
+   }
+   dbprint(pp-1,"// lifting of kbase of Ext^p(M,N)",
+   "//  - the columns present the lifting of kbase elements in Hom(F(p+q),G(q))",lift1);
+//-------  compute: cup-products of base-elements -----------------------------
+   f0 = nrows(P);
+   f1 = ncols(N1);
+   f2 = ncols(resM[p+q]);
+   for (i=1;i<=e1;i=i+1)
+   {
+     for (j=1;j<=e2;j=j+1)
+     {
+       mA = matrix(ideal(lift1[j]),f1,f2);
+       mB = matrix(ideal(kbMP[i]),f0,f1);
+       mC = mB*mA;
+       mC = matrix(ideal(mC),f0*f2,1);
+       cup= concat(cup,mC);
+     }
+  }
+  else { extMN,extMN0 = Hom(M,N,1); }
+  if (size(#)>0) { return(extMN,extMN0); }
+  return(extMN );
+   }
+   dbprint(pp-1,"// matrix of cup-products (in Ext^p+q)",cup,"////// end level 2 //////");
+//------- comptute: presentation of base-elements -----------------------------
+   cup = lift(ker,cup);
+   cup = lift_kbase(cup,extMP0);
+//------- special output ------------------------------------------------------
+   if (pp>2)
+   {
+     "// the associated matrices of the bilinear mapping 'cup' ";
+     "// corresponding to the kbase elements of Ext^p+q(M,P) are shown,";
+     "//  i.e. the rows of the final matrix are written as matrix of";
+     "//  a bilinear form on Ext^p x Ext^q";
+     matrix BL[e1][e2];
+     for (k=1;k<=e3;k=k+1)
+     {
+       "//----"+string(k)+". component:";
+       for (i=1;i<=e1;i=i+1)
+       {
+         for (j=1;j<=e2;j=j+1)
+         {
+           BL[i,j]=cup[k,j+e1*(i-1)];
+         }
+        }
+        print(BL);
+      }
+     "////// end level 3 //////";
+   }
+   if (size(#)) { return(cup,kbMN,kbNP,kbMP);}
+   else         { return(cup); }
+}
 example
+{ "EXAMPLE:"; echo=2;
+   ring r=0,(x,y),(c,dp);
+   ideal i=x2-y3;
+   qring qr=std(i);
+   module M=[-x,y],[-y2,x];
+   list EXT1=Ext(1,M,M,"");
+   kbase(EXT1[2]);
+}
+////////////////////////////////////////////////////////////////////////////////
+proc Hom (module M, module N, list #)
+USAGE:   Hom(M,N);  M,N=modules
+DISPLAY: degree of Hom(M,N);
+RETURN:  Hom(M,N), presentation as a quotient of a free module.
+         If a third argument is given (of any type) return a list of two
+         modules:
+           [1] = module Hom(M,N)
+           [2] = SB of  Hom(M,N)
+EXAMPLE: example Hom; shows an example
+{"EXAMPLE";  echo=2;
+  int p      = printlevel;
+  ring  rr   = 32003,(x,y,z),(dp,C);
+  ideal  I   = x4+y3+z2;
+  qring  o   = std(I);
+  module M   = [x,y,0,z],[y2,-x3,z,0],[z,0,-y,-x3],[0,z,x,-y2];
+  print(cupproduct(M,M,M,1,3));
+  printlevel = 3;
+  list l     = (cupproduct(M,M,M,1,3,"any"));
+  show(l[1]);show(l[2]);
+  printlevel = p;
+}
+///////////////////////////////////////////////////////////////////////////////
+proc Ext_R (intvec v, module M, list #)
+USAGE:   Ext_R(v,M[,p]);  v=int/intvec , M=module, p=int
+COMPUTE: A presentation of Ext^k(M',R); for k=v[1],v[2],..., M'=coker(M).
+         Let ...--> F2 --> F1 --M-> F0-->M'-->0 be a free resolution of M'. If
+<-- F0* <-A1-- F1* <-A2-- F2* <-A3--... denotes the dual sequence,
+         Fi*=Hom(Fi,R), then Ext^k = ker(Ak)/im(Ak+1) is presented as in the
+         following exact sequences:
+                 Fk-1* <-Ak-- Fk* <-syz(Ak)-- R^p
+                 Fk*/im(Ak+1) <-syz(Ak)-- R^p <-Ext^k-- R^q
+         Hence Ext^k=modulo(syz(Ak),Ak+1) presents Ext^k(M',R);
+RETURN:  Ext^k, of type module, a presentation of Ext^k(M',R) if v is of type
+         int, resp. a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.
+         In case of a third argument of type int return a list:
+           [1] = module Ext^k/list of Ext^k
+           [2] = SB of Ext^k/list of SB of Ext^k
+           [3] = matrix/list of matrices, each representing kbase of Ext^k
+                 (if finite dimensional)
+DISPLAY: printlevel >=0: degree of Ext^k for each k  (default)
+         printlevel >=1: Ak, Ak+1 and kbase of Ext^k in Fk*
+NOTE:    In order to compute Ext^k(M,R) use the command Ext_R(k,syz(M));
+         or the 2 commands: list L=mres(M,2);  Ext_R(k,L[2]);
+EXAMPLE: example Ext_R; shows examples
+{
+//Bemerkung (nur fuer internen Gebrauch): Ext_R(v,M[,"sres",p]); bewirkt folgendes:
+// Es wird der input mit sres statt mres aufgeloest. Bei den bisherigen Versuchen
+// ist dies insgesamt langsamer, da sres zu gro§e Objekte zurueckliefert.
+// Der Versuch, auch spaeter sres statt mres zu verwenden, ist ganz aufgegeben,
+// da mit sres spaeter syz und modulo viel laenger dauern!
+// In case M is known to be a SB, set attrib(M,"isSB",1); in order to
+// avoid  unnecessary SB computations
+//------------ initialisation -------------------------------------------------
+  module m1,m2,ret,ret0;
+  matrix ker,kb;
+  list L1,L2,L3,L,resl,K;
+  int k,max,ii,t1,t2;
+  int s = size(v);
+  intvec v1 = sort(v)[1];
+  max = v1[s];                 // the maximum integer occuring in intvec v
+  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
+  // --------------- Variante mit sres
+  for( ii=1; ii<=size(#); ii++ )
+  {
+    t2=1; // return a list if t2=1
+    if( typeof(#[ii])=="string" )
+    {
+      if ( #[ii]=="sres" ) { t1=1; t2=0; } // use sres instead of mres if t1=1
+    }
+  }
+//----------------- compute resolution of coker(M) ----------------------------
+  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
+  if( t1==1 )
+  {
+    if( attrib(M,"isSB")==0 ) { M=std(M); }
+    resl = sres(M,max+1);
+  }
+  else { resl = mres(M,max+1); }
+  for( ii=1; ii<=s; ii++ )
+  {
+//-----------------  apply Hom(_,R) at k-th place -----------------------------
+    k=v[ii];
+    if( k<0 )                   // Ext^k=0 for negative k
+    {
+      dbprint(p-1,"// Ext^i=0 for i<0!");
+      ret    = gen(1);
+      ret0   = std(ret);
+      L1[ii] = ret;
+      L2[ii] = ret0;
+      L3[ii] = matrix(kbase(ret0));
+      dbprint(p,"// degree of Ext^"+string(k)+":");
+      if( p>=0 ) { degree(ret0);"";}
+    }
+    else
+    {
+      m2 = transpose(resl[k+1]);
+      if( k==0 ) { m1=0*gen(nrows(m2)); }
+      else { m1 = transpose(resl[k]); }
+//----------------- presentation of ker(m2)/im(m1) ----------------------------
+      ker = syz(m2);
+      ret = modulo(ker,m1);
+      dbprint(p-1,"// Computing Ext^"+string(k)+":",
+         "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of M,",
+         "// then F"+string(k)+"*-->F"+string(k+1)+"* is given by:",m2,
+         "// and F"+string(k-1)+"*-->F"+string(k)+"* is given by:",m1,"");
+      ret0 = std(ret);
+      dbprint(p,"// degree of Ext^"+string(k)+":");
+      if( p>0 ) { degree(ret0);"";}
+      if( t2 )
+      {
+         if( vdim(ret0)>=0 )
+         {
+            kb = kbase(ret0);
+            if ( size(ker)!=0 ) { kb = matrix(ker)*kb; }
+            dbprint(p-1,
+            "// columns of matrix are kbase of Ext^"+string(k)+" in F"+string(k)+"*:",kb,"");
+            L3[ii] = kb;
+         }
+         L2[ii] = ret0;
+      }
+      L1[ii] = ret;
+    }
+  }
+  if( t2 )
+  {
+     if( s>1 ) { L = L1,L2,L3; return(L); }
+     else { L = ret,ret0,kb; return(L); }
+  }
+  else
+  {
+     if( s>1 ) { return(L1); }
+     else { return(ret); }
+  }
+}
+example
+{"EXAMPLE:";     echo=2;
+  int p      = printlevel;
+  printlevel = 1;
+  ring r     = 0,(x,y,z),dp;
+  ideal i    = x2y,y2z,z3x;
+  module E   = Ext_R(1,i);       //computes Ext^1(r/i,r)
+  is_zero(E);
+  module m   = [x],[0,y];
+  list L     = Ext_R(2..3,m);    //computes Ext^i(r^2/m,r), i=2,3
+  show(L);"";
+  qring R    = std(x2+yz);
+  intvec v   = 0,2,4;
+  printlevel = 2;                //shows what is going on
+  ideal i    = x,y,z;            //computes Ext^i(r/(x,y,z),r/(x2+yz)), i=0,2,4
+  list L     = Ext_R(v,i,1);     //over the qring R=r/(x2+yz), std and kbase
+  printlevel = p;
+}
+///////////////////////////////////////////////////////////////////////////////
+proc Ext (intvec v, module M, module N, list #)
+USAGE:   Ext(v,M,N[,any]);  v=int/intvec, M,N=modules
+COMPUTE: A presentation of Ext^k(M',N'); for k=v[1],v[2],... where
+         M'=coker(M) and N'=coker(N). Let
+<--M'<-- F0 <-M-- F1 <-- F2 <--...  resp. 0<--N'<-- G0 <--N- G1 be
+         a free resolution of M' resp. a presentations of N'. Consider
+                  0                  0
+                      |^                 |^                 |^
+               --> Hom(Fk-1,N') -Ak-> Hom(Fk,N') -Ak+1-> Hom(Fk+1,N')
+                      |^                 |^                 |^
+               --> Hom(Fk-1,G0) -Ak-> Hom(Fk,G0) -Ak+1-> Hom(Fk+1,G0)
+                                         |^                 |^
+                                         |C                 |B
+                                      Hom(Fk,G1) -----> Hom(Fk+1,G1)
+         (Ak,Ak+1 induced by M and B,C induced by N).
+         Let K=modulo(Ak+1,B), J=module(Ak)+module(C) and Ext=modulo(K,J),
+         then we have exact sequences
+                  R^p  --K-> Hom(Fk,G0) --Ak+1-> Hom(Fk+1,G0)/im(B)
+         R^q -Ext-> R^p --K->Hom(Fk,G0)/im(Ak)+im(C) --Ak+1->Hom(Fk+1,G0)/im(B)
+         Hence Ext presents Ext^k(M',N')
+RETURN:  Ext, of type module, a presentation of Ext^k(M',N') if v is of type
+         int, resp. a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.
+         In case of a third argument of any type return a list:
+           [1] = module Ext/list of Ext^k
+           [2] = SB of Ext/list of SB of Ext^k
+           [3] = matrix/list of matrices, each representing a kbase of Ext^k
+                 (if finite dimensional)
+DISPLAY: printlevel >=0: degree of Ext^k for each k  (default)
+         printlevel >=1: Ak, Ak+1 and kbase of Ext^k in Hom(Fk,G0)
+                            (if finite dimensional)
+NOTE:    In order to compute Ext^k(M,N) use the command Ext(k,syz(M),syz(N));
+         or: list P=mres(M,2); list Q=mres(N,2); Ext(k,P[2],Q[2]);
+EXAMPLE: example Ext;   shows examples
+{
 //---------- initialisation ---------------------------------------------------
+   module A,B,C,ker,imag,homMN,homMN0;
+   ideal  tes;
+   tes = simplify(ideal(M),10);
+   if (tes[1]==0)
+   {
+      "// Hom(ring,N)=N";
+      if (size(#)>0) { return(N,std(N)); }
+      return(N);
+   }
+   tes = simplify(ideal(N),10);
+   if (tes[1]==0)
+   {
+      "// Hom(M,ring)=M^*";
+      homMN = transpose(M);
+      if (size(#)>0) { return(homMN,std(homMN));}
+      return(homMN);
+   }
+   B = kohom(N,ncols(M));
+   A = kontrahom(M,nrows(N));
+   C = intersect(A,B);
+   C = reduce(C,std(ideal(0)));C=simplify(C,10);
+   ker   = lift(A,C)+syz(A);                              //ker(Hom(m,N))
+   imag  = kohom(N,ncols(M));                             //im(Hom(M,n))
+   homMN = modulo(ker,imag);
+   homMN0= std(homMN);
+   "// degree of Hom(M,N):";degree(homMN0);
+   if (size(#)>0) { return(homMN,homMN0); }
+   return(homMN);
+  int k,max,ii,l,row,col;
+  module A,B,C,D,M1,M2,N1,ker,imag,extMN,extMN0;
+  matrix kb;
+  list L1,L2,L3,L,resM,K;
+  ideal  test1;
+  intmat Be;
+  int s = size(v);
+  intvec v1 = sort(v)[1];
+  max = v1[s];                      // the maximum integer occuring in intvec v
+  int p = printlevel-voice+3;       // p=printlevel+1 (default: p=1)
+//---------- test: coker(N)=basering, coker(N)=0 ? ----------------------------
+  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
+  N1 = std(N);
+  if( size(N1)==0 )      //coker(N)=basering, in this case proc Ext_R is faster
+  { printlevel=printlevel+1;
+    if( size(#)==0 )
+    { def E = Ext_R(v,M);
+      printlevel=printlevel-1;
+      return(E);
+    }
+    else
+    { def E = Ext_R(v,M,#[1]);
+       printlevel=printlevel-1;
+       return(E);
+     }
+  }
+  if( dim(N1)==-1 )                          //coker(N)=0, all Ext-groups are 0
+  { dbprint(p-1,"2nd module presents 0, hence Ext^k=0, for all k");
+    for( ii=1; ii<=s; ii++ )
+    { k=v[ii];
+      extMN    = gen(1);
+      extMN0   = std(extMN);
+      L1[ii] = extMN;
+      L2[ii] = extMN0;
+      L3[ii] = matrix(kbase(extMN0));
+      if( p>0 ) { "// degree of Ext^"+string(k)+":"; degree(extMN0);""; }
+    }
+  }
+  else
+  {
+    if( size(N1) < size(N) ) { N=N1;}
+    row = nrows(N);
+//---------- resolution of M -------------------------------------------------
+    resM = mres(M,max+1);
+    for( ii=1; ii<=s; ii++ )
+    { k=v[ii];
+      if( k<0  )                                    // Ext^k is 0 for negative k
+      { dbprint(p-1,"// Ext^k=0 for k<0!");
+        extMN    = gen(1);
+        extMN0   = std(extMN);
+        L1[ii] = extMN;
+        L2[ii] = extMN0;
+        L3[ii] = matrix(kbase(extMN0));
+        if( p>0 ) { "// degree of Ext^"+string(k)+":"; degree(extMN0);""; }
+      }
+      else
+      { M2 = resM[k+1];
+        if( k==0 ) { M1=0*gen(nrows(M2)); }
+        else { M1 = resM[k]; }
+        col = ncols(M1);
+        D = kohom(N,col);
+//---------- computing homology ----------------------------------------------
+        imag  = kontrahom(M1,row);
+        A     = kontrahom(M2,row);
+        B     = kohom(N,ncols(M2));
+        ker   = modulo(A,B);
+        imag  = imag,D;
+        extMN = modulo(ker,imag);
+        dbprint(p-1,"// Computing Ext^"+string(k)+":",
+         "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of coker(M),",
+         "// and 0<--coker(N)<--G0<--G1 a presentation of coker(N),",
+         "// then Hom(F"+string(k)+",G0)-->Hom(F"+string(k+1)+",G0) is given by:",A,
+         "// and Hom(F"+string(k-1)+",G0) + Hom(F"+string(k)+",G1)-->Hom(F"+string(k)+",G0) is given by:",imag,"");
+        extMN0 = std(extMN);
+        if( p>0 ) { "// degree of Ext^"+string(k)+":"; degree(extMN0);""; }
+//---------- more information -------------------------------------------------
+        if( size(#)>0 )
+        { if( vdim(extMN0) >= 0 )
+          { kb = kbase(extMN0);
+            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
+            dbpri(p-1,"// columns of matrix are kbase of Ext^"+
+                     string(k)+" in Hom(F"+string(k)+",G0)",kb,"");
+            if( p>0 )
+            { for (l=1;l<=ncols(kb);l=l+1)
+              {
+                "// element",l,"of kbase of Ext^"+string(k)+" in Hom(F"+string(k)+",G0)";
+                "// as matrix: F"+string(k)+"-->G0";
+                print(matrix(ideal(kb[l]),row,col));
+              }
+              "";
+            }
+            L3[ii] = matrix(kb);
+          }
+          L2[ii] = extMN0;
+        }
+        L1[ii] = extMN;
+      }
+    }
+  }
+  if( size(#) )
+  {  if( s>1 ) { L = L1,L2,L3; return(L); }
+     else { L = extMN,extMN0,matrix(kb); return(L); }
+  }
+  else
+  {  if( s>1 ) { return(L1); }
+     else { return(extMN); }
+  }
+}
 example
+{ "EXAMPLE:";  echo = 2;
+   ring r=0,(x,y),(c,dp);
+   ideal i=x2-y3;
+   qring q=std(i);
+   module M=[-x,y],[-y2,x];
+   module hom=Hom(M,M);
+   print(hom); newline;
+   ring s=3,(x,y,z),(c,dp);
+   ideal i=x2+y5+z4;i=jacob(i);
+   qring rq=std(i);
+   matrix M[2][2]=xy,x3,5y,z2,x2;
+   matrix N[4][4]=x,y,z,x2,xyx2y,y3,xz2,x2z,z3;
+   print(M);
+   print(N);
+   print(Hom(M,N));
+{"EXAMPLE:";   echo=2;
+  int p      = printlevel;
+  printlevel = 1;
+  ring r     = 0,(x,y),dp;
+  ideal i    = x2-y3;
+  ideal j    = x2-y5;
+  list E     = Ext(0..2,i,j);    // Ext^k(r/i,r/j) for k=0,1,2 over r
+  qring R    = std(i);
+  ideal j    = fetch(r,j);
+  module M   = [-x,y],[-y2,x];
+  printlevel = 2;
+  module E1  = Ext(1,M,j);       // Ext^1(R^2/M,R/j) over R=r/i
+  list l     = Ext(4,M,M,1);     // Ext^4(R^2/M,R^2/M) over R=r/i
+  printlevel = p;
+}
 ////////////////////////////////////////////////////////////////////////////////
+proc qmod (module M, module N)
+USAGE:   qmod(M,N); M,N=modules, N a submodule of M
+COMPUTE: presentation S of M/N, i.e. M/N<<--F<--[S], F free of rank = size(M)
+RETURN:  module  S
+{
+  return(lift(M,N)+syz(M));
+proc Hom (module M, module N, list #)
+USAGE:   Hom(M,N,[any]);  M,N=modules
+COMPUTE: A presentation of Hom(M',N'), M'=coker(M), N'=coker(N) as follows:
+         Let ...-->F1 --M-> F0-->M'-->0 and ...-->G1 --N-> G0-->N'-->0  be
+         presentations of M' and N'. Consider
+               0
+                                         |^              |^
+--> Hom(M',N') ----> Hom(F0,N') ----> Hom(F1,N')
+                                         |^              |^
+         (A:  induced by M)          Hom(F0,G0) --A-> Hom(F1,G0)
+                                         |^              |^
+         (B,C:induced by N)              |C              |B
+                                     Hom(F0,G1) ----> Hom(F1,G1)
+         Let D=modulo(A,B) and Hom=modulo(D,C), then we have exact sequences
+              R^p  --D-> Hom(F0,G0) --A-> Hom(F1,G0)/im(B)
+              R^q -Hom-> R^p --D-> Hom(F0,G0)/im(C) --A-> Hom(F1,G0)/im(B).
+         Hence Hom presents Hom(M',N')
+RETURN:  Hom, of type module, presentation of Hom(M',N') or,
+         in case of 3 arguments, a list:
+           [1] = Hom
+           [2] = SB of Hom
+           [3] = kbase of coker(Hom) (if finite dimensional), represented by
+                 elements in Hom(F0,G0) via mapping D
+DISPLAY: printlevel >=0: degree of Hom  (default)
+         printlevel >=1: D and C and kbase of coker(Hom) in Hom(F0,G0)
+         printlevel >=2: elements of kbase of coker(Hom) as matrix :F0-->G0
+NOTE:    DISPLAY is as described only for a direct call of 'Hom'. Calling 'Hom'
+         from another proc has the same effect as decreasing printlevel by 1.
+EXAMPLE: example Hom;  shows examples
+{
+//---------- initialisation ---------------------------------------------------
+  int l,p;
+  matrix kb;
+  module A,B,C,D,homMN,homMN0;
+  list L;
+//---------- computation of Hom -----------------------------------------------
+  B = kohom(N,ncols(M));
+  A = kontrahom(M,nrows(N));
+  C = kohom(N,nrows(M));
+  D = modulo(A,B);
+  homMN = modulo(D,C);
+  homMN0= std(homMN);
+  p = printlevel-voice+3;       // p=printlevel+1 (default: p=1)
+  if( p>=0 ) { "// degree of Hom:"; degree(homMN0); ""; }
+  dbprint(p-1,"// given ...--> F1 --M-> F0 -->M'--> 0 and ...--> G1 --N-> G0 -->N'--> 0,",
+         "// show D=ker(Hom(F0,G0) --> Hom(F1,G0)/im(Hom(F1,G1) --> Hom(F1,G0)))",D,
+         "// show C=im(Hom(F0,G1) --> Hom(F0,G0))",C,"");
+//---------- extra output if size(#)>0 ----------------------------------------
+  if( size(#)>0 )
+  {  if( vdim(homMN0)>0 )
+     {  kb = kbase(homMN0);
+        kb = matrix(D)*kb;
+        if( p>2 )
+        {  for (l=1;l<=ncols(kb);l=l+1)
+           {
+             "// element",l,"of kbase of Hom in Hom(F0,G0) as matrix: F0-->G0:";
+             print(matrix(ideal(kb[l]),nrows(N),nrows(M)));
+           }
+        }
+        else
+       { dbprint(p-1,"// columns of matrix are kbase of Hom in Hom(F0,G0)",kb); }
+        L=homMN,homMN0,kb; return(L);
+     }
+     L=homMN,homMN0; return(L);
+  }
+  return(homMN);
+}
+example
+{"EXAMPLE:";  echo = 2;
+  int p     = printlevel;
+  printlevel= 1;   //in 'example proc' printlevel has to be increased by 1
+  ring r    = 0,(x,y),dp;
+  ideal i   = x2-y3,xy;
+  qring q   = std(i);
+  ideal i   = fetch(r,i);
+  module M  = [-x,y],[-y2,x],[x3];
+  module H  = Hom(M,i);
+  print(H);
+  printlevel= 2;
+  list L    = Hom(M,i,1);"";
+  ring s    = 3,(x,y,z),(c,dp);
+  ideal i   = jacob(ideal(x2+y5+z4));
+  qring rq=std(i);
+  matrix M[2][2]=xy,x3,5y,4z,x2;
+  matrix N[3][2]=x2,x,y3,3xz,x2z,z;
+  print(M);
+  print(N);
+  list l=Hom(M,N,1);
+  printlevel = p;
+}
+////////////////////////////////////////////////////////////////////////////////
+proc homology (matrix A,matrix B,module M,module N,list #)
+USAGE:   homology(A,B,M,N);
+COMPUTE: Let M and N be submodules of R^m and R^n presenting M'=R^m/M, N'=R^n/N
+         (R=basering) and let A,B matrices inducing maps R^k--A-->R^m--B-->R^n.
+         Compute a presentation R^q --H-> R^m of the module
+              ker(B)/im(A) := ker(M'/im(A) --B--> N'/im(BM)+im(BA)).
+         If B induces a map M'--B-->N' (i.e BM=0) and if R^k--A-->M'--B-->N' is
+         a complex (i.e. BA=0) then ker(B)/im(A) is the homology of this complex
+RETURN:  module H, a presentation of ker(B)/im(A)
+NOTE:    homology returns a free module of rank m if ker(B)=im(A)
+EXAMPLE: example homology; shows examples
+{
+  module ker,ima;
+  ker = modulo(B,N);
+  ima = A,M;
+  return(modulo(ker,ima));
+}
+example
+{"EXAMPLE";    echo=2;
+  ring r;
+  ideal id=maxideal(4);
+  qring qr=std(id);
+  module N=maxideal(3)*freemodule(2);
+  module M=maxideal(2)*freemodule(2);
+  module B=[2x,0],[x,y],[z2,y];
+  module A=M;
+  degree(std(homology(A,B,M,N)));"";
+  ring s=0,x,ds;
+  qring qs=std(x4);
+  module A=[x];module B=A;
+  module M=[x3];module N=M;
+  homology(A,B,M,N);
+}
 //////////////////////////////////////////////////////////////////////////////
+proc kernel (matrix A,module M,module N)
+USAGE:   kernel(A,M,N);
+COMPUTE: Let M and N be submodules of R^m and R^n presenting M'=R^m/M, N'=R^n/N
+         (R=basering) and let A:R^m-->R^n a matrix inducing a map A':M'-->N'.
+         Compute a presentation K of ker(A') as in the commutative diagram:
+                       ker(A') --->  M' --A'--> N'
+                           |^        |^         |^
+                           |         |          |
+                          R^r  ---> R^m --A--> R^n
+                           |^        |^         |^
+                           |K        |M         |N
+                           |         |          |
+                          R^s  ---> R^p -----> R^q
+RETURN:  module K, a presentation of ker(A')
+EXAMPLE: example kernel; shows examples
+{
+  module M1 = modulo(A,N);
+  return(modulo(M1,M));
+}
+example
+{"EXAMPLE";    echo=2;
+  ring r;
+  module N=[2x,x],[0,y];
+  module M=maxideal(1)*freemodule(2);
+  matrix A[2][2]=2x,0,x,y,z2,y;
+  module K=kernel(A,M,N);
+  degree(std(K));
+  print(K);
+}
+////////////////////////////////////////////////////////////////////////////////
+proc kohom (matrix M, int j)
+USAGE:   kohom(A,k); A=matrix, k=integer
+RETURN:  matrix Hom(R^k,A) i.e. let A be a matrix defining a map: F1 --> F2 of
+         free R-modules, the matrix of Hom(R^k,F1)-->Hom(R^k,F2) is computed
+EXAMPLE: example kohom; shows an example
+{
+  if (j==1)
+  { return(M);}
+  if (j>1)
+  { return(outer(M,diag(1,j))); }
+  else { return(0);}
+}
+example
+{"EXAMPLE:";   echo=2;
+  ring r;
+  matrix n[2][3]=x,y,5,z,77,33;
+  print(kohom(n,3));
+}
+/////////////////////////////////////////////////////////////////////////////////
 proc kontrahom (matrix M, int j)
-USAGE:   qmod(M,N); M,N=modules, N a submodule of M
 USAGE:   kontrahom(A,k); A=matrix, k=integer
 RETURN:  Hom(A,P^k), i.e. let A be a matrix defining a map: F1 --> F2 of free
          P-modules, the matrix of Hom(F2,P^k)-->Hom(F1,P^k) is computed
+RETURN:  matrix Hom(A,R^k), i.e. let A be a matrix defining a map: F1 --> F2 of
+         free  R-modules, the matrix of Hom(F2,R^k)-->Hom(F1,R^k) is computed
 EXAMPLE: example kontrahom; shows an example
+{
+  return(transpose(outer(diag(1,j),M)));
+if (j==1)
+  { return(transpose(M));}
+  if (j>1)
+  { return(transpose(outer(diag(1,j),M)));}
+  else { return(0);}
+}
 example
 { "EXAMPLE:";  echo=2;
    ring r;
    matrix n[2][3]=x,y,5,z,77,33;
    print(kontrahom(n,3));
+{"EXAMPLE:";  echo=2;
+  ring r;
+  matrix n[2][3]=x,y,5,z,77,33;
+  print(kontrahom(n,3));
+}
 ///////////////////////////////////////////////////////////////////////////////
-proc kohom (matrix M, int j)
-USAGE:   kohom(A,k); A=matrix, k=integer
-RETURN:  Hom(P^k,A) i.e. let A be a matrix defining a map: F1 --> F2 of free
-         P-modules, the matrix of Hom(P^k,F1)-->Hom(P^k,F2) is computed
-EXAMPLE: example kohom; shows an example
+{
-    return(outer(M,diag(1,j)));
+}
-example
-{ "EXAMPLE:";   echo=2;
-   ring r;
-   matrix n[2][3]=x,y,5,z,77,33;
-   print(kohom(n,3));
+}
-///////////////////////////////////////////////////////////////////////////////

Singular/LIB/inout.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: inout.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
 //system("random",787422842);
 //(GMG,BM)
 ///////////////////////////////////////////////////////////////////////////////
 LIBRARY:  inout.lib     PROCEDURES FOR MANIPULATING IN- AND OUTPUT
+// $Id: inout.lib,v 1.2 1997-04-28 19:27:20 obachman Exp $
+// system("random",787422842);
+// (GMG/BM, last modified 22.06.96)
+///////////////////////////////////////////////////////////////////////////////
+LIBRARY:  inout.lib     PROCEDURES FOR MANIPULATING IN- AND OUTPUT
  allprint(list);        print list if ALLprint is defined, with pause if >0
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc allprint (list #)
+proc allprint (list #)
 USAGE:   allprint(L);  L list
 CREATE:  display L[1], L[2], ... if an integer with name ALLprint is defined,
          makes "pause",   if ALLprint > 0
+         makes "pause",   if ALLprint > 0
          listvar(matrix), if ALLprint = 2
 RETURN:  no return value
 …
+   {
       int i;
       for( i=1; i<=size(#); i++ ) { print(#[i]); }
+      for( i=1; i<=size(#); i=i+1 ) { print(#[i]); }
       if( ALLprint==2 ) { pause; listvar(matrix); }
       if( ALLprint >0 ) { pause; }
+   }
+   }
    return();
+}
+}
 example
 { "EXAMPLE:"; echo = 2;
 …
    allprint("M =",M);
    kill ALLprint;
+}
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 RETURN:  no return value
 NOTE:    this is uesful to control the printing of comments or partial results
          in a procedure, e.g. for debugging a procedure.
+         in a procedure, e.g. for debugging a procedure.
          It is similair but not the same as the internal function dbprint
 EXAMPLE: example dbpri; shows an example
+{
    int i;
+   int i;
    if( defined(printlevel)==0 ) { int printlevel; export printlevel; }
    if( n<=printlevel )
+   {
       for( i=1; i<=size(#); i++ ) {  print(#[i]); }
+      for( i=1; i<=size(#); i=i+1 ) {  print(#[i]); }
+   }
    return();
+}
+}
 example
 { "EXAMPLE:"; echo = 2;
    ring s;
    module M=freemodule(3);
+   module M=freemodule(3);
    dbpri(0,"M =",M);
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc lprint
+proc lprint
 USAGE:   lprint(id[,n]);  id poly/ideal/vector/module/matrix, n integer
 RETURN:  string of id in a format fitting into lines of size n; if only one
          argument is present, n = pagewidth
+RETURN:  string of id in a format fitting into lines of size n; if only one
+         argument is present, n = pagewidth
 NOTE:    id is printed columnwise, each column separated by a blank line;
          hence lprint(transpose(id)); displays a matrix id in a format which
 …
    matrix M = matrix(#[1]);
    poly p,h,L; string s1,s,S; int jj,ii,a;
    for (jj=1; jj<=ncols(M); jj++)
+   {
       for (ii=1; ii<=nrows(M); ii++)
+   for (jj=1; jj<=ncols(M); jj=jj+1)
+   {
+      for (ii=1; ii<=nrows(M); ii=ii+1)
+      {
          a=2;
 …
             while (p != 0)
+            {
                if (a+size(string(h+L)) > n)
+               {
+               if (a+size(string(h+L)) > n)
+               {
                   s = string(h);
                   if (a != 0) { s = "  "+s; }
 …
+         }
          if (ii != nrows(M)) { s = s+","; S=S+newline+s; }
          else
+         else
+         {
             if (jj != ncols(M)) { s = s+","; S=S+newline+s+newline;}
 …
    // above) as input in order to reproduce m (by defining m1):
    execute("matrix M[2][3]="+s+";");
    module m1 = transpose(M);
+   module m1 = transpose(M);
    m-m1;
+}
 …
          if n is given, only the first n characters of each colum are shown
 RETURN:  no return value
 EXAMPLE: example pmat; shows an example
+EXAMPLE: example pmat; shows an example
+{
 //------------- main case: input is a matrix, no second argument---------------
    if ( size(#)==0)
+   {
+   {
       int elems,mlen,slen,c,r;
    //-------------- count maximal size of each column, and sum up -------------
       for ( c=1; c<=ncols(m); c++)
+      for ( c=1; c<=ncols(m); c=c+1)
       {  int len(c);
          for (r=1; r<=nrows(m); r++)
+         for (r=1; r<=nrows(m); r=r+1)
+         {
             elems = elems+1;
 …
       string aus; string sep = " ";
       if (mlen >= pagewidth) { sep = newline; }
       for (r=1; r<nrows(m); r++)
+      for (r=1; r<nrows(m); r=r+1)
       {  elems = r; aus = "";
          for (c=1; c<=ncols(m); c++)
+         for (c=1; c<=ncols(m); c=c+1)
+         {
             aus = aus + s(elems)[1,len(c)] + sep;
 …
    //--------------- print last row (no comma after last entry) ---------------
       aus = ""; elems = nrows(m);
       for (c=1; c<ncols(m); c++)
+      for (c=1; c<ncols(m); c=c+1)
+      {
          aus = aus + s(elems)[1,len(c)] + sep;
 …
    if ( typeof(#[1])=="int" )
    {  string aus,tmp; int ll,c,r;
       for ( r=1; r<=nrows(m); r++)
+      for ( r=1; r<=nrows(m); r=r+1)
       {  aus = "";
          for (c=1; c<=ncols(m); c++)
+         {
             tmp=string(m[r,c]);
             aus=aus+tmp[1,#[1]]+" ";
+         for (c=1; c<=ncols(m); c=c+1)
+         {
+            tmp=string(m[r,c]);
+            aus=aus+tmp[1,#[1]]+" ";
+         }
          aus;
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc rMacaulay
+proc rMacaulay
 USAGE:   rMacaulay(s[,n]);  s string, n integer
 RETURN:  a string which should be readable by Singular if s is a string read
          by Singular from a file which was produced by Macaulay_1 (='Macaulay
          classic'). If a second argument is present the first n lines of the
          file are deleted (which is useful if the file was prodeuced e.g. by
          the putstd command of Macaulay)
+         classic'). If a second argument is present the first n lines of the
+         file are deleted (which is useful if the file was prodeuced e.g. by
+         the putstd command of Macaulay)
 NOTE:    This does not always work with 'cut and paste' since, coming
          from the screen, the character \ is treated differently
 EXAMPLE: example rMacaulay; shows an example
+         from the screen, the character \ is treated differently
+EXAMPLE: example rMacaulay; shows an example
+{
    int n;
 …
 //------------------------ delete first n=#[2] lines --------------------------
    int ii=find(s0,newline); int jj;
    for ( jj=1; jj<=n; jj++)
+   for ( jj=1; jj<=n; jj=jj+1)
+   {
       s0 = s0[ii+1,size(s0)-ii];
 …
       ii = find(s0,newline);
+   }
    jj=jj+1;
+   jj=jj+1;
    string s(jj) = s0;
 //------------ delete blanks and \ at end of each string and add , ------------
 …
       s1 = ""; s2 = s(ii);
       kk = find(s(ii)," "); ll=kk+1;
       while( kk!=0 )
+      while( kk!=0 )
+      {
          while( s2[ll]==" ") { ll=ll+1; }
 …
 example
 {  "EXAMPLE:"; echo = 2;
    // Assume there exists a file 'Macid' with the following ideal in Macaulay
+   // Assume there exists a file 'Macid' with the following ideal in Macaulay
    // format:"
    // x[0]3-101/74x[0]2x[1]+7371x[0]x[1]2-13/83x[1]3-x[0]2x[2] \
+   // x[0]3-101/74x[0]2x[1]+7371x[0]x[1]2-13/83x[1]3-x[0]2x[2] \
    //     -4/71x[0]x[1]x[2]-65/64x[1]2x[2]-49/111x[0]x[2]2-x[1]x[2]2 \
    //     -747x[2]3+6072x[0]2x[3]
    // You can read this file into Singular and assign it to the string s1 by:
+   //     -747x[2]3+6072x[0]2x[3]
+   // You can read this file into Singular and assign it to the string s1 by:
    // string s1 = read("Macid");
    // This is equivalent to";
 …
    // You may wish to assign s1 to a Singular ideal id:
    string sid = "ideal id =",rMacaulay(s1),";";
    ring r = 0,x(0..3),dp;
+   ring r = 0,x(0..3),dp;
    execute sid;
    id; "";
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc show
+proc show (id, list #)
 USAGE:   show(id);   id any object of basering or of type ring/qring
          show(R,s);  R=ring, s=string (`s` = name of an object belonging to R)
 CREATE:  display id resp.`s` in a compact format together with some information
+         show(R,s);  R=ring, s=string (s = name of an object belonging to R)
+CREATE:  display id/s in a compact format together with some information
 RETURN:  no return value
 NOTE:    objects of type string, int, intvec, intmat belong to any ring.
          id may be a ring or a qring. In this case the minimal polynomial is
          displayed, and, for a qring, also the defining ideal.
          id must not be of type list
 CAUTION: show(R,s) does not work inside a procedure
+         id may be a ring or a qring. In this case the minimal polynomial is
+         displayed, and, for a qring, also the defining ideal
+         id may be of type list but the list must not contain a ring
+CAUTION: show(R,s) does not work inside a procedure
 EXAMPLE: example show; shows an example
+{
+{
 //------------- use funny names in order to avoid name conflicts --------------
+   int @li@, @ii;
+   string @s@,@@s;
    int @short@=short; short=1;
-   string @s@;
 //----------------------------- check syntax ----------------------------------
+   if( size(#)==1 )
+   {
+       def @id@ = #[1];
+   }
+   if( size(#)==2 )
+   {
+      if( typeof(#[1])=="ring" or typeof(#[1])=="qring")
+      {
+         def @R@ = #[1];
+         setring @R@;
+         def @id@=`#[2]`;
+      }
+   }
+   if( defined(@id@)!=voice ) { "// wrong syntax, type help show;";  return(); }
+   if( size(#)!= 0 )
+   {
+      if( typeof(#[1])=="int" ) { @li@=#[1]; }
+   }
+   if ( typeof(id)!="list" )
+   {
+      if( size(#)==0 )
+      {
+          def @id@ = id;
+      }
+      if( size(#)==1 )
+      {
+         if( typeof(#[1])=="int" )
+         {
+             def @id@ = id;
+         }
+         if( typeof(#[1])=="string" )
+         {
+            if( typeof(id)=="ring" or typeof(id)=="qring")
+            {
+               def @R@ = id;
+               setring @R@;
+               def @id@=`#[1]`;
+            }
+         }
+      }
+   }
+//----------------------- case: id is of type list ----------------------------
+   if ( typeof(id)=="list" )
+   {
+      @@s = tab(@li@)+"// list, "+string(size(id))+" element(s):";
+      @@s;
+      for ( @ii=1; @ii<=size(id); @ii++ )
+      {
+         if( typeof(id[@ii])!="none" )
+         {
+            def @id(@ii) = id[@ii];
+            show(@id(@ii),@li@+3);
+         }
+         else { tab(@li@+2),"//",id[@ii]; }
+      }
+      short=@short@; return();
+    }
+   if( defined(@id@)!=voice ) { "// wrong syntax, type help show;";  return(); }
 //-------------------- case: @id@ belongs to any ring -------------------------
    if( typeof(@id@)=="string" or typeof(@id@)=="int" or typeof(@id@)=="intvec"
+       or typeof(@id@)=="intmat" or typeof(@id@)=="list" )
+   {
+      "//", typeof(@id@);
+       or typeof(@id@)=="intmat" or typeof(@id@)=="list" )
+   {
+      if( typeof(@id@)!="intmat" )
+      {
+         @@s = tab(@li@)+"// "+typeof(@id@)+", size "+string(size(@id@));
+         @@s;
+      }
+      if( typeof(@id@)=="intmat" )
+      {
+         @@s = tab(@li@)+"// "+typeof(@id@)+", "+string(nrows(@id@))+" rows, "
+               + string(ncols(@id@))+" columns";
+         @@s;
+      }
       @id@;
       short=@short@; return();
+      short=@short@; return();
+   }
 //-------------------- case: @id@ belongs to basering -------------------------
    if( typeof(@id@)=="poly" or typeof(@id@)=="ideal" or typeof(@id@)=="matrix" )
+   {
+      if( typeof(@id@)=="ideal" ) { @s@=", "+string(ncols(@id@))+" generators";}
+      if( typeof(@id@)=="matrix") { @s@=", "+string(nrows(@id@))+"x"+string(ncols(@id@));}
+      "// "+ typeof(@id@)+ @s@;
+      print(matrix(@id@)); short=@short@; return();
+      if( typeof(@id@)=="ideal" )
+      {
+         @s@=", "+string(ncols(@id@))+" generator(s)";
+      }
+      if( typeof(@id@)=="matrix")
+      {
+         @s@=", "+string(nrows(@id@))+"x"+string(ncols(@id@));
+      }
+      @@s = tab(@li@)+"// "+ typeof(@id@)+ @s@;
+      @@s;
+      print(matrix(@id@));
+      short=@short@; return();
+   }
    if( typeof(@id@)=="vector" )
+   {
+      "//", typeof(@id@);
+      print(@id@);  short=@short@; return();
+      @@s = tab(@li@)+"// "+typeof(@id@);
+      @@s;
+      print(@id@);
+      short=@short@; return();
+   }
    if( typeof(@id@)=="module" )
+   {
+      @s@=", "+string(ncols(@id@))+" generators";
+      "// "+ typeof(@id@)+ @s@;
+      @s@=", "+string(ncols(@id@))+" generator(s)";
+      @@s = tab(@li@)+"// "+ typeof(@id@)+ @s@;
+      @@s;
       int @n@;
       for( @n@=1; @n@<=ncols(@id@); @n@++ ) { print(@id@[@n@]); }
+      for( @n@=1; @n@<=ncols(@id@); @n@=@n@+1 ) { print(@id@[@n@]); }
       short=@short@; return();
+   }
+   if( typeof(@id@)=="number" ) { "//", typeof(@id@); @id@; short=@short@; return(); }
+   if( typeof(@id@)=="map" )
+   {
+   if( typeof(@id@)=="number" )
+   {
+      @@s = tab(@li@)+"//", typeof(@id@);
+      @@s;
+      @id@; short=@short@; return();
+   }
+   if( typeof(@id@)=="map" )
+   {
       def @map = @id@;
+      "// nth variable of preimage ring is mapped to nth component of map";
+      if( size(#)==1 ) { type @map; }
+      if( size(#)==2 ) { type `#[2]`; }
+      @@s = tab(@li@)+"// i-th variable of preimage ring is mapped to @map[i]";
+      @@s;
+      if( size(#)==0 ) { type @map; }
+      if( size(#)==1 )
+      {
+         if( typeof(#[1])=="int" )    { type @map; }
+         if( typeof(#[1])=="string" ) { type `#[1]`; }
+      }
       short=@short@; return();
+   }
 //---------------------- case: @id@ is a ring/qring ---------------------------
    if( typeof(@id@)=="ring" or typeof(@id@)=="qring" )
+   {
+   if( typeof(@id@)=="ring" or typeof(@id@)=="qring" )
+   {
       setring @id@;
       string s="("+charstr(@id@)+"),("+varstr(@id@)+"),("+ordstr(@id@)+");";
       if( typeof(@id@)=="ring" )
+      {
          "// ring:",s;
          "// minpoly =", minpoly;
+      if( typeof(@id@)=="ring" )
+      {
+         @@s = tab(@li@)+"// ring:"; @@s,s;
+         @@s = tab(@li@)+"// minpoly ="; @@s,minpoly;
+      }
       if( typeof(@id@)=="qring" )
+      {
+         "// qring:",s;
+         "// minpoly =", minpoly;
+         "// quotient ring from ideal:"; ideal(@id@);
+      }
+      short=@short@; return();
+      {
+         @@s = tab(@li@)+"// qring:"; @@s,s;
+         @@s = tab(@li@)+"// minpoly ="; @@s, minpoly;
+         @@s = tab(@li@)+"// quotient ring from ideal:"; @@s;
+         ideal(@id@);
+      }
+      short=@short@; //return();
+   }
+}
 …
     show(r);
     ideal i=x^3+y^5-6*z^3,xy,x3-y2;
     show(i);
+    show(i,3);
     vector v=x*gen(1)+y*gen(3);
     show(v);
     module m=v,2*v+gen(4);
     show(m);
     ring S=(0,T),(a,b,c,d),ws(1,2,3,4);
     minpoly = T2+1;
     ideal i=a2+b,c2+T2*d2; i=std(i);
     qring Q=i;
     show(Q);
+    module m=v,2*v+gen(4);
+    list L = i,v,m;
+    show(L);
+    ring S=(0,T),(a,b,c,d),ws(1,2,3,4);
+    minpoly = T^2+1;
+    ideal i=a2+b,c2+T^2*d2; i=std(i);
+    qring Q=i;
+    show(Q);
     map F=r,a2,b^2,3*c3;
     show(F);
 …
           (default: n=pagewidth)
 NOTE:     may be used in connection with lprint
 EXAMPLE:  example split; shows an example
+EXAMPLE:  example split; shows an example
+{
    string line,re; int p,l;
    if( size(#)==0 ) { int n=pagewidth; }
+   if( size(#)==0 ) { int n=pagewidth; }
    else { int n=#[1]; }
    if( s[size(s),1] != newline ) { s=s+newline; }
 …
+      }
       re=re+line[p,l-1]+"\\"+newline;
       l=size(line);
+      l=size(line);
       if( l>=size(s) ) break;
       s=s[l+1,size(s)-l];
 …
 proc writelist (string fil, string nam, list L)
 USAGE:   writelist(fil,nam,L);  fil,nam=strings (file name, list name), L=list
+USAGE:   writelist(fil,nam,L);  fil,nam=strings (file-name, list-name), L=list
 CREATE:  a file with name `fil`, write the content of the list L into it and
          call the list `nam`.
 RETURN:  no return value
 NOTE:    The syntax of writelist uses and is similar to the the syntax of the
+NOTE:    The syntax of writelist uses and is similar to the syntax of the
          write command of Singular which does not manage lists properly.
          If, say, (fil,nam) = ("listfile","L1"),  writelist creates (resp.
          appends if listfile exists) a file with name listfile and stores
          there the list L under the name L1. The Singular command
          execute(read("listfile")); assignes the content of L (stored in
+         If, say, (fil,nam) = ("listfile","L1"),  writelist creates (resp.
+         appends if listfile exists) a file with name listfile and stores
+         there the list L under the name L1. The Singular command
+         execute(read("listfile")); assignes the content of L (stored in
          listfile) to a list L1.
          On a UNIX system, overwrite an existing file if fil=">...", resp.
+         On a UNIX system, overwrite an existing file if fil=">...", resp.
          append if fil=">>...".
 EXAMPLE: example writelist; shows an example
+{
    int i;
    write fil,"list "+nam+";";
+   write(fil,"list "+nam+";");
    if( fil[1]==">" ) { fil=fil[2..size(fil)]; }
    if( fil[1]==">" ) { fil=fil[2..size(fil)]; }
    for( i=1;i<=size(L);i++ )
+   {
      write fil,"   "+nam+"["+string(i)+"]=",string(L[i])+";";
+   for( i=1;i<=size(L);i=i+1 )
+   {
+     write(fil,"   "+nam+"["+string(i)+"]=",string(L[i])+";");
+   }
    return();
 …
    writelist("L","L",L);
    read("L");
+   // You might wish to remove the just created files `zumSpass` and `L`!
+}
+///////////////////////////////////////////////////////////////////////////////
+   system("sh","/bin/rm L zumSpass");
+   // Under UNIX, this removes the files 'L' and 'zumSpass'
+   // Type help system; to get more information about the shell escape
+   // If your operating system does not accept the shell escape, you
+   // have to remove the just created files 'zumSpass' and 'L' directly
+}
+///////////////////////////////////////////////////////////////////////////////

Singular/LIB/matrix.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: matrix.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
 //(GMG+BM)
+// $Id: matrix.lib,v 1.2 1997-04-28 19:27:22 obachman Exp $
+// (GMG/BM, last modified 22.06.96)
 ///////////////////////////////////////////////////////////////////////////////
 LIBRARY:  matrix.lib    PROCEDURES FOR MATRIX OPERATIONS
+LIBRARY:  matrix.lib    PROCEDURES FOR MATRIX OPERATIONS
  compress(A);           matrix, zero columns from A deleted
 …
  dsum(A1,A2,..);        matrix, direct sum of matrices A1,A2,...
  flatten(A);            ideal, generated by entries of matrix A
+ genericmat(n,m[,id]);  generic nxm matrix [entries from id]
  is_complex(c);         1 if list c is a complex, 0 if not
  outer(A,B);            matrix, outer product of matrices A and B
  skewmat(n[,id]);       generic skew-symmetric nxn matrix [entries from id]
  submat(A,r,c);         submatrix of A with rows/cols specified by intvec r/c
+ submat(A,r,c);         submatrix of A with rows/cols specified by intvec r/c
  symmat(n[,id]);        generic symmetric nxn matrix [entries from id]
  tensor(A,B);           matrix, tensor product of matrices A nd B
 …
 proc compress (A)
+USAGE:   compress(A); A matrix/intmat/ideal/module
+RETURN:  matrix/intmat/ideal/module, zero columns/generators from A deleted
+USAGE:   compress(A); A matrix/ideal/module/intmat/intvec
+RETURN:  same type, zero columns/generators from A deleted
+         (in an intvec zero elements are deleted)
 EXAMPLE: example compress; shows an example
+{
    if( typeof(A)=="matrix" ) { return(matrix(simplify(A,2))); }
    if( typeof(A)=="intmat" )
+   if( typeof(A)=="matrix" ) { return(matrix(simplify(A,2))); }
+   if( typeof(A)=="intmat" or typeof(A)=="intvec" )
+   {
       ring r=0,x,lp;
       module m=module(matrix(A));
       int c= size(m);
       intmat B[nrows(A)][c];
+      if( typeof(A)=="intvec" ) { intmat C=transpose(A); kill A; intmat A=C; }
+      module m = matrix(A);
+      intmat B[nrows(A)][size(m)];
       int i,j;
       for( i=1; i<=ncols(A); i++ )
+      {
          if( m[i]!=[0] )
+         {
+      for( i=1; i<=ncols(A); i=i+1 )
+      {
+         if( m[i]!=[0] )
+         {
             j=j+1;
             B[1..nrows(A),j]=A[1..nrows(A),i];
+         }
+      }
+      if( defined(C) ) { return(intvec(B)); }
       return(B);
+    }
 …
    intmat B[3][4]=1,0,3,0,4,0,5,0,6,0,7,0;
    compress(B);
+   intvec C=0,0,1,2,0,3;
+   compress(C);
+}
 ////////////////////////////////////////////////////////////////////////////////
 proc concat (list #)
 USAGE:   concat(A1,A2,..); A1,A2,... matrices
 RETURN:  matrix, concatenation of A1,A2,... . Number of rows of result matrix is
+USAGE:   concat(A1,A2,..); A1,A2,... matrices
+RETURN:  matrix, concatenation of A1,A2,... . Number of rows of result matrix is
          max(nrows(A1),nrows(A2),...)
 EXAMPLE: example concat; shows an example
 …
    int i;
    module B=module(#[1]);
    for( i=2; i<=size(#); i++ ) { B=B,module(#[i]); }
+   for( i=2; i<=size(#); i=i+1 ) { B=B,module(#[i]); }
    return(matrix(B));
+}
 …
 proc diag (list #)
 USAGE:   diag(p,n); p poly, n integer
+USAGE:   diag(p,n); p poly, n integer
          diag(A);   A matrix
 RETURN:  diag(p,n): diagonal matrix, p times unitmatrix of size n
          diag(A)  : n*mxn*m diagonal matrix with entries all the entries of  the
+         diag(A)  : n*mxn*m diagonal matrix with entries all the entries of the
                     nxm matrix A, taken from the 1st row, 2nd row etc of A
 EXAMPLE: example diag; shows an example
 …
       int i; ideal id=#[1];
       int n=ncols(id); matrix A[n][n];
       for( i=1; i<=n; i++ ) { A[i,i]=id[i]; }
+      for( i=1; i<=n; i=i+1 ) { A[i,i]=id[i]; }
+   }
    return(A);
 …
 ////////////////////////////////////////////////////////////////////////////////
-proc flatten (matrix A)
-USAGE:   flatten(A); A matrix
-RETURN:  ideal, generated by all entries from A
-EXAMPLE: example flatten; shows an example
+{
-   return(ideal(A));
+}
-example
-{ "EXAMPLE:"; echo = 2;
-   ring r=0,(x,y,z),ds;
-   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
-   print(A);
-   flatten(A);
+}
-////////////////////////////////////////////////////////////////////////////////
 proc dsum (list #)
 USAGE:   dsum(A1,A2,..); A1,A2,... matrices
+USAGE:   dsum(A1,A2,..); A1,A2,... matrices
 RETURN:  matrix, direct sum of A1,A2,...
 EXAMPLE: example dsum; shows an example
 …
    int i,N,a;
    list L;
    for( i=1; i<=size(#); i++ ) { N=N+nrows(#[i]); }
    for( i=1; i<=size(#); i++ )
+   {
       matrix B[N][ncols(#[i])];
+   for( i=1; i<=size(#); i=i+1 ) { N=N+nrows(#[i]); }
+   for( i=1; i<=size(#); i=i+1 )
+   {
+      matrix B[N][ncols(#[i])];
       B[a+1..a+nrows(#[i]),1..ncols(#[i])]=#[i];
       a=a+nrows(#[i]);
 …
    print(dsum(A,B,C));
+}
+////////////////////////////////////////////////////////////////////////////////
+proc flatten (matrix A)
+USAGE:   flatten(A); A matrix
+RETURN:  ideal, generated by all entries from A
+EXAMPLE: example flatten; shows an example
+{
+   return(ideal(A));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=0,(x,y,z),ds;
+   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
+   print(A);
+   flatten(A);
+}
+////////////////////////////////////////////////////////////////////////////////
+proc genericmat (int n,int m,list #)
+USAGE:   genericmat(n,m[,id]);  n,m=integers, id=ideal
+RETURN:  nxm matrix, with entries from id (default: id=maxideal(1))
+NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
+         genericmat(n,m); creates the generic nxm matrix
+EXAMPLE: example genericmat; shows an example
+{
+   if( size(#)==0 ) { ideal id=maxideal(1); }
+   if( size(#)==1 ) { ideal id=#[1]; }
+   if( size(#)>=2 ) { "// give 3 arguments, 3-rd argument must be an ideal"; }
+   matrix B[n][m]=id;
+   return(B);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring R=0,x(1..16),lp;
+   print(genericmat(4,4));    // the generic 4x4 matrix
+   changevar("R1",A_Z("a",4),R);
+   matrix A=genericmat(4,5,maxideal(1)^3);
+   print(A);
+   int n,m=4,3;
+   ideal i = ideal(randommat(1,n*m,maxideal(1),9));
+   print(genericmat(n,m,i));  // matrix of generic linear forms
+   kill R1;
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc is_complex (list c)
 USAGE:   is_complex(c); c = list of size-compatible modules or matrices
 RETURN:  1 if c[i]*c[i+1]=0 for all i, 0 if not.
+RETURN:  1 if c[i]*c[i+1]=0 for all i, 0 if not.
 NOTE:    Ideals are treated internally as 1-line matrices
 EXAMPLE: example is_complex; shows an example
 …
+}
 example
 { "EXAMPLE:";   echo = 2;
    ring r=32003,(x,y,z),ds;
    ideal i=x4+y5+z6,xyz,yx2+xz2+zy7;
+{ "EXAMPLE:";   echo = 2;
+   ring r=32003,(x,y,z),ds;
+   ideal i=x4+y5+z6,xyz,yx2+xz2+zy7;
    list L=res(i,0);
    is_complex(L);
 …
+{
    int i,j; list L;
+   int N=nrows(A)*nrows(B);
+   matrix C[N][ncols(B)];
+   for( i=1; i<=ncols(A); i++ )
+   {
+      for( j=1; j<=nrows(A); j++ )
+      {
+         C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B;
+      }
+      L[i]=C;
+   }
+   return(concat(L));
+   int triv = nrows(B)*ncols(B);
+   if( triv==1 )
+   {
+     return(B[1,1]*A);
+   }
+   else
+   {
+     int N = nrows(A)*nrows(B);
+     matrix C[N][ncols(B)];
+     for( i=1; i<=ncols(A); i=i+1 )
+     {
+       for( j=1; j<=nrows(A); j=j+1 )
+       {
+          C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B;
+       }
+       L[i]=C;
+     }
+     return(concat(L));
+   }
+}
 example
 …
 proc skewmat (int n, list #)
 USAGE:   skewmat(n[,id]);  n integer, id ideal
 RETURN:  skew-symmetric nxn matrix, with entries from id
+RETURN:  skew-symmetric nxn matrix, with entries from id
          (default: id=maxideal(1))
 NOTE:    if id has less than n*(n-1)/2 elements, the matrix is filled with 0's,
 …
    id = id,B[1..n,1..n];
    int i,j;
    for( i=0; i<=n-2; i++ )
+   {
       B[i+1,i+2..n]=id[j+1..j+n-i-1];
       j=j+n-i-1;
+   for( i=0; i<=n-2; i=i+1 )
+   {
+      B[i+1,i+2..n]=id[j+1..j+n-i-1];
+      j=j+n-i-1;
+   }
    matrix A=transpose(B);
 …
 proc submat (matrix A, intvec r, intvec c)
 USAGE:   submat(A,r,c);  A=matrix, r,c=intvec
 RETURN:  matrix, submatrix of A with rows specified by intvec r and columns
+RETURN:  matrix, submatrix of A with rows specified by intvec r and columns
          specified by intvec c
 EXAMPLE: example submat; shows an example
 …
    id = id,B[1..n,1..n];
    int i,j;
    for( i=0; i<=n-1; i++ )
+   {
       B[i+1,i+1..n]=id[j+1..j+n-i];
       j=j+n-i;
+   for( i=0; i<=n-1; i=i+1 )
+   {
+      B[i+1,i+1..n]=id[j+1..j+n-i];
+      j=j+n-i;
+   }
    matrix A=transpose(B);
    for( i=1; i<=n; i++ ) {  A[i,i]=0; }
+   for( i=1; i<=n; i=i+1 ) {  A[i,i]=0; }
    B=A+B;
    return(B);
 …
 { "EXAMPLE:"; echo = 2;
    ring R=0,x(1..10),lp;
    print(symmat(4));    // the generic symmetric matrix
+   print(symmat(4));    // the generic symmetric matrix
    changevar("R1",A_Z("a",5),R);
    matrix A=symmat(5,maxideal(1)^2);
 …
    int i,j;
    matrix C=B;
    for( i=2; i<=nrows(A); i++ ) { C=dsum(C,B); }
+   for( i=2; i<=nrows(A); i=i+1 ) { C=dsum(C,B); }
    matrix D[nrows(C)][ncols(A)*nrows(B)];
    for( j=1; j<=nrows(B); j++ )
+   {
       for( i=1; i<=nrows(A); i++ )
+   for( j=1; j<=nrows(B); j=j+1 )
+   {
+      for( i=1; i<=nrows(A); i=i+1 )
+      {
          D[(i-1)*nrows(B)+j,(j-1)*ncols(A)+1..j*ncols(A)]=A[i,1..ncols(A)];
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc gauss_col (matrix m)
+proc gauss_col (matrix A)
 USAGE:   gauss_col(A); A=matrix with constant coefficients
+RETURN:  matrix = col-reduced normal form of A
+RETURN:  matrix = col-reduced lower-triagonal normal form of A
+NOTE:    the procedure sets the global option-command: option(noredSB);
 EXAMPLE: example gauss_col; shows an example
+{
    def R=basering;
    changeord("@R","ds,c",R);
+   changeord("@R","ds,c",R);
    option(redSB); option(nointStrategy);
    matrix m = imap(R,m);
    m = matrix(std(m),nrows(m),ncols(m));
    setring R;
    m=imap(@R,m);
+   matrix A = imap(R,A);
+   A = matrix(std(A),nrows(A),ncols(A));
+   setring R;
+   A=imap(@R,A);
    option(noredSB);
    kill @R;
+   return(m);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring S;
+   matrix M[3][4] = 1,3,2,4,2,6,4,8,1,3,4,4;
+   print(M);
+   print(gauss_col(M));
+   return(A);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring S=0,x,dp;
+   matrix A[5][4] =  3, 1,1,-1,
+, 8,6,-7,
+,10,6,-7,
+, 4,3,-3,
+, 1,0, 3;
+   print(gauss_col(A));
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc gauss_row (matrix m)
+proc gauss_row (matrix A)
 USAGE:   gauss_row(A); A=matrix with constant coefficients
+RETURN:  matrix = row-reduced normal form of A
+RETURN:  matrix = row-reduced upper-triangular normal form of A
+NOTE:    may be used to solve a system of linear equations
+         see proc 'linearsolve' from 'solve.lib' for a different method
+         the procedure sets the global option-command: option(noredSB);
 EXAMPLE: example gauss_row; shows an example
+{
+   m = gauss_col(transpose(m));
+   return(transpose(m));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring S;
+   matrix M[3][4] = 1,3,2,4,2,6,4,8,1,3,4,4;
+   print(M);
+   print(gauss_row(M));
+   A = gauss_col(transpose(A));
+   return(transpose(A));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring S=0,x,dp;
+   matrix A[4][5] =  3, 1,1,-1,2,
+, 8,6,-7,1,
+,10,6,-7,1,
+, 4,3,-3,3;
+   print(gauss_row(A));
+}
 ////////////////////////////////////////////////////////////////////////////////

Singular/LIB/poly.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: poly.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
+// $Id: poly.lib,v 1.2 1997-04-28 19:27:22 obachman Exp $
 //system("random",787422842);
+//(GMG)
+///////////////////////////////////////////////////////////////////////////////
+LIBRARY:  poly.lib      PROCEDURES FOR MANIPULATING POLYS AND IDEALS
+ ishomog(poly/...);     int, =1 resp. =0 if input is homogeneous resp. not
+ maxcoef(poly/...);     maximal length of coefficient occuring in poly/...
+//(GMG, last modified 22.06.96)
+///////////////////////////////////////////////////////////////////////////////
+LIBRARY:  poly.lib      PROCEDURES FOR MANIPULATING POLYS, IDEALS, MODULES
+ cyclic(int);           ideal of cyclic n-roots
+ freerank(poly/...)     rank of coker(input) if coker is free else -1
+ is_homog(poly/...);    int, =1 resp. =0 if input is homogeneous resp. not
+ is_zero(poly/...);     int, =1 resp. =0 if coker(input) is 0 resp. not
+ maxcoef(poly/...);     maximal length of coefficient occuring in poly/...
  maxdeg(poly/...);      int/intmat = degree/s of terms of maximal order
+ maxdeg1(poly/...);     int = [weighted] maximal degree of input
  mindeg(poly/...);      int/intmat = degree/s of terms of minimal order
+ mindeg1(poly/...);     int = [weighted] minimal degree of input
  normalize(poly/...);   normalize poly/... such that leading coefficient is 1
            (parameters in square brackets [] are optional)
 …
 ///////////////////////////////////////////////////////////////////////////////
+proc ishomog (id)
+USAGE:   ishomog(id);  id  poly/ideal/vector/module/matrix
+proc cyclic (int n)
+USAGE:   cyclic(n);  n integer
+RETURN:  ideal of cyclic n-roots from 1-st n variables of basering
+EXAMPLE: example cyclic; shows examples
+{
+//----------------------------- procedure body --------------------------------
+   ideal m = maxideal(1);
+   m = m[1..n],m[1..n];
+   int i,j;
+   ideal s; poly t;
+   for ( j=0; j<=n-2; j=j+1 )
+   {
+      t=0;
+      for( i=1;i<=n;i=i+1 ) { t=t+product(m,i..i+j); }
+      s=s+t;
+   }
+   s=s,product(m,1..n)-1;
+   return (s);
+}
+//-------------------------------- examples -----------------------------------
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=0,(u,v,w,x,y,z),lp;
+   cyclic(nvars(basering));
+   homog(cyclic(5),z);
+}
+///////////////////////////////////////////////////////////////////////////////
+proc freerank
+USAGE:   freerank(M[,any]);  M=poly/ideal/vector/module/matrix
+COMPUTE: rank of module presented by M in case it is free. By definition this
+         is vdim(coker(M)/m*coker(M)) if coker(M) is free, where m = maximal
+         ideal of basering and M is considered as matrix (the 0-module is
+         free of rank 0)
+RETURN:  rank of coker(M) if coker(M) is free and -1 else;
+         in case of a second argument return a list:
+                L[1] = rank of coker(M) or -1
+                L[2] = minbase(M)
+NOTE:    freerank(syz(M)); computes the rank of M if M is free (and -1 else)
+         //* Zur Zeit noch ein Bug, da erste Bettizahl falsch berechnet wird:
+         //betti(0) ist -1 statt 0
+EXAMPLE: example freerank; shows examples
+{
+  int rk;
+  def M = simplify(#[1],10);
+  list mre = mres(M,2);
+  intmat B = betti(mre);
+  if ( ncols(B)>1 ) { rk = -1; }
+  else { rk = sum(B[1..nrows(B),1]); }
+  if (size(#) == 2) { list L=rk,mre[1]; return(L);}
+  return(rk);
+}
+example
+{"EXAMPLE";   echo=2;
+  ring r;
+  ideal i=x;
+  module M=[x,0,1],[-x,0,-1];
+  freerank(M);           // should be -1, coker(M) is not free
+                         // [1] should be 1, coker(syz(M))=M is free of rank 1
+  freerank(syz (M),"");  // [2] should be gen(2)+gen(1) (minimal relation of M)
+  freerank(i);
+  freerank(syz(i));      //* bug, should be 1, coker(syz(i))=i is free of rank 1
+}
+///////////////////////////////////////////////////////////////////////////////
+proc is_homog (id)
+USAGE:   is_homog(id);  id  poly/ideal/vector/module/matrix
 RETURN:  integer which is 1 if input is homogeneous (resp. weighted homogeneous
          if the monomial ordering consists of one block of type ws,Ws,wp or Wp,
 …
          degree, a module/matrix is homogeneous if all column vectors are
          homogeneous
+         //*** ergaenzen, wenn Matrizen Spalten Gewichte haben
+EXAMPLE: example ishomog; shows an example
+{
+         //*** ergaenzen, wenn Matrizen-Spalten Gewichte haben
+EXAMPLE: example is_homog; shows examples
+{
+//----------------------------- procedure body --------------------------------
    module M = module(matrix(id));
    M = simplify(M,2);                        // remove 0-columns
    intvec v = ringweights(basering);
    int i,j=1,1;
    for (i=1; i<=ncols(M); i++)
+   for (i=1; i<=ncols(M); i=i+1)
+   {
       if( M[i]!=jet(M[i],deg(lead(M[i])),v)-jet(M[i],deg(lead(M[i]))-1,v))
 …
    return(1);
+}
+//-------------------------------- examples -----------------------------------
 example
 { "EXAMPLE:"; echo = 2;
    ring r = 0,(x,y,z),wp(1,2,3);
    ishomog(x5-yz+y3);
+   is_homog(x5-yz+y3);
    ideal i = x6+y3+z2, x9-z3;
    ishomog(i);
+   is_homog(i);
    ring s = 0,(a,b,c),ds;
    vector v = [a2,0,ac+bc];
    vector w = [a3,b3,c4];
+   ishomog(v);
+   ishomog(w);
+}
+///////////////////////////////////////////////////////////////////////////////
+   is_homog(v);
+   is_homog(w);
+}
+///////////////////////////////////////////////////////////////////////////////
+proc is_zero
+USAGE:   is_zero(M[,any]); M=poly/ideal/vector/module/matrix
+RETURN:  integer, 1 if coker(M)=0 resp. 0 if coker(M)!=0, where M is considered
+         as matrix
+         if a second argument is given, return a list:
+                L[1] = 1 if coker(M)=0 resp. 0 if coker(M)!=0
+                L[2] = dim(M)
+EXAMPLE: example is_zero; shows examples
+{
+  int d=dim(std(#[1]));
+  int a = ( d==-1 );
+  if( size(#) >1 ) { list L=a,d; return(L); }
+  return(a);
+}
+example
+{ "EXAMPLE:";   echo=2;
+  ring r;
+  module m = [x],[y],[1,z];
+  is_zero(m,1);
+  qring q = std(ideal(x2+y3+z2));
+  ideal j = x2+y3+z2-37;
+  is_zero(j);
+}
+////////////////////////////////////////////////////////////////////////////////
 proc maxcoef (f)
 USAGE:   maxcoef(f);  f  poly/ideal/vector/module/matrix
 RETURN:  maximal length of coefficient of f of type int (by counting the
+RETURN:  maximal length of coefficient of f of type int (by counting the
          length of the string of each coefficient)
+EXAMPLE: example maxcoef; shows an example
+{
+EXAMPLE: example maxcoef; shows examples
+{
+//----------------------------- procedure body --------------------------------
    int max,s,ii,jj; string t;
    ideal i = ideal(matrix(f));
    i = simplify(i,6);         //* delete 0's and keep first of equal elements
+   i = simplify(i,6);            // delete 0's and keep first of equal elements
    poly m = var(1); matrix C;
    for (ii=2;ii<=nvars(basering);ii++) { m = m*var(ii); }
    for (ii=1; ii<=size(i); ii++)
+   for (ii=2;ii<=nvars(basering);ii=ii+1) { m = m*var(ii); }
+   for (ii=1; ii<=size(i); ii=ii+1)
+   {
       C = coef(i[ii],m);
       for (jj=1; jj<=ncols(C); jj++)
+      for (jj=1; jj<=ncols(C); jj=jj+1)
+      {
          t = string(C[2,jj]);  s = size(t);
 …
    return(max);
+}
+//-------------------------------- examples -----------------------------------
 example
 { "EXAMPLE:"; echo = 2;
 …
    poly g = 345x2-1234567890y+7/4z;
    maxcoef(g);
    ideal i = g,10/1234567890;
    maxcoef(i);
    // since i[2]=1/123456789
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc maxdeg (id)
+   ideal i = g,10/1234567890;
+   maxcoef(i);
+   // since i[2]=1/123456789
+}
+///////////////////////////////////////////////////////////////////////////////
+proc maxdeg (id)
 USAGE:   maxdeg(id);  id  poly/ideal/vector/module/matrix
+RETURN:  maximal degree/s of monomials of id independent of ring ordering
+         (maxdeg of each variable is 1)
+RETURN:  int/intmat, each component equals maximal degree of monomials in the
+         corresponding component of id, independent of ring ordering
+         (maxdeg of each var is 1)
          of type int if id is of type poly, of type intmat else
+EXAMPLE: example maxdeg; shows an example
+{
+//------------------- find maximal degree of given component ------------------
+NOTE:    proc maxdeg1 returns 1 integer, the absolut maximum; moreover, it has
+         an option for computing weighted degrees
+EXAMPLE: example maxdeg; shows examples
+{
+   //-------- subprocedure to find maximal degree of given component ----------
    proc findmaxdeg
+   {
 …
    int r,c = nrows(M), ncols(M); int i,j;
    intmat m[r][c];
    for (i=r; i>0; i--)
+   {
       for (j=c; j>0; j--) { m[i,j] = findmaxdeg(M[i,j]); }
+   }
    if( typeof(id)=="poly" ) { return(m[1,1]); }
+   for (i=r; i>0; i=i-1)
+   {
+      for (j=c; j>0; j=j-1) { m[i,j] = findmaxdeg(M[i,j]); }
+   }
+   if (typeof(id)=="poly") { return(m[1,1]); }
    return(m);
+}
+//-------------------------------- examples -----------------------------------
 example
 { "EXAMPLE:"; echo = 2;
    ring r = 0,(x,y,z),wp(-1,-2,-3);
    poly f = x+y2+z3;
    deg(f);               //deg returns weighted degree (in case of 1 block)!
+   deg(f);               //deg; returns weighted degree (in case of 1 block)!
    maxdeg(f);
    matrix m[2][2]=f+x10,1,0,f^2;
 …
 ///////////////////////////////////////////////////////////////////////////////
+proc maxdeg1 (id,list #)
+USAGE:   maxdeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
+RETURN:  integer, maximal [weighted] degree of monomials of id independent of
+         ring ordering, maxdeg1 of i-th variable is v[i] (default: v=1..1).
+NOTE:    This proc returns one integer while maxdeg returns, in general,
+         a matrix of integers. For one polynomial and if no intvec v is given
+         maxdeg is faster
+EXAMPLE: example maxdeg1; shows examples
+{
+   //-------- subprocedure to find maximal degree of given component ----------
+   proc findmaxdeg
+   {
+      poly c = #[1];
+      if (c==0) { return(-1); }
+      intvec v = #[2];
+   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
+      int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c);
+      int i = d;
+      if ( c == jet(c,-1,v))      //case: maxdeg is negative
+      {
+         i = -d;
+         while ( c == jet(c,i,v) ) { i = 2*(i-1); }
+         int o = (d != -i)*((i/ 2)+2) - 1;
+         int u = i+1;
+         int e = -1;
+      }
+      else                        //case: maxdeg is nonnegative
+      {
+         while ( c != jet(c,i,v) ) { i = 2*(i+1); }
+         int o = i-1;
+         int u = (d != i)*((i/ 2)-1);
+         int e = 1;
+      }
+   //----------------------- "quick search" for maxdeg ------------------------
+      while ( ( c==jet(c,i,v) )*( c!=jet(c,i-1,v) ) == 0 )
+      {
+         i = (o+e+u)/ 2;
+         if ( c!=jet(c,i,v) ) { u = i+1; }
+         else { o = i-1; }
+      }
+      return(i);
+   }
+//------------------------------ main program ---------------------------------
+   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
+   int c = ncols(M);
+   int i,n;
+   if( size(#)==0 )
+   {
+      int m = maxdeg(M[c]);
+      for (i=c-1; i>0; i=i-1)
+      {
+          n = maxdeg(M[i]);
+          m = (m>=n)*m + (m<n)*n;             //let m be the maximum of m and n
+      }
+   }
+   else
+   {
+      intvec v=#[1];                          //weight vector for the variables
+      int m = findmaxdeg(M[c],v);
+      for (i=c-1; i>0; i--)
+      {
+         n = findmaxdeg(M[i],v);
+         if( n>m ) { m=n; }
+      }
+   }
+   return(m);
+}
+//-------------------------------- examples -----------------------------------
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r = 0,(x,y,z),wp(-1,-2,-3);
+   poly f = x+y2+z3;
+   deg(f);                  //deg returns weighted degree (in case of 1 block)!
+   maxdeg1(f);
+   intvec v = ringweights(r);
+   maxdeg1(f,v);                             //weighted maximal degree
+   matrix m[2][2]=f+x10,1,0,f^2;
+   maxdeg1(m,v);                             //absolut weighted maximal degree
+}
+///////////////////////////////////////////////////////////////////////////////
 proc mindeg (id)
 USAGE:   mindeg(id);  id  poly/ideal/vector/module/matrix
+RETURN:  minimal degree/s of monomials of id independent of ring ordering
+         (mindeg of each variable is 1)
+         of type int if id is of type poly, of type intmat else
+EXAMPLE: example mindeg; shows an example
+{
+//------------------- find minimal degree of given component ------------------
+RETURN:  minimal degree/s of monomials of id, independent of ring ordering
+         (mindeg of each variable is 1) of type int if id of type poly, else
+         of type intmat.
+NOTE:    proc mindeg1 returns one integer, the absolut minimum; moreover it
+         has an option for computing weighted degrees.
+EXAMPLE: example mindeg; shows examples
+{
+   //--------- subprocedure to find minimal degree of given component ---------
    proc findmindeg
+   {
 …
       int o = i-1;
       int u = (d != i)*((i/ 2)-1);
 //----------------------- "quick search" for mindeg --------------------------
+   //----------------------- "quick search" for mindeg ------------------------
       while ( (jet(c,u)==0)*(jet(c,o)!=0) )
+      {
 …
    int r,c = nrows(M), ncols(M); int i,j;
    intmat m[r][c];
    for (i=r; i>0; i--)
+   {
       for (j=c; j>0; j--) { m[i,j] = findmindeg(M[i,j]); }
+   for (i=r; i>0; i=i-1)
+   {
+      for (j=c; j>0; j=j-1) { m[i,j] = findmindeg(M[i,j]); }
+   }
    if (typeof(id)=="poly") { return(m[1,1]); }
    return(m);
+}
+//-------------------------------- examples -----------------------------------
 example
 { "EXAMPLE:"; echo = 2;
    ring r = 0,(x,y,z),ls;
    poly f = x5+y2+z3;
    ord(f);                // ord returns weighted order of leading term!
    mindeg(f);
+   ord(f);                      // ord returns weighted order of leading term!
+   mindeg(f);                   // computes minimal degree
    matrix m[2][2]=x10,1,0,f^2;
+   mindeg(m);
+   mindeg(m);                   // computes matrix of minimum degrees
+}
+///////////////////////////////////////////////////////////////////////////////
+proc mindeg1 (id, list #)
+USAGE:   mindeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
+RETURN:  integer, minimal [weighted] degree of monomials of id independent of
+         ring ordering, mindeg1 of i-th variable is v[i] (default v=1..1).
+NOTE:    This proc returns one integer while mindeg returns, in general,
+         a matrix of integers. For one polynomial and if no intvec v is given
+         mindeg is faster.
+EXAMPLE: example mindeg1; shows examples
+{
+   //--------- subprocedure to find minimal degree of given component ---------
+   proc findmindeg
+   {
+      poly c = #[1];
+      intvec v = #[2];
+      if (c==0) { return(-1); }
+   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
+      int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c);
+      int i = d;
+      if ( jet(c,-1,v) !=0 )      //case: mindeg is negative
+      {
+         i = -d;
+         while ( jet(c,i,v) != 0 ) { i = 2*(i-1); }
+         int o = (d != -i)*((i/ 2)+2) - 1;
+         int u = i+1;
+         int e = -1; i=u;
+      }
+      else                        //case: inded is nonnegative
+      {
+         while ( jet(c,i,v) == 0 ) { i = 2*(i+1); }
+         int o = i-1;
+         int u = (d != i)*((i/ 2)-1);
+         int e = 1; i=u;
+      }
+   //----------------------- "quick search" for mindeg ------------------------
+      while ( (jet(c,i-1,v)==0)*(jet(c,i,v)!=0) == 0 )
+      {
+         i = (o+e+u)/ 2;
+         if (jet(c,i,v)==0) { u = i+1; }
+         else { o = i-1; }
+      }
+      return(i);
+   }
+//------------------------------ main program ---------------------------------
+   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
+   int c = ncols(M);
+   int i,n;
+   if( size(#)==0 )
+   {
+      int m = mindeg(M[c]);
+      for (i=c-1; i>0; i=i-1)
+      {
+          n = mindeg(M[i]);
+          m = (m<=n)*m + (m>n)*n;             //let m be the maximum of m and n
+      }
+   }
+   else
+   {
+      intvec v=#[1];                          //weight vector for the variables
+      int m = findmindeg(M[c],v);
+      for (i=c-1; i>0; i=i-1)
+      {
+         n = findmindeg(M[i],v);
+         m = (m<=n)*m + (m>n)*n;              //let m be the maximum of m and n
+      }
+   }
+   return(m);
+}
+//-------------------------------- examples -----------------------------------
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r = 0,(x,y,z),ls;
+   poly f = x5+y2+z3;
+   ord(f);                      // ord returns weighted order of leading term!
+   intvec v = 1,-3,2;
+   mindeg1(f,v);                // computes minimal weighted degree
+   matrix m[2][2]=x10,1,0,f^2;
+   mindeg1(m,1..3);             // computes absolut minimum of weighted degrees
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 RETURN:  object of same type with leading coefficient equal to 1
 EXAMPLE: example normalize; shows an example
+{
+{
    return(simplify(id,1));
+}
+//-------------------------------- examples -----------------------------------
 example
 {  "EXAMPLE:"; echo = 2;
 …
    module m=[9xy,0,3z3],[4z,6y,2x];
    normalize(m);
+   normalize(matrix(m));  // by automatic type conversion to module!
+}
+///////////////////////////////////////////////////////////////////////////////
+   normalize(matrix(m));             // by automatic type conversion to module!
+}
+///////////////////////////////////////////////////////////////////////////////

Singular/LIB/random.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: random.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
+// $Id: random.lib,v 1.2 1997-04-28 19:27:23 obachman Exp $
 //system("random",787422842);
 //(GMG+BM)
+//(GMG/BM, last modified 22.06.96)
 ///////////////////////////////////////////////////////////////////////////////
 …
 proc genericid (id, list #)
 USAGE:   genericid(id,[,p,b]);  id ideal/module, p,b integers
 RETURN:  system of generators of id which are generic, sparse, trigonal linear
+USAGE:   genericid(id,[,p,b]);  id ideal/module, k,p,b integers
+RETURN:  system of generators of id which are generic, sparse, triagonal linear
          combinations of given generators with coefficients in [1,b] and
          sparsety p percent, bigger p being sparser (default: p=75, b=30000)
 …
    if( size(#)==0 ) { int p=75; int b=30000; }
 //---------------- use sparsetriag for creation of genericid ------------------
    def i = simplify(id,10);
+   def i = simplify(id,10);
    i = i*sparsetriag(ncols(i),ncols(i),p,b);
    return(i);
+}
 example
 { "EXAMPLE:"; echo = 2;
    ring r=0,(t,x,y,z),ds;
    ideal i= x3+y4,z4+yx,t+x+y+z;
    genericid(i,0,10);
    module m=[x,0,0,0],[0,y2,0,0],[0,0,z3,0],[0,0,0,t4];
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=0,(t,x,y,z),ds;
+   ideal i= x3+y4,z4+yx,t+x+y+z;
+   genericid(i,0,10);
+   module m=[x,0,0,0],[0,y2,0,0],[0,0,z3,0],[0,0,0,t4];
    print(genericid(m));
+}
 …
    if( size(#)==0 ) { int k=size(id); int b=30000; }
 //--------------------------- create randomid ---------------------------------
    def i = id;
+   def i = id;
    i = matrix(id)*random(b,ncols(id),k);
    return(i);
+}
 example
 { "EXAMPLE:"; echo = 2;
    ring r=0,(x,y,z),dp;
    randomid(maxideal(2),2,9);
    module m=[x,0,1],[0,y2,0],[y,0,z3];
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=0,(x,y,z),dp;
+   randomid(maxideal(2),2,9);
+   module m=[x,0,1],[0,y2,0],[y,0,z3];
    show(randomid(m));
+}
 …
    int g=ncols(id);
    matrix rand[n][m]; matrix ra[1][m];
    for (int k=1; k<=n; k++)
+   for (int k=1; k<=n; k=k+1)
+   {
       ra = id*random(b,g,m);
 …
    A=randommat(2,3);
    print(A);
+}
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc sparseid (int k, int u, list #)
 USAGE:   sparseid(k,u[,o,p,b]);  k,u,o,p,b integers
 RETURN:  ideal having k generators in each degree d, u<=d<=o, p percent of
          terms in degree d are 0, the remaining have random coefficients
+RETURN:  ideal having k generators in each degree d, u<=d<=o, p percent of
+         terms in degree d are 0, the remaining have random coefficients
          in the interval [1,b], (default: o=u=d, p=75, b=30000)
 EXAMPLE: example sparseid; shows an example
 …
 //------------------ use sparsemat for creation of sparseid -------------------
    int ii; ideal i; intmat m;
    for ( ii=u; ii<=o; ii++)
+   {
+   for ( ii=u; ii<=o; ii=ii+1)
+   {
        m = sparsemat(size(maxideal(ii)),k,p,b);
        i = i+ideal(matrix(maxideal(ii))*m);
 …
    sparseid(3,4);"";
    sparseid(2,2,5,90,9);
+}
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
    if( p<40 )
+   {
       for( ii=1; ii<=t; ii++ )
+      for( ii=1; ii<=t; ii=ii+1 )
       { r=( random(1,100)>p ); v[1,ii]=r*random(1,b); h=h+r; }
+   }
 …
 example
 { "EXAMPLE:"; echo = 2;
    sparsemat(5,5);
    sparsemat(5,5,95);
    sparsemat(5,5,5);
+   sparsemat(5,5);"";
+   sparsemat(5,5,95);"";
+   sparsemat(5,5,5);"";
    sparsemat(5,5,50,100);
+}
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 USAGE:   sparsepoly(u[,o,p,b]);  u,o,p,b integers
 RETURN:  poly having only terms in degree d, u<=d<=o, p percent of the terms
          in degree d are 0, the remaining have random coefficients in [1,b),
+         in degree d are 0, the remaining have random coefficients in [1,b),
          (defaults: o=u=d, p=75, b=30000)
 EXAMPLE:  example sparsepoly; shows an example
 …
    int ii; poly f;
 //----------------- use sparseid for creation of sparsepoly -------------------
    for( ii=u; ii<=o; ii++ ) { f=f+sparseid(1,ii,ii,p,b)[1]; }
+   for( ii=u; ii<=o; ii=ii+1 ) { f=f+sparseid(1,ii,ii,p,b)[1]; }
    return(f);
+}
 …
    sparsepoly(5);"";
    sparsepoly(3,5,90,9);
+}
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
    if( n<=m ) { min=n-1; M[n,n]=1; }
    else { min=m; }
    for( ii=1; ii<=min; ii++ )
+   for( ii=1; ii<=min; ii=ii+1 )
+   {
       l=r+1; r=r+n-ii;
 …
 example
 { "EXAMPLE:"; echo = 2;
    sparsetriag(5,7);
    sparsetriag(7,5,90);
    sparsetriag(5,5,0);
+   sparsetriag(5,7);"";
+   sparsetriag(7,5,90);"";
+   sparsetriag(5,5,0);"";
    sparsetriag(5,5,50,100);
+}
 ///////////////////////////////////////////////////////////////////////////////
+}
+///////////////////////////////////////////////////////////////////////////////

Singular/LIB/redfac.lib

r6d09c56	r6f2edc
1		// $Id: redfac.lib,v 1.~~1.1.1 1997-04-25 15:13:26~~ obachman Exp $
	1	// $Id: redfac.lib,v 1.2 1997-04-28 19:27:24 obachman Exp $
2	2	//(RS)
3	3	///////////////////////////////////////////////////////////////////////////////

Singular/LIB/ring.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: ring.lib,v 1.1.1.1 1997-04-25 15:13:27 obachman Exp $
 //(GMG, 95/11/3)
 ///////////////////////////////////////////////////////////////////////////////
 LIBRARY:  ring.lib      PROCEDURES FOR MANIPULATING RINGS AND MAPS
+// $Id: ring.lib,v 1.2 1997-04-28 19:27:25 obachman Exp $
+//(GMG, last modified 03.11.95)
+///////////////////////////////////////////////////////////////////////////////
+LIBRARY:  ring.lib      PROCEDURES FOR MANIPULATING RINGS AND MAPS
  changechar("R",c[,r]); make a copy R of basering [ring r] with new char c
 …
 USAGE:   changechar(newr,c[,r]);  newr,c=strings, r=ring
 CREATE:  create a new ring with name `newr` and make it the basering if r is
          an existing ring [default: r=basering].
+         an existing ring [default: r=basering].
          The new ring differs from the old ring only in the characteristic.
          If, say, (newr,c) = ("R","0,A") and the ring r exists, the new
+         If, say, (newr,c) = ("R","0,A") and the ring r exists, the new
          basering will have name R characteristic 0 and one parameter A.
 RETURN:  No return value
 NOTE:    Works for qrings if map from old_char to new_char is implemented
          This proc uses 'execute' or calls a procedure using 'execute'.
          If you use it in your own proc, let the local names of your proc
+         This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
          start with @ (see the file HelpForProc)
 EXAMPLE: example changechar; shows an example
 …
    setring @r;
    ideal i = ideal(@r); int @q = size(i);
    if( @q!=0 )
+   if( @q!=0 )
       { string @s = "@newr1"; }
    else
+   else
       { string @s = newr; }
    string @newring = @s+"=("+c+"),("+varstr(@r)+"),("+ordstr(@r)+");";
 …
       ideal i = phi(i);
       attrib(i,"isSB",1);         //*** attrib funktioniert ?
       execute("qring "+newr+"=i;");
+      execute("qring "+newr+"=i;");
+   }
    export(`newr`);
 …
 USAGE:   changeord(newr,o[,r]);  newr,o=strings, r=ring/qring
 CREATE:  create a new ring with name `newr` and make it the basering if r is
          an existing ring/qring [default: r=basering].
+         an existing ring/qring [default: r=basering].
          The new ring differs from the old ring only in the ordering. If, say,
          (newr,o) = ("R","wp(2,3),dp") and the ring r exists and has >=3
+         (newr,o) = ("R","wp(2,3),dp") and the ring r exists and has >=3
          variables, the new basering will have name R and ordering wp(2,3),dp.
 RETURN:  No return value
 NOTE:    This proc uses 'execute' or calls a procedure using 'execute'.
          If you use it in your own proc, let the local names of your proc
+NOTE:    This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
          start with @ (see the file HelpForProc)
 EXAMPLE: example changeord; shows an example
 …
    setring @r;
    ideal i = ideal(@r); int @q = size(i);
    if( @q!=0 )
+   if( @q!=0 )
       { string @s = "@newr1"; }
    else
+   else
       { string @s = newr; }
    string @newring = @s+"=("+charstr(@r)+"),("+varstr(@r)+"),("+o+");";
 …
       ideal i = phi(i);
       attrib(i,"isSB",1);         //*** attrib funktioniert ?
       execute("qring "+newr+"=i;");
+      execute("qring "+newr+"=i;");
+   }
    export(`newr`);
 …
 USAGE:   changevar(newr,vars[,r]);  newr,vars=strings, r=ring/qring
 CREATE:  creates a new ring with name `newr` and makes it the basering if r
          is an existing ring/qring [default: r=basering].
+         is an existing ring/qring [default: r=basering].
          The new ring differs from the old ring only in the variables. If,
          say, (newr,vars) = ("R","t()") and the ring r exists and has n
          variables, the new basering will have name R and variables
+         say, (newr,vars) = ("R","t()") and the ring r exists and has n
+         variables, the new basering will have name R and variables
          t(1),...,t(n).
          If vars = "a,b,c,d", the new ring will have the variables a,b,c,d.
+         If vars = "a,b,c,d", the new ring will have the variables a,b,c,d.
 RETURN:  No return value
 NOTE:    This procedure is useful in connection with the procedure ringtensor,
          when a conflict between variable names must be avoided.
          This proc uses 'execute' or calls a procedure using 'execute'.
          If you use it in your own proc, let the local names of your proc
+         This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
          start with @ (see the file HelpForProc)
 EXAMPLE: example changevar; shows an example
 …
    setring @r;
    ideal i = ideal(@r); int @q = size(i);
    if( @q!=0 )
+   if( @q!=0 )
       { string @s = "@newr1"; }
    else
+   else
       { string @s = newr; }
    string @newring = @s+"=("+charstr(@r)+"),(";
    if( vars[size(vars)-1]=="(" and vars[size(vars)]==")" )
+   {
       @newring = @newring+vars[1,size(vars)-2]+"(1.."+string(nvars(@r))+")";
+   {
+      @newring = @newring+vars[1,size(vars)-2]+"(1.."+string(nvars(@r))+")";
+   }
    else { @newring = @newring+vars; }
 …
       ideal i = phi(i);
       attrib(i,"isSB",1);         //*** attrib funktioniert ?
       execute("qring "+newr+"=i;");
+      execute("qring "+newr+"=i;");
+   }
    export(`newr`);
 …
 CREATE:  Define a ring with name 's1', characteristic 's2', ordering 's4' and
          n variables with names derived from s3 and make it the basering.
          If s3 is a single letter, say s3="a", and if n<=26 then a and the
          following n-1 letters from the alphabeth (cyclic order) are taken as
          variables. If n>26 or if s3 is a single letter followed by (, say
          s3="T(", the variables are T(1),...,T(n).
+         If s3 is a single letter, say s3="a", and if n<=26 then a and the
+         following n-1 letters from the alphabeth (cyclic order) are taken as
+         variables. If n>26 or if s3 is a single letter followed by (, say
+         s3="T(", the variables are T(1),...,T(n).
 RETURN:  No return value
 NOTE:    This proc is useful for defining a ring in a procedure.
          This proc uses 'execute' or calls a procedure using 'execute'.
          If you use it in your own proc, let the local names of your proc
+         This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
          start with @ (see the file HelpForProc)
 EXAMPLE: example defring; shows an example
 …
    execute("keepring P"+string(n)+";");
    //the next comment is only shown if defringp is not called by another proc
    if (voice==2) { "// basering is now:",s; }
+}
+   if (voice==2) { "// basering is now:",s; }
+}
 example
 { "EXAMPLE:"; echo = 2;
 …
 proc extendring (string na, int n, string va, string o, list #)
 USAGE:   extendring(na,n,va,o[iv,i,r]);  na,va,o=strings,
+USAGE:   extendring(na,n,va,o[iv,i,r]);  na,va,o=strings,
          n,i=integers, r=ring, iv=intvec of positive integers or iv=0
 CREATE:  Define a ring with name `na` which extends the ring r by adding n new
 …
          the basering [default: (i,r)=(0,basering)]
          -- The characteristic is the characteristic of r
          -- The new vars are derived from va. If va is a single letter, say
             va="T", and if n<=26 then T and the following n-1 letters from
             T..Z..T (resp. T(1..n) if n>26) are taken as additional variables.
             If va is a single letter followed by (, say va="x(", the new
+         -- The new vars are derived from va. If va is a single letter, say
+            va="T", and if n<=26 then T and the following n-1 letters from
+            T..Z..T (resp. T(1..n) if n>26) are taken as additional variables.
+            If va is a single letter followed by (, say va="x(", the new
             variables are x(1),...,x(n)
          -- The ordering is the product ordering between the ordering of r and
             an ordering derived from `o` [and iv]
          -  If o contains a 'c' or a 'C' in front resp. at the end this is
+         -  If o contains a 'c' or a 'C' in front resp. at the end this is
             taken for the whole ordering in front resp. at the end. If o does
             not contain a 'c' or a 'C' the same rule applies to ordstr(r).
          -  If no intvec iv is given, or if iv=0, o may be any allowed ordstr,
             like "ds" or "dp(2),wp(1,2,3),Ds(2)" or "ds(a),dp(b),ls" if a and b
             are globally (!) defined integers and if a+b+1<=n
             If, however, a and b are local to a proc calling extendring, the
             intvec iv must be used to let extendring know the values of a and b
+            are globally (!) defined integers and if a+b+1<=n
+            If, however, a and b are local to a proc calling extendring, the
+            intvec iv must be used to let extendring know the values of a and b
          -  If an intvec iv !=0 is given, iv[1],iv[2],... is taken for the 1st,
 nd,... block of o, if o contains no substring "w" or "W" i.e. no
             weighted ordering (in the above case o="ds,dp,ls" and iv=a,b).
             If o contains a weighted ordering (only one (!) weighted block is
             allowed) iv[1] is taken as size for the weight-vector, the next
+nd,... block of o, if o contains no substring "w" or "W" i.e. no
+            weighted ordering (in the above case o="ds,dp,ls" and iv=a,b).
+            If o contains a weighted ordering (only one (!) weighted block is
+            allowed) iv[1] is taken as size for the weight-vector, the next
             iv[1] values of iv are taken as weights and the remaining values of
             iv as block-size for the remaining non-weighted blocks.
             e.g. o="dp,ws,Dp,ds", iv=3,2,3,4,2,5 creates the ordering
+            iv as block-size for the remaining non-weighted blocks.
+            e.g. o="dp,ws,Dp,ds", iv=3,2,3,4,2,5 creates the ordering
             dp(2),ws(2,3,4),Dp(5),ds
 RETURN:  No return value
 NOTE:    This proc is useful for adding deformation parameters.
          This proc uses 'execute' or calls a procedure using 'execute'.
          If you use it in your own proc, let the local names of your proc
+         This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
          start with @ (see the file HelpForProc)
 EXAMPLE: example extendring; shows an example
+{
 //--------------- initialization and place c/C of ordering properly -----------
    string @o1,@o2,@ro,@wstr,@v,@newring;
+   string @o1,@o2,@ro,@wstr,@v,@newring;
    int @i,@w,@ii,@k;
    intvec @iv,@iw;
 …
+   }
    @ro=ordstr(@r);
    if( @ro[1]=="c" or @ro[1]=="C" )
+   if( @ro[1]=="c" or @ro[1]=="C" )
       { @v=@ro[1,2]; @ro=@ro[3..size(@ro)]; }
    else
+   else
       { @wstr=@ro[size(@ro)-1,2]; @ro=@ro[1..size(@ro)-2]; }
    if( @k==0) { @o1=@v; @o2=@wstr; }
 …
       @k=n;                               //@k counts no of vars not yet ordered
       @w=find(o,"w")+find(o,"W");o=o+" ";
       if( @w!=0 )
+      {
          @wstr=o[@w..@w+1];
+      if( @w!=0 )
+      {
+         @wstr=o[@w..@w+1];
          o=o[1,@w-1]+"@"+o[@w+2,size(o)];
          @iw=@iv[2..@iv[1]+1];
          @wstr=@wstr+"("+string(@iw)+")";
          @k=@k-@iv[1];
+         @k=@k-@iv[1];
          @iv=@iv[@iv[1]+2..size(@iv)];
          @w=0;
+      }
       for( @ii=1; @ii<=size(@iv); @ii++ )
+      for( @ii=1; @ii<=size(@iv); @ii=@ii+1 )
+      {
          if( find(o,",",@w+1)!=0 )
 …
             @w=find(o,",",@w+1);
             if( o[@w-1]!="@" )
+            {
                o=o[1,@w-1]+"("+string(@iv[@ii])+")"+o[@w,size(o)];
+            {
+               o=o[1,@w-1]+"("+string(@iv[@ii])+")"+o[@w,size(o)];
                @w=find(o,",",@w+1);
                @k=@k-@iv[@ii];
 …
    if( n>26 or va[2]=="(" ) { @v = va[1]+"(1.."+string(n)+")"; }
    else                     { @v = A_Z(va,n); }
    if( @i==0 )
+   {
       @v=@v+","+varstr(@r);
       o=@o1+o+","+@ro+@o2;
+   }
    else
+   {
       @v=varstr(@r)+","+@v;
       o=@o1+@ro+","+o+@o2;
+   if( @i==0 )
+   {
+      @v=@v+","+varstr(@r);
+      o=@o1+o+","+@ro+@o2;
+   }
+   else
+   {
+      @v=varstr(@r)+","+@v;
+      o=@o1+@ro+","+o+@o2;
+   }
    @newring=@newring+@v+"),("+o+");";
 …
 { "EXAMPLE:"; echo = 2;
    ring r=0,(x,y,z),ds;
    show(r);"";
+   show(r);"";
    //no intvec given, no blocksize given: blocksize is derived from no of vars
    int t=5;
+   int t=5;
    extendring("R1",t,"a","dp");         //t global: "dp" -> "dp(5)"
    show(R1); "";
 …
    //the remaining variables, if intvec has too many components, the last
    //ones are ignored
    intvec v=3,2,3,4,1,3;
    extendring("R4",10,"A","ds,ws,Dp,dp",v,0,r);
          //v covers 3 blocks: v[1] (=3) : no of components of ws
+   intvec v=3,2,3,4,1,3;
+   extendring("R4",10,"A","ds,ws,Dp,dp",v,0,r);
+         //v covers 3 blocks: v[1] (=3) : no of components of ws
          //next v[1] values (=v[2..4]) give weights
          //remaining components of v are used for the remaining blocks
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc fetchall (R, list #)
+proc fetchall (R, list #)
 USAGE:   fetchall(R[,s]);  R=ring/qring, s=string
 CREATE:  fetch all objects of ring R (of type poly/ideal/vector/module/number/
 …
          If no 3rd argument is present, the names are the same as in R. If,
          say, f is a poly in R and the 3rd argument is the string "R", then f
          is maped to f_R etc.
+         is maped to f_R etc.
 RETURN:  no return value
 NOTE:    As fetch, this procedure maps the 1st, 2nd, ... variable of R to the
 st, 2nd, ... variable of the basering.
+st, 2nd, ... variable of the basering.
          The 3rd argument is useful in order to avoid conflicts of names, the
          empty string is allowed
 CAUTION: fetchall does not work inside a procedure
          //***at the moment it does not work if R contains a map
+         //***at the moment it does not work if R contains a map
 EXAMPLE: example fetchall; shows an example
+{
+{
    list @L@=names(R);
    int @ii@; string @s@;
    if( size(#) > 0 ) { @s@=@s@+"_"+#[1]; }
    for( @ii@=size(@L@); @ii@>0; @ii@-- )
+   {
+   for( @ii@=size(@L@); @ii@>0; @ii@=@ii@-1 )
+   {
       execute("def "+@L@[@ii@]+@s@+"=fetch(R,`@L@[@ii@]`);");
       execute("export "+@L@[@ii@]+@s@+";");
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc imapall (R, list #)
+proc imapall (R, list #)
 USAGE:   imapall(R[,s]);  R=ring/qring, s=string
 CREATE:  map all objects of ring R (of type poly/ideal/vector/module/number/
 …
          If no 3rd argument is present, the names are the same as in R. If,
          say, f is a poly in R and the 3rd argument is the string "R", then f
          is maped to f_R etc.
+         is maped to f_R etc.
 RETURN:  no return value
 NOTE:    As imap, this procedure maps the variables of R to the variables with
 …
          empty string is allowed
 CAUTION: imapall does not work inside a procedure
          //***at the moment it does not work if R contains a map
+         //***at the moment it does not work if R contains a map
 EXAMPLE: example imapall; shows an example
+{
+{
    list @L@=names(R);
    int @ii@; string @s@;
    if( size(#) > 0 ) { @s@=@s@+"_"+#[1]; }
    for( @ii@=size(@L@); @ii@>0; @ii@-- )
+   {
+   for( @ii@=size(@L@); @ii@>0; @ii@=@ii@-1 )
+   {
          execute("def "+@L@[@ii@]+@s@+"=imap(R,`@L@[@ii@]`);");
          execute("export "+@L@[@ii@]+@s@+";");
 …
 ///////////////////////////////////////////////////////////////////////////////
 proc mapall (R, ideal i, list #)
+proc mapall (R, ideal i, list #)
 USAGE:   mapall(R,i[,s]);  R=ring/qring, i=ideal of basering, s=string
 CREATE:  map all objects of ring R (of type poly/ideal/vector/module/number/
+CREATE:  map all objects of ring R (of type poly/ideal/vector/module/number/
          matrix, map) into the basering, by mapping the jth variable of R to
          the jth generator of the ideal i. If no 3rd argument is present, the
          names are the same as in R. If, say, f is a poly in R and the 3rd
          argument is the string "R", then f is maped to f_R etc.
+         names are the same as in R. If, say, f is a poly in R and the 3rd
+         argument is the string "R", then f is maped to f_R etc.
 RETURN:  no return value
 NOTE:    This procedure has the same effect as defining a map, say psi, by
          map psi=R,i; and then applying psi to all objects of R. In particular,
          maps from R to some ring S are composed with psi, creating thus a map
          from the basering to S.
+         from the basering to S.
          mapall may be combined with copyring to change vars for all objects.
          The 3rd argument is useful in order to avoid conflicts of names, the
          empty string is allowed
 CAUTION: mapall does not work inside a procedure
+CAUTION: mapall does not work inside a procedure
 EXAMPLE: example mapall; shows an example
+{
+{
    list @L@=names(R); map @psi@ = R,i;
    int @ii@; string @s@;
    if( size(#) > 0 ) { @s@=@s@+"_"+#[1]; }
    for( @ii@=size(@L@); @ii@>0; @ii@-- )
+   {
+   for( @ii@=size(@L@); @ii@>0; @ii@=@ii@-1 )
+   {
       execute("def "+@L@[@ii@]+@s@+"=@psi@(`@L@[@ii@]`);");
       execute("export "+@L@[@ii@]+@s@+";");
 …
 "   ideal j=x,y,z;";
 "   matrix M[2][3]=1,2,3,x,y,z;";
 "   map phi=R,x2,y2,z2; ";
+"   map phi=R,x2,y2,z2; ";
 "   ring S=0,(a,b,c),ds;";
 "   ideal i=c,a,b;";
 …
 USAGE:   ringtensor(s,r1,r2,...); s=string, r1,r2,...=rings
 CREATE:  A new base ring with name `s` if r1,r2,... are existing rings.
          If, say, s = "R" and the rings r1,r2,... exist, the new ring will
          have name R, variables from all rings r1,r2,... and as monomial
+         If, say, s = "R" and the rings r1,r2,... exist, the new ring will
+         have name R, variables from all rings r1,r2,... and as monomial
          ordering the block (product) ordering of r1,r2,... . Hence, R
          is the tensor product of the rings r1,r2,... with ordering matrix
 …
          variable with name x in r1, there is no access to x in r2).
          The procedure works also for quotient rings ri, if the characteristic
          of ri is compatible with the characteristic of r1 (i.e. if imap from
+         of ri is compatible with the characteristic of r1 (i.e. if imap from
          ri to r1 is implemented)
          This proc uses 'execute' or calls a procedure using 'execute'.
          If you use it in your own proc, let the local names of your proc
+         This proc uses 'execute' or calls a procedure using 'execute'.
+         If you use it in your own proc, let the local names of your proc
          start with @ (see the file HelpForProc)
 EXAMPLE: example ringtensor; shows an example
 …
    string @vars,@order,@oi,@s1;
 //---- gather variables, orderings and ideals (of qrings) from given rings ----
    for(@ii=1; @ii<=@n; @ii++)
+   {
       if( ordstr(#[@ii])[1]=="C" or ordstr(#[@ii])[1]=="c" )
+   for(@ii=1; @ii<=@n; @ii=@ii+1)
+   {
+      if( ordstr(#[@ii])[1]=="C" or ordstr(#[@ii])[1]=="c" )
            { @oi=ordstr(#[@ii])[3,size(ordstr(#[@ii]))-2]; }
       else { @oi=ordstr(#[@ii])[1,size(ordstr(#[@ii]))-2]; }
 …
       def @r(@ii)=#[@ii];
       setring @r(@ii);
       ideal i(@ii)=ideal(@r(@ii));
+      ideal i(@ii)=ideal(@r(@ii));
       int @q(@ii)=size(i(@ii));
       @q=@q+@q(@ii);
 …
+   {
       ideal i;
       for(@ii=1; @ii<=@n; @ii++)
+      for(@ii=1; @ii<=@n; @ii=@ii+1)
+      {
          if( @q(@ii)!=0 )
 …
+      }
       i=std(i);
       execute("qring "+s+"=i;");
+      execute("qring "+s+"=i;");
+   }
 //----------------------- export and keep created ring ------------------------
 …
    return();
+}
 example
+example
 { "EXAMPLE:"; echo = 2;
    ring r=32003,(x,y,u,v),dp;
    ring s=0,(a,b,c),wp(1,2,3);
    ring t=0,x(1..5),(c,ls);
    ringtensor("R",r,s,t);
+   ringtensor("R",r,s,t);
    type R;"";
    setring s;
 …
    changevar("T","d,e,f,g,h",t); //set T (change vars of t to d..h) the basering
    qring qT=std(d2+e2-f3);       //create qring of T mod d2+e2-f3
    ringtensor("Q",s,qS,t,qT);
+   ringtensor("Q",s,qS,t,qT);
    type Q;
    kill R,Q,S,T;

Singular/LIB/sing.lib

-                      r6d09c56
+                      r6f2edc
 // $Id: sing.lib,v 1.1.1.1 1997-04-25 15:13:27 obachman Exp $
+// $Id: sing.lib,v 1.2 1997-04-28 19:27:25 obachman Exp $
 //system("random",787422842);
 //(GMG+BM)
 ///////////////////////////////////////////////////////////////////////////////
 LIBRARY:  sing.lib      PROCEDURES FOR SINGULARITIES
+//(GMG/BM, last modified 26.06.96)
+///////////////////////////////////////////////////////////////////////////////
+LIBRARY:  sing.lib      PROCEDURES FOR SINGULARITIES
  deform(i);             infinitesimal deformations of ideal i
 …
  T12(i);                T1- and T2-module of ideal i
 LIB "inout.lib";
+LIB "inout.lib";
 LIB "random.lib";
 ///////////////////////////////////////////////////////////////////////////////
 proc deform (ideal id)
 USAGE:   deform(id);  id = ideal or poly
+USAGE:   deform(id); id=ideal or poly
 RETURN:  matrix, columns are kbase of infinitesimal deformations
 EXAMPLE: example deform; shows an example
+{
+EXAMPLE: example deform; shows an example
+{
    list L=T1(id,"");
+   return(L[2]*kbase(std(L[1])));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=32003,(x,y,z),ds;
+   ideal i=xy,xz,yz;
+   matrix T=deform(i);print(T);
+   print(deform(x3+y5+z2));
+   def K=L[1]; attrib(K,"isSB",1);
+   return(L[2]*kbase(K));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r   = 32003,(x,y,z),ds;
+   ideal i  = xy,xz,yz;
+   matrix T = deform(i);
+   print(T);
+   print(deform(x3+y5+z2));
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 example
 { "EXAMPLE:"; echo = 2;
    ring r=32003,(x,y,z),ds;
    ideal i= x5+y6+z6,x2+2y2+3z2;
+   ring r  = 32003,(x,y,z),ds;
+   ideal i = x5+y6+z6,x2+2y2+3z2;
    dim_slocus(i);
+}
 …
 proc is_active (poly f, id)
 USAGE:   is_active(f,id); f poly, id ideal or module
+USAGE:   is_active(f,id); f poly, id ideal or module
 RETURN:  1 if f is an active element modulo id (i.e. dim(id)=dim(id+f*R^n)+1,
          if id is a submodule of R^n) resp. 0 if f is not active.
          The basering may be a quotient ring
+         The basering may be a quotient ring
 NOTE:    regular parameters are active but not vice versa (id may have embedded
          components). proc is_reg tests whether f is a regular parameter
 EXAMPLE: example is_active; shows an example
+{
    if( size(id)==0 ) { return(1); }
+   if( size(id)==0 ) { return(1); }
    if( typeof(id)=="ideal" ) { ideal m=f; }
    if( typeof(id)=="module" ) { module m=f*freemodule(rank(id)); }
+   if( typeof(id)=="module" ) { module m=f*freemodule(nrows(id)); }
    return(dim(std(id))-dim(std(id+m)));
+}
 example
 { "EXAMPLE:"; echo = 2;
    ring r=32003,(x,y,z),ds;
    ideal i=yx3+y,yz3+y3z;
    poly f=x;
+   ring r   =32003,(x,y,z),ds;
+   ideal i  = yx3+y,yz3+y3z;
+   poly f   = x;
    is_active(f,i);
    qring q = std(x4y5);
    poly f=x;
    module m=[yx3+x,yx3+y3x];
+   qring q  = std(x4y5);
+   poly f   = x;
+   module m = [yx3+x,yx3+y3x];
    is_active(f,m);
+}
 …
 USAGE:   is_ci(i); i ideal
 RETURN:  intvec = sequence of dimensions of ideals (j[1],...,j[k]), for
+         k=1,...,size(j), where j is minimal base of i. i is a complete
+         intersection if last number equals nvars-size(i)
+NOTE:    dim(0-ideal) = -1. You may first apply simplify(i,10); in order to
+         delete zeroes and multiples from set of generators
+         k=1,...,size(j), where j is minimal base of i. i is a complete
+         intersection if last number equals nvars-size(i)
+NOTE:    dim(0-ideal) = -1. You may first apply simplify(i,10); in order to
+         delete zeroes and multiples from set of generators
+         printlevel >=0: display comments (default)
 EXAMPLE: example is_ci; shows an example
+{
 …
    i=minbase(i);
    int s = ncols(i);
+   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
 //--------------------------- compute dimensions ------------------------------
    for( n=1; n<=s; n++ )
+   {
+   for( n=1; n<=s; n=n+1 )
+   {
       id = i[1..n];
       dimvec[n] = dim(std(id));
+   }
    n = dimvec[s];
+//--------------------------- output ------------------------------------------
+   if( defined(printlevel) )
+   {
+      if( n+s !=nvars(basering) )
+      { dbprint(printlevel,"// no complete intersection"); }
+      if( n+s ==nvars(basering) )
+      { dbprint(printlevel,"// complete intersection of dim "+string(n)); }
+      dbprint(printlevel,"// dim-sequence:");
+   }
+   if( voice==2 )
+   {
+      if( n+s !=nvars(basering)) {"// no complete intersection"; }
+      if( n+s ==nvars(basering)) {"// complete intersection of dim",n;}
+      "// dim-sequence:";
+   }
+//--------------------------- output ------------------------------------------
+   if( n+s != nvars(basering) )
+   { dbprint(p,"// no complete intersection"); }
+   if( n+s == nvars(basering) )
+   { dbprint(p,"// complete intersection of dim "+string(n)); }
+   dbprint(p,"// dim-sequence:");
    return(dimvec);
+}
 example
 { "EXAMPLE:";   echo = 2;
    int printlevel=2;                // this forces the proc to display comments
    export printlevel;
    ring r=32003,(x,y,z),ds;
    ideal i=x4+y5+z6,xyz,yx2+xz2+zy7;
+{ "EXAMPLE:"; echo = 2;
+   int p      = printlevel;
+   printlevel = 1;                // display comments
+   ring r     = 32003,(x,y,z),ds;
+   ideal i    = x4+y5+z6,xyz,yx2+xz2+zy7;
    is_ci(i);
    i=xy,yz;
    is_ci(i);
    kill printlevel;
+   i          = xy,yz;
+   is_ci(i);
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
          generated by id[1]..id[i], k = 1..size(id); dim(0-ideal) = -1;
          id defines an isolated singularity if last number is 0
+NOTE:    printlevel >=0: display comments (default)
 EXAMPLE: example is_is; shows an example
+{
   int l; intvec dims; ideal j;
+  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
 //--------------------------- compute dimensions ------------------------------
    for( l=1; l<=ncols(i); l++ )
+   {
      j = i[1..l];
+   for( l=1; l<=ncols(i); l=l+1 )
+   {
+     j = i[1..l];
      dims[l] = dim(std(slocus(j)));
+   }
    if( voice==2 ) {"// dim of singular locus =",dims[size(dims)],"dim-sequence:";
                    "// isolated singularity if last number is 0"; }
+   dbprint(p,"// dim of singular locus = "+string(dims[size(dims)]),
+             "// isolated singularity if last number is 0 in dim-sequence:");
    return(dims);
+}
 example
 { "EXAMPLE:"; echo = 2;
+   ring r=32003,(x,y,z),ds;
+   ideal i=x2y,x4+y5+z6,yx2+xz2+zy7;
+   int p      = printlevel;
+   printlevel = 1;
+   ring r     = 32003,(x,y,z),ds;
+   ideal i    = x2y,x4+y5+z6,yx2+xz2+zy7;
    is_is(i);
+// isolated singularity if last number is 0
+   poly f=xy+yz;
+   poly f     = xy+yz;
    is_is(f);
+// isolated singularity if last number is 0
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
    id=std(id);
    d=size(q);
    for( ii=1; ii<=d; ii++ )
+   for( ii=1; ii<=d; ii=ii+1 )
+   {
       if( reduce(q[ii],id)!=0 )
 …
    return(1);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=32003,(x,y),ds;
+   ideal i=x8,y8;ideal j=(x+y)^4;
+   i=intersect(i,j); poly f=xy;
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r  = 32003,(x,y),ds;
+   ideal i = x8,y8;
+   ideal j = (x+y)^4;
+   i       = intersect(i,j);
+   poly f  = xy;
    is_reg(f,i);
+}
 …
 NOTE:    let R be the basering and id a submodule of R^n. The procedure checks
          injectivity of multiplication with i[k] on R^n/id+i[1..k-1].
+         The basering may be a quotient ring
+         The basering may be a quotient ring
+         printlevel >=0: display comments (default)
+         printlevel >=1: display comments during computation
 EXAMPLE: example is_regs; shows an example
+{
+   int d,ii,r;
+   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
    if( size(#)==0 ) { ideal id; }
    else { def id=#[1]; }
    if( size(i)==0 ) { return(0); }
+   int d,ii,r;
+   d=size(i);
+   d=size(i);
    if( typeof(id)=="ideal" ) { ideal m=1; }
    if( typeof(id)=="module" ) { module m=freemodule(rank(id)); }
    for( ii=1; ii<=d; ii++ )
+   {
       if( voice==2 )
+   if( typeof(id)=="module" ) { module m=freemodule(nrows(id)); }
+   for( ii=1; ii<=d; ii=ii+1 )
+   {
+      if( p>=2 )
       { "// checking whether element",ii,"is regular mod 1 ..",ii-1; }
+      if( is_reg(i[ii],id)==0 )
+      {
+         if( voice==2 )
+         {
+            "// elements 1 ..",ii-1,"are regular,",
+            ii,"is not regular mod 1 ..",ii-1;
+         }
+         return(0);
+      if( is_reg(i[ii],id)==0 )
+      {
+        dbprint(p,"// elements 1.."+string(ii-1)+" are regular, " +
+                string(ii)+" is not regular mod 1.."+string(ii-1));
+         return(0);
+      }
       id=id+i[ii]*m;
+   }
    if( voice==2 ) { "// elements are a regular sequence of length",d; }
+      id=id+i[ii]*m;
+   }
+   if( p>=1 ) { "// elements are a regular sequence of length",d; }
    return(1);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r1=32003,(x,y,z),ds;
+   ideal i=x8,y8,(x+y)^4;
+example
+{ "EXAMPLE:"; echo = 2;
+   int p      = printlevel;
+   printlevel = 1;
+   ring r1    = 32003,(x,y,z),ds;
+   ideal i    = x8,y8,(x+y)^4;
    is_regs(i);
    module m=[x,0,y];
    i=x8,(x+z)^4;;
+   module m   = [x,0,y];
+   i          = x8,(x+z)^4;;
    is_regs(i,m);
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 USAGE:   milnor(i); i ideal or poly
 RETURN:  Milnor number of i, if i is ICIS (isolated complete intersection
          singularity) in generic form, resp. -1 if not
+         singularity) in generic form, resp. -1 if not
 NOTE:    use proc nf_icis to put generators in generic form
+         printlevel >=0: display comments (default)
 EXAMPLE: example milnor; shows an example
+{
   i = simplify(i,10);     //delete zeroes and multiples from set of generators
+{
+  i = simplify(i,10);     //delete zeroes and multiples from set of generators
   int n = size(i);
   int l,q,m_nr;  ideal t;  intvec disc;
+//---------------------------- hypersurface case ------------------------------
+  if( n==1 )
+  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
+//---------------------------- hypersurface case ------------------------------
+  if( n==1 )
+  {
      i = std(jacob(i[1]));
      m_nr = vdim(i);
      if( m_nr<0 and voice==2 ) { "// no isolated singularity"; }
      return(m_nr);
+     m_nr = vdim(i);
+     if( m_nr<0 and p>=1 ) { "// no isolated singularity"; }
+     return(m_nr);
+  }
 //------------ isolated complete intersection singularity (ICIS) --------------
   for( l=n; l>0; l=l-1)
+  {
+      t      = minor(jacob(i),l);
+      i[l]   = 0;
+  {   t      = minor(jacob(i),l);
+      i[l]   = 0;
       q      = vdim(std(i+t));
       disc[l]= q;
       if( q ==-1 )
+      {
+         if( voice==2 )
+      {  if( p>=1 )
             {  "// not in generic form or no ICIS; use proc nf_icis to put";
             "// generators in generic form and then try milnor again!";  }
+            "// generators in generic form and then try milnor again!";  }
          return(q);
+      }
       m_nr = q-m_nr;
+      m_nr = q-m_nr;
+  }
 //---------------------------- change sign ------------------------------------
   if (m_nr < 0) { m_nr=-m_nr; }
   if( voice==2 ) { "//sequence of discriminant numbers:",disc; }
+//---------------------------- change sign ------------------------------------
+  if (m_nr < 0) { m_nr=-m_nr; }
+  if( p>=1 ) { "//sequence of discriminant numbers:",disc; }
   return(m_nr);
+}
 example
 { "EXAMPLE:"; echo = 2;
+   ring r=32003,(x,y,z),ds;
+   ideal j=x5+y6+z6,x2+2y2+3z2,xyz+yx;
+   int p      = printlevel;
+   printlevel = 1;
+   ring r     = 32003,(x,y,z),ds;
+   ideal j    = x5+y6+z6,x2+2y2+3z2,xyz+yx;
    milnor(j);
+   poly f=x7+y7+(x-y)^2*x2y2+z2;
+   milnor(f);
+   poly f     = x7+y7+(x-y)^2*x2y2+z2;
+   milnor(f);
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
          (isolated complete intersection singularity), return i if not
 NOTE:    this proc is useful in connection with proc milnor
+         printlevel >=0: display comments (default)
 EXAMPLE: example nf_icis; shows an example
+{
    i = simplify(i,10);  //delete zeroes and multiples from set of generators
    int p,b = 100,0;
+   int p,b = 100,0;
    int n = size(i);
    matrix mat=freemodule(n);
+//---------------------------- test: complete intersection? -------------------
+   int P = printlevel-voice+3;  // P=printlevel+1 (default: P=1)
+//---------------------------- test: complete intersection? -------------------
    intvec sl = is_ci(i);
    if( n+sl[n] != nvars(basering) )
+   {
     if( voice==2 ) { "// no complete intersection"; }
     return(i);
+   }
 //--------------- test: isolated singularity in generic form? -----------------
+   if( n+sl[n] != nvars(basering) )
+   {
+      dbprint(P,"// no complete intersection");
+      return(i);
+   }
+//--------------- test: isolated singularity in generic form? -----------------
    sl = is_is(i);
    if ( sl[n] != 0 )
+   {
       if( voice==2 ) { "// no isolated singularity"; }
+      dbprint(P,"// no isolated singularity");
       return(i);
+   }
 //------------ produce generic linear combinations of generators --------------
+//------------ produce generic linear combinations of generators --------------
    int prob;
    while ( sum(sl) != 0 )
+   while ( sum(sl) != 0 )
    {  prob=prob+1;
       p=p-25; b=b+10;
+      p=p-25; b=b+10;
       i = genericid(i,p,b);          // proc genericid from random.lib
       sl = is_is(i);
+   }
+   if( voice==2 ) { "// ICIS in generic form after",prob,"genericity loop(s)";}
+   return(i);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=32003,(x,y,z),ds;
+   ideal i=x3+y4,z4+yx;
+   nf_icis(i);
+   ideal j=x3+y4,xy,yz;
+   dbprint(P,"// ICIS in generic form after "+string(prob)+" genericity loop(s)");
+   return(i);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   int p      = printlevel;
+   printlevel = 1;
+   ring r     = 32003,(x,y,z),ds;
+   ideal i    = x3+y4,z4+yx;
+   nf_icis(i);
+   ideal j    = x3+y4,xy,yz;
    nf_icis(j);
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
   int cod  = nvars(basering)-dim(std(i));
   i        = i+minor(jacob(i),cod);
   return(i);
+}
 example
 { "EXAMPLE:"; echo = 2;
    ring r=32003,(x,y,z),ds;
    ideal i= x5+y6+z6,x2+2y2+3z2;
+  return(i);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r  = 32003,(x,y,z),ds;
+   ideal i = x5+y6+z6,x2+2y2+3z2;
    dim(std(slocus(i)));
+}
 …
 proc Tjurina (id, list #)
+USAGE:   Tjurina(id[,<any>]);  id ideal or poly (assume: id=ICIS)
+RETURN:  Tjurina(id): standard basis of Tjurina-module of id, displays Tjurina
+         number
+         Tjurina(id,...): If a second argument is present (of any type) return
+         a list of 4 objects: [1]=Tjurina number (int), [2]=basis of miniversal
+         deformation (module), [3]=SB of Tjurina module (module), [4]=Tjurina
+         module (module)
+NOTE:    if id is a poly the output will be of type ideal rather than module
+EXAMPLE: example Tjurina; shows an example
+USAGE:   Tjurina(id[,<any>]);  id=ideal or poly
+ASSUME:  id=ICIS (isolated complete intersection singularity)
+RETURN:  standard basis of Tjurina-module of id,
+         of type module if id=ideal, resp. of type ideal if id=poly.
+         If a second argument is present (of any type) return a list:
+           [1] = Tjurina number,
+           [2] = k-basis of miniversal deformation,
+           [3] = SB of Tjurina module,
+           [4] = Tjurina module
+DISPLAY: Tjurina number if printlevel >= 0 (default)
+NOTE:    Tjurina number = -1 implies that id is not an ICIS
+EXAMPLE: example Tjurina; shows examples
+{
 //---------------------------- initialisation ---------------------------------
   def i = simplify(id,10);
+  def i = simplify(id,10);
   int tau,n = 0,size(i);
   if( size(ideal(i))==1 ) { def m=i; }  // hypersurface case
 …
   def st1 = std(t1);                    // SB of Tjurina module/ideal
   tau = vdim(st1);                      // Tjurina number
+  def kB = kbase(st1);                  // basis of miniversal deformation
+  if( voice==2 )  { "// Tjurina number =",tau; }
+  if( size(#)>0 ) { return(tau,kB,st1,t1); }
+  dbprint(printlevel-voice+3,"// Tjurina number = "+string(tau));
+  if( size(#)>0 )
+  {
+     def kB = kbase(st1);               // basis of miniversal deformation
+     return(tau,kB,st1,t1);
+  }
   return(st1);
+}
 example
 { "EXAMPLE:"; echo = 2;
+   ring r=0,(x,y,z),ds;
+   poly f = x5+y6+z7+xyz;        // singularity T[5,6,7]
+   list T = Tjurina(f,"");
+   show(T[1]);                   // Tjurina number, should be 16
+   show(T[2]);                   // basis of miniversal deformation
+   show(T[3]);                   // SB of Tjurina ideal
+   show(T[4]); "";               // Tjurina ideal
+   ideal j=x2+y2+z2,x2+2y2+3z2;
+   show(kbase(Tjurina(j)));      // basis of miniversal deformation
+   hilb(Tjurina(j));             // Hilbert series of Tjurina module
+   int p      = printlevel;
+   printlevel = 1;
+   ring r     = 0,(x,y,z),ds;
+   poly f     = x5+y6+z7+xyz;        // singularity T[5,6,7]
+   list T     = Tjurina(f,"");
+   show(T[1]);                       // Tjurina number, should be 16
+   show(T[2]);                       // basis of miniversal deformation
+   show(T[3]);                       // SB of Tjurina ideal
+   show(T[4]); "";                   // Tjurina ideal
+   ideal j    = x2+y2+z2,x2+2y2+3z2;
+   show(kbase(Tjurina(j)));          // basis of miniversal deformation
+   hilb(Tjurina(j));                 // Hilbert series of Tjurina module
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc tjurina (ideal i)
+USAGE:   tjurina(id);  id ideal or poly (assume: id=ICIS)
+USAGE:   tjurina(id);  id=ideal or poly
+ASSUME:  id=ICIS (isolated complete intersection singularity)
 RETURN:  int = Tjurina number of id
+NOTE:    Tjurina number = -1 implies that id is not an ICIS
 EXAMPLE: example tjurina; shows an example
+{
    return(vdim(Tjurina(i)));
+   return(vdim(Tjurina(i)));
+}
 example
 …
    ring r=32003,(x,y,z),(c,ds);
    ideal j=x2+y2+z2,x2+2y2+3z2;
    tjurina(j);
+   tjurina(j);
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 proc T1 (ideal id, list #)
 USAGE:   T1(id[,<any>]);  id = ideal or poly
+RETURN:  T1(id): T1-module of id or T1-ideal if id is a poly. This is a
+RETURN:  T1(id): of type module/ideal if id is of type ideal/poly.
+         We call T1(id) the T1-module of id. It is a std basis of the
          presentation of 1st order deformations of P/id, if P is the basering.
          T1(id,...): If a second argument is present (of any type) return a
          list of 3 modules:
             [1]= presentation of infinitesimal deformations of id (=T1(id))
+         If a second argument is present (of any type) return a list of
+modules:
+            [1]= T1(id)
             [2]= generators of normal bundle of id, lifted to P
+            [3]= module of relations of [2], lifted to P ([2]*[3]=0 mod id)
+         The list contains all non-easy objects which must be computed anyway
+            [3]= module of relations of [2], lifted to P
+                 (note: transpose[3]*[2]=0 mod id)
+         The list contains all non-easy objects which must be computed
          to get T1(id).
+         The situation is described in detail in the procedure T1_expl
+         from library explain.lib
+NOTE:    T1(id) itself is usually of minor importance, nevertheless, from it
+         all relevant information can be obtained. Since no standard basis is
+         computed, the user has first to compute a standard basis before
+         applying vdim or hilb etc.. For example, matrix([2])*(kbase(std([1])))
+         represents a basis of 1st order semiuniversal deformation of id
+         (use proc 'deform', to get this in a direct and convenient way).
+         If the input is weighted homogeneous with weights w1,...,wn, use
+         ordering wp(w1..wn), even in the local case, which is equivalent but
+         faster than ws(w1..wn).
+DISPLAY: k-dimension of T1(id) if printlevel >= 0 (default)
+NOTE:    T1(id) itself is usually of minor importance. Nevertheless, from it
+         all relevant information can be obtained. The most important are
+         probably vdim(T1(id)); (which computes the Tjurina number),
+         hilb(T1(id)); and kbase(T1(id));
+         If T1 is called with two argument, then matrix([2])*(kbase([1]))
+         represents a basis of 1st order semiuniversal deformation of id
+         (use proc 'deform', to get this in a direct way).
          For a complete intersection the proc Tjurina is faster
 EXAMPLE: example T1; shows an example
 …
    ideal J=simplify(id,10);
 //--------------------------- hypersurface case -------------------------------
+  if( size(J)<2 )
+  {
+     ideal t1=J[1],jacob(J[1]); module nb=[1]; module pnb;
+  if( size(J)<2 )
+  {
+     ideal t1  = std(J+jacob(J[1]));
+     module nb = [1]; module pnb;
+     dbprint(printlevel-voice+3,"// dim T1 = "+string(vdim(t1)));
      if( size(#)>0 ) { return(t1*gen(1),nb,pnb); }
      return(t1);
 …
 //--------------------------- presentation of J -------------------------------
    int rk;
   def P=basering;
+  def P = basering;
    module jac, t1;
    jac=jacob(J);                 // jacobian matrix of J converted to module
    res(J,2,A);                   // compute presentation of J
    t1=transpose(A(2));           // transposed 1st syzygy module of J
+   jac  = jacob(J);                 // jacobian matrix of J converted to module
+   res(J,2,A);                      // compute presentation of J
+                                    // A(2) = 1st syzygy module of J
 //---------- go to quotient ring mod J and compute normal bundle --------------
   qring  R=std(J);
    module jac=fetch(P,jac);
    module t1=fetch(P,t1);
    res(t1,3,B);                  // resolve t1, B(2)=(J/J^2)*=normal_bdl
    t1=lift(B(2),jac)+B(3);       // pres. of normal_bdl/trivial_deformations
    rk=rank(t1);
+  qring  R    = std(J);
+   module jac = fetch(P,jac);
+   module t1  = transpose(fetch(P,A(2)));
+   res(t1,2,B);                     // resolve t1, B(2)=(J/J^2)*=normal_bdl
+   t1         = modulo(B(2),jac);   // pres. of normal_bdl/trivial_deformations
+   rk=nrows(t1);
 //-------------------------- pull back to basering ----------------------------
   setring P;
    t1 = fetch(R,t1)+J*freemodule(rk);  // T1-module, presentation of T1
+   if( size(#)>0 )
+   {
+      module B2 = fetch(R,B(2));        // (generators of) normal bundle
+      module B3 = fetch(R,B(3));        // presentation of normal bundle
+      return(t1,B2,B3);
+   }
+   return(t1);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=32003,(x,y,z),(c,ds);
+   ideal i=xy,xz,yz;
+   module T=T1(i);
+   vdim(std(T));                 // Tjurina number = dim_K(T1), should be 3
+   t1 = std(t1);
+   dbprint(printlevel-voice+3,"// dim T1 = "+string(vdim(t1)));
+   if( size(#)>0 )
+   {
+      module B2 = fetch(R,B(2));        // presentation of normal bundle
+      list L = t1,B2,A(2);
+      attrib(L[1],"isSB",1);
+      return(L);
+   }
+   return(t1);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   int p      = printlevel;
+   printlevel = 1;
+   ring r     = 32003,(x,y,z),(c,ds);
+   ideal i    = xy,xz,yz;
+   module T   = T1(i);
+   vdim(T);                      // Tjurina number = dim_K(T1), should be 3
    list L=T1(i,"");
    module kB = kbase(std(L[1]));
+   module kB  = kbase(L[1]);
    print(L[2]*kB);               // basis of 1st order miniversal deformation
    size(L[1]);                   // number of generators of T1-module
    show(L[2]);                   // (generators of) normal bundle
    print(L[3]);                  // relation matrix of normal bundle (mod i)
    print(L[2]*L[3]);             // should be 0 (mod i)
+   show(L[2]);                   // presentation of normal bundle
+   print(L[3]);                  // relations of i
+   print(transpose(L[3])*L[2]);  // should be 0 (mod i)
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 proc T2 (ideal id, list #)
 USAGE:   T2(id[,<any>]);  id = ideal
 RETURN:  T2(id): T2-module of id . This is a presentation of the module of
             obstructions of R=P/id, if P is the basering.
          T2(id,...): If a second argument is present (of any type) return a
          list of 6 modules and 1 ideal:
             [1]= presentation of module of obstructions (=T2(id))
+RETURN:  T2(id): T2-module of id . This is a std basis of a presentation of
+         the module of obstructions of R=P/id, if P is the basering.
+         If a second argument is present (of any type) return a list of
+modules and 1 ideal:
+            [1]= T2(id)
             [2]= standard basis of id (ideal)
             [3]= module of relations of id (=1st syzygy module of id)
+            [4]= presentation of [3] (=2nd syzygy module of id)
+            [5]= lifting of Koszul relations kos, kos=module([3]*matrix([5]))
+            [6]= generators of Hom_P([3]/kos,R), lifted to P
+            [7]= relations of Hom_P([3]/kos,R), lifted to P
+         The list contains all non-easy objects which must be computed anyway
+         to get T2(id). The situation is described in detail in the procedure
+         T2_expl from library explain.lib
+NOTE:    Since no standard basis is computed, the user has first to compute a
+         standard basis before applying vdim or hilb etc..
+         If the input is weighted homogeneous with weights w1,...,wn, use
+         ordering wp(w1..wn), even in the local case, which is equivalent but
+         faster than ws(w1..wn).
+         Use proc miniversal to get equations of miniversal deformation;
+            [4]= presentation of syz/kos
+            [5]= relations of Hom_P([3]/kos,R), lifted to P
+         The list contains all non-easy objects which must be computed
+         to get T2(id).
+DISPLAY: k-dimension of T2(id) if printlevel >= 0 (default)
+NOTE:    The most important information is probably vdim(T2(id)).
+         Use proc miniversal to get equations of miniversal deformation.
 EXAMPLE: example T2; shows an example
+{
 //--------------------------- initialisation ----------------------------------
   def P = basering;
+   ideal J = simplify(id,10);
+   module kos,L0,t2;
+   ideal J = id;
+   module kos,SK,B2,t2;
+   list L;
    int n,rk;
 //------------------- presentation of non-trivial syzygies --------------------
    res(J,3,A);                            // resolve J, A(2)=syz
+//------------------- presentation of non-trivial syzygies --------------------
+   res(J,2,A);                            // resolve J, A(2)=syz
    kos  = koszul(2,J);                    // module of Koszul relations
+   L0 = lift(A(2),kos);                   // lift Koszul relations to syz
+   t2   = L0+A(3);                        // presentation of syz/kos
+   SK   = modulo(A(2),kos);               // presentation of syz/kos
    ideal J0 = std(J);                     // standard basis of J
 //*** sollte bei der Berechnung von res mit anfallen, zu aendern!!
+//?*** sollte bei der Berechnung von res mit anfallen, zu aendern!!
 //---------------------- fetch to quotient ring mod J -------------------------
   qring R = J0;                           // make P/J the basering
    module A2' = transpose(fetch(P,A(2))); // dual of syz
    module t2  = transpose(fetch(P,t2));   // dual of syz/kos
    res(t2,3,B);                           // resolve t2
    t2 = lift(B(2),A2')+B(3);              // presentation of T2
    rk = rank(t2);
+   module A2' = transpose(fetch(P,A(2))); // dual of syz
+   module t2  = transpose(fetch(P,SK));   // dual of syz/kos
+   res(t2,2,B);                           // resolve (syz/kos)*
+   t2 = modulo(B(2),A2');                 // presentation of T2
+   rk = nrows(t2);
 //---------------------  fetch back to basering -------------------------------
   setring P;
    t2 = fetch(R,t2)+J*freemodule(rk);
+   if( size(#)>0 )
+   {
+      module B2 = fetch(R,B(2));        // generators of Hom_P(syz/kos,R)
+      module B3 = fetch(R,B(3));        // relations of Hom_P(syz/kos,R)
+      return(t2,J0,A(2),A(3),L0,B2,B3);
+   }
+   return(t2);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring  r = 32003,(x,y),(c,dp);
+   ideal j = x6-y4,x6y6,x2y4-x5y2;
+   module T= std(T2(j));
+   vdim(T);hilb(T);
+   ring r1=0,(x,y,z),dp;
+   ideal id=xy,xz,yz;
+   list L=T2(id,"");
+   vdim(std(L[1]));                      // vdim of T2
+   L[4];                                 // 2nd syzygy module of ideal
+   t2 = std(t2);
+   dbprint(printlevel-voice+3,"// dim T2 = "+string(vdim(t2)));
+   if( size(#)>0 )
+   {
+      B2 = fetch(R,B(2));        // generators of Hom_P(syz/kos,R)
+      L  = t2,J0,A(2),SK,B2;
+      return(L);
+   }
+   return(t2);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   int p      = printlevel;
+   printlevel = 1;
+   ring  r    = 32003,(x,y),(c,dp);
+   ideal j    = x6-y4,x6y6,x2y4-x5y2;
+   module T   = T2(j);
+   vdim(T);
+   hilb(T);"";
+   ring r1    = 0,(x,y,z),dp;
+   ideal id   = xy,xz,yz;
+   list L     = T2(id,"");
+   vdim(L[1]);                           // vdim of T2
+   print(L[3]);                          // syzygy module of id
+   printlevel = p;
+}
 ///////////////////////////////////////////////////////////////////////////////
 …
 proc T12 (ideal i, list #)
 USAGE:   T12(i[,any]);  i = ideal
+DISPLAY: dim T1 and dim T2  of i
+RETURN:  T12(i): list of 2 modules:
              presentation of T1-module =T1(i) , 1st order deformations
              presentation of T2-module =T2(i) , obstructions of R=P/i
          T12(i,...): If a second argument is present (of any type) return
          a list of 9 modules, matrices, integers:
              [1]= presentation of T1 (module)
              [2]= presentation of T2 (module)
              [3]= matrix, whose cols present infinitesimal deformations
              [4]= matrix, whose cols are generators of relations of i
              [5]= matrix, presenting Hom_P([4]/kos,R), lifted to P
              [6]= standard basis of T1-module
              [7]= standard basis of T2-module
              [8]= vdim of T1
+             [9]= vdim of T2
+RETURN:  T12(i): list of 2 modules:
+             std basis of T1-module =T1(i), 1st order deformations
+             std basid of T2-module =T2(i), obstructions of R=P/i
+         If a second argument is present (of any type) return a list of
+modules, matrices, integers:
+             [1]= standard basis of T1-module
+             [2]= standard basis of T2-module
+             [3]= vdim of T1
+             [4]= vdim of T2
+             [5]= matrix, whose cols present infinitesimal deformations
+             [6]= matrix, whose cols are generators of relations of i (=syz(i))
+             [7]= matrix, presenting Hom_P(syz/kos,R), lifted to P
+             [8]= presentation of T1-module, no std basis
+             [9]= presentation of T2-module, no std basis
+DISPLAY: k-dimension of T1 and T2 if printlevel >= 0 (default)
 NOTE:    Use proc miniversal from deform.lib to get miniversal deformation of i,
          the list contains all objects used by proc miniversal
 …
   def P = basering;
    i = simplify(i,10);
+   if (size(i)<2)
+   {
+     "// hypersurface, use proc 'Tjurina'";
+     //return([1]);
+   }
+   module jac,t1,t2,kos,sbt1,sbt2;
+   matrix L3,L4,L5;
+   module jac,t1,t2,sbt1,sbt2;
+   matrix Kos,Syz,SK,kbT1,Sx;
+   list L;
    ideal  i0 = std(i);
 //-------------------- presentation of non-trivial syzygies -------------------
    list I= res(i,3);                            // resolve i
    L4  = matrix(I[2]);                          // syz(i)
+   list I= res(i,2);                            // resolve i
+   Syz  = matrix(I[2]);                         // syz(i)
    jac = jacob(i);                              // jacobi ideal
+   t1  = transpose(I[2]);                       // dual of syzygies
+   kos = koszul(2,i);                           // koszul-relations
+   t2  = lift(I[2],kos)+I[3];                   // presentation of syz/kos
+   Kos = koszul(2,i);                           // koszul-relations
+   SK  = modulo(Syz,Kos);                       // presentation of syz/kos
 //--------------------- fetch to quotient ring  mod i -------------------------
   qring   Ox  = i0;                             // make P/i the basering
    module jac = fetch(P,jac);
    module t1  = fetch(P,t1);                    // Hom(syz,R)
    module t2  = transpose(fetch(P,t2));         // Hom(syz/kos,R)
    list resS  = res(t1,3);
    list resR  = res(t2,3);
    t2 = lift(resR[2],t1)+resR[3];               // presentation of T2
    r2 = rank(t2);
    t1 = lift(resS[2],jac)+resS[3];              // presentation of T1
    r1 = rank(t1);
    matrix L3  = resS[2];
    matrix L5  = resR[2];
+   module Jac = fetch(P,jac);
+   matrix No  = transpose(fetch(P,Syz));        // ker(No) = Hom(syz,Ox)
+   module So  = transpose(fetch(P,SK));         // Hom(syz/kos,R)
+   list resS  = res(So,2);
+   matrix Sx  = resS[2];
+   list resN  = res(No,2);
+   matrix Nx  = resN[2];
+   module T2  = modulo(Sx,No);                  // presentation of T2
+   r2         = nrows(T2);
+   module T1  = modulo(Nx,Jac);                 // presentation of T1
+   r1         = nrows(T1);
 //------------------------ pull back to basering ------------------------------
   setring P;
    t1   = fetch(Ox,t1)+i*freemodule(r1);
    t2   = fetch(Ox,t2)+i*freemodule(r2);
+   t1   = fetch(Ox,T1)+i*freemodule(r1);
+   t2   = fetch(Ox,T2)+i*freemodule(r2);
    sbt1 = std(t1);
    d1   = vdim(sbt1);
    sbt2=std(t2);
+   sbt2 = std(t2);
    d2   = vdim(sbt2);
+   "// dim T1  = ",d1;
+   "// dim T2  = ",d2;
+   dbprint(printlevel-voice+3,"// dim T1 = "+string(d1),"// dim T2 = "+string(d2));
    if  ( size(#)>0)
+   {
+     L3 = fetch(Ox,L3)*kbase(sbt1);
+     L5 = fetch(Ox,L5);
+     return(t1,t2,L3,L4,L5,sbt1,sbt2,d1,d2);
+   }
+   return(t1,t2);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=200,(x,y,z,u,v),(c,ws(4,3,2,3,4));
+   ideal i=xz-y2,yz2-xu,xv-yzu,yu-z3,z2u-yv,zv-u2;
+   //a cyclic quotient singularity
+   list L = T12(i,1);
+   print(L[3]);
+}
+///////////////////////////////////////////////////////////////////////////////
+      kbT1 = fetch(Ox,Nx)*kbase(sbt1);
+      Sx   = fetch(Ox,Sx);
+      L = sbt1,sbt2,d1,d2,kbT1,Syz,Sx,t1,t2;
+      return(L);
+   }
+   L = sbt1,sbt2;
+   return(L);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   int p      = printlevel;
+   printlevel = 1;
+   ring r     = 200,(x,y,z,u,v),(c,ws(4,3,2,3,4));
+   ideal i    = xz-y2,yz2-xu,xv-yzu,yu-z3,z2u-yv,zv-u2;
+                            //a cyclic quotient singularity
+   list L     = T12(i,1);
+   print(L[5]);             //matrix of infin. deformations
+   printlevel = p;
+}
+///////////////////////////////////////////////////////////////////////////////

Context Navigation

Changeset 6f2edc in git for Singular/LIB

Legend:

Singular/LIB/deform.lib

Singular/LIB/elim.lib

Singular/LIB/factor.lib

Singular/LIB/general.lib

Singular/LIB/homog.lib

Singular/LIB/homolog.lib

Singular/LIB/inout.lib

Singular/LIB/matrix.lib

Singular/LIB/poly.lib

Singular/LIB/random.lib

Singular/LIB/redfac.lib

Singular/LIB/ring.lib

Singular/LIB/sing.lib

Download in other formats: