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

//GP, last modified 23.10.06
///////////////////////////////////////////////////////////////////////////////
version="$Id$";
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);  intvec of n primes  <= 2134567879 in decreasing order
";

LIB "poly.lib";
LIB "crypto.lib";
///////////////////////////////////////////////////////////////////////////////
proc pTestSB(ideal I, ideal J, list L)
"USAGE:  pTestSB(I,J,L); I,J ideals, L intvec of primes
RETURN: 1 (resp. 0) if for a randomly chosen prime p not in L
        J mod p is (resp. is not) a standard basis of I mod p
EXAMPLE:example primList; shows an example
"
{
  int i,j,k,p;
  def R=basering;
  list r= ringlist(R);

  while(!j)
  {
    j=1;
    p=prime(random(1000000000,2134567879));
    for(i=1;i<=size(L);i++)
    {
      if(p==L[i]){j=0;break}
    }
    if(j)
    {
      for(i=1;i<=ncols(J);i++)
      {
        for(k=2;k<=size(J[i]);k++)
        {
          if((denominator(leadcoef(J[i][k])) mod
p)==0){j=0;break}
        }
        if(!j){break;}
      }
    }
  }
  r[1]=p;
  def @R=ring(r);
  setring @R;
  ideal I=imap(R,I);
  ideal J=imap(R,J);
  attrib(J,"isSB",1);
  j=1;
  if(size(reduce(I,J))!=0){j=0;}
  if(j)
  {
    ideal K=std(J);
    if(size(reduce(K,J))!=0){j=0;}
  }
  setring R;
  return(j);
}
example
{ "EXAMPLE:"; echo = 2;
   intvec L=2,3,5;
   ring r=0,(x,y,z),dp;
   ideal i=x+1,x+y+1;
   ideal j=x+1,y;
   pTestSB(i,i,L);
   pTestSB(i,j,L);
}

proc primeList(int n, list #)
"USAGE:  primeList(n); (resp. primeList(n,L); )
RETURN: the intvec of n greatest primes  <= 2147483647 (resp. n greatest
        primes < L[size(L)] union with L)
EXAMPLE:example primList; shows an example
"
{
  intvec L;
  int i,p;
  if(size(#)==0)
  {
      p=2147483647;
      L[1]=p;
   }
   else
  {
     L=#[1];
     p=prime(L[size(L)]-1);
     L[size(L)+1]=p;
  }
  if(p==2){ERROR("no more primes");}
  for(i=2;i<=n;i++)
  {
    p=prime(p-1);
    L[size(L)+1]=p;
  }
  return(L);
}
example
{ "EXAMPLE:"; echo = 2;
   intvec L=primeList(10);
   size(L);
   L[size(L)];
   L=primeList(5,L);
   size(L);
   L[size(L)];
}


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
         containing I and with high probability a standard basis of I.
         For further experiments see procedure modS.
EXAMPLE: example modStd; shows an example
"
{
  def R0=basering;
  list rl=ringlist(R0);
  if((npars(R0)>0)||(rl[1]>0))
  {
     ERROR("characteristic of basering should be zero, basering should have
no parameters");
  }

  int k,c;
  int pd=printlevel-voice+2;
  int j=1;
  int h=homog(I);
  int en=2134567879;
  int an=1000000000;

  intvec opt = option(get);     // Save current options
  intvec L=primeList(10);
  L[5]=prime(random(an,en));
  list T,TT,TH,LL;


  ideal J,K;

  option(redSB);

  while(1)
  {
     while(j<=size(L))
     {
        if(pd>2){c++;c;}
        rl[1]=L[j];
        def @r=ring(rl);
        setring @r;
        ideal i=fetch(R0,I);
        i=groebner(i);
        setring R0;
        T[size(T)+1]=fetch(@r,i);
        kill @r;
        j++;
     }
     if(pd>2){"lifting";}
  //================= delete unlucky primes ====================
  // unlucky if and only if the leading ideal is wrong
     LL=deleteUnluckyPrimes(T,L);
     TH=LL[1];
     L=LL[2];
  //============ now all leading ideals are the same ============
     for(j=1;j<=ncols(TH[1]);j++)
     {
         for(k=1;k<=size(L);k++)
         {
             TT[k]=TH[k][j];
         }
         J[j]=liftPoly(TT,L);
     }
     if(pd>2){"list of primes:";L;"pTest";}
     if(pTestSB(I,J,L))
     {
        attrib(J,"isSB",1);
        K=std(J);
        if(size(reduce(I,J))==0)
        {
           if(size(reduce(K,J))==0)
           {
               if(!((h)||(ord_test(R0)==-1)))
              {

"==================================================================";
                 "    The input is not homogeneous and the ordering is not
local.   ";
                 "WARNING: ideal generated by output may be greater then
input ideal";

"==================================================================";
              }
              option(set, opt);
              return(J);
           }
           if(pd>2){"pTest o.k. but result wrong";}
         }
         if(pd>2){"pTest o.k. but result wrong";}
      }
      j=size(L)+1;
      L=primeList(10,L);
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z,t),dp;
   ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
   ideal J=modStd(I);
   J;
   I=homog(I,t);
   J=modStd(I);
   J;
   ring s=0,(x,y,z),ds;
   ideal I=jacob(x5+y6+z7+xyz);
   ideal J=modStd(I);
   J;
}
///////////////////////////////////////////////////////////////////////////////
proc modS(ideal I, intvec L, list #)
"USAGE:  modS(I,L); I ideal, L intvec 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 if and only if 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;
   intvec 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,intvec L)
"USAGE:  deleteUnluckyPrimes(T,L);T list of polys, L intvec of primes
RETURN:  list L,T with T list of polys, L intvec 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);
        if(j==1) { L=L[2..size(L)]; }
        else
        {
          if (j==size(L)) { L=L[1..size(L)-1]; }
          else { L=L[1..j-1],L[j+1..size(L)]; }
        }
        j--;
     }
     j++;
  }
  return(list(T,L,lT));
}
example
{ "EXAMPLE:"; echo = 2;
   intvec 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, intvec L)
