source: git/Singular/LIB/noether.lib

/////////////////////////////////////////////////////////////////////////////
version="version noether.lib 4.1.2.0 Feb_2019 "; // $Id$
category="Commutative Algebra";
info="
LIBRARY: noether.lib   Noether normalization of an ideal (not necessary
                       homogeneous)
AUTHORS: A. Hashemi,  Amir.Hashemi@lip6.fr


OVERVIEW:
A library for computing the Noether normalization of an ideal that DOES NOT
require the computation of the dimension of the ideal.
It checks whether an ideal is in Noether position.  A modular version of
these algorithms is also provided.
The procedures are based on a paper of Amir Hashemi 'Efficient Algorithms for
Computing Noether Normalization' (presented in ASCM 2007)

This library computes also Castelnuovo-Mumford regularity and satiety of an
ideal.  A modular version of these algorithms is also provided.
The procedures are based on a paper of Amir Hashemi 'Computation of
Castelnuovo-Mumford regularity and satiety' (preprint 2008)


PROCEDURES:
 NPos_test(id);  checks whether monomial ideal id is in Noether position
 modNpos_test(id); the same as above using modular methods
 NPos(id);       Noether normalization of ideal id
 modNPos(id);      Noether normalization of ideal id by modular methods
 nsatiety(id); Satiety of ideal id
 modsatiety(id)  Satiety of ideal id by modular methods
 regCM(id);    Castelnuovo-Mumford regularity of ideal id
 modregCM(id); Castelnuovo-Mumford regularity of ideal id by modular methods
";
LIB "elim.lib";
LIB "algebra.lib";
LIB "polylib.lib";
LIB "ring.lib";
LIB "presolve.lib";

///////////////////////////////////////////////////////////////////////////////

proc NPos_test (ideal I)
"
USAGE:  NPos_test (I); I monomial ideal
RETURN: A list whose first element is 1, if i is in Noether position,
        0 otherwise. The second element of this list is a list of variables ordered
        such that those variables are listed first, of which a power belongs to the
        initial ideal of i. If i is in Noether position, the method returns furthermore
        the dimension of i.
