singular-data 4.0.3+ds-1 / usr / share / singular / LIB / assprimeszerodim.lib

///////////////////////////////////////////////////////////////////////////////
version="version assprimeszerodim.lib 4.0.0.0 Jun_2014 "; // $Id: 4c9703d9d449016b92ecb5f8786539f6d04d29ed $
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
          A. Steenpass    steenpass@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[, exactness]); I ideal, exactness int
ASSUME:  I is zero-dimensional in Q[variables]
RETURN:  the radical of I
NOTE:    A final test is applied to the result if exactness != 0 (default),
         otherwise no final test is done.
EXAMPLE: example zeroRadical; shows an example
"
{
   /* read optional parameter */
   int exactness = 1;
   if(size(#) > 0)
   {
      if(size(#) > 1 || typeof(#[1]) != "int")
      {
         ERROR("wrong optional parameter");
      }
      exactness = #[1];
   }

   /* compute a standard basis if necessary */
   if (!attrib(I, "isSB"))
   {
      I = modStd(I, exactness);
   }

   /* call modular() */
   // TODO: write deleteUnluckyPrimes_zeroRadical()
   if(exactness)
   {
      ideal F = modular("Assprimeszerodim::zeroRadP", list(I),
         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
         pTest_zeroRadical, finalTest_zeroRadical);
   }
   else
   {
      ideal F = modular("Assprimeszerodim::zeroRadP", list(I),
         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
         pTest_zeroRadical);
   }

   /* compute the squarefree parts */
   poly f;
   int k;
   int i;
   for(i = nvars(basering); i > 0; i--)
   {
      f = gcd(F[i], diff(F[i], var(i)));
      k = k + deg(f);
      F[i] = F[i]/f;
   }

   /* return the result */
   if(k == 0)
   {
      return(I);
   }
   else
   {
      return(modStd(I + F, exactness));
   }
}
example
{ "EXAMPLE:";  echo = 2;
   ring R = 0, (x,y), dp;
   ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
   zeroRadical(I);
}

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

/* The pTest for zeroRadical(), to be used in modular(). */
static proc pTest_zeroRadical(string command, list args, ideal result, int p)
{
   /* change to characteristic p */
   def br = basering;
   list lbr = ringlist(br);
   if(typeof(lbr[1]) == "int")
   {
      lbr[1] = p;
   }
   else
   {
      lbr[1][1] = p;
   }
   def rp = ring(lbr);
   setring(rp);
   ideal Ip = fetch(br, args)[1];
   ideal Fp_result = fetch(br, result);

   /* run the command and compare with given result */
   execute("ideal Fp = "+command+"(Ip);");
   int i;
   for(i = nvars(br); i > 0; i--)
   {
      if(Fp[i] != Fp_result[i])
      {
         setring(br);
         return(0);
      }
   }
   setring(br);
   return(1);
}

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

/* The finalTest for zeroRadical, to be used in modular(). */
static proc finalTest_zeroRadical(string command, list args, ideal F)
{
   int i;
   for(i = nvars(basering); i > 0; i--)
   {
      if(reduce(F[i], args[1]) != 0) { return(0); }
   }
   return(1);
}

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

proc assPrimes(def I, list #)
"USAGE:  assPrimes(I[, alg, exactness]); I ideal or module,
         alg string (optional), exactness int (optional)
         - alg = "GTZ":    method of Gianni/Trager/Zacharias (default)
         - alg = "EHV":    method of Eisenbud/Hunecke/Vasconcelos
         - alg = "Monico": method of Monico
ASSUME:  I is zero-dimensional over Q[variables]
RETURN:  a list of the associated primes of I
NOTE:    A final test is applied to the result if exactness != 0 (default),
         otherwise no final test is done.
EXAMPLE: example assPrimes; shows an example
"
{
   /* read input */
   if(typeof(I) != "ideal")
   {
      if(typeof(I) != "module")
      {
         ERROR("The first argument must be of type 'ideal' or 'module'.");
      }
      module M = I;
      kill I;
      ideal I = Ann(M);
      kill M;
   }

   /* read optional parameters */
   list defaults = list("GTZ", 1);
   int i;
   for(i = 1; i <= size(defaults); i++)
   {
      if(typeof(#[i]) != typeof(defaults[i]))
      {
         # = insert(#, defaults[i], i-1);
      }
   }
   if(size(#) != size(defaults))
   {
      ERROR("wrong optional parameters");
   }
   string alg = #[1];
   int exactness = #[2];
   int a;
   if(alg == "GTZ")
   {
      a = 1;
   }
   if(alg == "EHV")
   {
      a = 2;
   }
   if(alg == "Monico")
   {
      a = 3;
   }
   if(a == 0)   // alg != "GTZ" && alg != "EHV" && alg != "Monico"
   {
      ERROR("unknown method");
   }

   /* compute a standard basis if necessary */
   if(printlevel >= 10) { "========== Start modStd =========="; }
   if (!attrib(I, "isSB")) {
       I = modStd(I, exactness);
   }
   if(printlevel >= 10) { "=========== End modStd ==========="; }
   int d = vdim(I);
   if(d == -1) { ERROR("Ideal is not zero-dimensional."); }
   if(homog(I) == 1) { return(list(maxideal(1))); }

   /* compute the radical if necessary */
   ideal J = I;
   int isRad;
   poly f;
   isRad, f = pTestRad(I, d);
   while(!isRad)
   {
      J = zeroRadical(I, exactness);
      J = modStd(J, exactness);
      d = vdim(J);
      isRad, f = pTestRad(J, d);
   }
   I = J;
   kill J;

   /* call modular() */
   if(exactness)
   {
      ideal F = modular("Assprimeszerodim::modpSpecialAlgDep",
         list(I, f, d, a),
         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
         pTest_assPrimes, finalTest_assPrimes);
   }
   else
   {
      ideal F = modular("Assprimeszerodim::modpSpecialAlgDep",
         list(I, f, d, a),
         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
         pTest_assPrimes);
   }

   /* compute the components */
   list result;
   list H = factorize(F[1]);
   for(i = size(H[1])-1; i > 0; i--)
   {
      result[i] = I + ideal(quickSubst(H[1][i+1], f, I));
   }

   /* return the result */
   return(result);
}
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);
}

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

/* Computes a poly F in Q[T] such that
 * <F> = kernel(Q[T] --> basering, T |-> f),
 * T := last variable in the basering.
 */
static proc specialAlgDepEHV(ideal I, poly f)
{
   def R = basering;
   execute("ring QT = ("+charstr(R)+"), "+varstr(R, nvars(R))+", dp;");
   setring(R);
   map phi = QT, f;
   setring QT;
   ideal F = preimage(R, phi, I); // corresponds to std(I, f-T) in dp(n),dp(1)
   setring(R);
   ideal F = imap(QT, F);
   return(F);
}

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

/* Assume I is zero-dimensional.
 * Computes a poly F in Q[T] such that
 * <F> = kernel(Q[T] --> basering, T |-> f),
 * T := last variable in the basering.
 */
static proc specialAlgDepGTZ(ideal I, poly f)
{
   def R = basering;
   if(nvars(R) > 1)
   {
      def Rlp = changeord(list(list("dp", 1:(nvars(R)-1)), list("dp", 1:1)));
      setring(Rlp);
      poly f = imap(R, f);
      ideal I;
   }
   ideal K = maxideal(1);
   K[nvars(R)] = 2*var(nvars(R))-f;
   map phi = R, K;
   I = phi(I);
   I = std(I);
   ideal F = I[1];
   if(nvars(R) > 1)
   {
      setring(R);
      ideal F = imap(Rlp, F);
   }
   return(F);
}

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

/* Assume I is zero-dimensional.
 * Computes a poly F in Q[T], the characteristic polynomial of the map
 * basering/I ---> basering/I  defined by the multiplication with f,
 * T := last variable in the basering.
 * In case I is radical, it is the same polynomial as in specialAlgDepEHV.
 */
static proc specialAlgDepMonico(ideal I, poly f, int d)
{
   def R = basering;
   int j;
   matrix M[d][d];
   ideal J = std(I);
   ideal basis = kbase(J);
   poly vars = var(nvars(R));
   for(j = nvars(R)-1; j > 0; j--)
   {
     vars = var(j)*vars;
   }
   for(j = 1; j <= d; j++)
   {
     M[1..d, j] = coeffs(reduce(f*basis[j], J), basis, vars);
   }
   execute("ring QT = ("+charstr(R)+"), "+varstr(R, nvars(R))+", dp;");
   matrix M = imap(R, M);
   ideal F = det(M-var(1)*freemodule(d));
   setring(R);
   ideal F = imap(QT, F);
   return(F);
}

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

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), fetch(Rhelp, F)[1]);
}

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

/* Assume d = vector space dim(basering/J).
 * Tries to find a (sparse) linear form r such that d = deg(F), where
 * F is a poly in Q[T] such that <F> = kernel(Q[T]-->basering) mapping T to r.
 * If found, returns (1, r, F). If not found, returns (0, 0, 0).
 */
static proc findGen(ideal J, int d)
{
   def R = basering;
   int n = nvars(R);
   int okay;
   poly F;

   /* try trivial transformation */
   poly r = var(n);
   okay, F = specialTest(d, r, J);
   if(okay)
   {
      return(1, r, F);
   }

   /* try transformations of the form var(n) + var(i) */
   int i;
   for(i = 1; i < n; i++)
   {
      okay, F = specialTest(d, r+var(i), J);
      if(okay)
      {
         return(1, r+var(i), F);
      }
   }

   /* try transformations of the form var(n) + \sum var(i) */
   if(n > 2)
   {
      for(i = 1; i < n; i++)
      {
         r = r + var(i);
         okay, F = specialTest(d, r, J);
         if(okay)
         {
            return(1, r, F);
         }
      }
   }

   /* try random transformations */
   int N = 2;   // arbitrarily chosen
   for(i = N; i > 0; i--)
   {
      r = randomLast(100)[n];
      okay, F = specialTest(d, r, J);
      if(okay)
      {
         return(1, r, F);
      }
   }

   /* not found */
   return(0, 0, 0);
}

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

/* Assume d = vector space dim(basering/I).
 * Tests if I is radical over F_p, where p is some randomly chosen prime.
 * If yes, chooses a linear form r such that d = deg(squarefreepart(F)), where
 * F is a poly in Z/p[T] such that <F> = kernel(Z/p[T]-->Z/p[vars(basering)])
 * mapping T to r.
 * Returns (1, r), if I is radical over F_p, and (0, 0) otherwise.
 */
static proc pTestRad(ideal I, int d)
{
   int N = 2;   // Try N random primes. Value of N can be chosen arbitrarily.
   def R = basering;
   list rl = ringlist(R);
   int p;
   int okay;
   int i;
   for(i = N; i > 0; i--)
   {
      p = prime(random(100000000,536870912));

      // change to characteristic p
      if(typeof(rl[1]) == "int")
      {
         rl[1] = p;
      }
      else
      {
         rl[1][1] = p;
      }
      def Rp(i) = ring(rl);
      setring Rp(i);
      ideal I = imap(R, I);

      // find and test transformation
      poly r;
      poly F;
      okay, r, F = findGen(I, d);
      if(okay)
      {
         poly f = gcd(F, diff(F, var(1)));
         if(d == deg(F/f))   // F squarefree?
         {
            setring(R);
            return(1, imap(Rp(i), r));
         }
      }
      setring(R);
   }
   return(0, 0);
}

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

/* Computes an ideal F such that ncols(F) = nvars(basering),
 * < F[i] > = (I intersected with K[var(i)]), and F[i] is monic.
 */
static proc zeroRadP(ideal I)
{
   intvec opt = option(get);
   option(redSB);
   I = std(I);
   ideal F = finduni(I);   // F[i] generates I intersected with K[var(i)]
   F = simplify(F, 1);
   option(set, opt);
   return(F);
}

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

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);
}

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

/* Simple switch for specialAlgDepEHV, specialAlgDepGTZ, and
 * specialAlgDepMonico.
 */
static proc modpSpecialAlgDep(ideal I, poly f, int d, int alg)
{
   ideal F;
   if(alg == 1) { F = specialAlgDepEHV(I, f); }
   if(alg == 2) { F = specialAlgDepGTZ(I, f); }
   if(alg == 3) { F = specialAlgDepMonico(I, f, d); }
   F = simplify(F, 1);
   return(F);
}

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

/* The pTest for assPrimes(), to be used in modular(). */
static proc pTest_assPrimes(string command, list args, ideal F, int p)
{
   def br = basering;
   list lbr = ringlist(br);
   if(typeof(lbr[1]) == "int")
   {
      lbr[1] = p;
   }
   else
   {
      lbr[1][1] = p;
   }
   def rp = ring(lbr);
   setring(rp);
   list args_p = fetch(br, args);
   ideal F = fetch(br, F);
   execute("ideal Fp = "+command+"("
      +Tasks::argsToString("args_p", size(args_p))+");");
   int k = (Fp[1] == F[1]);
   setring br;
   return(k);
}

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

/* The finalTest for assPrimes(), to be used in modular(). */
static proc finalTest_assPrimes(string command, alias list args, ideal F)
{
   F = cleardenom(F[1]);
   if(deg(F[1]) != args[3]) { return(0); }
   if(gcd(F[1], diff(F[1], var(nvars(basering)))) != 1) { return(0); };
   if(quickSubst(F[1], args[2], args[1]) != 0) { return(0); }
   return(1);
}

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

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