Changeset a2e51b in git


Ignore:
Timestamp:
May 13, 2014, 4:08:04 PM (10 years ago)
Author:
Andreas Steenpass <steenpass@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
c6b5f62a2d93b0c3cf0e5bd7cb3a8c2eee553e1b
Parents:
a93f95f272cf2777251e0e930411dedb45d02ba0
git-author:
Andreas Steenpass <steenpass@mathematik.uni-kl.de>2014-05-13 16:08:04+02:00
git-committer:
Andreas Steenpass <steenpass@mathematik.uni-kl.de>2014-05-20 08:03:42+02:00
Message:
base modstd.lib on modular.lib

(cherry picked from commit 5ec123fda2041b4cace52475654a231bc4c8adaa)
Signed-off-by: Andreas Steenpass <steenpass@mathematik.uni-kl.de>

Conflicts:
	Singular/LIB/modstd.lib
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/modstd.lib

    ra93f95f ra2e51b  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="version modstd.lib 4.0.0.0 Dec_2013 "; // $Id$
    3 category = "Commutative Algebra";
     2version="version modstd.lib 4.0.0.0 May_2014 "; // $Id$
     3category="Commutative Algebra";
    44info="
    5 LIBRARY:  modstd.lib      Groebner basis of ideals
     5LIBRARY:  modstd.lib      Groebner bases of ideals using modular methods
    66
    77AUTHORS:  A. Hashemi      Amir.Hashemi@lip6.fr
    8 @*        G. Pfister      pfister@mathematik.uni-kl.de
    9 @*        H. Schoenemann  hannes@mathematik.uni-kl.de
    10 @*        A. Steenpass    steenpass@mathematik.uni-kl.de
    11 @*        S. Steidel      steidel@mathematik.uni-kl.de
     8          G. Pfister      pfister@mathematik.uni-kl.de
     9          H. Schoenemann  hannes@mathematik.uni-kl.de
     10          A. Steenpass    steenpass@mathematik.uni-kl.de
     11          S. Steidel      steidel@mathematik.uni-kl.de
    1212
    1313OVERVIEW:
    14 
    15   A library for computing the Groebner basis of an ideal in the polynomial
    16   ring over the rational numbers using modular methods. The procedures are
    17   inspired by the following paper:
    18   Elizabeth A. Arnold: Modular algorithms for computing Groebner bases.
    19   Journal of Symbolic Computation 35, 403-419 (2003).
     14  A library for computing Groebner bases of ideals in the polynomial ring over
     15  the rational numbers using modular methods.
     16
     17REFERENCES:
     18  E. A. Arnold: Modular algorithms for computing Groebner bases.
     19  J. Symb. Comp. 35, 403-419 (2003).
     20
     21  N. Idrees, G. Pfister, S. Steidel: Parallelization of Modular Algorithms.
     22  J. Symb. Comp. 46, 672-684 (2011).
    2023
    2124PROCEDURES:
    22  modStd(I);        standard basis of I using modular methods (chinese remainder)
    23  modS(I,L);        liftings to Q of standard bases of I mod p for p in L
    24  modHenselStd(I);  standard basis of I using modular methods (hensel lifting)
     25  modStd(I);    standard basis of I using modular methods
    2526";
    2627
    2728LIB "poly.lib";
    28 LIB "ring.lib";
    29 LIB "parallel.lib";
    30 
    31 ////////////////////////////////////////////////////////////////////////////////
    32 
    33 static proc mod_init()
    34 {
    35    newstruct("idealPrimeTest", "ideal Ideal");
    36 }
    37 
    38 ////////////////////////////////////////////////////////////////////////////////
    39 
    40 static proc redFork(ideal I, ideal J, int n)
    41 {
    42    attrib(J,"isSB",1);
    43    return(reduce(I,J,1));
    44 }
    45 
    46 ////////////////////////////////////////////////////////////////////////////////
    47 
    48 proc isIncluded(ideal I, ideal J, list #)
    49 "USAGE:  isIncluded(I,J); I,J ideals
    50 RETURN:  1 if J includes I,
    51 @*       0 if there is an element f in I which does not reduce to 0 w.r.t. J.
    52 EXAMPLE: example isIncluded; shows an example
    53 "
    54 {
    55    def R = basering;
    56    setring R;
    57 
    58    attrib(J,"isSB",1);
    59    int i,j,k;
    60 
    61    if(size(#) > 0)
    62    {
    63       int n = #[1];
    64       if(n >= ncols(I)) { n = ncols(I); }
    65       if(n > 1)
    66       {
    67          for(i = 1; i <= n - 1; i++)
    68          {
    69             //link l(i) = "MPtcp:fork";
    70             link l(i) = "ssi:fork";
    71             open(l(i));
    72 
    73             write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1)));
    74          }
    75 
    76          int t = timer;
    77          if(reduce(I[ncols(I)], J, 1) != 0)
    78          {
    79             for(i = 1; i <= n - 1; i++)
    80             {
    81                close(l(i));
     29LIB "modular.lib";
     30
     31proc modStd(ideal I, list #)
     32"USAGE:   modStd(I[, exactness]); I ideal, exactness int
     33RETURN:   a standard basis of I
     34NOTE:     The procedure computes a standard basis of I (over the rational
     35          numbers) by using modular methods.
     36       @* An optional parameter 'exactness' can be provided.
     37          If exactness = 1, the procedure computes a standard basis of I for
     38          sure; if exactness = 0, it computes a standard basis of I
     39          with high probability.
     40SEE ALSO: modular
     41EXAMPLE:  example modStd; shows an example"
     42{
     43    /* read optional parameter */
     44    int exactness = 1;
     45    if (size(#) > 0) {
     46        /* For compatibility, we only test size(#) > 4. This can be changed to
     47         * size(#) > 1 in the future. */
     48        if (size(#) > 4 || typeof(#[1]) != "int") {
     49            ERROR("wrong optional parameter");
     50        }
     51        exactness = #[1];
     52    }
     53
     54    /* save options */
     55    intvec opt = option(get);
     56    option(redSB);
     57
     58    /* choose the right command */
     59    string command = "groebner";
     60    if (npars(basering) > 0) {
     61        command = "Modstd::groebner_norm";
     62    }
     63
     64    /* call modular() */
     65    if (exactness) {
     66        I = modular(command, list(I), primeTest_std,
     67            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
     68    }
     69    else {
     70        I = modular(command, list(I), primeTest_std,
     71            deleteUnluckyPrimes_std, pTest_std);
     72    }
     73
     74    /* return the result */
     75    attrib(I, "isSB", 1);
     76    option(set, opt);
     77    return(I);
     78}
     79example
     80{
     81    "EXAMPLE:";
     82    echo = 2;
     83    ring R1 = 0, (x,y,z,t), dp;
     84    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     85    ideal J = modStd(I);
     86    J;
     87    I = homog(I, t);
     88    J = modStd(I);
     89    J;
     90
     91    ring R2 = 0, (x,y,z), ds;
     92    ideal I = jacob(x5+y6+z7+xyz);
     93    ideal J = modStd(I, 0);
     94    J;
     95
     96    ring R3 = 0, x(1..4), lp;
     97    ideal I = cyclic(4);
     98    ideal J1 = modStd(I, 1);   // default
     99    ideal J2 = modStd(I, 0);
     100    size(reduce(J1, J2));
     101    size(reduce(J2, J1));
     102}
     103
     104/* compute a normalized GB via groebner() */
     105static proc groebner_norm(ideal I)
     106{
     107    I = simplify(groebner(I), 1);
     108    attrib(I, "isSB", 1);
     109    return(I);
     110}
     111
     112/* test if the prime p is suitable for the input, i.e. it does not divide
     113 * the numerator or denominator of any of the coefficients */
     114static proc primeTest_std(int p, alias list args)
     115{
     116    /* erase zero generators */
     117    ideal I = simplify(args[1], 2);
     118
     119    /* clear denominators and count the terms */
     120    ideal J;
     121    ideal K;
     122    int n = ncols(I);
     123    intvec sizes;
     124    number cnt;
     125    int i;
     126    for(i = n; i > 0; i--) {
     127        J[i] = cleardenom(I[i]);
     128        cnt = leadcoef(J[i])/leadcoef(I[i]);
     129        K[i] = numerator(cnt)*var(1)+denominator(cnt);
     130    }
     131    sizes = size(J[1..n]);
     132
     133    /* change to characteristic p */
     134    def br = basering;
     135    list lbr = ringlist(br);
     136    if (typeof(lbr[1]) == "int") {
     137        lbr[1] = p;
     138    }
     139    else {
     140        lbr[1][1] = p;
     141    }
     142    def rp = ring(lbr);
     143    setring(rp);
     144    ideal Jp = fetch(br, J);
     145    ideal Kp = fetch(br, K);
     146
     147    /* test if any coefficient is missing */
     148    if (intvec(size(Kp[1..n])) != 2:n) {
     149        setring(br);
     150        return(0);
     151    }
     152    if (intvec(size(Jp[1..n])) != sizes) {
     153        setring(br);
     154        return(0);
     155    }
     156    setring(br);
     157    return(1);
     158}
     159
     160/* find entries in modresults which come from unlucky primes.
     161 * For this, sort the entries into categories depending on their leading
     162 * ideal and return the indices in all but the biggest category. */
     163static proc deleteUnluckyPrimes_std(alias list modresults)
     164{
     165    int size_modresults = size(modresults);
     166
     167    /* sort results into categories.
     168     * each category is represented by three entries:
     169     * - the corresponding leading ideal
     170     * - the number of elements
     171     * - the indices of the elements
     172     */
     173    list cat;
     174    int size_cat;
     175    ideal L;
     176    int i;
     177    int j;
     178    for (i = 1; i <= size_modresults; i++) {
     179        L = lead(modresults[i]);
     180        attrib(L, "isSB", 1);
     181        for (j = 1; j <= size_cat; j++) {
     182            if (size(L) == size(cat[j][1])
     183                && size(reduce(L, cat[j][1])) == 0
     184                && size(reduce(cat[j][1], L)) == 0) {
     185                cat[j][2] = cat[j][2]+1;
     186                cat[j][3][cat[j][2]] = i;
     187                break;
    82188            }
     189        }
     190        if (j > size_cat) {
     191            size_cat++;
     192            cat[size_cat] = list();
     193            cat[size_cat][1] = L;
     194            cat[size_cat][2] = 1;
     195            cat[size_cat][3] = list(i);
     196        }
     197    }
     198
     199    /* find the biggest categories */
     200    int cat_max = 1;
     201    int max = cat[1][2];
     202    for (i = 2; i <= size_cat; i++) {
     203        if (cat[i][2] > max) {
     204            cat_max = i;
     205            max = cat[i][2];
     206        }
     207    }
     208
     209    /* return all other indices */
     210    list unluckyIndices;
     211    for (i = 1; i <= size_cat; i++) {
     212        if (i != cat_max) {
     213            unluckyIndices = unluckyIndices + cat[i][3];
     214        }
     215    }
     216    return(unluckyIndices);
     217}
     218
     219/* test if 'command' applied to 'args' in characteristic p is the same as
     220   'result' mapped to characteristic p */
     221static proc pTest_std(string command, list args, ideal result, int p)
     222{
     223    /* change to characteristic p */
     224    def br = basering;
     225    list lbr = ringlist(br);
     226    if (typeof(lbr[1]) == "int") {
     227        lbr[1] = p;
     228    }
     229    else {
     230        lbr[1][1] = p;
     231    }
     232    def rp = ring(lbr);
     233    setring(rp);
     234    ideal Ip = fetch(br, args)[1];
     235    ideal Gp = fetch(br, result);
     236    attrib(Gp, "isSB", 1);
     237
     238    /* test if Ip is in Gp */
     239    int i;
     240    for (i = ncols(Ip); i > 0; i--) {
     241        if (reduce(Ip[i], Gp, 1) != 0) {
     242            setring(br);
    83243            return(0);
    84          }
    85          t = timer - t;
    86          if(t > 60) { t = 60; }
    87          int i_sleep = system("sh", "sleep "+string(t));
    88 
    89          j = ncols(I) - n;
    90 
    91          while(j >= 0)
    92          {
    93             for(i = 1; i <= n - 1; i++)
    94             {
    95                if(status(l(i), "read", "ready"))
    96                {
    97                   if(read(l(i)) != 0)
    98                   {
    99                      for(i = 1; i <= n - 1; i++)
    100                      {
    101                         close(l(i));
    102                      }
    103                      return(0);
    104                   }
    105                   else
    106                   {
    107                      if(j >= 1)
    108                      {
    109                         write(l(i), quote(redFork(eval(I[j]), eval(J), 1)));
    110                         j--;
    111                      }
    112                      else
    113                      {
    114                         k++;
    115                         close(l(i));
    116                      }
    117                   }
    118                }
    119             }
    120             if(k == n - 1)
    121             {
    122                j--;
    123             }
    124             i_sleep = system("sh", "sleep "+string(t));
    125          }
    126          return(1);
    127       }
    128    }
    129 
    130    for(i = ncols(I); i >= 1; i--)
    131    {
    132       if(reduce(I[i],J,1) != 0){ return(0); }
    133    }
    134    return(1);
    135 }
    136 example
    137 { "EXAMPLE:"; echo = 2;
    138    ring r=0,(x,y,z),dp;
    139    ideal I = x+1,x+y+1;
    140    ideal J = x+1,y;
    141    isIncluded(I,J);
    142    isIncluded(J,I);
    143    isIncluded(I,J,4);
    144 
    145    ring R = 0, x(1..5), dp;
    146    ideal I1 = cyclic(4);
    147    ideal I2 = I1,x(5)^2;
    148    isIncluded(I1,I2,4);
    149 }
    150 
    151 ////////////////////////////////////////////////////////////////////////////////
    152 
    153 proc pTestSB(ideal I, ideal J, list L, int variant, list #)
    154 "USAGE:  pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
    155 RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
    156          J mod p is (resp. is not) a standard basis of I mod p
    157 EXAMPLE: example pTestSB; shows an example
    158 "
    159 {
    160    int i,j,k,p;
    161    def R = basering;
    162    list r = ringlist(R);
    163 
    164    while(!j)
    165    {
    166       j = 1;
    167       p = prime(random(1000000000,2134567879));
    168       for(i = 1; i <= size(L); i++)
    169       {
    170          if(p == L[i]) { j = 0; break; }
    171       }
    172       if(j)
    173       {
    174          for(i = 1; i <= ncols(I); i++)
    175          {
    176             for(k = 2; k <= size(I[i]); k++)
    177             {
    178                if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; }
    179             }
    180             if(!j){ break; }
    181          }
    182       }
    183       if(j)
    184       {
    185          if(!primeTest(I,p)) { j = 0; }
    186       }
    187    }
    188    r[1] = p;
    189    def @R = ring(r);
    190    setring @R;
    191    ideal I = imap(R,I);
    192    ideal J = imap(R,J);
    193    attrib(J,"isSB",1);
    194 
    195    int t = timer;
    196    j = 1;
    197    if(isIncluded(I,J) == 0) { j = 0; }
    198 
    199    if(printlevel >= 11)
    200    {
    201       "isIncluded(I,J) takes "+string(timer - t)+" seconds";
    202       "j = "+string(j);
    203    }
    204 
    205    t = timer;
    206    if(j)
    207    {
    208       if(size(#) > 0)
    209       {
    210          ideal K = modpStd(I,p,variant,#[1])[1];
    211       }
    212       else
    213       {
    214          ideal K = groebner(I);
    215       }
    216       t = timer;
    217       if(isIncluded(J,K) == 0) { j = 0; }
    218 
    219       if(printlevel >= 11)
    220       {
    221          "isIncluded(J,K) takes "+string(timer - t)+" seconds";
    222          "j = "+string(j);
    223       }
    224    }
    225    setring R;
    226    return(j);
    227 }
    228 example
    229 { "EXAMPLE:"; echo = 2;
    230    intvec L = 2,3,5;
    231    ring r = 0,(x,y,z),dp;
    232    ideal I = x+1,x+y+1;
    233    ideal J = x+1,y;
    234    pTestSB(I,I,L,2);
    235    pTestSB(I,J,L,2);
    236 }
    237 
    238 ////////////////////////////////////////////////////////////////////////////////
    239 
    240 proc deleteUnluckyPrimes(list T, list L, int ho, list #)
    241 "USAGE:  deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer
    242 RETURN:  lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #),
    243          lT ideal
    244 NOTE:    - if ho = 1, the polynomials in T are homogeneous, else ho = 0,
    245 @*       - lT is prevalent, i.e. the most appearing leading ideal in T
    246 EXAMPLE: example deleteUnluckyPrimes; shows an example
    247 "
    248 {
    249    ho = ((ho)||(ord_test(basering) == -1));
    250    int j,k,c;
    251    intvec hl,hc;
    252    ideal cT,lT,cK;
    253    lT = lead(T[size(T)]);
    254    attrib(lT,"isSB",1);
    255    if(!ho)
    256    {
    257       for(j = 1; j < size(T); j++)
    258       {
    259          cT = lead(T[j]);
    260          attrib(cT,"isSB",1);
    261          if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
    262          {
    263             cK = cT;
    264             c++;
    265          }
    266       }
    267       if(c > size(T) div 2){ lT = cK; }
    268    }
    269    else
    270    {
    271       hl = hilb(lT,1);
    272       for(j = 1; j < size(T); j++)
    273       {
    274          cT = lead(T[j]);
    275          attrib(cT,"isSB",1);
    276          hc = hilb(cT,1);
    277          if(hl == hc)
    278          {
    279             for(k = 1; k <= size(lT); k++)
    280             {
    281                if(lT[k] < cT[k]) { lT = cT; c++; break; }
    282                if(lT[k] > cT[k]) { c++; break; }
    283             }
    284          }
    285          else
    286          {
    287             if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; }
    288          }
    289       }
    290    }
    291 
    292    int addList;
    293    if(size(#) > 0) { list M = #; addList = 1; }
    294    j = 1;
    295    attrib(lT,"isSB",1);
    296    while((j <= size(T))&&(c > 0))
    297    {
    298       cT = lead(T[j]);
    299       attrib(cT,"isSB",1);
    300       if((size(reduce(cT,lT)) != 0)||(size(reduce(lT,cT)) != 0))
    301       {
    302          T = delete(T,j);
    303          if(j == 1)
    304          {
    305             L = L[2..size(L)];
    306             if(addList == 1) { M = M[2..size(M)]; }
    307          }
    308          else
    309          {
    310             if(j == size(L))
    311             {
    312                L = L[1..size(L)-1];
    313                if(addList == 1) { M = M[1..size(M)-1]; }
    314             }
    315             else
    316             {
    317                L = L[1..j-1],L[j+1..size(L)];
    318                if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; }
    319             }
    320          }
    321          j--;
    322       }
    323       j++;
    324    }
    325 
    326    for(j = 1; j <= size(L); j++)
    327    {
    328       L[j] = bigint(L[j]);
    329    }
    330 
    331    if(addList == 0) { return(list(T,L,lT)); }
    332    if(addList == 1) { return(list(T,L,M,lT)); }
    333 }
    334 example
    335 { "EXAMPLE:"; echo = 2;
    336    list L = 2,3,5,7,11;
    337    ring r = 0,(y,x),Dp;
    338    ideal I1 = 2y2x,y6;
    339    ideal I2 = yx2,y3x,x5,y6;
    340    ideal I3 = y2x,x3y,x5,y6;
    341    ideal I4 = y2x,11x3y,x5;
    342    ideal I5 = y2x,yx3,x5,7y6;
    343    list T = I1,I2,I3,I4,I5;
    344    deleteUnluckyPrimes(T,L,1);
    345    list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5);
    346    deleteUnluckyPrimes(T,L,1,P);
    347 }
    348 
    349 ////////////////////////////////////////////////////////////////////////////////
    350 
    351 proc primeTest(def II, bigint p)
    352 {
    353    if(typeof(II) == "string")
    354    {
    355       ideal I = `II`.Ideal;
    356    }
    357    else
    358    {
    359       ideal I = II;
    360    }
    361 
    362    I = simplify(I, 2);   // erase zero generators
    363 
    364    int i,j;
    365    poly f;
    366    number cnt;
    367    for(i = 1; i <= size(I); i++)
    368    {
    369       f = cleardenom(I[i]);
    370       if(f == 0) { return(0); }
    371       cnt = leadcoef(I[i])/leadcoef(f);
    372       if((numerator(cnt) mod p) == 0) { return(0); }
    373       if((denominator(cnt) mod p) == 0) { return(0); }
    374       for(j = size(f); j > 0; j--)
    375       {
    376          if((leadcoef(f[j]) mod p) == 0) { return(0); }
    377       }
    378    }
    379    return(1);
    380 }
    381 
    382 ////////////////////////////////////////////////////////////////////////////////
    383 
    384 proc primeList(ideal I, int n, list #)
    385 "USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
    386          n integer
    387 RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
    388          < L[size(L)] union with L) such that none of these primes divides any
    389          coefficient occuring in I
    390 NOTE:    The number of cores to use can be defined by ncores, default is 1.
    391 EXAMPLE: example primeList; shows an example
    392 "
    393 {
    394    intvec L;
    395    int i,p;
    396    int ncores = 1;
    397 
    398 //-----------------  Initialize optional parameter ncores  ---------------------
    399    if(size(#) > 0)
    400    {
    401       if(size(#) == 1)
    402       {
    403          if(typeof(#[1]) == "int")
    404          {
    405             ncores = #[1];
    406             # = list();
    407          }
    408       }
    409       else
    410       {
    411          ncores = #[2];
    412       }
    413    }
    414 
    415    if(size(#) == 0)
    416    {
    417       p = 2147483647;
    418       while(!primeTest(I,p))
    419       {
    420          p = prime(p-1);
    421          if(p == 2) { ERROR("no more primes"); }
    422       }
    423       L[1] = p;
    424    }
    425    else
    426    {
    427       L = #[1];
    428       p = prime(L[size(L)]-1);
    429       while(!primeTest(I,p))
    430       {
    431          p = prime(p-1);
    432          if(p == 2) { ERROR("no more primes"); }
    433       }
    434       L[size(L)+1] = p;
    435    }
    436    if(p == 2) { ERROR("no more primes"); }
    437    if(ncores == 1)
    438    {
    439       for(i = 2; i <= n; i++)
    440       {
    441          p = prime(p-1);
    442          while(!primeTest(I,p))
    443          {
    444             p = prime(p-1);
    445             if(p == 2) { ERROR("no more primes"); }
    446          }
    447          L[size(L)+1] = p;
    448       }
    449    }
    450    else
    451    {
    452       int neededSize = size(L)+n-1;;
    453       list parallelResults;
    454       list arguments;
    455       int neededPrimes = neededSize-size(L);
    456       idealPrimeTest Id;
    457       Id.Ideal = I;
    458       export(Id);
    459       while(neededPrimes > 0)
    460       {
    461          arguments = list();
    462          for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
    463             *ncores; i > 0; i--)
    464          {
    465             p = prime(p-1);
    466             if(p == 2) { ERROR("no more primes"); }
    467             arguments[i] = list("Id", p);
    468          }
    469          parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores);
    470          for(i = size(arguments); i > 0; i--)
    471          {
    472             if(parallelResults[i])
    473             {
    474                L[size(L)+1] = arguments[i][2];
    475             }
    476          }
    477          neededPrimes = neededSize-size(L);
    478       }
    479       kill Id;
    480       if(size(L) > neededSize)
    481       {
    482          L = L[1..neededSize];
    483       }
    484    }
    485    return(L);
    486 }
    487 example
    488 { "EXAMPLE:"; echo = 2;
    489    ring r = 0,(x,y,z),dp;
    490    ideal I = 2147483647x+y, z-181;
    491    intvec L = primeList(I,10);
    492    size(L);
    493    L[1];
    494    L[size(L)];
    495    L = primeList(I,5,L);
    496    size(L);
    497    L[size(L)];
    498 }
    499 
    500 ////////////////////////////////////////////////////////////////////////////////
    501 
    502 static proc liftstd1(ideal I)
    503 {
    504    def R = basering;
    505    list rl = ringlist(R);
    506    list ordl = rl[3];
    507 
    508    int i;
    509    for(i = 1; i <= size(ordl); i++)
    510    {
    511       if((ordl[i][1] == "C") || (ordl[i][1] == "c"))
    512       {
    513          ordl = delete(ordl, i);
    514          break;
    515       }
    516    }
    517 
    518    ordl = insert(ordl, list("c", 0));
    519    rl[3] = ordl;
    520    def newR = ring(rl);
    521    setring newR;
    522    ideal I = imap(R,I);
    523 
    524    intvec opt = option(get);
    525    option(none);
    526    option(prompt);
    527 
    528    module M;
    529    for(i = 1; i <= size(I); i++)
    530    {
    531       M = M + module(I[i]*gen(1) + gen(i+1));
    532       M = M + module(gen(i+1));
    533    }
    534 
    535    module sM = std(M);
    536 
    537    ideal sI;
    538    if(attrib(R,"global"))
    539    {
    540       for(i = size(I)+1; i <= size(sM); i++)
    541       {
    542          sI[size(sI)+1] = sM[i][1];
    543       }
    544       matrix T = submat(sM,2..nrows(sM),size(I)+1..ncols(sM));
    545    }
    546    else
    547    {
    548       //"==========================================================";
    549       //"WARNING: Algorithm is not applicable if ordering is mixed.";
    550       //"==========================================================";
    551       for(i = 1; i <= size(sM)-size(I); i++)
    552       {
    553          sI[size(sI)+1] = sM[i][1];
    554       }
    555       matrix T = submat(sM,2..nrows(sM),1..ncols(sM)-size(I));
    556    }
    557 
    558    setring R;
    559    option(set, opt);
    560    return(imap(newR,sI),imap(newR,T));
    561 }
    562 example
    563 { "EXAMPLE:"; echo = 2;
    564    ring R = 0,(x,y,z),dp;
    565    poly f = x3+y7+z2+xyz;
    566    ideal i = jacob(f);
    567    matrix T;
    568    ideal sm = liftstd(i,T);
    569    sm;
    570    print(T);
    571    matrix(sm) - matrix(i)*T;
    572 
    573 
    574    ring S = 32003, x(1..5), lp;
    575    ideal I = cyclic(5);
    576    ideal sI;
    577    matrix T;
    578    sI,T = liftstd1(I);
    579    matrix(sI) - matrix(I)*T;
    580 }
    581 
    582 ////////////////////////////////////////////////////////////////////////////////
    583 
    584 proc modpStd(ideal I, int p, int variant, list #)
    585 "USAGE:  modpStd(I,p,variant,#); I ideal, p integer, variant integer
    586 ASSUME:  If size(#) > 0, then #[1] is an intvec describing the Hilbert series.
    587 RETURN:  ideal - a standard basis of I mod p, integer - p
    588 NOTE:    The procedure computes a standard basis of the ideal I modulo p and
    589          fetches the result to the basering. If size(#) > 0 the Hilbert driven
    590          standard basis computation std(.,#[1]) is used instead of groebner.
    591          The standard basis computation modulo p does also vary depending on the
    592          integer variant, namely
    593 @*       - variant = 1: std(.,#[1]) resp. groebner,
    594 @*       - variant = 2: groebner,
    595 @*       - variant = 3: homog. - std(.,#[1]) resp. groebner - dehomog.,
    596 @*       - variant = 4: fglm.
    597 EXAMPLE: example modpStd; shows an example
    598 "
    599 {
    600    def R0 = basering;
    601    list rl = ringlist(R0);
    602    rl[1] = p;
    603    def @r = ring(rl);
    604    setring @r;
    605    ideal i = fetch(R0,I);
    606 
    607    option(redSB);
    608 
    609    if(variant == 1)
    610    {
    611       if(size(#) > 0)
    612       {
    613          i = std(i, #[1]);
    614       }
    615       else
    616       {
    617          i = groebner(i);
    618       }
    619    }
    620 
    621    if(variant == 2)
    622    {
    623       i = groebner(i);
    624    }
    625 
    626    if(variant == 3)
    627    {
    628       list rl = ringlist(@r);
    629       int nvar@r = nvars(@r);
    630 
    631       int k;
    632       intvec w;
    633       for(k = 1; k <= nvar@r; k++)
    634       {
    635          w[k] = deg(var(k));
    636       }
    637       w[nvar@r + 1] = 1;
    638 
    639       rl[2][nvar@r + 1] = "homvar";
    640       rl[3][2][2] = w;
    641 
    642       def HomR = ring(rl);
    643       setring HomR;
    644       ideal i = imap(@r, i);
    645       i = homog(i, homvar);
    646 
    647       if(size(#) > 0)
    648       {
    649          if(w == 1)
    650          {
    651             i = std(i, #[1]);
    652          }
    653          else
    654          {
    655             i = std(i, #[1], w);
    656          }
    657       }
    658       else
    659       {
    660          i = groebner(i);
    661       }
    662 
    663       i = subst(i, homvar, 1);
    664       i = simplify(i, 34);
    665 
    666       setring @r;
    667       i = imap(HomR, i);
    668       i = interred(i);
    669       kill HomR;
    670    }
    671 
    672    if(variant == 4)
    673    {
    674       def R1 = changeord(list(list("dp",1:nvars(basering))));
    675       setring R1;
    676       ideal i = fetch(@r,i);
    677       i = std(i);
    678       setring @r;
    679       i = fglm(R1,i);
    680    }
    681 
    682    setring R0;
    683    return(list(fetch(@r,i),p));
    684 }
    685 example
    686 { "EXAMPLE:"; echo = 2;
    687    ring r = 0, x(1..4), dp;
    688    ideal I = cyclic(4);
    689    int p = 181;
    690    list P = modpStd(I,p,2);
    691    P;
    692 
    693    ring r2 = 0, x(1..5), lp;
    694    ideal I = cyclic(5);
    695    int q = 32003;
    696    list Q = modpStd(I,q,4);
    697    Q;
    698 }
    699 
    700 static proc algeModStd(ideal i,list #)
    701 {
    702 //reduces modStd over algebraic extensions to the one over a polynomial ring
    703    list L=#;
    704    def R=basering;
    705    int n=nvars(R);
    706    list rl=ringlist(R);
    707    poly p=minpoly;
    708    rl[2][n+1]=rl[1][2][1];
    709    rl[1]=rl[1][1];
    710    rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
    711    rl[3][size(rl[3])-1]=list("dp",1);
    712    def S=ring(rl);
    713    setring S;
    714    poly p=imap(R,p);
    715    ideal i=imap(R,i);
    716    i=i,p;
    717    ideal j=modStd(i,L);
    718    if(j[1]==p)
    719    {
    720       j[1]=0;
    721    }
    722    j=simplify(j,2);
    723    setring R;
    724    ideal j=imap(S,j);
    725    return(j);
    726 }
    727 ////////////////////////////// main procedures /////////////////////////////////
    728 
    729 proc modStd(ideal I, list #)
    730 "USAGE:  modStd(I); I ideal
    731 ASSUME:  If size(#) > 0, then # contains either 1, 2 or 4 integers such that
    732 @*       - #[1] is the number of available processors for the computation,
    733 @*       - #[2] is an optional parameter for the exactness of the computation,
    734                 if #[2] = 1, the procedure computes a standard basis for sure,
    735 @*       - #[3] is the number of primes until the first lifting,
    736 @*       - #[4] is the constant number of primes between two liftings until
    737            the computation stops.
    738 RETURN:  a standard basis of I if no warning appears;
    739 NOTE:    The procedure computes a standard basis of I (over the rational
    740          numbers) by using modular methods.
    741          By default the procedure computes a standard basis of I for sure, but
    742          if the optional parameter #[2] = 0, it computes a standard basis of I
    743          with high probability.
    744          The procedure distinguishes between different variants for the standard
    745          basis computation in positive characteristic depending on the ordering
    746          of the basering, the parameter #[2] and if the ideal I is homogeneous.
    747 @*       - variant = 1, if I is homogeneous,
    748 @*       - variant = 2, if I is not homogeneous, 1-block-ordering,
    749 @*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
    750                         > 1 block),
    751 @*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
    752 EXAMPLE: example modStd; shows an example
    753 "
    754 {
    755    int TT = timer;
    756    int RT = rtimer;
    757 
    758    def R0 = basering;
    759    list rl = ringlist(R0);
    760 
    761    int algebraic;
    762    if(size(#)>0)
    763    {
    764       if(#[1]<=0)
    765       {
    766          algebraic=1;
    767          #[1]=-#[1];
    768          if(#[1]==0){list LK;#=LK;}
    769       }
    770    }
    771 
    772    if((npars(R0) > 0) || (rl[1][1] > 0))
    773    {
    774       if(npars(R0)==1)
    775       {
    776          if(minpoly!=0)
    777          {
    778             list LM=#;
    779             if(size(LM)==0){LM[1]=0;}
    780             LM[1]=-LM[1];
    781             return(algeModStd(I,LM));
    782          }
    783       }
    784 
    785       ERROR("Characteristic of basering should be zero, basering should
    786              have no parameters.");
    787    }
    788 
    789    int index = 1;
    790    int i,k,c;
    791    int j = 1;
    792    int pTest, sizeTest;
    793    int en = 2134567879;
    794    int an = 1000000000;
    795    bigint N;
    796 
    797 //--------------------  Initialize optional parameters  ------------------------
    798    if(size(#) > 0)
    799    {
    800       if(size(#) == 1)
    801       {
    802          int n1 = #[1];
    803          int exactness = 1;
    804          if(n1 >= 10)
    805          {
    806             int n2 = n1 + 1;
    807             int n3 = n1;
    808          }
    809          else
    810          {
    811             int n2 = 10;
    812             int n3 = 10;
    813          }
    814       }
    815       if(size(#) == 2)
    816       {
    817          int n1 = #[1];
    818          int exactness = #[2];
    819          if(n1 >= 10)
    820          {
    821             int n2 = n1 + 1;
    822             int n3 = n1;
    823          }
    824          else
    825          {
    826             int n2 = 10;
    827             int n3 = 10;
    828          }
    829       }
    830       if(size(#) == 4)
    831       {
    832          int n1 = #[1];
    833          int exactness = #[2];
    834          if(n1 >= #[3])
    835          {
    836             int n2 = n1 + 1;
    837          }
    838          else
    839          {
    840             int n2 = #[3];
    841          }
    842          if(n1 >= #[4])
    843          {
    844             int n3 = n1;
    845          }
    846          else
    847          {
    848             int n3 = #[4];
    849          }
    850       }
    851    }
    852    else
    853    {
    854       int n1 = 1;
    855       int exactness = 1;
    856       int n2 = 10;
    857       int n3 = 10;
    858    }
    859 
    860    if(printlevel >= 10)
    861    {
    862       "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
    863        +", exactness = "+string(exactness);
    864    }
    865 
    866 //-------------------------  Save current options  -----------------------------
    867    intvec opt = option(get);
    868 
    869    option(redSB);
    870 
    871 //--------------------  Initialize the list of primes  -------------------------
    872    int tt = timer;
    873    int rt = rtimer;
    874    intvec L = primeList(I,n2,n1);
    875    if(printlevel >= 10)
    876    {
    877       "CPU-time for primeList: "+string(timer-tt)+" seconds.";
    878       "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
    879    }
    880    L[5] = prime(random(an,en));
    881 
    882 //---------------------  Decide which variant to take  -------------------------
    883    int variant;
    884    int h = homog(I);
    885 
    886    tt = timer;
    887    rt = rtimer;
    888 
    889    if(!hasMixedOrdering())
    890    {
    891       if(h)
    892       {
    893          variant = 1;
    894          if(printlevel >= 10) { "variant = 1"; }
    895 
    896          rl[1] = L[5];
    897          def @r = ring(rl);
    898          setring @r;
    899          def @s = changeord(list(list("dp",1:nvars(basering))));
    900          setring @s;
    901          ideal I = std(fetch(R0,I));
    902          intvec hi = hilb(I,1);
    903          setring R0;
    904          kill @r,@s;
    905       }
    906       else
    907       {
    908          string ordstr_R0 = ordstr(R0);
    909          int neg = 1 - attrib(R0,"global");
    910 
    911          if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
    912          {
    913             variant = 2;
    914             if(printlevel >= 10) { "variant = 2"; }
    915          }
    916          else
    917          {
    918             string order;
    919             if(system("nblocks") <= 2)
    920             {
    921                if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
    922                                        + find(ordstr_R0, "rp") <= 0)
    923                {
    924                   order = "simple";
    925                }
    926             }
    927 
    928             if((order == "simple") || (size(rl) > 4))
    929             {
    930                variant = 2;
    931                if(printlevel >= 10) { "variant = 2"; }
    932             }
    933             else
    934             {
    935                rl[1] = L[5];
    936                def @r = ring(rl);
    937                setring @r;
    938 
    939                def @s = changeord(list(list("dp",1:nvars(basering))));
    940                setring @s;
    941                ideal I = std(fetch(R0,I));
    942                if(dim(I) == 0)
    943                {
    944                   variant = 4;
    945                   if(printlevel >= 10) { "variant = 4"; }
    946                }
    947                else
    948                {
    949                   variant = 3;
    950                   if(printlevel >= 10) { "variant = 3"; }
    951 
    952                   int nvar@r = nvars(@r);
    953                   intvec w;
    954                   for(i = 1; i <= nvar@r; i++)
    955                   {
    956                      w[i] = deg(var(i));
    957                   }
    958                   w[nvar@r + 1] = 1;
    959 
    960                   list hiRi = hilbRing(fetch(R0,I),w);
    961                   intvec W = hiRi[2];
    962                   @s = hiRi[1];
    963                   setring @s;
    964 
    965                   Id(1) = std(Id(1));
    966                   intvec hi = hilb(Id(1), 1, W);
    967                }
    968 
    969                setring R0;
    970                kill @r,@s;
    971             }
    972          }
    973       }
    974    }
    975    else
    976    {
    977       if(exactness == 1) { return(groebner(I)); }
    978       if(h)
    979       {
    980          variant = 1;
    981          if(printlevel >= 10) { "variant = 1"; }
    982          rl[1] = L[5];
    983          def @r = ring(rl);
    984          setring @r;
    985          def @s = changeord(list(list("dp",1:nvars(basering))));
    986          setring @s;
    987          ideal I = std(fetch(R0,I));
    988          intvec hi = hilb(I,1);
    989          setring R0;
    990          kill @r,@s;
    991       }
    992       else
    993       {
    994          string ordstr_R0 = ordstr(R0);
    995          int neg = 1 - attrib(R0,"global");
    996 
    997          if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
    998          {
    999             variant = 2;
    1000             if(printlevel >= 10) { "variant = 2"; }
    1001          }
    1002          else
    1003          {
    1004             string order;
    1005             if(system("nblocks") <= 2)
    1006             {
    1007                if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
    1008                                        + find(ordstr_R0, "rp") <= 0)
    1009                {
    1010                   order = "simple";
    1011                }
    1012             }
    1013 
    1014             if((order == "simple") || (size(rl) > 4))
    1015             {
    1016                variant = 2;
    1017                if(printlevel >= 10) { "variant = 2"; }
    1018             }
    1019             else
    1020             {
    1021                variant = 3;
    1022                if(printlevel >= 10) { "variant = 3"; }
    1023 
    1024                rl[1] = L[5];
    1025                def @r = ring(rl);
    1026                setring @r;
    1027                int nvar@r = nvars(@r);
    1028                intvec w;
    1029                for(i = 1; i <= nvar@r; i++)
    1030                {
    1031                   w[i] = deg(var(i));
    1032                }
    1033                w[nvar@r + 1] = 1;
    1034 
    1035                list hiRi = hilbRing(fetch(R0,I),w);
    1036                intvec W = hiRi[2];
    1037                def @s = hiRi[1];
    1038                setring @s;
    1039 
    1040                Id(1) = std(Id(1));
    1041                intvec hi = hilb(Id(1), 1, W);
    1042 
    1043                setring R0;
    1044                kill @r,@s;
    1045             }
    1046          }
    1047       }
    1048    }
    1049    if(algebraic){variant=2;}
    1050 
    1051    list P,T1,T2,T3,LL;
    1052 
    1053    ideal J,K,H;
    1054 
    1055 //-----  If there is more than one processor available, we parallelize the  ----
    1056 //-----  main standard basis computations in positive characteristic        ----
    1057 
    1058    if(n1 > 1)
    1059    {
    1060       ideal I_for_fork = I;
    1061       export(I_for_fork);           // I available for each link
    1062 
    1063 //-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
    1064 //-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
    1065 
    1066       for(i = 1; i <= n1; i++)
    1067       {
    1068          //link l(i) = "MPtcp:fork";
    1069          link l(i) = "ssi:fork";
    1070          open(l(i));
    1071          if((variant == 1) || (variant == 3))
    1072          {
    1073             write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
    1074                                                   eval(variant), eval(hi))));
    1075          }
    1076          if((variant == 2) || (variant == 4))
    1077          {
    1078             write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
    1079                                                   eval(variant))));
    1080          }
    1081       }
    1082 
    1083       int t = timer;
    1084       if((variant == 1) || (variant == 3))
    1085       {
    1086          P = modpStd(I_for_fork, L[1], variant, hi);
    1087       }
    1088       if((variant == 2) || (variant == 4))
    1089       {
    1090          P = modpStd(I_for_fork, L[1], variant);
    1091       }
    1092       t = timer - t;
    1093       if(t > 60) { t = 60; }
    1094       int i_sleep = system("sh", "sleep "+string(t));
    1095       T1[1] = P[1];
    1096       T2[1] = bigint(P[2]);
    1097       index++;
    1098 
    1099       j = j + n1 + 1;
    1100    }
    1101 
    1102 //--------------  Main standard basis computations in positive  ----------------
    1103 //----------------------  characteristic start here  ---------------------------
    1104 
    1105    list arguments_farey, results_farey;
    1106 
    1107    while(1)
    1108    {
    1109       tt = timer; rt = rtimer;
    1110 
    1111       if(printlevel >= 10) { "size(L) = "+string(size(L)); }
    1112 
    1113       if(n1 > 1)
    1114       {
    1115          while(j <= size(L) + 1)
    1116          {
    1117             for(i = 1; i <= n1; i++)
    1118             {
    1119                //--- ask if link l(i) is ready otherwise sleep for t seconds ---
    1120                if(status(l(i), "read", "ready"))
    1121                {
    1122                   //--- read the result from l(i) ---
    1123                   P = read(l(i));
    1124                   T1[index] = P[1];
    1125                   T2[index] = bigint(P[2]);
    1126                   index++;
    1127 
    1128                   if(j <= size(L))
    1129                   {
    1130                      if((variant == 1) || (variant == 3))
    1131                      {
    1132                         write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
    1133                                                   eval(variant), eval(hi))));
    1134                         j++;
    1135                      }
    1136                      if((variant == 2) || (variant == 4))
    1137                      {
    1138                         write(l(i), quote(modpStd(I_for_fork,
    1139                                                   eval(L[j]), eval(variant))));
    1140                         j++;
    1141                      }
    1142                   }
    1143                   else
    1144                   {
    1145                      k++;
    1146                      close(l(i));
    1147                   }
    1148                }
    1149             }
    1150             //--- k describes the number of closed links ---
    1151             if(k == n1)
    1152             {
    1153                j++;
    1154             }
    1155             i_sleep = system("sh", "sleep "+string(t));
    1156          }
    1157       }
    1158       else
    1159       {
    1160          while(j <= size(L))
    1161          {
    1162             if((variant == 1) || (variant == 3))
    1163             {
    1164                P = modpStd(I, L[j], variant, hi);
    1165             }
    1166             if((variant == 2) || (variant == 4))
    1167             {
    1168                P = modpStd(I, L[j], variant);
    1169             }
    1170 
    1171             T1[index] = P[1];
    1172             T2[index] = bigint(P[2]);
    1173             index++;
    1174             j++;
    1175          }
    1176       }
    1177 
    1178       if(printlevel >= 10)
    1179       {
    1180          "CPU-time for computing list is "+string(timer - tt)+" seconds.";
    1181          "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
    1182       }
    1183 
    1184 //------------------------  Delete unlucky primes  -----------------------------
    1185 //-------------  unlucky if and only if the leading ideal is wrong  ------------
    1186 
    1187       LL = deleteUnluckyPrimes(T1,T2,h);
    1188       T1 = LL[1];
    1189       T2 = LL[2];
    1190 
    1191 //-------------------  Now all leading ideals are the same  --------------------
    1192 //-------------------  Lift results to basering via farey  ---------------------
    1193 
    1194       tt = timer; rt = rtimer;
    1195       N = T2[1];
    1196       for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
    1197       H = chinrem(T1,T2);
    1198       if(n1 == 1)
    1199       {
    1200          J = farey(H,N);
    1201       }
    1202       else
    1203       {
    1204          for(i = ncols(H); i > 0; i--)
    1205          {
    1206             arguments_farey[i] = list(ideal(H[i]), N);
    1207          }
    1208          results_farey = parallelWaitAll("farey", arguments_farey, 0, n1);
    1209          for(i = ncols(H); i > 0; i--)
    1210          {
    1211             J[i] = results_farey[i][1];
    1212          }
    1213       }
    1214       if(printlevel >= 10)
    1215       {
    1216          "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
    1217          "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
    1218       }
    1219 
    1220 //----------------  Test if we already have a standard basis of I --------------
    1221 
    1222       tt = timer; rt = rtimer;
    1223       if((variant == 1) || (variant == 3))
    1224       {
    1225          pTest = pTestSB(I,J,L,variant,hi);
    1226       }
    1227       if((variant == 2) || (variant == 4))
    1228       {
    1229          pTest = pTestSB(I,J,L,variant);
    1230       }
    1231 
    1232       if(printlevel >= 10)
    1233       {
    1234          "CPU-time for pTest is "+string(timer - tt)+" seconds.";
    1235          "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
    1236       }
    1237 
    1238       if(pTest)
    1239       {
    1240          if(printlevel >= 10)
    1241          {
    1242             "CPU-time for computation without final tests is "
    1243             +string(timer - TT)+" seconds.";
    1244             "Real-time for computation without final tests is "
    1245             +string(rtimer - RT)+" seconds.";
    1246          }
    1247 
    1248          attrib(J,"isSB",1);
    1249 
    1250          if(exactness == 0)
    1251          {
    1252             option(set, opt);
    1253             if(n1 > 1) { kill I_for_fork; }
    1254             return(J);
    1255          }
    1256 
    1257          if(exactness == 1)
    1258          {
    1259             tt = timer; rt = rtimer;
    1260             sizeTest = 1 - isIncluded(I,J,n1);
    1261 
    1262             if(printlevel >= 10)
    1263             {
    1264                "CPU-time for checking if I subset <G> is "
    1265                +string(timer - tt)+" seconds.";
    1266                "Real-time for checking if I subset <G> is "
    1267                +string(rtimer - rt)+" seconds.";
    1268             }
    1269 
    1270             if(sizeTest == 0)
    1271             {
    1272                tt = timer; rt = rtimer;
    1273                K = std(J);
    1274 
    1275                if(printlevel >= 10)
    1276                {
    1277                   "CPU-time for last std-computation is "
    1278                   +string(timer - tt)+" seconds.";
    1279                   "Real-time for last std-computation is "
    1280                   +string(rtimer - rt)+" seconds.";
    1281                }
    1282 
    1283                if(size(reduce(K,J)) == 0)
    1284                {
    1285                   option(set, opt);
    1286                   if(n1 > 1) { kill I_for_fork; }
    1287                   return(J);
    1288                }
    1289             }
    1290          }
    1291       }
    1292 
    1293 //--------------  We do not already have a standard basis of I  ----------------
    1294 //-----------  Therefore do the main computation for more primes  --------------
    1295 
    1296       T1 = H;
    1297       T2 = N;
    1298       index = 2;
    1299 
    1300       j = size(L) + 1;
    1301       tt = timer; rt = rtimer;
    1302       L = primeList(I,n3,L,n1);
    1303       if(printlevel >= 10)
    1304       {
    1305          "CPU-time for primeList: "+string(timer-tt)+" seconds.";
    1306          "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
    1307       }
    1308 
    1309       if(n1 > 1)
    1310       {
    1311          for(i = 1; i <= n1; i++)
    1312          {
    1313             open(l(i));
    1314             if((variant == 1) || (variant == 3))
    1315             {
    1316                write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
    1317                                                      eval(variant), eval(hi))));
    1318             }
    1319             if((variant == 2) || (variant == 4))
    1320             {
    1321                write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
    1322                                                      eval(variant))));
    1323             }
    1324          }
    1325          j = j + n1;
    1326          k = 0;
    1327       }
    1328    }
    1329 }
    1330 example
    1331 { "EXAMPLE:"; echo = 2;
    1332    ring R1 = 0,(x,y,z,t),dp;
    1333    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
    1334    ideal J = modStd(I);
    1335    J;
    1336    I = homog(I,t);
    1337    J = modStd(I);
    1338    J;
    1339 
    1340    ring R2 = 0,(x,y,z),ds;
    1341    ideal I = jacob(x5+y6+z7+xyz);
    1342    ideal J1 = modStd(I,1,0);
    1343    J1;
    1344 
    1345    ring R3 = 0,x(1..4),lp;
    1346    ideal I = cyclic(4);
    1347    ideal J1 = modStd(I,1);
    1348    ideal J2 = modStd(I,1,0);
    1349    size(reduce(J1,J2));
    1350    size(reduce(J2,J1));
    1351 }
    1352 
    1353 ////////////////////////////////////////////////////////////////////////////////
    1354 
    1355 proc modS(ideal I, list L, list #)
    1356 "USAGE:  modS(I,L); I ideal, L intvec of primes
    1357          if size(#)>0 std is used instead of groebner
    1358 RETURN:  an ideal which is with high probability a standard basis
    1359 NOTE:    This procedure is designed for fast experiments.
    1360          It is not tested whether the result is a standard basis.
    1361          It is not tested whether the result generates I.
    1362 EXAMPLE: example modS; shows an example
    1363 "
    1364 {
    1365    int j;
    1366    bigint N = 1;
    1367    def R0 = basering;
    1368    ideal J;
    1369    list T;
    1370    list rl = ringlist(R0);
    1371    if((npars(R0)>0) || (rl[1]>0))
    1372    {
    1373       ERROR("Characteristic of basering should be zero.");
    1374    }
    1375    for(j = 1; j <= size(L); j++)
    1376    {
    1377       N = N*L[j];
    1378       rl[1] = L[j];
    1379       def @r = ring(rl);
    1380       setring @r;
    1381       ideal I = fetch(R0,I);
    1382       if(size(#) > 0)
    1383       {
    1384          I = std(I);
    1385       }
    1386       else
    1387       {
    1388          I = groebner(I);
    1389       }
    1390       setring R0;
    1391       T[j] = fetch(@r,I);
    1392       kill @r;
    1393    }
    1394    L = deleteUnluckyPrimes(T,L,homog(I));
    1395    // unlucky if and only if the leading ideal is wrong
    1396    J = farey(chinrem(L[1],L[2]),N);
    1397    attrib(J,"isSB",1);
    1398    return(J);
    1399 }
    1400 example
    1401 { "EXAMPLE:"; echo = 2;
    1402    list L = 3,5,11,13,181,32003;
    1403    ring r = 0,(x,y,z,t),dp;
    1404    ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
    1405    I = homog(I,t);
    1406    ideal J = modS(I,L);
    1407    J;
    1408 }
    1409 
    1410 ////////////////////////////////////////////////////////////////////////////////
    1411 
    1412 proc modHenselStd(ideal I, list #)
    1413 "USAGE:  modHenselStd(I);
    1414 RETURN:  a standard basis of I;
    1415 NOTE:    The procedure computes a standard basis of I (over the rational
    1416          numbers) by using  modular computations and Hensellifting.
    1417          For further experiments see procedure modS.
    1418 EXAMPLE: example modHenselStd; shows an example
    1419 "
    1420 {
    1421    int i,j;
    1422 
    1423    bigint p = 2134567879;
    1424    if(size(#)!=0) { p=#[1]; }
    1425    while(!primeTest(I,p))
    1426    {
    1427       p = prime(random(2000000000,2134567879));
    1428    }
    1429 
    1430    def R = basering;
    1431    module F,PrevG,PrevZ,Z2;
    1432    ideal testG,testG1,G1,G2,G3,Gp;
    1433    list L = p;
    1434    list rl = ringlist(R);
    1435    rl[1] = int(p);
    1436 
    1437    def S = ring(rl);
    1438    setring S;
    1439    intvec opt = option(get);
    1440    option(redSB);
    1441    module Z,M,Z2;
    1442    ideal I = imap(R,I);
    1443    ideal Gp,G1,G2,G3;
    1444    Gp,Z = liftstd1(I);
    1445    attrib(Gp,"isSB",1);
    1446    module ZZ = syz(I);
    1447    attrib(ZZ,"isSB",1);
    1448    Z = reduce(Z,ZZ);
    1449 
    1450    setring R;
    1451    Gp = imap(S,Gp);
    1452    PrevZ = imap(S,Z);
    1453    PrevG = module(Gp);
    1454    F = module(I);
    1455    testG = farey(Gp,p);
    1456    attrib(testG,"isSB",1);
    1457    while(1)
    1458    {
    1459       i++;
    1460       G1 = ideal(1/(p^i) * sum_id(F*PrevZ,(-1)*PrevG));
    1461       setring S;
    1462       G1 = imap(R,G1);
    1463       G2 = reduce(G1,Gp);
    1464       G3 = sum_id(G1,(-1)*G2);
    1465       M = lift(Gp,G3);
    1466       Z2 = (-1)*Z*M;
    1467 
    1468       setring R;
    1469       G2 = imap(S,G2);
    1470       Z2 = imap(S,Z2);
    1471       PrevG = sum_id(PrevG, module(p^i*G2));
    1472       PrevZ = sum_id(PrevZ, multiply(poly(p^i),Z2));
    1473       testG1 = farey(ideal(PrevG),p^(i+1));
    1474       attrib(testG1,"isSB",1);
    1475       if(size(reduce(testG1,testG)) == 0)
    1476       {
    1477          if(size(reduce(I,testG1)) == 0)       // I is in testG1
    1478          {
    1479             if(pTestSB(I,testG1,L,2))
    1480             {
    1481                G3 = std(testG1);               // testG1 is SB
    1482                if(size(reduce(G3,testG1)) == 0)
    1483                {
    1484                   option(set, opt);
    1485                   return(G3);
    1486                }
    1487             }
    1488          }
    1489       }
    1490       testG = testG1;
    1491       attrib(testG,"isSB",1);
    1492    }
    1493 }
    1494 example
    1495 { "EXAMPLE:"; echo = 2;
    1496    ring r = 0,(x,y,z),dp;
    1497    ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
    1498    ideal J = modHenselStd(I);
    1499    J;
    1500 }
    1501 
    1502 ////////////////////////////////////////////////////////////////////////////////
    1503 
    1504 static proc sum_id(list #)
    1505 {
    1506    if(typeof(#[1])=="ideal")
    1507    {
    1508       ideal M;
    1509    }
    1510    else
    1511    {
    1512       module M;
    1513    }
    1514 
    1515    int i;
    1516    for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; }
    1517    return(M);
    1518 }
    1519 
    1520 ////////////////////////////////////////////////////////////////////////////////
    1521 
    1522 static proc multiply(poly p, list #)
    1523 {
    1524    if(typeof(#[1])=="ideal")
    1525    {
    1526       ideal M;
    1527    }
    1528    else
    1529    {
    1530       module M;
    1531    }
    1532 
    1533    int i;
    1534    for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; }
    1535    return(M);
     244        }
     245    }
     246
     247    /* compute command(args) */
     248    execute("Ip = "+command+"(Ip);");
     249
     250    /* test if Gp is in Ip */
     251    for (i = ncols(Gp); i > 0; i--) {
     252        if (reduce(Gp[i], Ip, 1) != 0) {
     253            setring(br);
     254            return(0);
     255        }
     256    }
     257    setring(br);
     258    return(1);
     259}
     260
     261/* test if 'result' is a GB of the input ideal */
     262static proc finalTest_std(string command, alias list args, ideal result)
     263{
     264    /* test if args[1] is in result */
     265    attrib(result, "isSB", 1);
     266    int i;
     267    for (i = ncols(args[1]); i > 0; i--) {
     268        if (reduce(args[1][i], result, 1) != 0) {
     269            return(0);
     270        }
     271    }
     272
     273    /* test if result is a GB */
     274    ideal G = std(result);
     275    if (reduce_parallel(G, result)) {
     276        return(0);
     277    }
     278    return(1);
     279}
     280
     281/* return 1, if I_reduce is _not_ in G_reduce,
     282 *        0, otherwise
     283 * (same as size(reduce(I_reduce, G_reduce))).
     284 * Uses parallelization. */
     285static proc reduce_parallel(ideal I_reduce, ideal G_reduce)
     286{
     287    exportto(Modstd, I_reduce);
     288    exportto(Modstd, G_reduce);
     289    int size_I = ncols(I_reduce);
     290    int chunks = Modular::par_range(size_I);
     291    intvec range;
     292    int i;
     293    for (i = chunks; i > 0; i--) {
     294        range = Modular::par_range(size_I, i);
     295        task t(i) = "Modstd::reduce_task", list(range);
     296    }
     297    startTasks(t(1..chunks));
     298    waitAllTasks(t(1..chunks));
     299    int result = 0;
     300    for (i = chunks; i > 0; i--) {
     301        if (getResult(t(i))) {
     302            result = 1;
     303            break;
     304        }
     305    }
     306    kill I_reduce;
     307    kill G_reduce;
     308    return(result);
     309}
     310
     311/* compute a chunk of reductions for reduce_parallel */
     312static proc reduce_task(intvec range)
     313{
     314    int result = 0;
     315    int i;
     316    for (i = range[1]; i <= range[2]; i++) {
     317        if (reduce(I_reduce[i], G_reduce, 1) != 0) {
     318            result = 1;
     319            break;
     320        }
     321    }
     322    return(result);
    1536323}
    1537324
Note: See TracChangeset for help on using the changeset viewer.