singular-data 4.0.3+ds-1 / usr / share / singular / LIB / mregular.lib

//////////////////////////////////////////////////////////////////////////////
version="version mregular.lib 4.0.0.0 Jun_2013 "; // $Id: 65c9ff28baa31013620d93b4f3033bcb6d5632ad $
category="Commutative Algebra";
info="
LIBRARY: mregular.lib   Castelnuovo-Mumford regularity of homogeneous ideals
AUTHORS: I.Bermejo,     ibermejo@ull.es
@*       Ph.Gimenez,    pgimenez@agt.uva.es
@*       G.-M.Greuel,   greuel@mathematik.uni-kl.de

OVERVIEW:
 A library for computing the Castelnuovo-Mumford regularity of a homogeneous
 ideal that DOES NOT require the computation of a minimal graded free
 resolution of the ideal.
 It also determines depth(basering/ideal) and satiety(ideal).
 The procedures are based on 3 papers by Isabel Bermejo and Philippe Gimenez:
 'On Castelnuovo-Mumford regularity of projective curves' Proc.Amer.Math.Soc.
 128(5) (2000), 'Computing the Castelnuovo-Mumford regularity of some
 subschemes of Pn using quotients of monomial ideals', Proceedings of
 MEGA-2000, J. Pure Appl. Algebra 164 (2001), and 'Saturation and
 Castelnuovo-Mumford regularity', Preprint (2004).

PROCEDURES:
 regIdeal(id,[,e]);    regularity of homogeneous ideal id
 depthIdeal(id,[,e]);  depth of S/id with S=basering, id homogeneous ideal
 satiety(id,[,e]);     saturation index of homogeneous ideal id
 regMonCurve(li);      regularity of projective monomial curve defined by li
 NoetherPosition(id);  Noether normalization of ideal id
 is_NP(id);            checks whether variables are in Noether position
 is_nested(id);        checks whether monomial ideal id is of nested type
";