"USAGE:  liftPoly(T,L); T list of polys, L intvec of primes
RETURN:  poly p in Q[x] such that p mod L[i]=T[i]
EXAMPLE: example liftPoly; shows an example
"
{
   int i;
   list TT;
   for(i=size(T);i>0;i--)
   { TT[i]=ideal(T[i]); }
   T=TT;
   ideal hh=chinrem(T,L);
   poly h=hh[1];
   poly p=lead(h);
   poly result;
   number n;
   number N=L[1];
   for(i=size(L);i>1;i--)
   {
      N=N*L[i];
   }
   while(h!=0)
   {
     n=Farey(N,leadcoef(h));
     result=result+n*p;
     h=h-lead(h);
     p=leadmonom(h);
   }
   return(result);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,(x,y),dp;
   intvec L=32003,181,241,499;
   list T=ideal(x2+7000x+13000),ideal(x2+100x+147y+40),ideal(x2+120x+191y+10),ideal(x2+x+67y+100);
   liftPoly(T,L);
}
///////////////////////////////////////////////////////////////////////////
proc liftPoly1(list T, intvec L)
"USAGE:  liftPoly1(T,L); T list of polys, L intvec of primes
RETURN:  poly p in Q[x] such that p mod L[i]=T[i]
EXAMPLE: example liftPoly1; 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;
   intvec L=32003,181,241,499;
   list T=x2+7000x+13000,x2+100x+147y+40,x2+120x+191y+10,x2+x+67y+100;
   liftPoly1(T,L);
}
 ///////////////////////////////////////////////////////////////////////////////
proc fareyIdeal(ideal I,intvec 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,intvec L,number N)
"USAGE:  chineseR(T,L,N);
RETURN: x such that x = T[i] mod L[i], N=product(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),intvec(2,3,5),30);
}

///////////////////////////////////////////////////////////////////////////////
proc pStd(int p,ideal i)
"USAGE:  pStd(p,i);p integer, i ideal;
RETURN:  an ideal G which is the groebner base for i
EXAMPLE: example pStd; shows an example
"
{
  def r=basering;
  list rl=ringlist(r);
  rl[1]=p;
  def r1=ring(rl);
  setring r1;
  option(redSB);
  ideal j=fetch(r,i);
  ideal GP=groebner(j);
  setring r;
  ideal G=fetch(r1,GP);
  attrib(G,"isSB",1);
  matrix Z=transmat(p,i,G);
  matrix G1=gstrich1(p,Z,i,G);
  ideal g1=G1;
  ideal g22=reduce(g1,G);
  matrix G22=transpose(matrix(g22));
  matrix M=redmat(G,G1,G22);
  matrix Z2=-M*Z;
  kill r1;
  number c=p;
  matrix G0=transpose(matrix(G));
  G0= MmodN(G0+ (c)* G22,c^2);
  matrix GF=fareyMatrix(G0,c^2);
  Z=MmodN(Z+(c)*Z2,c^2);
  matrix C=transpose(G);
  int n=3;
  while(GF<>C)
  {
    C=GF;
    G1= gstrich2(c,Z,i,G0,n);
    g1=G1;
    g22=reduce(g1,G);
    G22=transpose(matrix(g22));
    M=redmat(G,G1,G22);
    Z2=-M*Z;
    Z=MmodN(Z+(c^(n-1))*Z2,c^n);
    G0= MmodN(G0+ (c^(n-1))* G22,c^n);
    GF=fareyMatrix(G0,c^n);
    n++;
  }
  return(ideal(GF));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
   ideal J=pStd(32003,I);
   J;
}
///////////////////////////////////////////////////////////////////////////
proc transmat(int p,ideal i,ideal G)
"USAGE:  transmat(p,I,G); p integer, I,G ideal;
RETURN:  the transformationmatrix Z for the ideal i mod p and the groebner base for i mod p
EXAMPLE: example transmat; shows an example
"
{
  def r=basering;
  int n=nvars(r);
  list rl=ringlist(r);
  rl[1]=p;
  def r1=ring(rl);
  setring r1;
  ideal i=fetch(r,i);
  ideal G=fetch(r,G);
  attrib(G,"isSB",1);
  ring rhelp=p,x(1..n),dp;
  list lhelp=ringlist(rhelp);
  list l=lhelp[3];
  setring r;
  rl[3]=l;
  def r2=ring(rl);
  setring r2;
  ideal i=fetch(r,i);
  option(redSB);
  ideal j=std(i);
  matrix T=lift(i,j);
  setring r1;
  matrix T=fetch(r2,T);
  ideal j=fetch(r2,j);
  matrix M=lift(j,G);
  matrix Z=transpose(T*M);
  setring r;
  matrix Z=fetch(r1,Z);
  return(Z);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal i=3x3+x2+1,11y5+y3+2,5z4+z2+4;
   ideal g=x3-60x2-60, z4-36z2+37, y5+33y3+66;
   int p=181;
   matrix Z=transmat(p,i,g);
   Z;
}

