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

///////////////////////////////////////////////////////////////////////////////
version="version modstd.lib 4.0.2.0 Aug_2015 "; // $Id: 4b6b912051bdc20f10ea65b7c6b57406b17d96fd $
category="Commutative Algebra";
info="
LIBRARY:  modstd.lib      Groebner bases of ideals using modular methods

AUTHORS:  A. Hashemi      Amir.Hashemi@lip6.fr
          G. Pfister      pfister@mathematik.uni-kl.de
          H. Schoenemann  hannes@mathematik.uni-kl.de
          A. Steenpass    steenpass@mathematik.uni-kl.de
          S. Steidel      steidel@mathematik.uni-kl.de

OVERVIEW:
  A library for computing Groebner bases of ideals in the polynomial ring over
  the rational numbers using modular methods.

REFERENCES:
  E. A. Arnold: Modular algorithms for computing Groebner bases.
  J. Symb. Comp. 35, 403-419 (2003).

  N. Idrees, G. Pfister, S. Steidel: Parallelization of Modular Algorithms.
  J. Symb. Comp. 46, 672-684 (2011).

PROCEDURES:
  modStd(I);    standard basis of I using modular methods
";

LIB "poly.lib";
LIB "modular.lib";

proc modStd(ideal I, list #)
"USAGE:   modStd(I[, exactness]); I ideal, exactness int
RETURN:   a standard basis of I
NOTE:     The procedure computes a standard basis of I (over the rational
          numbers) by using modular methods.
       @* An optional parameter 'exactness' can be provided.
          If exactness = 1, the procedure computes a standard basis of I for
          sure; if exactness = 0, it computes a standard basis of I
          with high probability.
SEE ALSO: modular
EXAMPLE:  example modStd; shows an example"
{
    /* read optional parameter */
    int exactness = 1;
    if (size(#) > 0) {
        /* For compatibility, we only test size(#) > 4. This can be changed to
         * size(#) > 1 in the future. */
        if (size(#) > 4 || typeof(#[1]) != "int") {
            ERROR("wrong optional parameter");
        }
        exactness = #[1];
    }

    /* save options */
    intvec opt = option(get);
    option(redSB);

    /* choose the right command */
    string command = "groebner";
    if (npars(basering) > 0) {
        command = "Modstd::groebner_norm";
    }

    /* call modular() */
    if (exactness) {
        I = modular(command, list(I), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }
    else {
        I = modular(command, list(I), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std);
    }

    /* return the result */
    attrib(I, "isSB", 1);
    option(set, opt);
    return(I);
}
example
{
    "EXAMPLE:";
    echo = 2;
    ring R1 = 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 R2 = 0, (x,y,z), ds;
    ideal I = jacob(x5+y6+z7+xyz);
    ideal J = modStd(I, 0);
    J;

    ring R3 = 0, x(1..4), lp;
    ideal I = cyclic(4);
    ideal J1 = modStd(I, 1);   // default
    ideal J2 = modStd(I, 0);
    size(reduce(J1, J2));
    size(reduce(J2, J1));
}

/* compute a normalized GB via groebner() */
static proc groebner_norm(ideal I)
{
    I = simplify(groebner(I), 1);
    attrib(I, "isSB", 1);
    return(I);
}

/* test if the prime p is suitable for the input, i.e. it does not divide
 * the numerator or denominator of any of the coefficients */
static proc primeTest_std(int p, alias list args)
{
    /* erase zero generators */
    ideal I = simplify(args[1], 2);

    /* clear denominators and count the terms */
    ideal J;
    ideal K;
    int n = ncols(I);
    intvec sizes;
    number cnt;
    int i;
    for(i = n; i > 0; i--) {
        J[i] = cleardenom(I[i]);
        cnt = leadcoef(J[i])/leadcoef(I[i]);
        K[i] = numerator(cnt)*var(1)+denominator(cnt);
    }
    sizes = size(J[1..n]);

    /* 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 Jp = fetch(br, J);
    ideal Kp = fetch(br, K);

    /* test if any coefficient is missing */
    if (intvec(size(Kp[1..n])) != 2:n) {
        setring(br);
        return(0);
    }
    if (intvec(size(Jp[1..n])) != sizes) {
        setring(br);
        return(0);
    }
    setring(br);
    return(1);
}

/* find entries in modresults which come from unlucky primes.
 * For this, sort the entries into categories depending on their leading
 * ideal and return the indices in all but the biggest category. */
static proc deleteUnluckyPrimes_std(alias list modresults)
{
    int size_modresults = size(modresults);

    /* sort results into categories.
     * each category is represented by three entries:
     * - the corresponding leading ideal
     * - the number of elements
     * - the indices of the elements
     */
    list cat;
    int size_cat;
    ideal L;
    int i;
    int j;
    for (i = 1; i <= size_modresults; i++) {
        L = lead(modresults[i]);
        attrib(L, "isSB", 1);
        for (j = 1; j <= size_cat; j++) {
            if (size(L) == size(cat[j][1])
                && size(reduce(L, cat[j][1])) == 0
                && size(reduce(cat[j][1], L)) == 0) {
                cat[j][2] = cat[j][2]+1;
                cat[j][3][cat[j][2]] = i;
                break;
            }
        }
        if (j > size_cat) {
            size_cat++;
            cat[size_cat] = list();
            cat[size_cat][1] = L;
            cat[size_cat][2] = 1;
            cat[size_cat][3] = list(i);
        }
    }

    /* find the biggest categories */
    int cat_max = 1;
    int max = cat[1][2];
    for (i = 2; i <= size_cat; i++) {
        if (cat[i][2] > max) {
            cat_max = i;
            max = cat[i][2];
        }
    }

    /* return all other indices */
    list unluckyIndices;
    for (i = 1; i <= size_cat; i++) {
        if (i != cat_max) {
            unluckyIndices = unluckyIndices + cat[i][3];
        }
    }
    return(unluckyIndices);
}

/* test if 'command' applied to 'args' in characteristic p is the same as
   'result' mapped to characteristic p */
static proc pTest_std(string command, alias list args, alias 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 Gp = fetch(br, result);
    attrib(Gp, "isSB", 1);

    /* test if Ip is in Gp */
    int i;
    for (i = ncols(Ip); i > 0; i--) {
        if (reduce(Ip[i], Gp, 1) != 0) {
            setring(br);
            return(0);
        }
    }

    /* compute command(args) */
    execute("Ip = "+command+"(Ip);");

    /* test if Gp is in Ip */
    for (i = ncols(Gp); i > 0; i--) {
        if (reduce(Gp[i], Ip, 1) != 0) {
            setring(br);
            return(0);
        }
    }
    setring(br);
    return(1);
}

/* test if 'result' is a GB of the input ideal */
static proc finalTest_std(string command, alias list args, ideal result)
{
    /* test if args[1] is in result */
    attrib(result, "isSB", 1);
    int i;
    for (i = ncols(args[1]); i > 0; i--) {
        if (reduce(args[1][i], result, 1) != 0) {
            return(0);
        }
    }

    /* test if result is a GB */
    ideal G = std(result);
    if (reduce_parallel(G, result)) {
        return(0);
    }
    return(1);
}

/* return 1, if I_reduce is _not_ in G_reduce,
 *        0, otherwise
 * (same as size(reduce(I_reduce, G_reduce))).
 * Uses parallelization. */
static proc reduce_parallel(def I_reduce, def G_reduce)
{
    exportto(Modstd, I_reduce);
    exportto(Modstd, G_reduce);
    int size_I = ncols(I_reduce);
    int chunks = Modular::par_range(size_I);
    intvec range;
    int i;
    for (i = chunks; i > 0; i--) {
        range = Modular::par_range(size_I, i);
        task t(i) = "Modstd::reduce_task", list(range);
    }
    startTasks(t(1..chunks));
    waitAllTasks(t(1..chunks));
    int result = 0;
    for (i = chunks; i > 0; i--) {
        if (getResult(t(i))) {
            result = 1;
            break;
        }
    }
    kill I_reduce;
    kill G_reduce;
    return(result);
}

/* compute a chunk of reductions for reduce_parallel */
static proc reduce_task(intvec range)
{
    int result = 0;
    int i;
    for (i = range[1]; i <= range[2]; i++) {
        if (reduce(I_reduce[i], G_reduce, 1) != 0) {
            result = 1;
            break;
        }
    }
    return(result);
}

////////////////////////////////////////////////////////////////////////////////
/*
 * The following procedures are kept for backward compatibility with the old
 * version of modstd.lib. As of now (May 2014), they are still needed in
 * modnormal.lib, modwalk.lib, and symodstd.lib. They can be removed here as
 * soon as they are not longer needed in these libraries.
 */

LIB "parallel.lib";

static proc mod_init()
{
   newstruct("idealPrimeTest", "ideal Ideal");
}

static proc redFork(ideal I, ideal J, int n)
{
   attrib(J,"isSB",1);
   return(reduce(I,J,1));
}

proc isIncluded(ideal I, ideal J, list #)
"USAGE:  isIncluded(I,J); I,J ideals
RETURN:  1 if J includes I,
@*       0 if there is an element f in I which does not reduce to 0 w.r.t. J.
EXAMPLE: example isIncluded; shows an example
"
{
   def R = basering;
   setring R;

   attrib(J,"isSB",1);
   int i,j,k;

   if(size(#) > 0)
   {
      int n = #[1];
      if(n >= ncols(I)) { n = ncols(I); }
      if(n > 1)
      {
         for(i = 1; i <= n - 1; i++)
         {
            //link l(i) = "MPtcp:fork";
            link l(i) = "ssi:fork";
            open(l(i));

            write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1)));
         }

         int t = timer;
         if(reduce(I[ncols(I)], J, 1) != 0)
         {
            for(i = 1; i <= n - 1; i++)
            {
               close(l(i));
            }
            return(0);
         }
         t = timer - t;
         if(t > 60) { t = 60; }
         int i_sleep = system("sh", "sleep "+string(t));

         j = ncols(I) - n;

         while(j >= 0)
         {
            for(i = 1; i <= n - 1; i++)
            {
               if(status(l(i), "read", "ready"))
               {
                  if(read(l(i)) != 0)
                  {
                     for(i = 1; i <= n - 1; i++)
                     {
                        close(l(i));
                     }
                     return(0);
                  }
                  else
                  {
                     if(j >= 1)
                     {
                        write(l(i), quote(redFork(eval(I[j]), eval(J), 1)));
                        j--;
                     }
                     else
                     {
                        k++;
                        close(l(i));
                     }
                  }
               }
            }
            if(k == n - 1)
            {
               j--;
            }
            i_sleep = system("sh", "sleep "+string(t));
         }
         return(1);
      }
   }

   for(i = ncols(I); i >= 1; i--)
   {
      if(reduce(I[i],J,1) != 0){ return(0); }
   }
   return(1);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal I = x+1,x+y+1;
   ideal J = x+1,y;
   isIncluded(I,J);
   isIncluded(J,I);
   isIncluded(I,J,4);

   ring R = 0, x(1..5), dp;
   ideal I1 = cyclic(4);
   ideal I2 = I1,x(5)^2;
   isIncluded(I1,I2,4);
}

proc deleteUnluckyPrimes(list T, list L, int ho, list #)
"USAGE:  deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer
RETURN:  lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #),
         lT ideal
NOTE:    - if ho = 1, the polynomials in T are homogeneous, else ho = 0,
@*       - lT is prevalent, i.e. the most appearing leading ideal in T
EXAMPLE: example deleteUnluckyPrimes; shows an example
"
{
   ho = ((ho)||(ord_test(basering) == -1));
   int j,k,c;
   intvec hl,hc;
   ideal cT,lT,cK;
   lT = lead(T[size(T)]);
   attrib(lT,"isSB",1);
   if(!ho)
   {
      for(j = 1; j < size(T); j++)
      {
         cT = lead(T[j]);
         attrib(cT,"isSB",1);
         if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
         {
            cK = cT;
            c++;
         }
      }
      if(c > size(T) div 2){ lT = cK; }
   }
   else
   {
      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; c++; break; }
               if(lT[k] > cT[k]) { c++; break; }
            }
         }
         else
         {
            if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; }
         }
      }
   }

   int addList;
   if(size(#) > 0) { list M = #; addList = 1; }
   j = 1;
   attrib(lT,"isSB",1);
   while((j <= size(T))&&(c > 0))
   {
      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)];
            if(addList == 1) { M = M[2..size(M)]; }
         }
         else
         {
            if(j == size(L))
            {
               L = L[1..size(L)-1];
               if(addList == 1) { M = M[1..size(M)-1]; }
            }
            else
            {
               L = L[1..j-1],L[j+1..size(L)];
               if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; }
            }
         }
         j--;
      }
      j++;
   }

   for(j = 1; j <= size(L); j++)
   {
      L[j] = bigint(L[j]);
   }

   if(addList == 0) { return(list(T,L,lT)); }
   if(addList == 1) { return(list(T,L,M,lT)); }
}
example
{ "EXAMPLE:"; echo = 2;
   list L = 2,3,5,7,11;
   ring r = 0,(y,x),Dp;
   ideal I1 = 2y2x,y6;
   ideal I2 = yx2,y3x,x5,y6;
   ideal I3 = y2x,x3y,x5,y6;
   ideal I4 = y2x,11x3y,x5;
   ideal I5 = y2x,yx3,x5,7y6;
   list T = I1,I2,I3,I4,I5;
   deleteUnluckyPrimes(T,L,1);
   list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5);
   deleteUnluckyPrimes(T,L,1,P);
}