LIB "general.lib";
LIB "algebra.lib";
LIB "sing.lib";
LIB "poly.lib";
//////////////////////////////////////////////////////////////////////////////
//
proc regIdeal (ideal i, list #)
"
USAGE:   regIdeal (i[,e]); i ideal, e integer
RETURN:  an integer, the Castelnuovo-Mumford regularity of i.
         (returns -1 if i is not homogeneous)
ASSUME:  i is a homogeneous ideal of the basering S=K[x(0)..x(n)].
         e=0:  (default)
               If K is an infinite field, makes random changes of coordinates.
               If K is a finite field, works over a transcendental extension.
         e=1:  Makes random changes of coordinates even when K is finite.
               It works if it terminates, but may result in an infinite
               loop. After 30 loops, a warning message is displayed and
               -1 is returned.
NOTE:    If printlevel > 0 (default = 0), additional info is displayed:
         dim(S/i), depth(S/i) and end(H^(depth(S/i))(S/i)) are computed,
         and an upper bound for the a-invariant of S/i is given.
         The algorithm also determines whether the regularity is attained
         or not at the last step of a minimal graded free resolution of i,
         and if the answer is positive, the regularity of the Hilbert
         function of S/i is given.
EXAMPLE: example regIdeal; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
   int e,ii,jj,H,h,d,time,lastv,sat,firstind;
   int lastind,ch,nesttest,NPtest,nl,N,acc;
   intmat ran;
   def r0 = basering;
   int n = nvars(r0)-1;
   if ( size(#) > 0 )
   {
      e = #[1];
   }
   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
   execute(s);
   ideal i,sbi,I,J,K,chcoord,m;
   poly P;
   map phi;
   i = fetch(r0,i);
   time=rtimer;
   sbi=std(i);
   ch=char(r1);
//----- Check ideal homogeneous
   if ( homog(sbi) == 0 )
   {
       "// WARNING from proc regIdeal from lib mregular.lib:
// The ideal is not homogeneous!";
       return (-1);
   }
   I=simplify(lead(sbi),1);
   attrib(I,"isSB",1);
   d=dim(I);
//----- If the ideal i is not proper:
   if ( d == -1 )
   {
     dbprint(printlevel-voice+2,
             "// The ideal i is (1)!
// Its Castelnuovo-Mumford regularity is:");
     return (0);
   }
//----- If the ideal i is 0:
   if ( size(I) == 0 )
   {
     dbprint(printlevel-voice+2,
             "// The ideal i is (0)!
// Its Castelnuovo-Mumford regularity is:");
     return (0);
   }
//----- When the ideal i is 0-dimensional:
   if ( d == 0 )
   {
     H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
     time=rtimer-time;
     // Additional information:
     dbprint(printlevel-voice+2,
             "// Dimension of S/i : 0");
     dbprint(printlevel-voice+2,
             "// Time for computing regularity: " + string(time) + " sec.");
     dbprint(printlevel-voice+2,
"// The Castelnuovo-Mumford regularity of i coincides with its satiety, and
 // with the regularity of the Hilbert function of S/i. Its value is:");
     return (H);
   }
//----- Determine the situation: NT, or NP, or nothing.
//----- Choose the method depending on the situation, on the
//----- characteristic of the ground field, and on the option argument
//----- in order to get the mon. ideal of nested type associated to i
   if ( e == 1 )
   { ch=0; }
   NPtest=is_NP(I);
   if ( NPtest == 1 )
   {
     nesttest=is_nested(I);
   }
   if ( ch != 0 )
   {
     if ( NPtest == 0 )
     {
       N=d*n-d*(d-1)/2;
       s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
       execute(s);
       ideal chcoord,m,i,I;
       poly P;
       map phi;
       i=imap(r1,i);
       chcoord=select1(maxideal(1),1..(n-d+1));
       acc=0;
       for ( ii = 1; ii<=d; ii++ )
       {
         matrix trex[1][n-d+ii+1]=t((1+acc)..(n-d+ii+acc)),1;
         m=select1(maxideal(1),1..(n-d+1+ii));
         for ( jj = 1; jj<=n-d+ii+1; jj++ )
         {
           P=P+trex[1,jj]*m[jj];
         }
         chcoord[n-d+1+ii]=P;
         P=0;
         acc=acc+n-d+ii;
         kill trex;
       }
       phi=rtr,chcoord;
       I=simplify(lead(std(phi(i))),1);
       setring r1;
       I=imap(rtr,I);
       attrib(I,"isSB",1);
     }
     else
     {
       if ( nesttest == 0 )
       {
         N=d*(d-1)/2;
         s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
         execute(s);
         ideal chcoord,m,i,I;
         poly P;
         map phi;
         i=imap(r1,i);
         chcoord=select1(maxideal(1),1..(n-d+2));
         acc=0;
         for ( ii = 1; ii<=d-1; ii++ )
         {
           matrix trex[1][ii+1]=t((1+acc)..(ii+acc)),1;
           m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
           for ( jj = 1; jj<=ii+1; jj++ )
           {
             P=P+trex[1,jj]*m[jj];
           }
           chcoord[n-d+2+ii]=P;
           P=0;
           acc=acc+ii;
           kill trex;
         }
         phi=rtr,chcoord;
         I=simplify(lead(std(phi(i))),1);
         setring r1;
         I=imap(rtr,I);
         attrib(I,"isSB",1);
       }
     }
   }
   else
   {
     if ( NPtest == 0 )
     {
       while ( nl < 30 )
       {
         chcoord=select1(maxideal(1),1..(n-d+1));
         nl=nl+1;
         for ( ii = 1; ii<=d; ii++ )
         {
           ran=random(100,1,n-d+ii);
           ran=intmat(ran,1,n-d+ii+1);
           ran[1,n-d+ii+1]=1;
           m=select1(maxideal(1),1..(n-d+1+ii));
           for ( jj = 1; jj<=n-d+ii+1; jj++ )
           {
             P=P+ran[1,jj]*m[jj];
           }
           chcoord[n-d+1+ii]=P;
           P=0;
         }
         phi=r1,chcoord;
         dbprint(printlevel-voice+2,"// (1 random change of coord.)");
         I=simplify(lead(std(phi(i))),1);
         attrib(I,"isSB",1);
         NPtest=is_NP(I);
         if ( NPtest == 1 )
         {
           break;
         }
       }
       if ( NPtest == 0 )
       {
       "// WARNING from proc regIdeal from lib mregular.lib:
// The procedure has entered in 30 loops and could not put the variables
// in Noether position: in your example the method using random changes
// of coordinates may enter an infinite loop when the field is finite.
// Try removing this optional argument.";
         return (-1);
       }
       i=phi(i);
       nesttest=is_nested(I);
     }
     if ( nesttest == 0 )
     {
       while ( nl < 30 )
       {
         chcoord=select1(maxideal(1),1..(n-d+2));
         nl=nl+1;
         for ( ii = 1; ii<=d-1; ii++ )
         {
           ran=random(100,1,ii);
           ran=intmat(ran,1,ii+1);
           ran[1,ii+1]=1;
           m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
           for ( jj = 1; jj<=ii+1; jj++ )
           {
             P=P+ran[1,jj]*m[jj];
           }
           chcoord[n-d+2+ii]=P;
           P=0;
         }
         phi=r1,chcoord;
         dbprint(printlevel-voice+2,"// (1 random change of coord.)");
         I=simplify(lead(std(phi(i))),1);
         attrib(I,"isSB",1);
         nesttest=is_nested(I);
         if ( nesttest == 1 )
         {
           break;
         }
       }
       if ( nesttest == 0 )
       {
       "// WARNING from proc regIdeal from lib mregular.lib:
// The procedure has entered in 30 loops and could not find a monomial
// ideal of nested type with the same regularity as your ideal: in your
// example the method using random changes of coordinates may enter an
// infinite loop when the field is finite.
// Try removing this optional argument.";
         return (-1);
       }
     }
   }
//
// At this stage, we have obtained a monomial ideal I of nested type
// such that reg(i)=reg(I). We now compute reg(I).
//
//----- When S/i is Cohen-Macaulay:
   for ( ii = n-d+2; ii <= n+1; ii++ )
   {
     K=K+select(I,ii);
   }
   if ( size(K) == 0 )
   {
     s="ring nr = ",charstr(r0),",x(0..n-d),dp;";
     execute(s);
     ideal I;
     I = imap(r1,I);
     H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
     time=rtimer-time;
     // Additional information:
     dbprint(printlevel-voice+2,
               "// S/i is Cohen-Macaulay");
     dbprint(printlevel-voice+2,
               "// Dimension of S/i ( = depth(S/i) ): "+string(d));
     dbprint(printlevel-voice+2,
               "// Regularity attained at the last step of m.g.f.r. of i: YES");
     dbprint(printlevel-voice+2,
               "// Regularity of the Hilbert function of S/i: " + string(H-d));
     dbprint(printlevel-voice+2,
               "// Time for computing regularity: " + string(time) + " sec.");
     dbprint(printlevel-voice+2,
               "// The Castelnuovo-Mumford regularity of i is:");
     return(H);
   }
//----- When d=1:
   if ( d == 1 )
   {
     H=maxdeg1(simplify(reduce(quotient(I,maxideal(1)),I),2))+1;
     sat=H;
     J=subst(I,x(n),1);
     s = "ring nr = ",charstr(r0),",x(0..n-1),dp;";
     execute(s);
     ideal J=imap(r1,J);
     attrib(J,"isSB",1);
     h=maxdeg1(minbase(quotient(J,maxideal(1))))+1;
     time=rtimer-time;
     if ( h > H )
     {
       H=h;
     }
     // Additional information:
     dbprint(printlevel-voice+2,
               "// Dimension of S/i: 1");
     dbprint(printlevel-voice+2,
               "// Depth of S/i: 0");
     dbprint(printlevel-voice+2,
               "// Satiety of i: "+string(sat));
     dbprint(printlevel-voice+2,
               "// Upper bound for the a-invariant of S/i: end(H^1(S/i)) <= "+
               string(h-2));
     if ( H == sat )
     {
       dbprint(printlevel-voice+2,
                   "// Regularity attained at the last step of m.g.f.r. of i: YES");
       dbprint(printlevel-voice+2,
                   "// Regularity of the Hilbert function of S/i: "+string(H));
     }
     else
     {
       dbprint(printlevel-voice+2,
                   "// Regularity attained at the last step of m.g.f.r. of i: NO");
     }
     dbprint(printlevel-voice+2,
                   "// Time for computing regularity: "+ string(time) + " sec.");
     dbprint(printlevel-voice+2,
                   "// The Castelnuovo-Mumford regularity of i is:");
     return(H);
   }
//----- Now d>1 and S/i is not Cohen-Macaulay:
//
//----- First, determine the last variable really occuring
   lastv=n-d;
   h=n;
   while ( lastv == n-d and h > n-d )
   {
     K=select(I,h+1);
     if ( size(K) == 0 )
     {
       h=h-1;
     }
     else
     {
       lastv=h;
     }
   }
//----- and compute Castelnuovo-Mumford regularity:
   s = "ring nr = ",charstr(r0),",x(0..lastv),dp;";
   execute(s);
   ideal I,K,KK,LL;
   I=imap(r1,I);
   attrib(I,"isSB",1);
   K=simplify(reduce(quotient(I,maxideal(1)),I),2);
   H=maxdeg1(K)+1;
   firstind=H;
   KK=minbase(subst(I,x(lastv),1));
   for ( ii = n-lastv; ii<=d-2; ii++ )
   {
     LL=minbase(subst(I,x(n-ii-1),1));
     attrib(LL,"isSB",1);
     s = "ring mr = ",charstr(r0),",x(0..n-ii-1),dp;";
     execute(s);
     ideal K,KK;
     KK=imap(nr,KK);
     attrib(KK,"isSB",1);
     K=simplify(reduce(quotient(KK,maxideal(1)),KK),2);
     h=maxdeg1(K)+1;
     if ( h > H )
     {
       H=h;
     }
     setring nr;
     kill mr;
     KK=LL;
   }
   // We must determine one more sat. index:
   s = "ring mr = ",charstr(r0),",x(0..n-d),dp;";
   execute(s);
   ideal KK,K;
   KK=imap(nr,KK);
   attrib(KK,"isSB",1);
   K=simplify(reduce(quotient(KK,maxideal(1)),KK),2);
   h=maxdeg1(K)+1;
   lastind=h;
   if ( h > H )
   {
     H=h;
   }
   setring nr;
   kill mr;
   time=rtimer-time;
   // Additional information:
   dbprint(printlevel-voice+2,
               "// Dimension of S/i: "+string(d));
   dbprint(printlevel-voice+2,
               "// Depth of S/i: "+string(n-lastv));
   dbprint(printlevel-voice+2,
               "// end(H^"+string(n-lastv)+"(S/i)) = "
               +string(firstind-n+lastv-1));
   dbprint(printlevel-voice+2,
               "// Upper bound for the a-invariant of S/i: end(H^"
               +string(d)+"(S/i)) <= "+string(lastind-d-1));
   if ( H == firstind )
   {
     dbprint(printlevel-voice+2,
                   "// Regularity attained at the last step of m.g.f.r. of i: YES");
     dbprint(printlevel-voice+2,
                   "// Regularity of the Hilbert function of S/i: "
                   +string(H-n+lastv));
   }
   else
   {
     dbprint(printlevel-voice+2,
                   "// Regularity attained at the last step of m.g.f.r. of i: NO");
   }
   dbprint(printlevel-voice+2,
               "// Time for computing regularity: "+ string(time) + " sec.");
   dbprint(printlevel-voice+2,
               "// The Castelnuovo-Mumford regularity of i is:");
   return(H);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z,t,w),dp;
   ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
   regIdeal(i);
   regIdeal(lead(std(i)));
// Additional information is displayed if you change printlevel (=1);
}
////////////////////////////////////////////////////////////////////////////////
/*
Out-commented examples:
//
   ring s=0,x(0..5),dp;
   ideal i=x(2)^2-x(4)*x(5),x(1)*x(2)-x(0)*x(5),x(0)*x(2)-x(1)*x(4),
           x(1)^2-x(3)*x(5),x(0)*x(1)-x(2)*x(3),x(0)^2-x(3)*x(4);
   regIdeal(i);
   // Our procedure works when a min. graded free resol. can
   // not be computed. In this easy example, regularity can also
   // be obtained using a m.g.f.r.:
   nrows(betti(mres(i,0)));
   ring r1=0,(x,y,z,t),dp;
// Ex.2.5 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
   ideal i  = x17y14-y31, x20y13, x60-y36z24-x20z20t20;
   regIdeal(i);
// Ex.2.9 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
   int k=43;
   ideal j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
   regIdeal(j);
   k=14;
   j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
   regIdeal(j);
   k=22;
   j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
   regIdeal(j);
   k=315;
   j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
   regIdeal(j);
// Example in Rk.2.10 in [Bermejo-Gimenez], ProcAMS 128(5):
   ideal h=x2-3xy+5xt,xy-3y2+5yt,xz-3yz,2xt-yt,y2-yz-2yt;
   regIdeal(h);
// The initial ideal is not saturated
   regIdeal(lead(std(h)));
// More examples:
   i=y4-t3z, x3t-y2z2, x3y2-t2z3, x6-tz5;
   regIdeal(i);
//
   regIdeal(maxideal(4));
//
   ring r2=0,(x,y,z,t,w),dp;
   ideal i = xy-zw,x3-yw2,x2z-y2w,y3-xz2,-y2z3+xw4+tw4+w5,-yz4+x2w3+xtw3+xw4,
            -z5+x2tw2+x2w3+yw4;
   regIdeal(i);
//
   ring r3=0,(x,y,z,t,w,u),dp;
   ideal i=imap(r2,i);
   regIdeal(i);
// Next example is the defining ideal of the 2nd. Veronesean of P3, a variety
// in P8 which is arithmetically Cohen-Macaulay:
   ring r4=0,(a,b,c,d,x(0..9)),dp;
   ideal i= x(0)-ab,x(1)-ac,x(2)-ad,x(3)-bc,x(4)-bd,x(5)-cd,
            x(6)-a2,x(7)-b2,x(8)-c2,x(9)-d2;
   ideal ei=eliminate(i,abcd);
   ring r5=0,x(0..9),dp;
   ideal i=imap(r4,ei);
   regIdeal(i);
// Here is an example where the computation of a m.g.f.r. of I costs:
   ring r8=0,(x,y,z,t,u,a,b),dp;
   ideal i=u-b40,t-a40,x-a23b17,y-a22b18+ab39,z-a25b15;
   ideal ei=eliminate(i,ab); // It takes a few seconds to compute the ideal
   ring r9=0,(x,y,z,t,u),dp;
   ideal i=imap(r8,ei);
   regIdeal(i);   // This is very fast.
// Now you can use mres(i,0) to compute a m.g.f.r. of the ideal!
//
// The computation of the m.g.f.r. of the following example did not succeed
// using the command mres:
   ring r10=0,(x(0..8),s,t),dp;
   ideal i=x(0)-st24,x(1)-s2t23,x(2)-s3t22,x(3)-s9t16,x(4)-s11t14,x(5)-s18t7,
           x(6)-s24t,x(7)-t25,x(8)-s25;
   ideal ei=eliminate(i,st);
   ring r11=0,x(0..8),dp;
   ideal i=imap(r10,ei);
   regIdeal(i);
// More examples where not even sres works:
// Be careful: elimination takes some time here, but it succeeds!
   ring r12=0,(s,t,u,x(0..14)),dp;
   ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,x(4)-s2t7u6,x(5)-s7t7u,
           x(6)-s10t5,x(7)-s4t6u5,x(8)-s13tu,x(9)-s14u,x(10)-st2u12,x(11)-s3t9u3,
           x(12)-s15,x(13)-t15,x(14)-u15;
   ideal ei=eliminate(i,stu);
   size(ei);
   ring r13=0,x(0..14),dp;
   ideal i=imap(r12,ei);
   size(i);
   regIdeal(i);
*/
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

proc depthIdeal (ideal i, list #)
"
USAGE:   depthIdeal (i[,e]); i ideal, e integer
RETURN:  an integer, the depth of S/i where S=K[x(0)..x(n)] is the basering.
         (returns -1 if i is not homogeneous or if i=(1))
ASSUME:  i is a proper homogeneous ideal.
         e=0:  (default)
               If K is an infinite field, makes random changes of coordinates.
               If K is a finite field, works over a transcendental extension.
         e=1:  Makes random changes of coordinates even when K is finite.
               It works if it terminates, but may result in an infinite
               loop. After 30 loops, a warning message is displayed and
               -1 is returned.
NOTE:    If printlevel > 0 (default = 0), dim(S/i) is also displayed.
EXAMPLE: example depthIdeal; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
   int e,ii,jj,h,d,time,lastv,ch,nesttest,NPtest,nl,N,acc;
   intmat ran;
   def r0 = basering;
   int n = nvars(r0)-1;
   if ( size(#) > 0 )
   {
      e = #[1];
   }
   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
   execute(s);
   ideal i,sbi,I,J,K,chcoord,m;
   poly P;
   map phi;
   i = fetch(r0,i);
   time=rtimer;
   sbi=std(i);
   ch=char(r1);
//----- Check ideal homogeneous
   if ( homog(sbi) == 0 )
   {
       "// WARNING from proc depthIdeal from lib mregular.lib:
// The ideal is not homogeneous!";
       return (-1);
   }
   I=simplify(lead(sbi),1);
   attrib(I,"isSB",1);
   d=dim(I);
//----- If the ideal i is not proper:
   if ( d == -1 )
     {
       "// WARNING from proc depthIdeal from lib mregular.lib:
// The ideal i is (1)!";
       return (-1);
     }
//----- If the ideal i is 0:
   if ( size(I) == 0 )
     {
       dbprint(printlevel-voice+2,
               "// The ideal i is (0)!
// The depth of S/i is:");
       return (d);
     }
//----- When the ideal i is 0-dimensional:
   if ( d == 0 )
     {
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i : 0 (S/i is Cohen-Macaulay)");
       dbprint(printlevel-voice+2,
               "// Time for computing the depth: " + string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The depth of S/i is:");
       return (0);
     }
//----- Determine the situation: NT, or NP, or nothing.
//----- Choose the method depending on the situation, on the
//----- characteristic of the ground field, and on the option argument
//----- in order to get the mon. ideal of nested type associated to i
   if ( e == 1 )
     {
       ch=0;
     }
   NPtest=is_NP(I);
   if ( NPtest == 1 )
     {
       nesttest=is_nested(I);
     }
   if ( ch != 0 )
     {
       if ( NPtest == 0 )
         {
           N=d*n-d*(d-1)/2;
           s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
           execute(s);
           ideal chcoord,m,i,I;
           poly P;
           map phi;
           i=imap(r1,i);
           chcoord=select1(maxideal(1),1..(n-d+1));
           acc=0;
           for ( ii = 1; ii<=d; ii++ )
             {
               matrix trex[1][n-d+ii+1]=t((1+acc)..(n-d+ii+acc)),1;
               m=select1(maxideal(1),1..(n-d+1+ii));
               for ( jj = 1; jj<=n-d+ii+1; jj++ )
                 {
                   P=P+trex[1,jj]*m[jj];
                 }
               chcoord[n-d+1+ii]=P;
               P=0;
               acc=acc+n-d+ii;
               kill trex;
             }
               phi=rtr,chcoord;
               I=simplify(lead(std(phi(i))),1);
               setring r1;
               I=imap(rtr,I);
               attrib(I,"isSB",1);
         }
       else
         {
           if ( nesttest == 0 )
             {
               N=d*(d-1)/2;
               s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
               execute(s);
               ideal chcoord,m,i,I;
               poly P;
               map phi;
               i=imap(r1,i);
               chcoord=select1(maxideal(1),1..(n-d+2));
               acc=0;
               for ( ii = 1; ii<=d-1; ii++ )
                 {
                   matrix trex[1][ii+1]=t((1+acc)..(ii+acc)),1;
                    m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
                   for ( jj = 1; jj<=ii+1; jj++ )
                     {
                       P=P+trex[1,jj]*m[jj];
                     }
                   chcoord[n-d+2+ii]=P;
                   P=0;
                   acc=acc+ii;
                   kill trex;
                 }
               phi=rtr,chcoord;
               I=simplify(lead(std(phi(i))),1);
               setring r1;
               I=imap(rtr,I);
               attrib(I,"isSB",1);
             }
         }
     }
   else
     {
       if ( NPtest == 0 )
         {
           while ( nl < 30 )
             {
               chcoord=select1(maxideal(1),1..(n-d+1));
               nl=nl+1;
               for ( ii = 1; ii<=d; ii++ )
                 {
                   ran=random(100,1,n-d+ii);
                   ran=intmat(ran,1,n-d+ii+1);
                   ran[1,n-d+ii+1]=1;
                   m=select1(maxideal(1),1..(n-d+1+ii));
                   for ( jj = 1; jj<=n-d+ii+1; jj++ )
                     {
                       P=P+ran[1,jj]*m[jj];
                     }
                   chcoord[n-d+1+ii]=P;
                   P=0;
                 }
               phi=r1,chcoord;
               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
               I=simplify(lead(std(phi(i))),1);
               attrib(I,"isSB",1);
               NPtest=is_NP(I);
               if ( NPtest == 1 )
                 {
                   break;
                 }
             }
           if ( NPtest == 0 )
             {
       "// WARNING from proc depthIdeal from lib mregular.lib:
// The procedure has entered in 30 loops and could not put the variables
// in Noether position: in your example the method using random changes
// of coordinates may enter an infinite loop when the field is finite.
// Try removing this optional argument.";
       return (-1);
             }
           i=phi(i);
           nesttest=is_nested(I);
         }
       if ( nesttest == 0 )
         {
           while ( nl < 30 )
             {
               chcoord=select1(maxideal(1),1..(n-d+2));
               nl=nl+1;
               for ( ii = 1; ii<=d-1; ii++ )
                 {
                   ran=random(100,1,ii);
                   ran=intmat(ran,1,ii+1);
                   ran[1,ii+1]=1;
                   m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
                   for ( jj = 1; jj<=ii+1; jj++ )
                     {
                       P=P+ran[1,jj]*m[jj];
                     }
                   chcoord[n-d+2+ii]=P;
                   P=0;
                 }
               phi=r1,chcoord;
               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
               I=simplify(lead(std(phi(i))),1);
               attrib(I,"isSB",1);
               nesttest=is_nested(I);
               if ( nesttest == 1 )
                 {
                   break;
                 }
             }
           if ( nesttest == 0 )
             {
       "// WARNING from proc depthIdeal from lib mregular.lib:
// The procedure has entered in 30 loops and could not find a monomial
// ideal of nested type with the same depth as your ideal: in your
// example the method using random changes of coordinates may enter an
// infinite loop when the field is finite.
// Try removing this optional argument.";
       return (-1);
             }
         }
     }
//
// At this stage, we have obtained a monomial ideal I of nested type
// such that depth(S/i)=depth(S/I). We now compute depth(I).
//
//----- When S/i is Cohen-Macaulay:
   for ( ii = n-d+2; ii <= n+1; ii++ )
     {
       K=K+select(I,ii);
     }
   if ( size(K) == 0 )
     {
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: "+string(d)+" (S/i is Cohen-Macaulay)");
       dbprint(printlevel-voice+2,
               "// Time for computing depth: " + string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The depth of S/i is:");
       return(d);
     }
//----- When d=1 (and S/i is not Cohen-Macaulay) ==> depth =0:
   if ( d == 1 )
     {
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: 1");
       dbprint(printlevel-voice+2,
               "// Time for computing depth: "+ string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The depth of S/i is:");
       return(0);

     }
//----- Now d>1 and S/i is not Cohen-Macaulay:
//
//----- First, determine the last variable really occuring
       lastv=n-d;
       h=n;
       while ( lastv == n-d and h > n-d )
         {
           K=select(I,h+1);
           if ( size(K) == 0 )
             {
               h=h-1;
             }
           else
             {
               lastv=h;
             }
         }
//----- and compute the depth:
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: "+string(d));
       dbprint(printlevel-voice+2,
               "// Time for computing depth: "+ string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The depth of S/i is:");
       return(n-lastv);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z,t,w),dp;
   ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
   depthIdeal(i);
   depthIdeal(lead(std(i)));
// Additional information is displayed if you change printlevel (=1);
}
////////////////////////////////////////////////////////////////////////////////
/*
Out-commented examples:
   ring s=0,x(0..5),dp;
   ideal i=x(2)^2-x(4)*x(5),x(1)*x(2)-x(0)*x(5),x(0)*x(2)-x(1)*x(4),
           x(1)^2-x(3)*x(5),x(0)*x(1)-x(2)*x(3),x(0)^2-x(3)*x(4);
   depthIdeal(i);
   // Our procedure works when a min. graded free resol. can
   // not be computed. In this easy example, depth can also
   // be obtained using a m.g.f.r. (Auslander-Buchsbaum formula):
   nvars(s)-ncols(betti(mres(i,0)))+1;
   ring r1=0,(x,y,z,t),dp;
// Ex.2.5 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
   ideal i  = x17y14-y31, x20y13, x60-y36z24-x20z20t20;
   depthIdeal(i);
// Ex.2.9 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
   int k=43;
   ideal j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
   depthIdeal(j);
// Example in Rk.2.10 in [Bermejo-Gimenez], ProcAMS 128(5):
   ideal h=x2-3xy+5xt,xy-3y2+5yt,xz-3yz,2xt-yt,y2-yz-2yt;
   depthIdeal(h);
// The initial ideal is not saturated
   depthIdeal(lead(std(h)));
// More examples:
   i=y4-t3z, x3t-y2z2, x3y2-t2z3, x6-tz5;
   depthIdeal(i);
//
   depthIdeal(maxideal(4));
//
   ring r2=0,(x,y,z,t,w),dp;
   ideal i = xy-zw,x3-yw2,x2z-y2w,y3-xz2,-y2z3+xw4+tw4+w5,-yz4+x2w3+xtw3+xw4,
            -z5+x2tw2+x2w3+yw4;
   depthIdeal(i);
//
   ring r3=0,(x,y,z,t,w,u),dp;
   ideal i=imap(r2,i);
   depthIdeal(i);
// Next example is the defining ideal of the 2nd. Veronesean of P3, a variety
// in P8 which is arithmetically Cohen-Macaulay:
   ring r4=0,(a,b,c,d,x(0..9)),dp;
   ideal i= x(0)-ab,x(1)-ac,x(2)-ad,x(3)-bc,x(4)-bd,x(5)-cd,
            x(6)-a2,x(7)-b2,x(8)-c2,x(9)-d2;
   ideal ei=eliminate(i,abcd);
   ring r5=0,x(0..9),dp;
   ideal i=imap(r4,ei);
   depthIdeal(i);
// Here is an example where the computation of a m.g.f.r. of I costs:
   ring r8=0,(x,y,z,t,u,a,b),dp;
   ideal i=u-b40,t-a40,x-a23b17,y-a22b18+ab39,z-a25b15;
   ideal ei=eliminate(i,ab); // It takes a few seconds to compute the ideal
   ring r9=0,(x,y,z,t,u),dp;
   ideal i=imap(r8,ei);
   depthIdeal(i);   // This is very fast.
// Now you can use mres(i,0) to compute a m.g.f.r. of the ideal!
//
// Another one:
   ring r10=0,(x(0..8),s,t),dp;
   ideal i=x(0)-st24,x(1)-s2t23,x(2)-s3t22,x(3)-s9t16,x(4)-s11t14,x(5)-s18t7,
           x(6)-s24t,x(7)-t25,x(8)-s25;
   ideal ei=eliminate(i,st);
   ring r11=0,x(0..8),dp;
   ideal i=imap(r10,ei);
   depthIdeal(i);
// More examples where not even sres works:
// Be careful: elimination takes some time here, but it succeeds!
ring r12=0,(s,t,u,x(0..14)),dp;
ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,x(4)-s2t7u6,x(5)-s7t7u,
        x(6)-s10t5,x(7)-s4t6u5,x(8)-s13tu,x(9)-s14u,x(10)-st2u12,x(11)-s3t9u3,
        x(12)-s15,x(13)-t15,x(14)-u15;
ideal ei=eliminate(i,stu);
size(ei);
ring r13=0,x(0..14),dp;
ideal i=imap(r12,ei);
size(i);
depthIdeal(i);
//
*/
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

proc satiety (ideal i, list #)
"
USAGE:   satiety (i[,e]); i ideal, e integer
RETURN:  an integer, the satiety of i.
         (returns -1 if i is not homogeneous)
ASSUME:  i is a homogeneous ideal of the basering S=K[x(0)..x(n)].
         e=0:  (default)
               The satiety is computed determining the fresh elements in the
               socle of i. It works over arbitrary fields.
         e=1:  Makes random changes of coordinates to find a monomial ideal
               with same satiety. It works over infinite fields only. If K
               is finite, it works if it terminates, but may result in an
               infinite loop. After 30 loops, a warning message is displayed
               and -1 is returned.
THEORY:  The satiety, or saturation index, of a homogeneous ideal i is the
         least integer s such that, for all d>=s, the degree d part of the
         ideals i and isat=sat(i,maxideal(1))[1] coincide.
NOTE:    If printlevel > 0 (default = 0), dim(S/i) is also displayed.
EXAMPLE: example satiety; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
   int e,ii,jj,h,d,time,lastv,nesttest,NPtest,nl,sat;
   intmat ran;
   def r0 = basering;
   int n = nvars(r0)-1;
   if ( size(#) > 0 )
   {
      e = #[1];
   }
   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
   execute(s);
   ideal i,sbi,I,K,chcoord,m,KK;
   poly P;
   map phi;
   i = fetch(r0,i);
   time=rtimer;
   sbi=std(i);
//----- Check ideal homogeneous
   if ( homog(sbi) == 0 )
   {
       "// WARNING from proc satiety from lib mregular.lib:
// The ideal is not homogeneous!";
       return (-1);
   }
   I=simplify(lead(sbi),1);
   attrib(I,"isSB",1);
   d=dim(I);
//----- If the ideal i is not proper:
   if ( d == -1 )
     {
       dbprint(printlevel-voice+2,
               "// The ideal i is (1)!
// Its satiety is:");
       return (0);
     }
//----- If the ideal i is 0:
   if ( size(I) == 0 )
     {
       dbprint(printlevel-voice+2,
               "// The ideal i is (0)!
// Its satiety is:");
       return (0);
     }
//----- When the ideal i is 0-dimensional:
   if ( d == 0 )
     {
       sat=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: 0");
       dbprint(printlevel-voice+2,
               "// Time for computing the satiety: " + string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The satiety of i is:");
       return (sat);
     }
//----- When one has option e=1:
//
//----- Determine the situation: NT, or NP, or nothing.
//----- Choose the method depending on the situation in order to
//----- get the mon. ideal of nested type associated to i
   if ( e == 1 )
     {
       NPtest=is_NP(I);
       if ( NPtest == 0 )
         {
           while ( nl < 30 )
             {
               chcoord=select1(maxideal(1),1..(n-d+1));
               nl=nl+1;
               for ( ii = 1; ii<=d; ii++ )
                 {
                   ran=random(100,1,n-d+ii);
                   ran=intmat(ran,1,n-d+ii+1);
                   ran[1,n-d+ii+1]=1;
                   m=select1(maxideal(1),1..(n-d+1+ii));
                   for ( jj = 1; jj<=n-d+ii+1; jj++ )
                     {
                       P=P+ran[1,jj]*m[jj];
                     }
                   chcoord[n-d+1+ii]=P;
                   P=0;
                 }
               phi=r1,chcoord;
               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
               I=simplify(lead(std(phi(i))),1);
               attrib(I,"isSB",1);
               NPtest=is_NP(I);
               if ( NPtest == 1 )
                 {
                   break;
                 }
             }
           if ( NPtest == 0 )
             {
       "// WARNING from proc satiety from lib mregular.lib:
// The procedure has entered in 30 loops and could not put the variables
// in Noether position: in your example the method using random changes
// of coordinates may enter an infinite loop when the field is finite.
// Try removing the optional argument.";
       return (-1);
             }
           i=phi(i);
         }
       nesttest=is_nested(I);
       if ( nesttest == 0 )
         {
           while ( nl < 30 )
             {
               chcoord=select1(maxideal(1),1..(n-d+2));
               nl=nl+1;
               for ( ii = 1; ii<=d-1; ii++ )
                 {
                   ran=random(100,1,ii);
                   ran=intmat(ran,1,ii+1);
                   ran[1,ii+1]=1;
                   m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
                   for ( jj = 1; jj<=ii+1; jj++ )
                     {
                       P=P+ran[1,jj]*m[jj];
                     }
                   chcoord[n-d+2+ii]=P;
                   P=0;
                 }
               phi=r1,chcoord;
               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
               I=simplify(lead(std(phi(i))),1);
               attrib(I,"isSB",1);
               nesttest=is_nested(I);
               if ( nesttest == 1 )
                 {
                   break;
                 }
             }
           if ( nesttest == 0 )
             {
       "// WARNING from proc satiety from lib mregular.lib:
// The procedure has entered in 30 loops and could not find a monomial
// ideal of nested type with the same satiety as your ideal: in your
// example the method using random changes of coordinates may enter an
// infinite loop when the field is finite.
// Try removing the optional argument.";
       return (-1);
             }
         }
//
// At this stage, we have obtained a monomial ideal I of nested type
// such that depth(S/i)=depth(S/I). We now compute depth(I).
//
//----- When S/i is Cohen-Macaulay:
//
       for ( ii = n-d+2; ii <= n+1; ii++ )
         {
           K=K+select(I,ii);
         }
       if ( size(K) == 0 )
         {
           time=rtimer-time;
           // Additional information:
           dbprint(printlevel-voice+2,
                   "// Dimension of S/i: "+string(d));
           dbprint(printlevel-voice+2,
                   "// Time for computing satiety: " + string(time) + " sec.");
           dbprint(printlevel-voice+2,
                   "// The satiety of i is:");
           return(0);
         }
//----- When d=1 (and S/i is not Cohen-Macaulay) ==> depth =0:
       if ( d == 1 )
         {
           KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
           sat=maxdeg1(KK)+1;
           time=rtimer-time;
           // Additional information:
           dbprint(printlevel-voice+2,
                   "// Dimension of S/i: 1");
           dbprint(printlevel-voice+2,
                   "// Time for computing satiety: "+ string(time) + " sec.");
           dbprint(printlevel-voice+2,
                   "// The satiety of i is:");
           return(sat);
         }
//----- Now d>1 and S/i is not Cohen-Macaulay:
//
//----- First, determine the last variable really occuring
       lastv=n-d;
       h=n;
       while ( lastv == n-d and h > n-d )
         {
           K=select(I,h+1);
           if ( size(K) == 0 )
             {
               h=h-1;
             }
           else
             {
               lastv=h;
             }
         }
//----- and compute the satiety:
       sat=0;
       if ( lastv == n )
         {
           KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
           sat=maxdeg1(KK)+1;
         }
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: "+string(d));
       dbprint(printlevel-voice+2,
               "// Time for computing satiety: "+ string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The satiety of i is:");
       return(sat);
     }
//---- If no option: direct computation
   sat=maxdeg1(reduce(quotient(i,maxideal(1)),sbi))+1;
   time=rtimer-time;
       // Additional information:
   dbprint(printlevel-voice+2,
           "// Dimension of S/i: "+string(d)+";");
   dbprint(printlevel-voice+2,
           "// Time for computing satiety: "+ string(time) + " sec.");
   dbprint(printlevel-voice+2,
           "// The satiety of i is:");
   return(sat);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z,t,w),dp;
   ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
   satiety(i);
   ideal I=lead(std(i));
   satiety(I);   // First  method: direct computation
   satiety(I,1); // Second method: doing changes of coordinates
// Additional information is displayed if you change printlevel (=1);
}
////////////////////////////////////////////////////////////////////////////////
/*
Out-commented examples:
   ring s1=0,(x,y,z,t),dp;
   ideal I=zt3,z2t2,yz2t,xz2t,xy2t,x3y;
   satiety(I);
   satiety(I,1);
// Another example:
   ring s2=0,(z,y,x),dp;
   ideal I=z38,z26y2,z14y4,z12x,z10x5,z8x9,z6x16,z4x23,z2y6,y32;
   satiety(I);
   satiety(I,1);
// One more:
   ring s3=0,(s,t,u,x(0..8)),dp;
   ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,
           x(4)-s2t7u6,x(5)-s7t7u,x(6)-s15,x(7)-t15,x(8)-u15;
   ideal ei=eliminate(i,stu);
   size(ei);
   ring s4=0,x(0..8),dp;
   ideal i=imap(s3,ei);
   ideal m=maxideal(1);
   m[8]=m[8]+m[7];
   map phi=m;
   ideal phii=phi(i);
   ideal nI=lead(std(phii));
   ring s5=0,x(0..7),dp;
   ideal nI=imap(s4,nI);
   satiety(nI);
   satiety(nI,1);
   ideal I1=subst(nI,x(7),1);
   ring s6=0,x(0..6),dp;
   ideal I1=imap(s5,I1);
   satiety(I1);
   satiety(I1,1);
   ideal I2=subst(I1,x(6),1);
   ring s7=0,x(0..5),dp;
   ideal I2=imap(s6,I2);
   satiety(I2);
   satiety(I2,1);
//
*/
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

proc regMonCurve (list #)
"
USAGE:   regMonCurve (a0,...,an) ; ai integers with a0=0 < a1 < ... < an=:d
RETURN:  an integer, the Castelnuovo-Mumford regularity of the projective
         monomial curve C in Pn(K) parametrically defined by
              x(0) = t^d , x(1) = s^(a1)t^(d-a1) , ..... , x(n) = s^d
         where K is the field of complex numbers.
         (returns -1 if a0=0 < a1 < ... < an is not satisfied)
ASSUME:  a0=0 < a1 < ... < an are integers.
NOTES:   1. The defining ideal of the curve C, I in S=K[x(0),...,x(n)], is
            determined by elimination.
         2. The procedure regIdeal has been improved in this case since one
            knows beforehand that the monomial ideal J=lead(std(I)) is of
            nested type if the monomial ordering is dp, and that
            reg(C)=reg(J) (see preprint 'Saturation and Castelnuovo-Mumford
            regularity' by Bermejo-Gimenez, 2004).
         3. If printlevel > 0 (default = 0) additional info is displayed:
            - It says whether C is arithmetically Cohen-Macaulay or not.
            - If C is not arith. Cohen-Macaulay, end(H^1(S/I)) is computed
              and an upper bound for the a-invariant of S/I is given.
            - It also determines one step of the minimal graded free
              resolution (m.g.f.r.) of I where the regularity is attained
              and gives the value of the regularity of the Hilbert function
              of S/I when reg(I) is attained at the last step of a m.g.f.r.
EXAMPLE: example regMonCurve; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
   int ii,H,h,hh,time,ttime,firstind,lastind;
   int n = size(#)-1;
//------------------  Check assumptions on integers  -------------------------
   if ( #[1] != 0 )
   {"// WARNING from proc regMonCurve from lib mregular.lib:
// USAGE: your input must be a list of integers a0,a1,...,an such that
// a0=0 < a1 < a2 < ... < an";
      return(-1);
   }
   for ( ii=1; ii<= n; ii++ )
   {
      if ( #[ii] >= #[ii+1] )
      {
      "// WARNING from proc regMonCurve from lib mregular.lib:
// USAGE: your input must be a list of integers a0,a1,...,an such that
// a0=0 < a1 < a2 < ... < an";
      return(-1);
      }
   }
   ring R=0,(x(0..n),s,t),dp;
   ideal param,m,i;
   poly f(0..n);
   for (ii=0;ii<=n;ii++)
      {
      f(ii)=s^(#[n+1]-#[ii+1])*t^(#[ii+1]);
      param=param+f(ii);
      }
   m=subst(maxideal(1),s,0);
   m=simplify(subst(m,t,0),2);
   i=matrix(m)-matrix(param);
   ttime=rtimer;
   i=eliminate(i,st);
   ring r=0,(x(1..n),x(0)),dp;
   ideal i,I;
   i=imap(R,i);
   I=minbase(lead(std(i)));
   attrib(I,"isSB",1);
   ttime=rtimer-ttime;
   time=rtimer;
   ring nr=0,x(1..n),dp;
   ideal I,K,KK,J;
   I=imap(r,I);
   attrib(I,"isSB",1);
   K=select(I,n);
//------------------ Cohen-Macaulay case ------------
   if ( size(K) == 0 )
     {
       ring mr=0,x(1..n-1),dp;
       ideal I=imap(nr,I);
       H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// The sequence of integers defines a monomial curve C in P"
               + string(n));
       dbprint(printlevel-voice+2,
               "//    C is arithmetically Cohen-Macaulay");
       dbprint(printlevel-voice+2,
               "//    Regularity attained at the last step of a m.g.f.r. of I(C)");
       dbprint(printlevel-voice+2,
               "//    Regularity of the Hilbert function of S/I(C): "
               + string(H-2));
       dbprint(printlevel-voice+2,
               "//    Time for computing ideal I(C) (by elimination): "
               + string(ttime) + " sec.");
       dbprint(printlevel-voice+2,
               "//    Time for computing reg(C) once I(C) has been determined: "
               + string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The Castelnuovo-Mumford regularity of C is:");
       return(H);
     }
   else
     {
       KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
       firstind=maxdeg1(KK)+1;
       J=subst(I,x(n),1);
       ring mr=0,x(1..n-1),dp;
       ideal J=imap(nr,J);
       lastind=maxdeg1(minbase(quotient(J,maxideal(1))))+1;
       H=firstind;
       if ( lastind > H )
         {
           H=lastind;
         }
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// The sequence of integers defines a monomial curve C in P"
               + string(n));
       dbprint(printlevel-voice+2,
               "//    C is not arithmetically Cohen-Macaulay");
       dbprint(printlevel-voice+2,
               "//    end(H^1(S/I(C))) = "
               +string(firstind-2));
       dbprint(printlevel-voice+2,
               "//    Upper bound for the a-invariant of S/I(C): end(H^2(S/I(C))) <= "
               +string(lastind-3));
       if ( H == firstind )
         {
           dbprint(printlevel-voice+2,
                   "//    Regularity attained at the last step of a m.g.f.r. of I(C)");
           dbprint(printlevel-voice+2,
                   "//    Regularity of the Hilbert function of S/I(C): "
                   + string(H-1));
         }
       else
         {
           dbprint(printlevel-voice+2,
                   "//    Regularity attained at the second last step of a m.g.f.r. of I(C)");
           dbprint(printlevel-voice+2,
                   "//    (and not attained at the last step)");
         }
       dbprint(printlevel-voice+2,
               "//    Time for computing ideal I(C) (by elimination): "
               + string(ttime) + " sec.");
       dbprint(printlevel-voice+2,
               "//    Time for computing reg(C) once I(C) has been determined: "
               + string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// The Castelnuovo-Mumford regularity of C is:");
       return(H);
     }
}
example
{ "EXAMPLE:"; echo = 2;
// The 1st example is the twisted cubic:
   regMonCurve(0,1,2,3);
// The 2nd. example is the non arithm. Cohen-Macaulay monomial curve in P4
// parametrized by: x(0)-s6,x(1)-s5t,x(2)-s3t3,x(3)-st5,x(4)-t6:
   regMonCurve(0,1,3,5,6);
// Additional information is displayed if you change printlevel (=1);
}
////////////////////////////////////////////////////////////////////////////////
/*
Out-commented examples:
//
// The sequence of integers must be strictly increasing
   regMonCurve(1,4,6,9);
   regMonCurve(0,3,8,5,23);
   regMonCurve(0,4,7,7,9);
//
// A curve in P3 s.t. the regularity is attained at the last step:
   regMonCurve(0,2,12,15);
//
// A curve in P4 s.t. the regularity attained at the last but one
// but NOT at the last step (Ex. 3.3 Preprint 2004):
   regMonCurve(0,5,9,11,20);
//
// A curve in P8 s.t. the m.g.f.r. of the defining ideal is not easily
// obtained through m.g.f.r.:
   regMonCurve(0,1,2,3,9,11,18,24,25);
//
// A curve in P11 of degree 37:
   regMonCurve(0,1,2,7,16,17,25,27,28,30,36,37);
// It takes some time to compute the eliminated ideal; the computation of
// the regularity is then rather fast as one can check using proc regIdeal:
   ring q=0,(s,t,x(0..11)),dp;
   ideal i=x(0)-st36,x(1)-s2t35,x(2)-s7t30,x(3)-s16t21,x(4)-s17t20,x(5)-s25t12,
           x(6)-s27t10,x(7)-s28t9,x(8)-s30t7,x(9)-s36t,x(10)-s37,x(11)-t37;
   ideal ei=eliminate(i,st);
   ring qq=0,x(0..11),dp;
   ideal i=imap(q,ei);
   regIdeal(i);
//
// A curve in P14 of degree 55:
   regMonCurve(0,1,2,7,16,17,25,27,28,30,36,37,40,53,55);
//
*/
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

proc NoetherPosition (ideal i)
"
USAGE:   NoetherPosition (i); i ideal
RETURN:  ideal such that, for the homogeneous linear transformation
         map phi=S,NoetherPosition(i);
         one has that K[x(n-d+1),...,x(n)] is a Noether normalization of
         S/phi(i) where S=K[x(0),...x(n)] is the basering and d=dim(S/i).
         (returns -1 if i = (0) or (1)).
ASSUME:  The field K is infinite and i is a nonzero proper ideal.
NOTE:    1. It works also if K is a finite field if it terminates, but
            may result in an infinite loop. If the procedure enters more
            than 30 loops, -1 is returned and a warning message is displayed.@*
         2. If printlevel > 0 (default = 0), additional info is displayed:
            dim(S/i) and K[x(n-d+1),...,x(n)] are given.
EXAMPLE: example NoetherPosition; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
   int ii,jj,d,time,nl,NPtest;
   intmat ran;
   def r0 = basering;
   ideal K,chcoord;
   int n = nvars(r0)-1;
   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
   execute(s);
   ideal i,sbi,I,K,chcoord,m;
   poly P;
   map phi;
   i = fetch(r0,i);
   time=rtimer;
   sbi=std(i);
   I=simplify(lead(sbi),1);
   attrib(I,"isSB",1);
   d=dim(I);
//----- If the ideal i is not proper:
   if ( d == -1 )
     {
       "// WARNING from proc NoetherPosition from lib mregular.lib:
// The ideal i is (1)!";
       return (-1);
     }
//----- If the ideal i is 0:
   if ( size(I) == 0 )
     {
       "// WARNING from proc NoetherPosition from lib mregular.lib:
// The ideal i is (0)!";
       return (-1);
     }
//----- When the ideal i is 0-dimensional:
   if ( d == 0 )
     {
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: 0");
       dbprint(printlevel-voice+2,
               "// Time for computing a Noether normalization: "
               + string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// K is a Noether normalization of S/phi(i)");
       dbprint(printlevel-voice+2,
               "// where the map phi: S --> S is:");
       setring r0;
       return (maxideal(1));
     }
   NPtest=is_NP(I);
   if ( NPtest == 1 )
     {
       K=x(n-d+1..n);
       setring r0;
       K=fetch(r1,K);
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: " + string(d) );
       dbprint(printlevel-voice+2,
               "// Time for computing a Noether normalization: " +
               string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// K[" + string(K) +
               "] is a Noether normalization of S/phi(i)");
       dbprint(printlevel-voice+2,
               "// where the map phi: S --> S is:");
       return (maxideal(1));
     }
//---- Otherwise, random change of coordinates and
//---- test for Noether normalization.
//---- If we were unlucky, another change of coord. will be done:
   while ( nl < 30 )
   {
     chcoord=select1(maxideal(1),1..(n-d+1));
     nl=nl+1;
     for ( ii = 1; ii<=d; ii++ )
     {
       ran=random(100,1,n-d+ii);
       ran=intmat(ran,1,n-d+ii+1);
       ran[1,n-d+ii+1]=1;
       m=select1(maxideal(1),1..(n-d+1+ii));
       for ( jj = 1; jj<=n-d+ii+1; jj++ )
       {
       P=P+ran[1,jj]*m[jj];
       }
       chcoord[n-d+1+ii]=P;
       P=0;
     }
     phi=r1,chcoord;
     dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     I=simplify(lead(std(phi(i))),1);
     attrib(I,"isSB",1);
     NPtest=is_NP(I);
     if ( NPtest == 1 )
       {
       K=x(n-d+1..n);
       setring r0;
       K=fetch(r1,K);
       chcoord=fetch(r1,chcoord);
       time=rtimer-time;
       // Additional information:
       dbprint(printlevel-voice+2,
               "// Dimension of S/i: " + string(d) );
       dbprint(printlevel-voice+2,
               "// Time for computing a Noether normalization: " +
               string(time) + " sec.");
       dbprint(printlevel-voice+2,
               "// K[" + string(K) +
               "] is a Noether normalization of S/phi(i)");
       dbprint(printlevel-voice+2,
               "// where the map phi: S --> S is:");
       return (chcoord);
       }
   }
       "// WARNING from proc NoetherPosition from lib mregular.lib:
// The procedure has entered in more than 30 loops: in your example
// the method may enter an infinite loop over a finite field!";
       return (-1);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z,t,u),dp;
   ideal i1=y,z,t,u; ideal i2=x,z,t,u; ideal i3=x,y,t,u; ideal i4=x,y,z,u;
   ideal i5=x,y,z,t; ideal i=intersect(i1,i2,i3,i4,i5);
   map phi=r,NoetherPosition(i);
   phi;
   ring r5=5,(x,y,z,t,u),dp;
   ideal i=imap(r,i);
   map phi=r5,NoetherPosition(i);
   phi;
// Additional information is displayed if you change printlevel (=1);
}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

proc is_NP (ideal i)
"
USAGE:   is_NP (i); i ideal
RETURN:  1  if K[x(n-d+1),...,x(n)] is a Noether normalization of
            S/i where S=K[x(0),...x(n)] is the basering, and d=dim(S/i),
         0  otherwise.
         (returns -1 if i=(0) or i=(1)).
ASSUME:  i is a nonzero proper homogeneous ideal.
NOTE:    1. If i is not homogeneous and is_NP(i)=1 then K[x(n-d+1),...,x(n)]
            is a Noether normalization of S/i. The converse may be wrong if
            the ideal is not homogeneous.
         2. is_NP is used in the procedures regIdeal, depthIdeal, satiety,
            and NoetherPosition.
EXAMPLE: example is_NP; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
   int ii,d,dz;
   def r0 = basering;
   int n = nvars(r0)-1;
   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
   execute(s);
   ideal i,sbi,I,J;
   i = fetch(r0,i);
   sbi=std(i);
   I=simplify(lead(sbi),1);
   attrib(I,"isSB",1);
   d=dim(I);
//----- If the ideal i is not proper:
   if ( d == -1 )
     {
       "// WARNING from proc is_NP from lib mregular.lib:
// The ideal i is (1)!";
       return (-1);
     }
//----- If the ideal i is 0:
   if ( size(I) == 0 )
     {
       "// WARNING from proc is_NP from lib mregular.lib:
// The ideal i is (0)!";
       return (-1);
     }
//----- When the ideal i is 0-dimensional:
   if ( d == 0 )
     {
       return (1);
     }
//----- Check Noether position
   J=I;
   for ( ii = n-d+1; ii <= n; ii++ )
     {
       J=subst(J,x(ii),0);
     }
   attrib(J,"isSB",1);
   dz=dim(J);
   if ( dz == d )
     {
       return (1);
     }
   else
     {
       return(0);
     }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z,t,u),dp;
   ideal i1=y,z,t,u; ideal i2=x,z,t,u; ideal i3=x,y,t,u; ideal i4=x,y,z,u;
   ideal i5=x,y,z,t; ideal i=intersect(i1,i2,i3,i4,i5);
   is_NP(i);
   ideal ch=x,y,z,t,x+y+z+t+u;
   map phi=ch;
   is_NP(phi(i));
}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

proc is_nested (ideal i)
"
USAGE:   is_nested (i); i monomial ideal
RETURN:  1 if i is of nested type, 0 otherwise.
         (returns -1 if i=(0) or i=(1)).
ASSUME:  i is a nonzero proper monomial ideal.
NOTES:   1. The ideal must be monomial, otherwise the result has no meaning
            (so check this before using this procedure).@*
         2. is_nested is used in procedures depthIdeal, regIdeal and satiety.@*
         3. When i is a monomial ideal of nested type of S=K[x(0)..x(n)],
            the a-invariant of S/i coincides with the upper bound obtained
            using the procedure regIdeal with printlevel > 0.
THEORY:  A monomial ideal is of nested type if its associated primes are all
         of the form (x(0),...,x(i)) for some i<=n.
         (see definition and effective criterion to check this property in
         the preprint 'Saturation and Castelnuovo-Mumford regularity' by
         Bermejo-Gimenez, 2004).
EXAMPLE: example is_nested; shows some examples
"
{
//--------------------------- initialisation ---------------------------------
   int ii,d,tev,lastv,h,NPtest;
   def r0 = basering;
   int n = nvars(r0)-1;
   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
   execute(s);
   ideal I,K,KK,LL;
   I = fetch(r0,i);
   I=minbase(I);
   attrib(I,"isSB",1);
   d=dim(I);
//----- If the ideal i is not proper:
   if ( d == -1 )
   {
       "// WARNING from proc is_nested from lib mregular.lib:
// The ideal i is (1)!";
       return (-1);
   }
//----- If the ideal i is 0:
   if ( size(I) == 0 )
     {
       "// WARNING from proc is_nested from lib mregular.lib:
// The ideal i is (0)!";
       return (-1);
     }
//----- When the ideal i is 0-dimensional:
   if ( d == 0 )
     {
       return (1);
     }
//----- Check Noether position
   NPtest=is_NP(I);
   if ( NPtest != 1 )
   {
       return (0);
   }
//----- When ideal is 1-dim. + var. in Noether position -> Nested Type
   if ( d == 1 )
     {
       return (1);
     }
//----- Determ. of the last variable really occuring
   lastv=n-d;
   h=n;
   while ( lastv == n-d and h > n-d )
     {
       K=select(I,h+1);
       if ( size(K) == 0 )
         {
           h=h-1;
         }
       else
         {
           lastv=h;
         }
     }
//----- Check the second property by evaluation when NP + d>1
   KK=subst(I,x(lastv),1);
   for ( ii = n-lastv; ii<=d-2; ii++ )
   {
     LL=minbase(subst(I,x(n-ii-1),1));
     attrib(LL,"isSB",1);
     tev=size(reduce(KK,LL));
     if ( tev > 0 )
     {
       return(0);
     }
     KK=LL;
   }
   return(1);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s=0,(x,y,z,t),dp;
   ideal i1=x2,y3; ideal i2=x3,y2,z2; ideal i3=x3,y2,t2;
   ideal i=intersect(i1,i2,i3);
   is_nested(i);
   ideal ch=x,y,z,z+t;
   map phi=ch;
   ideal I=lead(std(phi(i)));
   is_nested(I);
}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////