////////////////////////////////////////////////////////////////////////////////
version="$Id$";
category = "Commutative Algebra";
info="
LIBRARY:  assPrimes.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.

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

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

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

proc zeroR(ideal I, list #)
"USAGE:  zeroR(I,[n]); I ideal, optional: n number of processors (for parallel computing)
ASSUME:  I is zero-dimensional in Q[variables]
RETURN:  the radical of I
EXAMPLE: example zeroR; shows an example
"
{
   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  -------------------------
// nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10)
   intvec L = primeList(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"))
               {
                  P = read(l(i));                                        // read the result from 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));
                  }
               }
            }
            if(k == n)                                                   // k describes the number of closed links
            {
               j++;
            }
            i_sleep = system("sh", "sleep "+string(t));                  // sleep for t seconds
         }
      }
      else
      {
         while(j<=size(L))
         {
            P = zeroRadP(I, L[j]);
            CO1[index] = P[1];
            CO2[index] = bigint(P[2]);
            index++;
            j++;
         }
      }

      //hier deleteUnlucky einbauen, streichen, wenn der Grad nicht stimmt
      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;

// nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10,L)
         L = primeList(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(std(I + F)); }
}
example
{ "EXAMPLE:";  echo = 2;
   ring R = 0, (x,y), dp;
   ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
   zeroR(I);
}

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

proc assPrimes(list #)
"USAGE:  assPrimes(I,[n],[a]); I ideal or module,
         optional: n number of processors (for parallel computing), a
         - a=1: method of Eisenbud/Hunecke/Vasconcelos
         - a=2: method of Gianni/Trager/Zacharias
         - a=3: mathod of Monico
         assPrimes(I)  chooses n=a=1
ASSUME:  I is zero-dimensional over Q[variables]
RETURN:  a list pr of associated primes of I:
EXAMPLE: example assPrimes; shows an example
"

{
   int T = timer;
   int RT = rtimer;

   ideal I;
   if(typeof(#[1])=="ideal")
   {
      I = #[1];
   }
   else
   {
      module M = #[1];
      I = Ann(M);
   }

   def SPR = basering;
   I = std(I);
   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);
   if(size(f) == nvars(SPR))
   {
      int spT = pTestRad(d,f,I);
      if(printlevel>=10)
      {
         "pTestRad(d,f,I) = "+string(spT);
      }
      if(!spT)
      {
         if(typeof(attrib(#[1],"isRad")) == "int")
         {
            if(attrib(#[1],"isRad") == 0)
            {
               I = zeroR(I);
               I = std(I);
               d = vdim(I);
               f = findGen(I);
            }
         }
         else
         {
            I = zeroR(I);
            I = std(I);
            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.";
   }
//---------------------  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;
   }

   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

   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;


//--------------------  Initialize the list of primes  -------------------------
// nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10)
   intvec L = primeList(10);
   L[5] = prime(random(1000000000,2134567879));

//-----  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++)
            {
               if(status(l(i), "read", "ready"))                         // ask if link l(i) is ready otherwise sleep for t seconds
               {
                  P = read(l(i));                                        // read the result from 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));
                  }
               }
            }
            if(k == n)                                                  // k describes the number of closed links
            {
               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.";
      }
      //hier deleteUnlucky einbauen, streichen, wenn der Grad nicht stimmt
      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);

                     if(typeof(#[1])=="ideal")
                     {
                        Re[i-1] = #[1] + ideal(H[1][i]);
                     }
                     else
                     {
                        Re[i-1] = M + H[1][i]*freemodule(nrows(M));
                     }
                  }

                  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;

// nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10,L)
      L = primeList(10,L);

      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(baserng/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; }
      }
   }
   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 choosen
//=== 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 primeList(int n, list #)
"USAGE:  primeList(n); (resp. primeList(n,L); )
RETURN: the intvec of n greatest primes  <= 536870912 (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 = prime(536870912);
      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)];
}

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

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: mathod 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);

//Amrhei
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;

//4 bodies with equal masses, before symmetrisation.
//We are looking for the real positive solutions
ring R=0,(B,D,F,b,d,f),dp;

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;

ring R=0,(B,b,D,d,F,f),dp;
ideal I=(b-d)*(B-D)-F+1,
(b-d)*(B+D-F)+(B-D),
(b-d)^2-(b+d)+f+1,
B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;


//prime
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;

*/