ASSUME: i is a nonzero monomial ideal.
"
{
//--------------------------- initialisation ---------------------------------
   int  time,ii,j,k,l,d,t,jj;
   intvec v;
   def r0 = basering;
   int n = nvars(r0)-1;
   list L,Y,P1,P2,P3;
   if (I[1]==1)
   {
     print("The ideal is 1");return(1);
   }
   for ( ii = 1; ii <= n+1; ii++ )
   {
     L[ii]=0;
   }
   for ( ii = 1; ii <= size(I); ii++ )
   {
     Y=variables(I[ii]);
     l=rvar(Y[1][1]);
     if (size(Y[1])==1)
     {
       L[l]=1;
       P1=insert(P1,Y[1][1]);
     }
     if (L[l]==0)
     {
       L[l]=-1;
     }
   }
   t=size(P1);
   if (t==0)
   {
     for ( jj = 1; jj <= n+1; jj++ )
     {
       P3=insert(P3,varstr(jj));
     }
   }
   else
   {
     P2=findvars(ideal(P1[1..t]))[3];
     for ( jj = 1; jj <= size(P2[1]); jj++ )
     {
       P3=insert(P3,P2[1][jj]);
     }
   }
   if (L[n+1]==-1)
   {
     return(list(0,P1+P3));
   }
   for ( ii = 1; ii <= n; ii++ )
   {
     if (L[ii]==-1)
     {
       return(list(0,P1+P3));
     }
     if (L[ii]==0 and L[ii+1]==1)
     {
       return(list(0,P1+P3));
     }
   }
   d=n+1-sum(L);
   print("The dimension of the ideal is:");print(d);
   return(list(1,P1+P3));
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(X,Y,a,b),dp;
poly f=X^8+a*Y^4-Y;
poly g=Y^8+b*X^4-X;
poly h=diff(f,X)*diff(g,Y)-diff(f,Y)*diff(g,X);
ideal i=f,g,h;
NPos_test(i);
}
//////////////////////////////////////////
proc modNpos_test (ideal i)
"USAGE: modNpos_test(i); i an ideal
RETURN: 1 if i is in Noether position 0  otherwise.
NOTE:   This test is a probabilistic test, and it computes the initial of the ideal modulo the prime number 2147483647 (the biggest prime less than 2^31).
"
{
  "// WARNING:
// The procedure is probabilistic and  it computes the initial of the ideal modulo the prime number 2147483647";
  int p;
  def br=basering;
  setring br;
  ideal I;
  list #;
  option(redSB);
  p=2147483647;
  #=ringlist(br);
  #[1]=p;
  def oro=ring(#);
  setring oro;
  ideal sbi,lsbi;
  sbi=fetch(br,i);
  lsbi=lead(std(sbi));
  setring br;
  I=fetch(oro,lsbi);
  I=simplify(I,1);
  attrib(I,"isSB",1);
  return(NPos_test(I));
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(X,Y,a,b),dp;
poly f=X^8+a*Y^4-Y;
poly g=Y^8+b*X^4-X;
poly h=diff(f,X)*diff(g,Y)-diff(f,Y)*diff(g,X);
ideal i=f,g,h;
modNpos_test(i);
}


///////////////////////////////////////////////////////////////////////////////
proc NPos (ideal i)
"USAGE:  NPos(i); i ideal
RETURN:  A linear map phi such that  phi(i) is in Noether position
KEYWORDS: Noether position
"
{
//--------------------------- initialisation ---------------------------------
  int ii,jj,d,time,n,nl,zz;
  intmat ran;
  def r0 = basering;
  ideal K,chcoord;
  n = nvars(r0)-1;
  list l1;
  for (zz = 0; zz<= n; zz++)
  {
   l1[zz+1] = "x("+string(zz)+")";
  }
  ring r1 = create_ring(ring_list(r0)[1], l1, "dp", "no_minpoly");
  ideal i,sbi,I,K,chcoord,m,L;
  list #;
  poly P;
  map phi;
  i = fetch(r0,i);
  time=rtimer;
  system("--ticks-per-sec",10);
  i=std(i);
  sbi=sort(lead(i))[1];
  #=NPos_test(sbi);
  if ( #[1]== 1 )
  {
    return ("The ideal is in Noether position and the time of this computation is:",rtimer-time,"/10 sec.");
  }
  else
  {
    L=maxideal(1);
    chcoord=maxideal(1);
    for ( ii = 1; ii<=n+1; ii++ )
    {
      chcoord[rvar(#[2][ii])]=L[ii];
    }
    phi=r1,chcoord;
    sbi=phi(sbi);
    if ( NPos_test(sbi)[1] == 1 )
    {
      setring r0;
      chcoord=fetch(r1,chcoord);
      return (chcoord,"and the time of this computation is:",rtimer-time,"/10 sec.");
    }
  }
  while ( nl < 30 )
  {
    nl=nl+1;
    I=i;
    L=maxideal(1);
    for ( ii = n; ii>=0; ii-- )
    {
      chcoord=select1(maxideal(1),1..ii);
      ran=random(100,1,ii);
      ran=intmat(ran,1,ii+1);
      ran[1,ii+1]=1;
      m=select1(maxideal(1),1..(ii+1));
      for ( jj = 1; jj<=ii+1; jj++ )
      {
        P=P+ran[1,jj]*m[jj];
      }
      chcoord[ii+1]=P;
      L[ii+1]=P;
      P=0;
      phi=r1,chcoord;
      I=phi(I);
      if ( NPos_test(sort(lead(std(I)))[1])[1] == 1 )
      {
        K=x(ii..n);
        setring r0;
        K=fetch(r1,K);
        ideal L=fetch(r1,L);
        return (L,"and the time of this computation is:",rtimer-time,"/10 sec.");
      }
    }
  }
  "// WARNING:
// 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,a,b),dp;
poly f=X^8+a*Y^4-Y;
poly g=Y^8+b*X^4-X;
poly h=diff(f,X)*diff(g,Y)-diff(f,Y)*diff(g,X);
ideal i=f,g,h;
NPos(i);
}
///////////////////////////////////////////////////////////////////////////////
proc modNPos (ideal i)
"USAGE:  modNPos(i); i ideal
RETURN:  A linear map phi such that  phi(i) is in Noether position
NOTE:    It uses the procedure  modNpos_test to test Noether position.
"
{
//--------------------------- initialisation ---------------------------------
   int ii,jj,d,time,n,nl,zz;
   intmat ran;
   def r0 = basering;
   ideal K,chcoord;
   n = nvars(r0)-1;
   list l2;
   for (zz = 0; zz<= n; zz++)
   {
     l2[zz+1] = "x("+string(zz)+")";
   }
   ring r1 = create_ring(ring_list(r0)[1], l2, "dp", "no_minpoly");
   ideal i,sbi,I,K,chcoord,m,L;
   poly P;
   list #;
   map phi;
   i = fetch(r0,i);
   time=rtimer;
   system("--ticks-per-sec",10);
   #=modNpos_test(i);
   if ( #[1]== 1 )
   {
     return ("The ideal is in Noether position and the time of this computation is:",rtimer-time,"/10 sec.");
   }
   else
   {
     L=maxideal(1);
     chcoord=maxideal(1);
     for ( ii = 1; ii<=n+1; ii++ )
     {
       chcoord[rvar(#[2][ii])]=L[ii];
     }
     phi=r1,chcoord;
     I=phi(i);
     if ( modNpos_test(I)[1] == 1 )
     {
       setring r0;
       chcoord=fetch(r1,chcoord);
       return (chcoord,"and the time of this computation is:",rtimer-time,"/10 sec.");
     }
   }
   while ( nl < 30 )
   {
     nl=nl+1;
     I=i;
     L=maxideal(1);
     for ( ii = n; ii>=0; ii-- )
     {
       chcoord=select1(maxideal(1),1..ii);
       ran=random(100,1,ii);
       ran=intmat(ran,1,ii+1);
       ran[1,ii+1]=1;
       m=select1(maxideal(1),1..(ii+1));
       for ( jj = 1; jj<=ii+1; jj++ )
       {
         P=P+ran[1,jj]*m[jj];
       }
       chcoord[ii+1]=P;
       L[ii+1]=P;
       P=0;
       phi=r1,chcoord;
       I=phi(I);
       if ( modNpos_test(I)[1] == 1 )
       {
         K=x(ii..n);
         setring r0;
         K=fetch(r1,K);
         ideal L=fetch(r1,L);
         return (L,"and the time of this computation is:",rtimer-time,"/10 sec.");
       }
     }
   }
   "// WARNING:
// 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,a,b),dp;
poly f=X^8+a*Y^4-Y;
poly g=Y^8+b*X^4-X;
poly h=diff(f,X)*diff(g,Y)-diff(f,Y)*diff(g,X);
ideal i=f,g,h;
modNPos(i);
}

////////////////////////////////////////////////////////////////////////////////////
static proc TestLastVarIsInGenericPos (ideal i)
"USAGE:   TestLastVarIsInGenericPos (i); i a monomial ideal,
RETURN:  1 if the last variable is in generic position for i and 0 otherwise.
THEORY:  The last variable is in generic position if the quotient of the ideal
         with respect to this variable is equal to the quotient of the ideal with respect to the maximal ideal.
"
{
//--------------------------- initialisation ---------------------------------
  int n,ret,zz;
  def r0 = basering;
  n = nvars(r0)-1;
  list l3;
  for (zz = 0; zz<= n; zz++)
  {
    l3[zz+1] = "x("+string(zz)+")";
  }
  ring r1 = create_ring(ring_list(r0)[1], l3, "dp", "no_minpoly");
  ideal I,i;
  i = fetch(r0,i);
  attrib(i,"isSB",1);
  I=quotient(select(i,n+1),x(n));
  I=I*maxideal(1);
  ret=1;
  if (size(reduce(I,i,5)) <> 0)
  {
    ret=0;
  }
  return(ret);
}

////////////////////////////////////////////////////////////////////////////////////
proc nsatiety (ideal i)
"USAGE:   nsatiety (i); i ideal,
RETURN:  an integer, the satiety of i.
         (returns -1 if i is not homogeneous)
ASSUME:  i is a homogeneous ideal of the basering R=K[x(0)..x(n)].
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)) coincide.
KEYWORDS: saturation
"
{
//--------------------------- initialisation ---------------------------------
  int e,ii,jj,h,d,time,lastv,nl,ret,zz;
  intmat ran;
  def r0 = basering;
  int n = nvars(r0)-1;
  list l4;
  for (zz = 0; zz<= n; zz++)
  {
    l4[zz+1] = "x("+string(zz)+")";
  }
  ring r1 = create_ring(ring_list(r0)[1], l4, "dp", "no_minpoly");
  ideal i,sbi,I,K,chcoord,m,L;
  poly P;
  map phi;
  i = fetch(r0,i);
  time=rtimer;
  system("--ticks-per-sec",100);
  sbi=std(i);
//----- Check ideal homogeneous
  if ( homog(sbi) == 0 )
  {
    dbprint(2,"The ideal is not homogeneous, and time for this test is: " + string(rtimer-time) + "/100sec.");
    return ();
  }
  I=simplify(lead(sbi),1);
  attrib(I,"isSB",1);
  K=select(I,n+1);
  if (size(K) == 0)
  {
    dbprint(2,"sat(i)=0 and the time of this computation: " + string(rtimer-time) + "/100sec.");
    return();
  }
  if (TestLastVarIsInGenericPos(I) == 1 )
  {
    dbprint(2,"sat(i)=" + string(maxdeg1(K)) + " and the time of this computation: " + string(rtimer-time) + "/100sec.");
    return();
  }
  while ( nl < 5 )
  {
    nl=nl+1;
    chcoord=select1(maxideal(1),1..n);
    ran=random(100,1,n);
    ran=intmat(ran,1,n+1);
    ran[1,n+1]=1;
    m=select1(maxideal(1),1..(n+1));
    for ( jj = 1; jj<=n+1; jj++ )
    {
      P=P+ran[1,jj]*m[jj];
    }
    chcoord[n+1]=P;
    P=0;
    phi=r1,chcoord;
    L=std(phi(i));
    I=simplify(lead(L),1);
    attrib(I,"isSB",1);
    K=select(I,n+1);
    if (size(K) == 0)
    {
      dbprint(2,"sat(i)=0 and the time of this computation: " + string(rtimer-time) + "/100sec.");
      return();
    }
    if (TestLastVarIsInGenericPos(I) == 1 )
    {
      dbprint(2,"sat(i)=" + string(maxdeg1(K)) + " and the time of this computation: " + string(rtimer-time) + "/100sec.");
      return();
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(t,a,b,c,d),dp;
ideal i=b4-a3d, ab3-a3c, bc4-ac3d-bcd3+ad4, c6-bc3d2-c3d3+bd5, ac5-b2c3d-ac2d3+b2d4, a2c4-a3d3+b3d3-a2cd3, b3c3-a3d3, ab2c3-a3cd2+b3cd2-ab2d3, a2bc3-a3c2d+b3c2d-a2bd3, a3c3-a3bd2, a4c2-a3b2d;
nsatiety(i);
}


//////////////////////////////////////////////////////////////////////////////
proc modsatiety (ideal i)
"USAGE:   modsatiety(i); i ideal,
RETURN:  an integer, the satiety of i.
         (returns -1 if i is not homogeneous)
ASSUME:  i is a homogeneous ideal of the basering R=K[x(0)..x(n)].
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)) coincide.
NOTE:    This is a probabilistic procedure, and it computes the initial of the ideal modulo the prime number 2147483647 (the biggest prime less than 2^31).
KEYWORDS: saturation
"
{
//--------------------------- initialisation ---------------------------------
  "// WARNING: The characteristic of base field must be zero.
// The procedure is probabilistic and  it computes the
//initial ideals modulo the prime number 2147483647.";
  int e,ii,jj,h,d,time,lastv,nl,ret,s1,d1,siz,j,si,u,k,p,zz;
  intvec v1;
  intmat ran;
  def r0 = basering;
  int n = nvars(r0)-1;
  list l5;
  for (zz = 0; zz<= n; zz++)
  {
    l5[zz+1] = "x("+string(zz)+")";
  }
  ring r1 = create_ring(ring_list(r0)[1], l5, "dp", "no_minpoly");
  ideal i,sbi,I,K,chcoord,m,L,sbi1,lsbi1,id1;
  vector V1;
  list #,LL,PL,Gb1,VGb1,Gb2,VGb2,Res1,Res2;
  poly P;
  map phi;
  time=rtimer;
  system("--ticks-per-sec",100);
  i = fetch(r0,i);
//----- Check ideal homogeneous
  if ( homog(i) == 0 )
  {
    "// WARNING: The ideal is not homogeneous.";
    dbprint(2,"Time for this test is: " + string(rtimer-time) + "/100sec.");
    return ();
  }
  option(redSB);
  p=2147483647;
  list r2=ringlist(r1);
  r2[1]=p;
  def oro=ring(r2);
  setring oro;
  ideal sbi=fetch(r1,i);
  sbi=std(sbi);
  setring r1;
  sbi=fetch(oro,sbi);
  kill oro;
  I=simplify(lead(sbi),1);
  attrib(I,"isSB",1);
  K=select(I,n+1);
  if (size(K) == 0)
  {
    dbprint(2,"msat(i)=0 and the time of this computation: " + string(rtimer-time) + "/100sec.");
    return();
  }
  if (TestLastVarIsInGenericPos(I) == 1 )
  {
    dbprint(2,"msat(i)=" + string(maxdeg1(K)) + " and the time of this computation: " + string(rtimer-time) + "/100sec.");
    return();
  }
  while ( nl < 30 )
  {
    nl=nl+1;
    chcoord=select1(maxideal(1),1..n);
    ran=random(100,1,n);
    ran=intmat(ran,1,n+1);
    ran[1,n+1]=1;
    m=select1(maxideal(1),1..(n+1));
    for ( jj = 1; jj<=n+1; jj++ )
    {
      P=P+ran[1,jj]*m[jj];
    }
    chcoord[n+1]=P;
    P=0;
    phi=r1,chcoord;
    sbi=phi(i);
    list r2=ringlist(r1);
    r2[1]=p;
    def oro=ring(r2);
    setring oro;
    ideal sbi=fetch(r1,sbi);
    sbi=std(sbi);
    setring r1;
    sbi=fetch(oro,sbi);
    kill oro;
    lsbi1=lead(sbi);
    attrib(lsbi1,"isSB",1);
    K=select(lsbi1,n+1);
    if (size(K) == 0)
    {
      dbprint(2,"msat(i)=0 and the time of this computation: " + string(rtimer-time) + "/100sec.");
      return();
    }
    if (TestLastVarIsInGenericPos(lsbi1) == 1 )
    {
      dbprint(2,"msat(i)=" + string(maxdeg1(K)) + " and the time of this computation: " + string(rtimer-time) + "/100sec.");
      return();
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(t,a,b,c,d),dp;
ideal i=b4-a3d, ab3-a3c, bc4-ac3d-bcd3+ad4, c6-bc3d2-c3d3+bd5, ac5-b2c3d-ac2d3+b2d4, a2c4-a3d3+b3d3-a2cd3, b3c3-a3d3, ab2c3-a3cd2+b3cd2-ab2d3, a2bc3-a3c2d+b3c2d-a2bd3, a3c3-a3bd2, a4c2-a3b2d;
modsatiety(i);
}

//////////////////////////////////////////////////////////////////////////////
//
proc regCM (ideal i)
"USAGE:  regCM (i); i ideal
RETURN:  the Castelnuovo-Mumford regularity of i.
         (returns -1 if i is not homogeneous)
ASSUME:  i is a homogeneous ideal.
KEYWORDS: regularity
"
{
//--------------------------- initialisation ---------------------------------
  int e,ii,jj,H,h,d,time,nl,zz;
  def r0 = basering;
  int n = nvars(r0)-1;
  list l6;
  for (zz = 0; zz<= n; zz++)
  {
    l6[zz+1] = "x("+string(zz)+")";
  }
  ring r1 = create_ring(ring_list(r0)[1], l6, "dp", "no_minpoly");
  ideal i,sbi,I,J,K,L;
  list #;
  poly P;
  map phi;
  i = fetch(r0,i);
  time=rtimer;
  system("--ticks-per-sec",100);
  sbi=std(i);
//----- Check ideal homogeneous
  if ( homog(sbi) == 0 )
  {
    "// The ideal is not homogeneous!";
    return (-1);
  }
  I=simplify(lead(sbi),1);
  attrib(I,"isSB",1);
  d=dim(I);
  if (char(r1) > 0 and d == 0)
  {
    def r2=changechar(0,r1);
    setring r2;
    ideal sbi,I,i,K,T;
    map phi;
    I = fetch(r1,I);
    i=I;
    attrib(I,"isSB",1);
  }
  else
  {
    def r2=changechar(ring_list(r1)[1],r1);
    setring r2;
    ideal sbi,I,i,K,T,ic,Ic;
    map phi;
    I = imap(r1,I);
    Ic=I;
    attrib(I,"isSB",1);
    i = imap(r1,i);
    ic=i;
  }
  K=select(I,n+1);
  if (size(K) == 0)
  {
    h=0;
  }
  else
  {
    if (TestLastVarIsInGenericPos(I) == 1)
    {
      h=maxdeg1(K);
    }
    else
    {
      while ( nl < 30 )
      {
        nl=nl+1;
        phi=r2,randomLast(100);
        T=phi(i);
        I=simplify(lead(std(T)),1);
        attrib(I,"isSB",1);
        K=select(I,n+1);
        if (size(K) == 0)
        {
          h=0;break;
        }
        if (TestLastVarIsInGenericPos(I) == 1 )
        {
          h=maxdeg1(K);break;
        }
      }
      i=T;
    }
  }
  list l7;
  for ( ii = n; ii>=n-d+1; ii-- )
  {
    l7=list();
    i=subst(i,x(ii),0);
    for (zz = 0; zz<= ii-1; zz++)
    {
      l7[zz+1] = "x("+string(zz)+")";
    }
    ring mr = create_ring(ring_list(r1)[1], l7, "dp", "no_minpoly");
    ideal i,sbi,I,J,K,L,T;
    poly P;
    map phi;
    i=imap(r2,i);
    I=simplify(lead(std(i)),1);
    attrib(I,"isSB",1);
    K=select(I,ii);
    if (size(K) == 0)
    {
      H=0;
    }
    else
    {
      if (TestLastVarIsInGenericPos(I) == 1)
      {
        H=maxdeg1(K);
      }
      else
      {
        while ( nl < 30 )
        {
          nl=nl+1;
          phi=mr,randomLast(100);
          T=phi(i);
          I=simplify(lead(std(T)),1);
          attrib(I,"isSB",1);
          K=select(I,ii);
          if (size(K) == 0)
          {
            H=0;break;
          }
          if (TestLastVarIsInGenericPos(I) == 1 )
          {
            H=maxdeg1(K);break;
          }
        }
        setring r2;
        i=imap(mr,T);
        kill mr;
      }
    }
    if (H > h)
    {
      h=H;
    }
  }
  if (nl < 30)
  {
    dbprint(2,"reg(i)=" + string(h) + " and the time of this computation: " + string(rtimer-time) + " sec./100");
    return();
  }
  else
  {
    I=Ic;
    attrib(I,"isSB",1);
    i=ic;
    K=subst(select(I,n+1),x(n),1);
    K=K*maxideal(maxdeg1(I));
    if (size(reduce(K,I,5)) <> 0)
    {
      nl=0;
      while ( nl < 30 )
      {
        nl=nl+1;
        phi=r1,randomLast(100);
        sbi=phi(i);
        I=simplify(lead(std(sbi)),1);
        attrib(I,"isSB",1);
        K=subst(select(I,n+1),x(n),1);
        K=K*maxideal(maxdeg1(I));
        if (size(reduce(K,I,5)) == 0)
        {
          break;
        }
      }
    }
    h=maxdeg1(simplify(reduce(quotient(I,maxideal(1)),I),2))+1;
    list l8;
    for ( ii = n; ii> n-d+1; ii-- )
    {
      l8=list();
      sbi=subst(sbi,x(ii),0);
      for (zz = 0; zz<= ii-1; zz++)
      {
        l8[zz+1] = "x("+string(zz)+")";
      }
      ring mr = create_ring(ring_list(r0)[1], l8, "dp", "no_minpoly");
      ideal sbi,I,L,K,T;
      map phi;
      sbi=imap(r1,sbi);
      I=simplify(lead(std(sbi)),1);
      attrib(I,"isSB",1);
      K=subst(select(I,ii),x(ii-1),1);
      K=K*maxideal(maxdeg1(I));
      if (size(reduce(K,I,5)) <> 0)
      {
        nl=0;
        while ( nl < 30 )
        {
          nl=nl+1;
          L=randomLast(100);
          phi=mr,L;
          T=phi(sbi);
          I=simplify(lead(std(T)),1);
          attrib(I,"isSB",1);
          K=subst(select(I,ii),x(ii-1),1);
          K=K*maxideal(maxdeg1(I));
          if (size(reduce(K,I,5)) == 0)
          {
            sbi=T;
            break;
          }
        }
      }
      H=maxdeg1(simplify(reduce(quotient(I,maxideal(1)),I),2))+1;
      if (H > h)
      {
        h=H;
      }
      setring r1;
      sbi=fetch(mr,sbi);
      kill mr;
    }
    sbi=subst(sbi,x(n-d+1),0);
    list l9;
    for (zz = 0; zz<= n-d; zz++)
    {
      l9[zz+1] = "x("+string(zz)+")";
    }
    ring mr = create_ring(ring_list(r0)[1], l9, "dp", "no_minpoly");
    ideal sbi,I,L,K,T;
    map phi;
    sbi=imap(r1,sbi);
    I=simplify(lead(std(sbi)),1);
    attrib(I,"isSB",1);
    H=maxdeg1(simplify(reduce(quotient(I,maxideal(1)),I),2))+1;
    if (H > h)
    {
      h=H;
    }
    dbprint(2,"reg(i)=" + string(h) + " and the time of this computation: " + string(rtimer-time) + " sec./100");
    return();
  }
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(t,a,b,c,d),dp;
ideal i=b4-a3d, ab3-a3c, bc4-ac3d-bcd3+ad4, c6-bc3d2-c3d3+bd5, ac5-b2c3d-ac2d3+b2d4, a2c4-a3d3+b3d3-a2cd3, b3c3-a3d3, ab2c3-a3cd2+b3cd2-ab2d3, a2bc3-a3c2d+b3c2d-a2bd3, a3c3-a3bd2, a4c2-a3b2d;
regCM(i);
}

//////////////////////////////////////////////////////////////////////////////
//
proc modregCM(ideal i)
"USAGE:  modregCM(i); i ideal
RETURN:  an integer, the Castelnuovo-Mumford regularity of i.
         (returns -1 if i is not homogeneous)
ASSUME:  i is a homogeneous ideal and the characteristic of base field is zero..
NOTE:    This is a probabilistic procedure, and it computes the initial of the ideal modulo the prime number 2147483647 (the biggest prime less than 2^31).
KEYWORDS: regularity
"
{
//--------------------------- initialisation ---------------------------------
  "// WARNING: The characteristic of base field must be zero.
// This procedure is probabilistic and  it computes the initial
//ideals modulo the prime number 2147483647";
  int e,ii,jj,H,h,d,time,p,nl;
  def r0 = basering;
  int n = nvars(r0)-1;
  list l12;
  for (int zz = 0; zz <= n; zz++)
  {
   l12[zz+1] = "x("+string(zz)+")";
  }
  ring r1 = create_ring(ring_list(r0)[1], l12, "dp", "no_minpoly");
  ideal i,sbi,I,J,K,L,lsbi1,lsbi2;
  list #;
  poly P;
  map phi;
  i = fetch(r0,i);
  time=rtimer;
  system("--ticks-per-sec",100);
//----- Check ideal homogeneous
  if ( homog(i) == 0 )
  {
    "// The ideal is not homogeneous!";
    return (-1);
  }
  option(redSB);
  p=2147483647;
  #=ringlist(r1);
  #[1]=p;
  def oro=ring(#);
  setring oro;
  ideal sbi,lsbi;
  sbi=fetch(r1,i);
  lsbi=lead(std(sbi));
  setring r1;
  lsbi1=fetch(oro,lsbi);
  lsbi1=simplify(lsbi1,1);
  attrib(lsbi1,"isSB",1);
  kill oro;
  I=lsbi1;
  d=dim(I);
  K=select(I,n+1);
  if (size(K) == 0)
  {
    h=0;
  }
  else
  {
    if (TestLastVarIsInGenericPos(I) == 1)
    {
      h=maxdeg1(K);
    }
    else
    {
      while ( nl < 30 )
      {
        nl=nl+1;
        phi=r1,randomLast(100);
        sbi=phi(i);
        #=ringlist(r1);
        #[1]=p;
        def oro=ring(#);
        setring oro;
        ideal sbi,lsbi;
        sbi=fetch(r1,sbi);
        lsbi=lead(std(sbi));
        setring r1;
        lsbi1=fetch(oro,lsbi);
        lsbi1=simplify(lsbi1,1);
        attrib(lsbi1,"isSB",1);
        kill oro;
        I=lsbi1;
        K=select(I,n+1);
        if (size(K) == 0)
        {
          h=0;break;
        }
        if (TestLastVarIsInGenericPos(I) == 1 )
        {
          h=maxdeg1(K);break;
        }
      }
      i=sbi;
    }
  }
  list l11;int zz;
  for ( ii = n; ii>=n-d+1; ii-- )
  {
    l11=list();
    i=subst(i,x(ii),0);
    for (zz = 0; zz <= ii-1; zz++)
    {
     l11[zz+1] = "x("+string(zz)+")";
    }
    ring mr = create_ring(0, l11, "dp");
    ideal i,sbi,I,J,K,L,lsbi1;
    poly P;
    list #;
    map phi;
    i=imap(r1,i);
    #=ringlist(mr);
    #[1]=p;
    def oro=ring(#);
    setring oro;
    ideal sbi,lsbi;
    sbi=fetch(mr,i);
    lsbi=lead(std(sbi));
    setring mr;
    lsbi1=fetch(oro,lsbi);
    lsbi1=simplify(lsbi1,1);
    attrib(lsbi1,"isSB",1);
    kill oro;
    I=lsbi1;
    K=select(I,ii);
    if (size(K) == 0)
    {
      H=0;
    }
    else
    {
      if (TestLastVarIsInGenericPos(I) == 1)
      {
        H=maxdeg1(K);
      }
      else
      {
        nl=0;
        while ( nl < 30 )
        {
          nl=nl+1;
          phi=mr,randomLast(100);
          sbi=phi(i);
          #=ringlist(mr);
          #[1]=p;
          def oro=ring(#);
          setring oro;
          ideal sbi,lsbi;
          sbi=fetch(mr,sbi);
          lsbi=lead(std(sbi));
          setring mr;
          lsbi1=fetch(oro,lsbi);
          lsbi1=simplify(lsbi1,1);
          kill oro;
          I=lsbi1;
          attrib(I,"isSB",1);
          K=select(I,ii);
          if (size(K) == 0)
          {
            H=0;break;
          }
          if (TestLastVarIsInGenericPos(I) == 1 )
          {
            H=maxdeg1(K);break;
          }
        }
        setring r1;
        i=imap(mr,sbi);
        kill mr;
      }
    }
    if (H > h)
    {
      h=H;
    }
  }
  dbprint(2,"mreg(i)=" + string(h) + " and the time of this computation: " + string(rtimer-time) + "sec./100");
  return();
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(t,a,b,c,d),dp;
ideal i=b4-a3d, ab3-a3c, bc4-ac3d-bcd3+ad4, c6-bc3d2-c3d3+bd5, ac5-b2c3d-ac2d3+b2d4, a2c4-a3d3+b3d3-a2cd3, b3c3-a3d3, ab2c3-a3cd2+b3cd2-ab2d3, a2bc3-a3c2d+b3c2d-a2bd3, a3c3-a3bd2, a4c2-a3b2d;
modregCM(i);
}
/*
//////////////////////////////////////////////////////////////
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(X,Y,a,b),dp;
poly f=X^8+a*Y^4-Y;
poly g=Y^8+b*X^4-X;
poly h=diff(f,X)*diff(g,Y)-diff(f,Y)*diff(g,X);
ideal i=f,g,h;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(x,y,z,a,b),dp;
ideal i=2*y^2*(y^2+x^2)+(b^2-3*a^2)*y^2-2*b*y^2*(x+y)+2*a^2*b*(y+x)-a^2*x^2+a^2*(a^2-b^2),4*y^3+4*y*(y^2+x^2)-2*b*y^2-4*b*y*(y+x)+2*(b^2-3*a^2)*y+2*a^2*b,4*x*y^2-2*b*y^2-2*a^2*x+2*a^2*b;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(t,a,b,c,d),dp;
ideal i=b4-a3d, ab3-a3c, bc4-ac3d-bcd3+ad4, c6-bc3d2-c3d3+bd5, ac5-b2c3d-ac2d3+b2d4, a2c4-a3d3+b3d3-a2cd3, b3c3-a3d3, ab2c3-a3cd2+b3cd2-ab2d3, a2bc3-a3c2d+b3c2d-a2bd3, a3c3-a3bd2, a4c2-a3b2d;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(a,b,c,d,e),dp;
ideal i=6*b4*c3+21*b4*c2*d+15b4cd2+9b4d3-8b2c2e-28b2cde+36b2d2e-144b2c-648b2d-120, 9b4c4+30b4c3d+39b4c2d2+18b4cd3-24b2c3e-16b2c2de+16b2cd2e+24b2d3e-432b2c2-720b2cd-432b2d2+16c2e2-32cde2+16d2e2+576ce-576de-240c+5184,-15b2c3e+15b2c2de-81b2c2+216b2cd-162b2d2+40c2e2-80cde2+40d2e2+1008ce-1008de+5184, -4b2c2+4b2cd-3b2d2+22ce-22de+261;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(c,b,d,p,q),dp;
ideal i=2*(b-1)^2+2*(q-p*q+p^2)+c^2*(q-1)^2-2*b*q+2*c*d*(1-q)*(q-p)+2*b*p*q*d*(d-c)+b^2*d^2*(1-2*p)+2*b*d^2*(p-q)+2*b*d*c*(p-1)+2*b*p*q*(c+1)+(b^2-2*b)*p^2*d^2+2*b^2*p^2+4*b*(1-b)*p+d^2*(p-q)^2,d*(2*p+1)*(q-p)+c*(p+2)*(1-q)+b*(b-2)*d+b*(1-2*b)*p*d+b*c*(q+p-p*q-1)+b*(b+1)*p^2*d, -b^2*(p-1)^2+2*p*(p-q)-2*(q-1),b^2+4*(p-q*q)+3*c^2*(q-1)*(q-1)-3*d^2*(p-q)^2+3*b^2*d^2*(p-1)^2+b^2*p*(p-2)+6*b*d*c*(p+q+q*p-1);
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(a,b,c,d,e,f),dp;
ideal i=2adef+3be2f-cef2,4ad2f+5bdef+cdf2,2abdf+3b2ef-bcf2,4a2df+5abef+acf2,4ad2e+3bde2+7cdef, 2acde+3bce2-c2ef, 4abde+3b2e2-4acdf+2bcef-c2f2, 4a2de+3abe2+7acef, 4acd2+5bcde+c2df, 4abd2+3b2de+7bcdf, 16a2d2-9b2e2+32acdf-18bcef+7c2f2, 2abcd+3b2ce-bc2f, 4a2cd+5abce+ac2f, 4a2bd+3ab2e+7abcf, abc2f-cdef2, ab2cf-bdef2, 2a2bcf+3be2f2-cef3, ab3f-3bdf3, 2a2b2f-4adf3+3bef3-cf4, a3bf+4aef3, 3ac3e-cde3, 3b2c2e-bc3f+2cd2ef, abc2e-cde2f, 6a2c2e-4ade3-3be4+ce3f, 3b3ce-b2c2f+2bd2ef, 2a2bce+3be3f-ce2f2, 3a3ce+4ae3f, 4bc3d+cd3e, 4ac3d-3bc3e-2cd2e2+c4f, 8b2c2d-4ad4-3bd3e-cd3f, 4b3cd+3bd3f, 4ab3d+3b4e-b3cf-6bd2f2, 4a4d+3a3be+a3cf-8ae2f2;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(x,y,z,t,u,v,w),dp;
ideal i=2tw+2wy-wz,2uw2-10vw2+20w3-7tu+35tv-70tw, 6tw2+2w2y-2w2z-21t2-7ty+7tz, 2v3-4uvw-5v2w+6uw2+7vw2-15w3-42vy, 6tw+9wy+2vz-3wz-21x, 9uw3-45vw3+135w4+14tv2-70tuw+196tvw-602tw2-14v2z+28uwz+14vwz-28w2z+147ux-735vx+2205wx-294ty+98tz+294yz-98z2, 36tw3+6w3y-9w3z-168t2w-14v2x+28uwx+14vwx-28w2x-28twy+42twz+588tx+392xy-245xz, 2uvw-6v2w-uw2+13vw2-5w3-28tw+14wy, u2w-3uvw+5uw2-28tw+14wy, tuw+tvw-11tw2-2vwy+8w2y+uwz-3vwz+5w2z-21wx, 5tuw-17tvw+33tw2-7uwy+22vwy-39w2y-2uwz+6vwz-10w2z+63wx, 20t2w-12uwx+30vwx-15w2x-10twy-8twz+4wyz, 4t2w-6uwx+12vwx-6w2x+2twy-2wy2-2twz+wyz, 8twx+8wxy-4wxz;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(a,b,c,d,x,w,u,v),dp;
ideal i=a+b+c+d,u+v+w+x, 3ab+3ac+3bc+3ad+3bd+3cd+2,bu+cu+du+av+cv+dv+aw+bw+dw+ax+bx+cx,bcu+bdu+cdu+acv+adv+cdv+abw+adw+bdw+abx+acx+bcx,abc+abd+acd+bcd,bcdu+acdv+abdw+abcx;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(b,x,y,z,s,t,u,v,w),dp;
ideal i=su+bv, tu+bw,tv+sw,sx+by,tx+bz,ty+sz,vx+uy,wx+uz,wy+vz;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(t,a,b,c,d,e,f,g,h),dp;
ideal i=a+c+d-e-h,2df+2cg+2eh-2h2-h-1,3df2+3cg2-3eh2+3h3+3h2-e+4h, 6bdg-6eh2+6h3-3eh+6h2-e+4h, 4df3+4cg3+4eh3-4h4-6h3+4eh-10h2-h-1, 8bdfg+8eh3-8h4+4eh2-12h3+4eh-14h2-3h-1, 12bdg2+12eh3-12h4+12eh2-18h3+8eh-14h2-h-1, -24eh3+24h4-24eh2+36h3-8eh+26h2+7h+1;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(a,b,c,d,e,f,g,h,k,l),dp;
ideal i=f2h-1,ek2-1,g2l-1, 2ef2g2hk2+f2g2h2k2+2ef2g2k2l+2f2g2hk2l+f2g2k2l2+ck2, 2e2fg2hk2+2efg2h2k2+2e2fg2k2l+4efg2hk2l+2fg2h2k2l+2efg2k2l2+2fg2hk2l2+2bfh, 2e2f2ghk2+2ef2gh2k2+2e2f2gk2l+4ef2ghk2l+2f2gh2k2l+2ef2gk2l2+2f2ghk2l2+2dgl, e2f2g2k2+2ef2g2hk2+2ef2g2k2l+2f2g2hk2l+f2g2k2l2+bf2, 2e2f2g2hk+2ef2g2h2k+2e2f2g2kl+4ef2g2hkl+2f2g2h2kl+2ef2g2kl2+2f2g2hkl2+2cek, e2f2g2k2+2ef2g2hk2+f2g2h2k2+2ef2g2k2l+2f2g2hk2l+dg2, -e2f2g2hk2-ef2g2h2k2-e2f2g2k2l-2ef2g2hk2l-f2g2h2k2l-ef2g2k2l2-f2g2hk2l2+a2;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(b,c,d,e,f,g,h,j,k,l),dp;
ideal i=-k9+9k8l-36k7l2+84k6l3-126k5l4+126k4l5-84k3l6+36k2l7-9kl8+l9, -bk8+8bk7l+k8l-28bk6l2-8k7l2+56bk5l3+28k6l3-70bk4l4-56k5l4+56bk3l5+70k4l5-28bk2l6-56k3l6+8bkl7+28k2l7-bl8-8kl8+l9, ck7-7ck6l-k7l+21ck5l2+7k6l2-35ck4l3-21k5l3+35ck3l4+35k4l4-21ck2l5-35k3l5+7ckl6+21k2l6-cl7-7kl7+l8, -dk6+6dk5l+k6l-15dk4l2-6k5l2+20dk3l3+15k4l3-15dk2l4-20k3l4+6dkl5+15k2l5-dl6-6kl6+l7, ek5-5ek4l-k5l+10ek3l2+5k4l2-10ek2l3-10k3l3+5ekl4+10k2l4-el5-5kl5+l6, -fk4+4fk3l+k4l-6fk2l2-4k3l2+4fkl3+6k2l3-fl4-4kl4+l5, gk3-3gk2l-k3l+3gkl2+3k2l2-gl3-3kl3+l4, -hk2+2hkl+k2l-hl2-2kl2+l3, jk-jl-kl+l2;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,x(0..10),dp;
ideal i=x(1)*x(0),x(1)*x(2),x(2)*x(3),x(3)*x(4),x(4)*x(5),x(5)*x(6),x(6)*x(7),x(7)*x(8),x(8)*x(9),x(9)*x(10),x(10)*x(0);
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(a,b,c,d,e,f,g,h,j,k,l,m,n,o,p,q,s),dp;
ideal i=ag,gj+am+np+q,bl,nq,bg+bk+al+lo+lp+b+c,ag+ak+jl+bm+bn+go+ko+gp+kp+lq+a+d+f+h+o+p,gj+jk+am+an+mo+no+mp+np+gq+kq+e+j+q+s-1,jm+jn+mq+nq,jn+mq+2nq,gj+am+2an+no+np+2gq+kq+q+s,2ag+ak+bn+go+gp+lq+a+d,bg+al, an+gq, 2jm+jn+mq, gj+jk+am+mo+2mp+np+e+2j+q, jl+bm+gp+kp+a+f+o+2p,lp+b,jn+mq,gp+a;
}
example
{ "EXAMPLE:"; echo = 2;
ring r=0,(a,b,c,d,e,f,g,h,v,w,k,l,m,n,o,p,q,s,t,u),dp;
ideal i=af+bg+ch+dv+ew-1/2, a2f+b2g+c2h+d2v+e2w-1/3,tdw+agk+ahl+bhm+avn+bvo+cvp+awq+bwu+cws-1/6, a3f+b3g+c3h+d3v+e3w-1/4, tdew+abgk+achl+bchm+advn+bdvo+cdvp+aewq+bewu+cews-1/8, td2w+a2gk+a2hl+b2hm+a2vn+b2vo+c2vp+a2wq+b2wu+c2ws-1/12, ahkm+tawn+tbwo+avko+tcwp+avlp+bvmp+awku+awls+bwms-1/24, a4f+b4g+c4h+d4v+e4w-1/5, tde2w+ab2gk+ac2hl+bc2hm+ad2vn+bd2vo+cd2vp+ae2wq+be2wu+ce2ws-1/10, td2ew+a2bgk+a2chl+b2chm+a2dvn+b2dvo+c2dvp+a2ewq+b2ewu+c2ews-1/15,achkm+taewn+tbewo+advko+tcewp+advlp+bdvmp+aewku+aewls+bewms-1/30,t2d2w+a2gk2+a2hl2+2abhlm+b2hm2+a2vn2+2abvno+b2vo2+2acvnp+2bcvop+c2vp2+2tadwq+a2wq2+2tbdwu+2abwqu+b2wu2+2tcdws+2acwqs+2bcwus+c2ws2-1/20,td3w+a3gk+a3hl+b3hm+a3vn+b3vo+c3vp+a3wq+b3wu+c3ws-1/20,abhkm+tadwn+tbdwo+abvko+tcdwp+acvlp+bcvmp+abwku+acwls+bcwms-1/40,a2hkm+ta2wn+tb2wo+a2vko+tc2wp+a2vlp+b2vmp+a2wku+a2wls+b2wms-1/60,tawko+tawlp+tbwmp+avkmp+awkms-1/20;
}
*/