proc primeTest(def II, bigint p)
{
   if(typeof(II) == "string")
   {
      ideal I = `II`.Ideal;
   }
   else
   {
      ideal I = II;
   }

   I = simplify(I, 2);   // erase zero generators

   int i,j;
   poly f;
   number cnt;
   for(i = 1; i <= size(I); i++)
   {
      f = cleardenom(I[i]);
      if(f == 0) { return(0); }
      cnt = leadcoef(I[i])/leadcoef(f);
      if((bigint(numerator(cnt)) mod p) == 0) { return(0); }
      if((bigint(denominator(cnt)) mod p) == 0) { return(0); }
      for(j = size(f); j > 0; j--)
      {
         if((bigint(leadcoef(f[j])) mod p) == 0) { return(0); }
      }
   }
   return(1);
}

proc primeList(ideal I, int n, list #)
"USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
         n integer
RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
         < L[size(L)] union with L) such that none of these primes divides any
         coefficient occuring in I
NOTE:    The number of cores to use can be defined by ncores, default is 1.
EXAMPLE: example primeList; shows an example
"
{
   intvec L;
   int i,p;
   int ncores = 1;

//-----------------  Initialize optional parameter ncores  ---------------------
   if(size(#) > 0)
   {
      if(size(#) == 1)
      {
         if(typeof(#[1]) == "int")
         {
            ncores = #[1];
            # = list();
         }
      }
      else
      {
         ncores = #[2];
      }
   }

   if(size(#) == 0)
   {
      p = 2147483647;
      while(!primeTest(I,p))
      {
         p = prime(p-1);
         if(p == 2) { ERROR("no more primes"); }
      }
      L[1] = p;
   }
   else
   {
      L = #[1];
      p = prime(L[size(L)]-1);
      while(!primeTest(I,p))
      {
         p = prime(p-1);
         if(p == 2) { ERROR("no more primes"); }
      }
      L[size(L)+1] = p;
   }
   if(p == 2) { ERROR("no more primes"); }
   if(ncores == 1)
   {
      for(i = 2; i <= n; i++)
      {
         p = prime(p-1);
         while(!primeTest(I,p))
         {
            p = prime(p-1);
            if(p == 2) { ERROR("no more primes"); }
         }
         L[size(L)+1] = p;
      }
   }
   else
   {
      int neededSize = size(L)+n-1;;
      list parallelResults;
      list arguments;
      int neededPrimes = neededSize-size(L);
      idealPrimeTest Id;
      Id.Ideal = I;
      export(Id);
      while(neededPrimes > 0)
      {
         arguments = list();
         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
            *ncores; i > 0; i--)
         {
            p = prime(p-1);
            if(p == 2) { ERROR("no more primes"); }
            arguments[i] = list("Id", p);
         }
         parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores);
         for(i = size(arguments); i > 0; i--)
         {
            if(parallelResults[i])
            {
               L[size(L)+1] = arguments[i][2];
            }
         }
         neededPrimes = neededSize-size(L);
      }
      kill Id;
      if(size(L) > neededSize)
      {
         L = L[1..neededSize];
      }
   }
   return(L);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),dp;
   ideal I = 2147483647x+y, z-181;
   intvec L = primeList(I,10);
   size(L);
   L[1];
   L[size(L)];
   L = primeList(I,5,L);
   size(L);
   L[size(L)];
}

////////////////////////////// further examples ////////////////////////////////

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