source: git/Singular/LIB/modstd.lib @ 967543

//GP, last modified 23.10.06
///////////////////////////////////////////////////////////////////////////////
version="$Id: modstd.lib,v 1.2 2007-01-09 12:44:33 pfister Exp $";
category="Commutative Algebra";
info="
LIBRARY: modstd.lib  Grobner basis of ideals
AUTHORS: A. Hashemi,     Amir.Hashemi@lip6.fr
@*       G. Pfister      pfister@mathematik.uni-kl.de
@*       H. Schoenemann  hannes@mathematik.uni-kl.de

NOTE:
 A library for computing the Grobner basis of an ideal in the polynomial
 ring over the rational numbers using modular methods.The procedures are
 inspired by the following paper:
 Elizabeth A. Arnold:
 Modular Algorithms for Computing Groebner Bases , Journal of Symbolic
 Computation , April 2003, Volume 35, (4), p. 403-419.



PROCEDURES:
modStd(I);     compute a standard basis of I using modular methods
modS(I,L);     liftings to Q of standard bases of I mod p for p in L 
primeList(n);  list of n primes  <= 2134567879 in decreasing order
";

LIB "poly.lib";
LIB "krypto.lib";
///////////////////////////////////////////////////////////////////////////////
proc modStd(ideal I)
"USAGE:  modStd(I); 
RETURN:  a standard basis of I if no warning appears; 
NOTE:    the procedure computes a standard basis of I (over the 
         rational numbers) by using  modular methods. If a 
         warning appears then the result is a standard basis with no defined
         relation to I; this is a sign that not enough prime numbers have 
         been used. For further experiments see procedure modS.
