////////////////////////////////////////////////////////////////////////////////
version="$Id$";
category = "Commutative Algebra";
info="
LIBRARY:  assprimeszerodim.lib   associated primes of a zero-dimensional ideal

AUTHORS:  N. Idrees       nazeranjawwad@gmail.com
@*        G. Pfister      pfister@mathematik.uni-kl.de
@*        S. Steidel      steidel@mathematik.uni-kl.de

OVERVIEW:

  A library for computing the associated primes and the radical of a
  zero-dimensional ideal in the polynomial ring over the rational numbers,
  Q[x_1,...,x_n], using modular computations.

SEE ALSO: primdec_lib

PROCEDURES:
 zeroRadical(I);       computes the radical of I
 assPrimes(I);         computes the associated primes of I
";

LIB "primdec.lib";
LIB "modstd.lib";

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

proc zeroRadical(ideal I, list #)
"USAGE:  zeroRadical(I,[n]); I ideal, optional: n number of processors (for
         parallel computing)
ASSUME:  I is zero-dimensional in Q[variables]
NOTE:    Parallelization is just applicable using 32-bit Singular version since
         MP-links are not compatible with 64-bit Singular version.
RETURN:  the radical of I
EXAMPLE: example zeroRadical; shows an example
"
{
   return(zeroR(modStd(I,#),#));
}
example
{ "EXAMPLE:";  echo = 2;
   ring R = 0, (x,y), dp;
   ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
   zeroRadical(I);
}

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

static proc zeroR(ideal I, list #)
// compute the radical of I provided that I is zero-dimensional in Q[variables]
// and a standard basis
{
   attrib(I,"isSB",1);
   int i, k;
   int j = 1;
   int index = 1;
   int crit;

   list CO1, CO2, P;
   ideal G, F;
   bigint N;
   poly f;

//---------------------  Initialize optional parameter  ------------------------
   if(size(#) > 0)
   {
      int n = #[1];
      if((n > 1) && (1 - system("with","MP")))
      {
     "========================================================================";
     "There is no MP available on your system. Since this is necessary to     ";
     "parallelize the algorithm, the computation will be done without forking.";
     "========================================================================";
         n = 1;
      }
   }
   else
   {
      int n = 1;
   }

//--------------------  Initialize the list of primes  -------------------------
   intvec L = primeList(I,10);
   L[5] = prime(random(100000000,536870912));

   if(n > 1)
   {

//-----  Create n links l(1),...,l(n), open all of them and compute  -----------
//-----  polynomial F for the primes L[2],...,L[n + 1].              -----------

      for(i = 1; i <= n; i++)
      {
         link l(i) = "MPtcp:fork";
         open(l(i));
         write(l(i), quote(zeroRadP(eval(I), eval(L[i + 1]))));
      }

      int t = timer;
      P = zeroRadP(I, L[1]);
      t = timer - t;
      if(t > 60) { t = 60; }
      int i_sleep = system("sh", "sleep "+string(t));
      CO1[index] = P[1];
      CO2[index] = bigint(P[2]);
      index++;

      j = j + n + 1;
   }

//---------  Main computations in positive characteristic start here  ----------

   while(!crit)
   {
      if(n > 1)
      {
         while(j <= size(L) + 1)
         {
            for(i = 1; i <= n; i++)
            {
               if(status(l(i), "read", "ready"))
               {
                  //--- read the result from l(i) ---
                  P = read(l(i));
                  CO1[index] = P[1];
                  CO2[index] = bigint(P[2]);
                  index++;

                  if(j <= size(L))
                  {
                     write(l(i), quote(zeroRadP(eval(I), eval(L[j]))));
                     j++;
                  }
                  else
                  {
                     k++;
                     close(l(i));
                  }
               }
            }
            //--- k describes the number of closed links ---
            if(k == n)
            {
               j++;
            }
            //--- sleep for t seconds ---
            i_sleep = system("sh", "sleep "+string(t));
         }
      }
      else
      {
         while(j<=size(L))
         {
            P = zeroRadP(I, L[j]);
            CO1[index] = P[1];
            CO2[index] = bigint(P[2]);
            index++;
            j++;
         }
      }

      // insert deleteUnluckyPrimes
      G = chinrem(CO1,CO2);
      N = CO2[1];
      for(i = 2; i <= size(CO2); i++){N = N*CO2[i];}
      F = farey(G,N);

      crit = 1;
      for(i = 1; i <= nvars(basering); i++)
      {
         if(reduce(F[i],I) != 0) { crit = 0; break; }
      }

      if(!crit)
      {
         CO1 = G;
         CO2 = N;
         index = 2;

         j = size(L) + 1;
         L = primeList(I,10,L);

         if(n > 1)
         {
            for(i = 1; i <= n; i++)
            {
               open(l(i));
               write(l(i), quote(zeroRadP(eval(I), eval(L[j+i-1]))));
            }
            j = j + n;
            k = 0;
         }
      }
   }

   k = 0;
   for(i = 1; i <= nvars(basering); i++)
   {
      f = gcd(F[i],diff(F[i],var(i)));
      k = k + deg(f);
      F[i] = F[i]/f;
   }

   if(k == 0) { return(I); }
   else { return(modStd(I + F)); }
}

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

proc assPrimes(list #)
"USAGE:  assPrimes(I,[a],[n]); I ideal or module,
         optional: int n: number of processors (for parallel computing), int a:
         - a = 1: method of Eisenbud/Hunecke/Vasconcelos
         - a = 2: method of Gianni/Trager/Zacharias
         - a = 3: method of Monico
         assPrimes(I) chooses n = a = 1
ASSUME:  I is zero-dimensional over Q[variables]
NOTE:    Parallelization is just applicable using 32-bit Singular version since
         MP-links are not compatible with 64-bit Singular version.
RETURN:  a list Re of associated primes of I:
EXAMPLE: example assPrimes; shows an example
"
{
   ideal I;
   if(typeof(#[1]) == "ideal")
   {
      I = #[1];
   }
   else
   {
      module M = #[1];
      I = Ann(M);
   }

//---------------------  Initialize optional parameter  ------------------------
   if(size(#) > 1)
   {
      if(size(#) == 2)
      {
         int alg = #[2];
         int n = 1;
      }
      if(size(#) == 3)
      {
         int alg = #[2];
         int n = #[3];
         if((n > 1) && (1 - system("with","MP")))
         {
     "========================================================================";
     "There is no MP available on your system. Since this is necessary to     ";
     "parallelize the algorithm, the computation will be done without forking.";
     "========================================================================";
            n = 1;
         }
      }
   }
   else
   {
      int alg = 1;
      int n = 1;
   }

   def SPR = basering;
   I = modStd(I,n);
   int d = vdim(I);
   if(d == -1) { ERROR("Ideal is not zero-dimensional."); }
   if(homog(I) == 1){ return(list(maxideal(1))); }
   poly f = findGen(I);

   int T = timer;
   int RT = rtimer;
   int TT;
   if(size(f) == nvars(SPR))
   {
      TT = timer;
      int spT = pTestRad(d,f,I);
      if(printlevel>=10)
      {
         "pTestRad(d,f,I) = "+string(spT)+" and takes
         "+string(timer-TT)+" seconds.";
      }
      if(!spT)
      {
         if(typeof(attrib(#[1],"isRad")) == "int")
         {
            if(attrib(#[1],"isRad") == 0)
            {
               TT = timer;
               I = zeroR(I,n);
               if(printlevel>=10)
               {
                  "zeroR(I,n) takes "+string(timer-TT)+" seconds.";
               }
               TT = timer;
               I = std(I);
               if(printlevel>=10)
               {
                  "std(I) takes "+string(timer-TT)+" seconds.";
               }
               d = vdim(I);
               f = findGen(I);
            }
         }
         else
         {
            TT = timer;
            I = zeroR(I,n);
            if(printlevel>=10)
            {
               "zeroR(I,n) takes "+string(timer-TT)+" seconds.";
            }
            TT = timer;
            I = std(I);
            if(printlevel>=10)
            {
               "std(I) takes "+string(timer-TT)+" seconds.";
            }
            d = vdim(I);
            f = findGen(I);
         }
      }
   }
   if(printlevel>=10)
   {
      "Real-time for radical-check is "+string(rtimer - RT)+" seconds.";
      "CPU-time for radical-check is "+string(timer - T)+" seconds.";
   }

   export(SPR);
   poly f_for_fork = f;
   export(f_for_fork);           // f available for each link
   ideal I_for_fork = I;
   export(I_for_fork);           // I available for each link

//--------------------  Initialize the list of primes  -------------------------
   intvec L = primeList(I,10);
   L[5] = prime(random(1000000000,2134567879));

   list Re;

   ring rHelp = 0,T,dp;
   list CO1,CO2,P;
   ideal F,G,testF;
   bigint N;

   list ringL = ringlist(SPR);
   int i,k,e,int_break;
   int j = 1;
   int index = 1;

//-----  If there is more than one processor available, we parallelize the  ----
//-----  main standard basis computations in positive characteristic        ----

   if(n > 1)
   {

//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
//-----  standard basis for the primes L[2],...,L[n + 1].              ---------

      for(i = 1; i <= n; i++)
      {
         link l(i) = "MPtcp:fork";
         open(l(i));
         write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[i + 1]),
                                                          eval(alg))));
      }

      int t = timer;
      P = modpSpecialAlgDep(ringL, L[1], alg);
      t = timer - t;
      if(t > 60) { t = 60; }
      int i_sleep = system("sh", "sleep "+string(t));
      CO1[index] = P[1];
      CO2[index] = bigint(P[2]);
      index++;

      j = j + n + 1;
   }

//---------  Main computations in positive characteristic start here  ----------

   int tt = timer;
   int rt = rtimer;

   while(1)
   {
      tt = timer;
      rt = rtimer;

      if(n > 1)
      {
         while(j <= size(L) + 1)
         {
            for(i = 1; i <= n; i++)
            {
               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
               if(status(l(i), "read", "ready"))
               {
                  //--- read the result from l(i) ---
                  P = read(l(i));
                  CO1[index] = P[1];
                  CO2[index] = bigint(P[2]);
                  index++;

                  if(j <= size(L))
                  {
                     write(l(i), quote(modpSpecialAlgDep(eval(ringL),
                                                         eval(L[j]),
                                                         eval(alg))));
                     j++;
                  }
                  else
                  {
                     k++;
                     close(l(i));
                  }
               }
            }
            //--- k describes the number of closed links ---
            if(k == n)
            {
               j++;
            }
            i_sleep = system("sh", "sleep "+string(t));
         }
      }
      else
      {
         while(j<=size(L))
         {
            P = modpSpecialAlgDep(ringL, L[j], alg);
            CO1[index] = P[1];
            CO2[index] = bigint(P[2]);
            index++;
            j++;
         }
      }
      if(printlevel>=10)
      {
         "Real-time for computing list in assPrimes is "+string(rtimer - rt)+
         " seconds.";
         "CPU-time for computing list in assPrimes is "+string(timer - tt)+
         " seconds.";
      }

      // insert deleteUnluckyPrimes
      G = chinrem(CO1,CO2);
      N = CO2[1];
      for(j = 2; j <= size(CO2); j++){N = N*CO2[j];}
      F = farey(G,N);
      if(F[1]-testF[1]==0)
      {
         if(printlevel>=10)
         {
            "size(L) = "+string(size(L));
         }
         F = cleardenom(F[1]);

         e = deg(F[1]);
         if(e == d)
         {
            list H = factorize(F[1]);

            int s = size(H[1]);
            for(i = 1; i <= s; i++)
            {
               if(H[2][i] != 1)
               {
                  int_break = 1;
               }
            }

            if(int_break == 0)
            {
               setring SPR;
               map phi = rHelp,var(nvars(SPR));
               list H = phi(H);
               if(printlevel>=10)
               {
                  "Real-time without test is "+string(rtimer - RT)+" seconds.";
                  "CPU-time without test is "+string(timer - T)+" seconds.";
               }
               T = timer;
               RT = rtimer;

               ideal F = phi(F);
               poly F1;

               if(n > 1)
               {
                  open(l(1));
                  write(l(1), quote(quickSubst(eval(F[1]), eval(f), eval(I))));
                  int t_sleep = timer;
               }
               else
               {
                  F1 = quickSubst(F[1],f,I);
                  if(F1 != 0) { int_break = 1; }
               }

               if(int_break == 0)
               {
                  for(i = 2; i <= s; i++)
                  {
                     H[1][i] = quickSubst(H[1][i],f,I);
                     Re[i-1] = I + ideal(H[1][i]);
                  }

                  if(n > 1)
                  {
                     t_sleep = timer - t_sleep;
                     if(t_sleep > 5) { t_sleep = 5; }

                     while(1)
                     {
                        if(status(l(1), "read", "ready"))
                        {
                           F1 = read(l(1));
                           close(l(1));
                           break;
                        }
                        i_sleep = system("sh", "sleep "+string(t_sleep));
                     }
                     if(F1 != 0) { int_break = 1; }
                  }
                  if(printlevel>=10)
                  {
                     "Real-time for test is "+string(rtimer - RT)+" seconds.";
                     "CPU-time for test is "+string(timer - T)+" seconds.";
                  }
                  if(int_break == 0)
                  {
                     kill f_for_fork;
                     kill I_for_fork;
                     kill SPR;
                     return(Re);
                  }
               }
            }
         }
      }

      int_break = 0;
      testF = F;
      CO1 = G;
      CO2 = N;
      index = 2;

      j = size(L) + 1;

      setring(SPR);
      L = primeList(I,10,L);
      setring rHelp;

      if(n > 1)
      {
         for(i = 1; i <= n; i++)
         {
            open(l(i));
            write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[j+i-1]),
                                                             eval(alg))));
         }
         j = j + n;
         k = 0;
      }
   }
}
example
{ "EXAMPLE:";  echo = 2;
   ring R=0,(a,b,c,d,e,f),dp;
   ideal I=
   2fb+2ec+d2+a2+a,
   2fc+2ed+2ba+b,
   2fd+e2+2ca+c+b2,
   2fe+2da+d+2cb,
   f2+2ea+e+2db+c2,
   2fa+f+2eb+2dc;
   assPrimes(I);
}
////////////////////////////////////////////////////////////////////////////////

static proc specialAlgDepEHV(poly p, ideal I)
{
//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
//=== mapping T to p
   def R = basering;
   execute("ring Rhelp="+charstr(R)+",T,dp;");
   setring R;
   map phi = Rhelp,p;
   setring Rhelp;
   ideal F = preimage(R,phi,I); //corresponds to std(I,p-T) in dp(n),dp(1)
   export(F);
   setring R;
   list L = Rhelp;
   return(L);
}

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

static proc specialAlgDepGTZ(poly p, ideal I)
{
//=== assume I is zero-dimensional
//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
//=== mapping T to p
   def R = basering;
   execute("ring Rhelp = "+charstr(R)+",T,dp;");
   setring R;
   map phi = Rhelp,p;
   def Rlp = changeord("dp("+string(nvars(R)-1)+"),dp(1)");
   setring Rlp;
   poly p = imap(R,p);
   ideal K = maxideal(1);
   K[nvars(R)] = 2*var(nvars(R))-p;
   map phi = R,K;
   ideal I = phi(I);
   I = std(I);
   poly q = subst(I[1],var(nvars(R)),var(1));
   setring Rhelp;
   map psi=Rlp,T;
   ideal F=psi(q);
   export(F);
   setring R;
   list L=Rhelp;
   return(L);
}

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

static proc specialAlgDepMonico(poly p, ideal I)
{
//=== assume I is zero-dimensional
//=== computes a poly F in Q[T], the characteristic polynomial of the map
//=== basering/I ---> baserng/I  defined by the multiplication with p
//=== in case I is radical it is the same poly as in specialAlgDepEHV
   def R = basering;
   execute("ring Rhelp = "+charstr(R)+",T,dp;");
   setring R;
   map phi = Rhelp,p;
   poly q;
   int j;
   matrix m ;
   poly va = var(1);
   ideal J = std(I);
   ideal ba = kbase(J);
   int d = vdim(J);
   matrix n[d][d];
   for(j = 2; j <= nvars(R); j++)
   {
     va = va*var(j);
   }
   for(j = 1; j <= d; j++)
   {
     q = reduce(p*ba[j],J);
     m = coeffs(q,ba,va);
     n[j,1..d] = m[1..d,1];
   }
   setring Rhelp;
   matrix n = imap(R,n);
   ideal F = det(n-T*freemodule(d));
   export(F);
   setring R;
   list L = Rhelp;
   return(L);
}

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

static proc specialTest(int d, poly p, ideal I)
{
//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
//=== mapping T to p and test if d=deg(F)
   def R = basering;
   execute("ring Rhelp="+charstr(R)+",T,dp;");
   setring R;
   map phi = Rhelp,p;
   setring Rhelp;
   ideal F = preimage(R,phi,I);
   int e=deg(F[1]);
   setring R;
   return((e==d));
}

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

static proc findGen(ideal J, list #)
{
//=== try to find a sparse linear form r such that
//=== vector space dim(basering/J)=deg(F),
//=== F a poly in Q[T] such that <F>=kernel(Q[T]--->basering) mapping T to r
//=== if not found returns a generic (randomly choosen) r
   int d = vdim(J);
   def R = basering;
   int n = nvars(R);
   list rl = ringlist(R);
   if(size(#) > 0) { int p = #[1]; }
   else { int p = prime(random(1000000000,2134567879)); }
   rl[1] = p;
   def @R = ring(rl);
   setring @R;
   ideal J = imap(R,J);
   poly r = var(n);
   int i,k;
   k = specialTest(d,r,J);
   if(!k)
   {
      for(i = 1; i < n; i++)
      {
         k = specialTest(d,r+var(i),J);
         if(k){ r = r + var(i); break; }
      }
   }
   if((!k) && (n > 2))
   {
      for(i = 1; i < n; i++)
      {
         r = r + var(i);
         k = specialTest(d,r,J);
         if(k){ break; }
      }
   }
   setring R;
   poly r = randomLast(100)[nvars(R)];
   if(k){ r = imap(@R,r); }
   return(r);
}

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

static proc pTestRad(int d, poly p1, ideal I)
{
//=== computes a poly F in Z/q1[T] such that
//===                    <F> = kernel(Z/q1[T]--->Z/q1[vars(basering)])
//=== mapping T to p1 and test if d=deg(squarefreepart(F)), q1 a prime randomly
//=== chosen
//=== If not choose randomly another prime q2 and another linear form p2 and
//=== computes a poly F in Z/q2[T] such that
//===                    <F> = kernel(Z/q2[T]--->Z/q2[vars(basering)])
//=== mapping T to p2 and test if d=deg(squarefreepart(F))
//=== if the test is positive then I is radical
   def R = basering;
   list rl = ringlist(R);
   int q1 = prime(random(100000000,536870912));
   rl[1] = q1;
   ring Shelp1 = q1,T,dp;
   setring R;
   def Rhelp1 = ring(rl);
   setring Rhelp1;
   poly p1 = imap(R,p1);
   ideal I = imap(R,I);
   map phi = Shelp1,p1;
   setring Shelp1;
   ideal F = preimage(Rhelp1,phi,I);
   poly f = gcd(F[1],diff(F[1],var(1)));
   int e = deg(F[1]/f);
   setring R;
   if(e != d)
   {
      poly p2 = findGen(I,q1);
      setring Rhelp1;
      poly p2 = imap(R,p2);
      phi = Shelp1,p2;
      setring Shelp1;
      F = preimage(Rhelp1,phi,I);
      f = gcd(F[1],diff(F[1],var(1)));
      e = deg(F[1]/f);
      setring R;
      if(e == d){ return(1); }
      if(e != d)
      {
         int q2 = prime(random(100000000,536870912));
         rl[1] = q2;
         ring Shelp2 = q2,T,dp;
         setring R;
         def Rhelp2 = ring(rl);
         setring Rhelp2;
         poly p1 = imap(R,p1);
         ideal I = imap(R,I);
         map phi = Shelp2,p1;
         setring Shelp2;
         ideal F = preimage(Rhelp2,phi,I);
         poly f = gcd(F[1],diff(F[1],var(1)));
         e = deg(F[1]/f);
         setring R;
         if(e == d){ return(1); }
      }
   }
   return((e==d));
}

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

static proc zeroRadP(ideal I, int p)
{
//=== computes F=(F_1,...,F_n) such that <F_i>=IZ/p[x_1,...,x_n] intersected
//=== with Z/p[x_i], F_i monic
   def R0 = basering;
   list ringL = ringlist(R0);
   ringL[1] = p;
   def @r = ring(ringL);
   setring @r;
   ideal I = fetch(R0,I);
   option(redSB);
   I = std(I);
   ideal F = finduni(I); //F[i] generates I intersected with K[var(i)]
   int i;
   for(i = 1; i <= size(F); i++){ F[i] = simplify(F[i],1); }
   setring R0;
   return(list(fetch(@r,F),p));
}

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

static proc quickSubst(poly h, poly r, ideal I)
{
//=== assume h is in Q[x_n], r is in Q[x_1,...,x_n], computes h(r) mod I
   attrib(I,"isSB",1);
   int n = nvars(basering);
   poly q = 1;
   int i,j,d;
   intvec v;
   list L;
   for(i = 1; i <= size(h); i++)
   {
      L[i] = list(leadcoef(h[i]),leadexp(h[i])[n]);
   }
   d = L[1][2];
   i = 0;
   h = 0;

   while(i <= d)
   {
      if(L[size(L)-j][2] == i)
      {
         h = reduce(h+L[size(L)-j][1]*q,I);
         j++;
      }
      q = reduce(q*r,I);
      i++;
   }
   return(h);
}

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

static proc modpSpecialAlgDep(list ringL, int p, list #)
{
//=== prepare parallel computing
//===  #=1: method of Eisenbud/Hunecke/Vasconcelos
//===  #=2: method of Gianni/Trager/Zacharias
//===  #=3: method of Monico

   def R0 = basering;

   ringL[1] = p;
   def @r = ring(ringL);
   setring @r;
   poly f = fetch(SPR,f_for_fork);
   ideal I = fetch(SPR,I_for_fork);
   if(size(#) > 0)
   {
      if(#[1] == 1) { list M = specialAlgDepEHV(f,I); }
      if(#[1] == 2) { list M = specialAlgDepGTZ(f,I); }
      if(#[1] == 3) { list M = specialAlgDepMonico(f,I); }
   }
   else
   {
      list M = specialAlgDepEHV(f,I);
   }
   def @S = M[1];

   setring R0;
   return(list(imap(@S,F),p));
}

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

/*
Examples:
=========

//=== Test for zeroR
ring R = 0,(x,y),dp;
ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;

//=== Cyclic_6
ring R = 0,x(1..6),dp;
ideal I = cyclic(6);

//=== Amrhein
ring R = 0,(a,b,c,d,e,f),dp;
ideal I = 2fb+2ec+d2+a2+a,
          2fc+2ed+2ba+b,
          2fd+e2+2ca+c+b2,
          2fe+2da+d+2cb,
          f2+2ea+e+2db+c2,
          2fa+f+2eb+2dc;

//=== Becker-Niermann
ring R = 0,(x,y,z),dp;
ideal I = x2+xy2z-2xy+y4+y2+z2,
          -x3y2+xy2z+xyz3-2xy+y4,
          -2x2y+xy4+yz4-3;

//=== Roczen
ring R = 0,(a,b,c,d,e,f,g,h,k,o),dp;
ideal I = o+1, k4+k, hk, h4+h, gk, gh, g3+h3+k3+1,
          fk, f4+f, eh, ef, f3h3+e3k3+e3+f3+h3+k3+1,
          e3g+f3g+g, e4+e, dh3+dk3+d, dg, df, de,
          d3+e3+f3+1, e2g2+d2h2+c, f2g2+d2k2+b,
          f2h2+e2k2+a;

//=== FourBodyProblem
//=== 4 bodies with equal masses, before symmetrisation.
//=== We are looking for the real positive solutions
ring R = 0,(B,b,D,d,F,f),dp;
ideal I = (b-d)*(B-D)-2*F+2,
          (b-d)*(B+D-2*F)+2*(B-D),
          (b-d)^2-2*(b+d)+f+1,
          B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;

//=== Reimer_5
ring R = 0,(x,y,z,t,u),dp;
ideal I = 2x2-2y2+2z2-2t2+2u2-1,
          2x3-2y3+2z3-2t3+2u3-1,
          2x4-2y4+2z4-2t4+2u4-1,
          2x5-2y5+2z5-2t5+2u5-1,
          2x6-2y6+2z6-2t6+2u6-1;

//=== ZeroDim.example_12
ring R = 0, (x, y, z), lp;
ideal I = 7xy+x+3yz+4y+2z+10,
          x3+x2y+3xyz+5xy+3x+2y3+6y2z+yz+1,
          3x4+2x2y2+3x2y+4x2z2+xyz+xz2+6y2z+5z4;

//=== ZeroDim.example_27
ring R = 0, (w, x, y, z), lp;
ideal I = -2w2+9wx+9wy-7wz-4w+8x2+9xy-3xz+8x+6y2-7yz+4y-6z2+8z+2,
          3w2-5wx-3wy-6wz+9w+4x2+2xy-2xz+7x+9y2+6yz+5y+7z2+7z+5,
          7w2+5wx+3wy-5wz-5w+2x2+9xy-7xz+4x-4y2-5yz+6y-4z2-9z+2,
          8w2+5wx-4wy+2wz+3w+5x2+2xy-7xz-7x+7y2-8yz-7y+7z2-8z+8;

//=== Cassou_1
ring R = 0, (b,c,d,e), dp;
ideal I = 6b4c3+21b4c2d+15b4cd2+9b4d3+36b4c2+84b4cd+30b4d2+72b4c+84b4d-8b2c2e
          -28b2cde+36b2d2e+48b4-32b2ce-56b2de-144b2c-648b2d-32b2e-288b2-120,
          9b4c4+30b4c3d+39b4c2d2+18b4cd3+72b4c3+180b4c2d+156b4cd2+36b4d3+216b4c2
          +360b4cd+156b4d2-24b2c3e-16b2c2de+16b2cd2e+24b2d3e+288b4c+240b4d
          -144b2c2e+32b2d2e+144b4-432b2c2-720b2cd-432b2d2-288b2ce+16c2e2-32cde2
          +16d2e2-1728b2c-1440b2d-192b2e+64ce2-64de2-1728b2+576ce-576de+64e2
          -240c+1152e+4704,
          -15b2c3e+15b2c2de-90b2c2e+60b2cde-81b2c2+216b2cd-162b2d2-180b2ce
          +60b2de+40c2e2-80cde2+40d2e2-324b2c+432b2d-120b2e+160ce2-160de2-324b2
          +1008ce-1008de+160e2+2016e+5184,
          -4b2c2+4b2cd-3b2d2-16b2c+8b2d-16b2+22ce-22de+44e+261;

================================================================================

The following timings are conducted on an Intel Xeon X5460 with 4 CPUs, 3.16 GHz
each, 64 GB RAM under the Gentoo Linux operating system by using the 32-bit
version of Singular 3-1-1.
The results of the timinings are summarized in the following table where
assPrimes* denotes the parallelized version of the algorithm on 4 CPUs and
 (1) approach of Eisenbud, Hunecke, Vasconcelos (cf. specialAlgDepEHV),
 (2) approach of Gianni, Trager, Zacharias      (cf. specialAlgDepGTZ),
 (3) approach of Monico                         (cf. specialAlgDepMonico).

Example             | minAssGTZ |     assPrimes         assPrimes*
                    |           |  (1)   (2)   (3)   (1)   (2)   (3)
--------------------------------------------------------------------
Cyclic_6            |      5    |    5    10     7     4     7     6
Amrhein             |      1    |    3     3     5     1     2    21
Becker-Niermann     |      -    |    0     0     1     0     0     0
Roczen              |      0    |    3     2     0     2     4     1
FourBodyProblem     |      -    |  139   139   148    96    83    96
Reimer_5            |      -    |  132   128   175    97    70   103
ZeroDim.example_12  |    170    |  125   125   125    67    68    63
ZeroDim.example_27  |     27    |  215   226   215   113   117   108
Cassou_1            |    525    |  112   112   112    56    56    57

minAssGTZ (cf. primdec_lib) runs out of memory for Becker-Niermann,
FourBodyProblem and Reimer_5.

================================================================================

//=== One component at the origin

ring R = 0, (x,y), dp;
poly f1 = (y5 + y4x7 + 2x8);
poly f2 = (y3 + 7x4);
poly f3 = (y7 + 2x12);
poly f = f1*f2*f3 + y19;
ideal I = f, diff(f, x), diff(f, y);

ring R = 0, (x,y), dp;
poly f1 = (y5 + y4x7 + 2x8);
poly f2 = (y3 + 7x4);
poly f3 = (y7 + 2x12);
poly f4 = (y11 + 2x18);
poly f = f1*f2*f3*f4 + y30;
ideal I = f, diff(f, x), diff(f, y);

ring R = 0, (x,y), dp;
poly f1 = (y15 + y14x7 + 2x18);
poly f2 = (y13 + 7x14);
poly f3 = (y17 + 2x22);
poly f = f1*f2*f3 + y49;
ideal I = f, diff(f, x), diff(f, y);

ring R = 0, (x,y), dp;
poly f1 = (y15 + y14x20 + 2x38);
poly f2 = (y19 + 3y17x50 + 7x52);
poly f = f1*f2 + y36;
ideal I = f, diff(f, x), diff(f, y);

//=== Several components

ring R = 0, (x,y), dp;
poly f1 = (y5 + y4x7 + 2x8);
poly f2 = (y13 + 7x14);
poly f = f1*subst(f2, x, x-3, y, y+5);
ideal I = f, diff(f, x), diff(f, y);

ring R = 0, (x,y), dp;
poly f1 = (y5  + y4x7 + 2x8);
poly f2 = (y3 + 7x4);
poly f3 = (y7 + 2x12);
poly f = f1*f2*subst(f3, y, y+5);
ideal I = f, diff(f, x), diff(f, y);

ring R = 0, (x,y), dp;
poly f1 = (y5  + 2x8);
poly f2 = (y3 + 7x4);
poly f3 = (y7 + 2x12);
poly f = f1*subst(f2,x,x-1)*subst(f3, y, y+5);
ideal I = f, diff(f, x), diff(f, y);
*/