///////////////////////////////////////////////////////////////////////////
proc gstrich1(int p, matrix Z, ideal i, ideal gp)
"USAGE:  gstrich1 (p,Z,i,gp); p integer, Z matrix, i,gp ideals;
RETURN:  a matrix G such that (Z*F-GP)/p, where F and GP are the matrices of the ideals i and gp
"
{
  matrix F=transpose(matrix(i));
  matrix GP=transpose(matrix(gp));
  matrix G=(Z*F-GP)/p;
  return(G);
}
///////////////////////////////////////////////////////////////////////////
proc gstrich2(number p, matrix Z, ideal i, ideal gp, int n)
"USAGE:  gstrich2 (p,Z,i,gp,n); p,n integer, Z matrix, i,gp ideals;
RETURN:  a matrix G such that (Z*F-GP)/(p^(n-1)), where F and GP are the matrices of the ideals i and gp
"
{
  matrix F=transpose(matrix(i));
  matrix GP=transpose(matrix(gp));
  matrix G=(Z*F-GP)/(p^(n-1));
  return(G);
}
///////////////////////////////////////////////////////////////////////////
proc redmat(ideal i, matrix h, matrix g)
"USAGE:  redmat(i,h,g); i ideal , h,g matrices;
RETURN:  a matrix M such that i=M*h+g
"
{
  matrix c=h-g;
  ideal f=transpose(c);
  matrix N=lift(i,f);
  matrix M=transpose(N);
  return(M);
}
///////////////////////////////////////////////////////////////////////////
proc fareyMatrix(matrix m,number N)
"USAGE:  fareyMatrix(m,y); m matrix, y integer;
RETURN:  a matrix k of the matrix m with Farey rational numbers a/b as coefficients
EXAMPLE: example fareyMatrix; shows an example
"
{
  ideal I=m;
  poly result,p;
  int i,j;
  number n;
  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;
  }
  matrix k=transpose(I);
  return(k);
}
example
{"EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   matrix m[3][1]=x3+682794673x2+682794673,z4+204838402z2+819353608,    y5+186216729y3+372433458;
   int p=32003;
   matrix b=fareyMatrix(m,p^2);
   b;
}
///////////////////////////////////////////////////////////////////////////
proc MmodN(matrix Z,number N)
"USAGE:  MmodN(Z,N);Z matrix, N number;
RETURN:  the matrix Z mod N
EXAMPLE: example MmodN;
"
{
  int i,j,k;
  poly m,p;
  number c;
  for(i=1;i<=nrows(Z);i++)
  {
    for(j=1;j<=ncols(Z);j++)
    {
      for(k=1;k<=size(Z[i,j]);k++)
      {
        m=leadmonom(Z[i,j][k]);
        c=leadcoef(Z[i,j][k]) mod N;
        p=p+c*m;
      }
      Z[i,j]=p;
      p=0;
    }
  }
  return(Z);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),dp;
   matrix m[3][1]= x3+10668x2+10668, z4-12801z2+12802, y5-8728y3+14547;
   number p=32003;
   matrix b=MmodN(m,p^2);
   b;
}
///////////////////////////////////////////////////////////////////////////////
/*
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);

*/

/*
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);

*/