EXAMPLE: example modStd; shows an example
"
{
  int aa=timer;
  def R0=basering;
  list rl=ringlist(R0);
  if((npars(R0)>0)||(rl[1]>0))
  {
     ERROR("characteristic of basering should be zero");
  }
  int l,j,k,q;
  int en=2134567879;
  int an=1000000000;
  intvec hi,hl,hc,hpl,hpc;
  list T,TT;
  list L=primeList(5);
  L[6]=prime(random(an,en));
  ideal J,cT,lT,K;
  ideal I0=I; 
  int h=homog(I);
  if((!h)&&(ord_test(R0)==0))
  {
     ERROR("input is not homogeneous and ordering is not local");
  }
  if(h)
  {
     execute("ring gn="+string(L[6])+",x(1.."+string(nvars(R0))+"),dp;");
     ideal I=fetch(R0,I);
     ideal J=std(I);
     hi=hilb(J,1);
     setring R0;
  }
  for (j=1;j<=size(L);j++)
  {
    rl[1]=L[j];
    def oro=ring(rl);
    setring oro;
    ideal I=fetch(R0,I);
    option(redSB);
    if(h)
    {
       ideal I1=std(I,hi);
    }
    else
    {
      if(ord_test(R0)==-1)
      {
         ideal I1=std(I);
      }
      else
      {
         matrix M;
         ideal I1=liftstd(I,M);
      }
    }
    setring R0;
    T[j]=fetch(oro,I1);
    kill oro;
  }
  //================= delete unlucky primes ====================
  // unlucky iff the leading ideal is wrong
  list LL=deleteUnluckyPrimes(T,L);
  T=LL[1];
  L=LL[2];
  lT=LL[3];
  //============ now all leading ideals are the same ============
  for(j=1;j<=ncols(T[1]);j++)
  {
    for(k=1;k<=size(L);k++)
    {
      TT[k]=T[k][j];
    }
    J[j]=liftPoly(TT,L); 
  }
  //=========== chooses more primes up to the moment the result becomes stable
  while(1)
  {
     k=0;
     q=prime(random(an,en));
     while(k<size(L))
     {
        k++;
        if(L[k]==q)
        {
           k=0;
           q=prime(random(an,en));
        }
     }
     L[size(L)+1]=q;
     rl[1]=L[size(L)];
     def @r=ring(rl);
     setring @r;
     ideal i=fetch(R0,I);
     option(redSB);
     if(h)
     {
       i=std(i,hi);
     }
     else
     {
       if(ord_test(R0)==-1)
       {
          i=std(i);
       }
       else
       {
          matrix M;
          i=liftstd(i,M);
       }
     }
     setring R0;
     T[size(T)+1]=fetch(@r,i);
     kill @r;
     cT=lead(T[size(T)]);
     attrib(cT,"isSB",1);
     if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
     {
        T=delete(T,size(T));
        L=delete(L,size(L));
        k=0;
     }
     else
     {
        for(j=1;j<=ncols(T[1]);j++)
        {
           for(k=1;k<=size(L);k++)
           {
              TT[k]=T[k][j];
           }
           K[j]=liftPoly(TT,L);
        }
        k=1;
        for(j=1;j<=size(K);j++)
        {
           if(K[j]-J[j]!=0)
           {
              k=0;
              J=K;
              break;
           }
        }
     }
     if(k){break;}
  }
  //============  test for standard basis and I=J =======
"computed";timer-aa;aa=timer;
  J=std(J);
  I0=reduce(I0,J);
  if(size(I0)>0){"WARNING: The input ideal is not contained 
                        in the ideal generated by the standardbasis";}
  attrib(J,"isSB",1);
"verification";timer-aa;
  return(J);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
   ideal J=modStd(I);
   J;
}   
///////////////////////////////////////////////////////////////////////////////
proc modS(ideal I, list L, list #)
"USAGE:  modS(I,L); I ideal, L list of primes
         if size(#)>0 std is used instead of groebner
RETURN:  an ideal which is with high probability a standard basis
NOTE:    This procedure is designed for fast experiments.
         It is not tested whether the result is a standard basis.
         It is not tested whether the result generates I.
EXAMPLE: example modS; shows an example
"
{
  int j,k;
  list T,TT;
  def R0=basering;
  ideal J,cT,lT,K;
  ideal I0=I;
  list rl=ringlist(R0);
  if((npars(R0)>0)||(rl[1]>0))
  {
     ERROR("characteristic of basering should be zero");
  }
  for (j=1;j<=size(L);j++)
  {
    rl[1]=L[j];
    def @r=ring(rl);
    setring @r;
    ideal i=fetch(R0,I);
    option(redSB);
    if(size(#)>0)
    {
       i=std(i);
    }
    else
    {
       i=groebner(i)
    }
    setring R0;
    T[j]=fetch(@r,i);
    kill @r;
  }
  //================= delete unlucky primes ====================
  // unlucky iff the leading ideal is wrong
  list LL=deleteUnluckyPrimes(T,L);
  T=LL[1];
  L=LL[2];
  //============ now all leading ideals are the same ============
  for(j=1;j<=ncols(T[1]);j++)
  {
    for(k=1;k<=size(L);k++)
    {
      TT[k]=T[k][j];
    }
    J[j]=liftPoly(TT,L);  
  }
  attrib(J,"isSB",1);
  return(J);
}
example
{ "EXAMPLE:"; echo = 2;
   list L=3,5,11,13,181;
   ring r=0,(x,y,z),dp;
   ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
   ideal J=modS(I,L);
   J;
}  
/////////////////////////////////////////////////////////////////////////////// 
proc deleteUnluckyPrimes(list T,list L)
"USAGE:  deleteUnluckyPrimes(T,L);T list of polys, L list of primes  
RETURN:  list L,T with T list of polys, L list of primes
EXAMPLE: example deleteUnluckyPrimes; shows an example
NOTE:    works only for homogeneous ideals with global orderings or
         for ideals with local orderings
"
{
  int j,k;
  intvec hl,hc;
  ideal cT,lT;

  lT=lead(T[size(T)]);
  attrib(lT,"isSB",1);
  hl=hilb(lT,1);
  for (j=1;j<size(T);j++)
  {
     cT=lead(T[j]);
     attrib(cT,"isSB",1);
     hc=hilb(cT,1);
     if(hl==hc)
     {
        for(k=1;k<=size(lT);k++)
        {
           if(lT[k]<cT[k]){lT=cT;break;}
           if(lT[k]>cT[k]){break;}
        }
     }
     else
     {
        if(hc<hl){lT=cT;hl=hilb(lT,1);}
     }
  }
  j=1;
  attrib(lT,"isSB",1);
  while(j<=size(T))
  {
     cT=lead(T[j]);
     attrib(cT,"isSB",1);
     if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
     {
        T=delete(T,j);
        L=delete(L,j);
        j--;
     }
     j++;
  }
  return(list(T,L,lT));
}
example
{ "EXAMPLE:"; echo = 2;
   list L=2,3,5,7,11;
   ring r=0,(y,x),Dp;
   ideal I1=y2x,y6;
   ideal I2=yx2,y3x,x5,y6;
   ideal I3=y2x,x3y,x5,y6;
   ideal I4=y2x,x3y,x5;
   ideal I5=y2x,yx3,x5,y6;
   list T=I1,I2,I3,I4,I5;
   list TT=deleteUnluckyPrimes(T,L);
   TT;
}  
///////////////////////////////////////////////////////////////////////////////
proc liftPoly(list T, list L)
"USAGE:  liftPoly(T,L); T list of polys, L list of primes
RETURN:  poly p in Q[x] such that p mod L[i]=T[i]
EXAMPLE: example liftPoly; shows an example
"
{
   poly result;
   int i;
   poly p;
   list TT;
   number n;
  
   number N=L[1];
   for(i=2;i<=size(L);i++)
   {
      N=N*L[i];
   }
   while(1)
   {
     p=leadmonom(T[1]);
     for(i=2;i<=size(T);i++)
     {
        if(leadmonom(T[i])>p)
        {
          p=leadmonom(T[i]);
        }
     }
     if (p==0) {return(result);}
     for(i=1;i<=size(T);i++)
     {
       if(p==leadmonom(T[i]))
       {
          TT[i]=leadcoef(T[i]);
          T[i]=T[i]-lead(T[i]);
       }
       else
       {
          TT[i]=0;
       }
     }
     n=chineseR(TT,L,N);
     n=Farey(N,n); 
     result=result+n*p;
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,(x,y),dp;
   list L=32003,181,241,499;
   list T=x2+7000x+13000,x2+100x+147y+40,x2+120x+191y+10,x2+x+67y+100;
   liftPoly(T,L);
}
 /////////////////////////////////////////////////////////////////////////////// 
proc fareyIdeal(ideal I,list L)
{
   poly result,p;
   int i,j;
   number n;
   number N=L[1];
   for(i=2;i<=size(L);i++)
   {
      N=N*L[i];
   }

   for(i=1;i<=size(I);i++)
   {
     p=I[i];
     result=lead(p);
     while(1)
     {
        if (p==0) {break;}
        p=p-lead(p);
        n=Farey(N,leadcoef(p)); 
        result=result+n*leadmonom(p);
     }
     I[i]=result;
   }
   return(I);
}
///////////////////////////////////////////////////////////////////////////////
proc Farey (number P, number N)
"USAGE:  Farey (P,N); P, N number;
RETURN:  a rational number a/b such that a/b=N mod P 
         and |a|,|b|<(P/2)^{1/2}
"
{
   if (P<0){P=-P;}
   if (N<0){N=N+P;}
   number A,B,C,D,E;
   E=P;
   B=1;
   while (N!=0)
   {
        if (2*N^2<P)
        {           
           return(N/B);
        }
        D=E mod N;
        C=A-(E-E mod N)/N*B;
        E=N;
        N=D;
        A=B;
        B=C;
   }
   return(0);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,x,dp;
   Farey(32003,12345);
}
 ///////////////////////////////////////////////////////////////////////////////
proc chineseR(list T,list L,number N)
"USAGE:  chineseR(T,L);
RETURN: x such that x = T[i] mod L[i]
NOTE:   chinese remainder theorem
EXAMPLE:example chineseR; shows an example
"
{
   number x;
   if(size(L)==1)
   {
      x=T[1] mod L[1];
      return(x);
   }
   int i;
   int n=size(L);
   list M;
   for(i=1;i<=n;i++)
   {
      M[i]=N/L[i];
   }
   list S=eexgcdN(M);
   for(i=1;i<=n;i++)
   {
      x=x+S[i]*M[i]*T[i];
   }
   x=x mod N;
   return(x);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,x,dp;
   chineseR(list(24,15,7),list(2,3,5));
}
  
///////////////////////////////////////////////////////////////////////////////
proc primeList(int n)
"USAGE:  primeList(n);
RETURN: the list of n greatest primes  <= 2134567879   
EXAMPLE:example primList; shows an example
"
{
  list L;
  int i;
  int p=2134567879;
  for(i=1;i<=n;i++)
  {
    L[i]=p;
    p=prime(p-1);
  }
  return(L);
}
example
{ "EXAMPLE:"; echo = 2;
   list L=primeList(10);
   size(L);
   L[size(L)];
}   
///////////////////////////////////////////////////////////////////////////////
/*
ring r=0,(x,y,z),lp;
poly s1 = 5x3y2z+3y3x2z+7xy2z2;
poly s2 = 3xy2z2+x5+11y2z2;
poly s3 = 4xyz+7x3+12y3+1;
poly s4 = 3x3-4y3+yz2;
ideal i =  s1, s2, s3, s4;

ring r=0,(x,y,z),lp;
poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
ideal i =  s1, s2, s3, s4;

ring r=0,(x,y,z),lp;
poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
poly s3 = 8x3 + 12y3 + xz2 + 3;
poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
ideal i =  s1, s2, s3, s4;

int n = 6;
ring r = 0,(x(1..n)),lp;
ideal i = cyclic(n);
ring s=0,(x(1..n),t),lp;
ideal i=imap(r,i);
i=homog(i,t);

ring r=0,(x(1..4),s),(dp(4),dp);
poly s1 =1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4);
poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
ideal i =  s1, s2, s3, s4, s5;

ring r=0,(x,y,z),ds;
int a =16;
int b =15;
int c =4;
int t =1;
poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2;
ideal i= jacob(f);

ring r=0,(x,y,z),ds;
int a =25;
int b =25;
int c =5;
int t =1;
poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2;
ideal i= jacob(f),f;

ring r=0,(x,y,z),ds;
int a=10;
poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
ideal i= jacob(f);

ring r=0,(x,y,z),ds;
int a =6;
int b =8;
int c =10;
int alpha =5;
int beta= 5;
int t= 1;
poly f =x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)+x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
ideal i= jacob(f);

*/