Changeset 604d8b in git for Singular/LIB/modstd.lib


Ignore:
Timestamp:
Jun 10, 2010, 3:02:22 PM (13 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
b41ecb2f02b86cc0b679c2ecb470330dae622a51
Parents:
6e0750067834212a8910aad022e669ca646b5786
Message:
checking in on behalf of Stefan Steidel

git-svn-id: file:///usr/local/Singular/svn/trunk@12850 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/modstd.lib

    r6e0750 r604d8b  
    11//GP, last modified 23.10.06
    2 ///////////////////////////////////////////////////////////////////////////////
     2////////////////////////////////////////////////////////////////////////////////
    33version="$Id$";
    4 category="Commutative Algebra";
     4category = "Commutative Algebra";
    55info="
    6 LIBRARY: modstd.lib  Grobner basis of ideals
    7 AUTHORS: A. Hashemi,     Amir.Hashemi@lip6.fr
    8 @*       G. Pfister      pfister@mathematik.uni-kl.de
    9 @*       H. Schoenemann  hannes@mathematik.uni-kl.de
    10 
    11 NOTE:
    12  A library for computing the Grobner basis of an ideal in the polynomial
    13  ring over the rational numbers using modular methods. The procedures are
    14  inspired by the following paper:
    15  Elizabeth A. Arnold:
    16  Modular Algorithms for Computing Groebner Bases , Journal of Symbolic
    17  Computation , April 2003, Volume 35, (4), p. 403-419.
    18 
    19 
     6LIBRARY:  modstd.lib  Groebner basis of ideals
     7
     8AUTHORS:  A. Hashemi      Amir.Hashemi@lip6.fr
     9@*        G. Pfister      pfister@mathematik.uni-kl.de
     10@*        H. Schoenemann  hannes@mathematik.uni-kl.de
     11@*        S. Steidel      steidel@mathematik.uni-kl.de
     12
     13OVERVIEW:
     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).
    2020
    2121PROCEDURES:
    22 modStd(I);     compute a standard basis of I using modular methods
    23 modS(I,L);     liftings to Q of standard bases of I mod p for p in L
    24 primeList(n);  intvec of n primes  <= 2134567879 in decreasing order
     22 modStd(I);        standard basis of I using modular methods (chinese remainder)
     23 modHenselStd(I);  standard basis of I using modular methods (hensel lifting)
     24 modS(I,L);        liftings to Q of standard bases of I mod p for p in L
    2525";
    2626
    2727LIB "poly.lib";
    28 LIB "crypto.lib";
    29 ///////////////////////////////////////////////////////////////////////////////
    30 proc pTestSB(ideal I, ideal J, list L)
    31 "USAGE:  pTestSB(I,J,L); I,J ideals, L intvec of primes
    32 RETURN: 1 (resp. 0) if for a randomly chosen prime p not in L
    33         J mod p is (resp. is not) a standard basis of I mod p
    34 EXAMPLE:example primList; shows an example
     28LIB "ring.lib";
     29
     30////////////////////////////////////////////////////////////////////////////////
     31
     32static proc redFork(ideal I, ideal J, int n)
     33{
     34   attrib(J,"isSB",1);
     35   return(reduce(I,J,1));
     36}
     37
     38////////////////////////////////////////////////////////////////////////////////
     39
     40proc isIncluded(ideal I, ideal J, list #)
     41"USAGE:  isIncluded(I,J); I,J ideals
     42RETURN:  1 if J includes I,
     43         0 if there is an element f in I which does not reduce to 0 w.r.t. J.
     44EXAMPLE: example isIncluded; shows an example
    3545"
    3646{
    37   int i,j,k,p;
    38   def R=basering;
    39   list r= ringlist(R);
    40 
    41   while(!j)
    42   {
    43     j=1;
    44     p=prime(random(1000000000,2134567879));
    45     for(i=1;i<=size(L);i++)
    46     {
    47       if(p==L[i]){j=0;break}
    48     }
    49     if(j)
    50     {
    51       for(i=1;i<=ncols(J);i++)
    52       {
    53         for(k=2;k<=size(J[i]);k++)
    54         {
    55           if((denominator(leadcoef(J[i][k])) mod
    56 p)==0){j=0;break}
    57         }
    58         if(!j){break;}
    59       }
    60     }
    61   }
    62   r[1]=p;
    63   def @R=ring(r);
    64   setring @R;
    65   ideal I=imap(R,I);
    66   ideal J=imap(R,J);
    67   attrib(J,"isSB",1);
    68   j=1;
    69   if(size(reduce(I,J))!=0){j=0;}
    70   if(j)
    71   {
    72     ideal K=std(J);
    73     if(size(reduce(K,J))!=0){j=0;}
    74   }
    75   setring R;
    76   return(j);
     47   attrib(J,"isSB",1);
     48   int i,j,k;
     49
     50   if(size(#) > 0)
     51   {
     52      int n = #[1];
     53      if((n > 1) && (n < ncols(I)) && system("with","MP"))
     54      {
     55         for(i = 1; i <= n - 1; i++)
     56         {
     57            link l(i) = "MPtcp:fork";
     58            open(l(i));
     59            write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1)));
     60         }
     61
     62         int t = timer;
     63         if(reduce(I[ncols(I)], J, 1) != 0)
     64         {
     65            for(i = 1; i <= n - 1; i++)
     66            {
     67               close(l(i));
     68            }
     69            return(0);
     70         }
     71         t = timer - t;
     72         if(t > 60) { t = 60; }
     73         int i_sleep = system("sh", "sleep "+string(t));
     74
     75         j = ncols(I) - n;
     76
     77         while(j >= 0)
     78         {
     79            for(i = 1; i <= n - 1; i++)
     80            {
     81               if(status(l(i), "read", "ready"))
     82               {
     83                  if(read(l(i)) != 0)
     84                  {
     85                     for(i = 1; i <= n - 1; i++)
     86                     {
     87                        close(l(i));
     88                     }
     89                     return(0);
     90                  }
     91                  else
     92                  {
     93                     if(j >= 1)
     94                     {
     95                        write(l(i), quote(redFork(eval(I[j]), eval(J), 1)));
     96                        j--;
     97                     }
     98                     else
     99                     {
     100                        k++;
     101                        close(l(i));
     102                     }
     103                  }
     104               }
     105            }
     106            if(k == n - 1)
     107            {
     108               j--;
     109            }
     110            i_sleep = system("sh", "sleep "+string(t));
     111         }
     112         return(1);
     113      }
     114   }
     115
     116   for(i = ncols(I); i >= 1; i--)
     117   {
     118      if(reduce(I[i],J,1) != 0){ return(0); }
     119   }
     120   return(1);
    77121}
    78122example
    79123{ "EXAMPLE:"; echo = 2;
    80    intvec L=2,3,5;
    81124   ring r=0,(x,y,z),dp;
    82    ideal i=x+1,x+y+1;
    83    ideal j=x+1,y;
    84    pTestSB(i,i,L);
    85    pTestSB(i,j,L);
    86 }
    87 
    88 proc primeList(int n, list #)
    89 "USAGE:  primeList(n); (resp. primeList(n,L); )
    90 RETURN: the intvec of n greatest primes  <= 2147483647 (resp. n greatest
    91         primes < L[size(L)] union with L)
    92 EXAMPLE:example primList; shows an example
     125   ideal I = x+1,x+y+1;
     126   ideal J = x+1,y;
     127   isIncluded(I,J);
     128   isIncluded(J,I);
     129   isIncluded(I,J,4);
     130
     131   ring R = 0, x(1..5), dp;
     132   ideal I1 = cyclic(4);
     133   ideal I2 = I1,x(5)^2;
     134   isIncluded(I1,I2,4);
     135}
     136
     137////////////////////////////////////////////////////////////////////////////////
     138
     139proc pTestSB(ideal I, ideal J, list L, int variant, list #)
     140"USAGE:  pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
     141RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
     142         J mod p is (resp. is not) a standard basis of I mod p
     143EXAMPLE: example pTestSB; shows an example
    93144"
    94145{
    95   intvec L;
    96   int i,p;
    97   if(size(#)==0)
    98   {
    99       p=2147483647;
    100       L[1]=p;
    101    }
    102    else
    103   {
    104      L=#[1];
    105      p=prime(L[size(L)]-1);
    106      L[size(L)+1]=p;
    107   }
    108   if(p==2){ERROR("no more primes");}
    109   for(i=2;i<=n;i++)
    110   {
    111     p=prime(p-1);
    112     L[size(L)+1]=p;
    113   }
    114   return(L);
     146   int i,j,k,p;
     147   def R = basering;
     148   list r = ringlist(R);
     149
     150   while(!j)
     151   {
     152      j = 1;
     153      p = prime(random(1000000000,2134567879));
     154      for(i = 1; i <= size(L); i++)
     155      {
     156         if(p == L[i]) { j = 0; break; }
     157      }
     158      if(j)
     159      {
     160         for(i = 1; i <= ncols(I); i++)
     161         {
     162            for(k = 2; k <= size(I[i]); k++)
     163            {
     164               if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; }
     165            }
     166            if(!j){ break; }
     167         }
     168      }
     169      if(j)
     170      {
     171         if(!primeTest(I,p)) { j = 0; }
     172      }
     173   }
     174   r[1] = p;
     175   def @R = ring(r);
     176   setring @R;
     177   ideal I = imap(R,I);
     178   ideal J = imap(R,J);
     179   attrib(J,"isSB",1);
     180
     181   int t = timer;
     182   j = 1;
     183   if(isIncluded(I,J) == 0) { j = 0; }
     184
     185   if(printlevel >= 10)
     186   {
     187      "isIncluded(I,J) takes "+string(timer - t)+" seconds";
     188      "j = "+string(j);
     189   }
     190
     191   t = timer;
     192   if(j)
     193   {
     194      if(size(#) > 0)
     195      {
     196         ideal K = modpStd(I,p,variant,#[1])[1];
     197      }
     198      else
     199      {
     200         ideal K = groebner(I);
     201      }
     202      t = timer;
     203      if(isIncluded(J,K) == 0) { j = 0; }
     204
     205      if(printlevel >= 10)
     206      {
     207         "isIncluded(K,J) takes "+string(timer - t)+" seconds";
     208         "j = "+string(j);
     209      }
     210   }
     211   setring R;
     212   return(j);
    115213}
    116214example
    117215{ "EXAMPLE:"; echo = 2;
    118    intvec L=primeList(10);
     216   intvec L = 2,3,5;
     217   ring r = 0,(x,y,z),dp;
     218   ideal I = x+1,x+y+1;
     219   ideal J = x+1,y;
     220   pTestSB(I,I,L,2);
     221   pTestSB(I,J,L,2);
     222}
     223
     224////////////////////////////////////////////////////////////////////////////////
     225
     226proc deleteUnluckyPrimes(list T, list L, int ho, list #)
     227"USAGE:  deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer
     228RETURN:  lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #), lT ideal
     229NOTE:    - if ho = 1, the polynomials in T are homogeneous, else ho = 0,
     230         - lT is prevalent, i.e. the most appearing leading ideal in T
     231EXAMPLE: example deleteUnluckyPrimes; shows an example
     232"
     233{
     234   ho = ((ho)||(ord_test(basering) == -1));
     235   int j,k,c;
     236   intvec hl,hc;
     237   ideal cT,lT,cK;
     238   lT = lead(T[size(T)]);
     239   attrib(lT,"isSB",1);
     240   if(!ho)
     241   {
     242      for(j = 1; j < size(T); j++)
     243      {
     244         cT = lead(T[j]);
     245         attrib(cT,"isSB",1);
     246         if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
     247         {
     248            cK = cT;
     249            c++;
     250         }
     251      }
     252      if(c > size(T)/2){ lT = cK; }
     253   }
     254   else
     255   {
     256      hl = hilb(lT,1);
     257      for(j = 1; j < size(T); j++)
     258      {
     259         cT = lead(T[j]);
     260         attrib(cT,"isSB",1);
     261         hc = hilb(cT,1);
     262         if(hl == hc)
     263         {
     264            for(k = 1; k <= size(lT); k++)
     265            {
     266               if(lT[k] < cT[k]) { lT = cT; c++; break; }
     267               if(lT[k] > cT[k]) { c++; break; }
     268            }
     269         }
     270         else
     271         {
     272            if(hc < hl){ lT = cT; hl = hilb(lT,1); c++ }
     273         }
     274      }
     275   }
     276
     277   int addList;
     278   if(size(#) > 0) { list M = #; addList = 1; }
     279   j = 1;
     280   attrib(lT,"isSB",1);
     281   while((j <= size(T))&&(c > 0))
     282   {
     283      cT = lead(T[j]);
     284      attrib(cT,"isSB",1);
     285      if((size(reduce(cT,lT)) != 0)||(size(reduce(lT,cT)) != 0))
     286      {
     287         T = delete(T,j);
     288         if(j == 1)
     289         {
     290            L = L[2..size(L)];
     291            if(addList == 1) { M = M[2..size(M)]; }
     292         }
     293         else
     294         {
     295            if(j == size(L))
     296            {
     297               L = L[1..size(L)-1];
     298               if(addList == 1) { M = M[1..size(M)-1]; }
     299            }
     300            else
     301            {
     302               L = L[1..j-1],L[j+1..size(L)];
     303               if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; }
     304            }
     305         }
     306         j--;
     307      }
     308      j++;
     309   }
     310
     311   for(j = 1; j <= size(L); j++)
     312   {
     313      L[j] = bigint(L[j]);
     314   }
     315
     316   if(addList == 0) { return(list(T,L,lT)); }
     317   if(addList == 1) { return(list(T,L,M,lT)); }
     318}
     319example
     320{ "EXAMPLE:"; echo = 2;
     321   list L = 2,3,5,7,11;
     322   ring r = 0,(y,x),Dp;
     323   ideal I1 = 2y2x,y6;
     324   ideal I2 = yx2,y3x,x5,y6;
     325   ideal I3 = y2x,x3y,x5,y6;
     326   ideal I4 = y2x,11x3y,x5;
     327   ideal I5 = y2x,yx3,x5,7y6;
     328   list T = I1,I2,I3,I4,I5;
     329   deleteUnluckyPrimes(T,L,1);
     330   list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5);
     331   deleteUnluckyPrimes(T,L,1,P);
     332}
     333
     334////////////////////////////////////////////////////////////////////////////////
     335
     336proc primeTest(ideal I, bigint p)
     337{
     338   int i,j;
     339   for(i = 1; i <= size(I); i++)
     340   {
     341      for(j = 1; j <= size(I[i]); j++)
     342      {
     343         if((leadcoef(I[i][j]) mod p) == 0) { return(0); }
     344      }
     345   }
     346   return(1);
     347}
     348
     349////////////////////////////////////////////////////////////////////////////////
     350
     351proc primeList(ideal I, int n, list #)
     352"USAGE:  primeList(I,n); ( resp. primeList(I,n,L); ) I ideal, n integer
     353RETURN:  the intvec of n greatest primes  <= 2147483647 (resp. n greatest primes
     354         < L[size(L)] union with L) such that none ot these primes divides any
     355         coefficient occuring in I
     356EXAMPLE: example primList; shows an example
     357"
     358{
     359   intvec L;
     360   int i,p;
     361   if(size(#) == 0)
     362   {
     363      p = 2147483647;
     364      while(!primeTest(I,p))
     365      {
     366         p = prime(p-1);
     367         if(p == 2) { ERROR("no more primes"); }
     368      }
     369      L[1] = p;
     370   }
     371   else
     372   {
     373      L = #[1];
     374      p = prime(L[size(L)]-1);
     375      while(!primeTest(I,p))
     376      {
     377         p = prime(p-1);
     378         if(p == 2) { ERROR("no more primes"); }
     379      }
     380      L[size(L)+1] = p;
     381   }
     382   if(p == 2) { ERROR("no more primes"); }
     383   for(i = 2; i <= n; i++)
     384   {
     385      p = prime(p-1);
     386      while(!primeTest(I,p))
     387      {
     388         p = prime(p-1);
     389         if(p == 2) { ERROR("no more primes"); }
     390      }
     391      L[size(L)+1] = p;
     392   }
     393   return(L);
     394}
     395example
     396{ "EXAMPLE:"; echo = 2;
     397   ring r = 0,(x,y,z),dp;
     398   ideal I = 2147483647x+y, z-181;
     399   intvec L = primeList(I,10);
     400   size(L);
     401   L[1];
     402   L[size(L)];
     403   L = primeList(I,5,L);
    119404   size(L);
    120405   L[size(L)];
    121    L=primeList(5,L);
    122    size(L);
    123    L[size(L)];
    124 }
    125 
    126 
    127 proc modStd(ideal I)
    128 "USAGE:  modStd(I);
     406}
     407
     408////////////////////////////////////////////////////////////////////////////////
     409
     410proc liftstd1(ideal I)
     411{
     412   def R = basering;
     413   list rl = ringlist(R);
     414   list ordl = rl[3];
     415
     416   int i;
     417   for(i = 1; i <= size(ordl); i++)
     418   {
     419      if((ordl[i][1] == "C") || (ordl[i][1] == "c"))
     420      {
     421         ordl = delete(ordl, i);
     422         break;
     423      }
     424   }
     425
     426   ordl = insert(ordl, list("c", 0));
     427   rl[3] = ordl;
     428   def newR = ring(rl);
     429   setring newR;
     430   ideal I = imap(R,I);
     431
     432   option(none);
     433   option(prompt);
     434
     435   module M;
     436   for(i = 1; i <= size(I); i++)
     437   {
     438      M = M + module(I[i]*gen(1) + gen(i+1));
     439      M = M + module(gen(i+1));
     440   }
     441
     442   module sM = std(M);
     443
     444   ideal sI;
     445   if(attrib(R,"global"))
     446   {
     447      for(i = size(I)+1; i <= size(sM); i++)
     448      {
     449         sI[size(sI)+1] = sM[i][1];
     450      }
     451      matrix T = submat(sM,2..nrows(sM),size(I)+1..ncols(sM));
     452   }
     453   else
     454   {
     455      "==========================================================";
     456      "WARNING: Algorithm is not applicable if ordering is mixed.";
     457      "==========================================================";
     458      for(i = 1; i <= size(sM)-size(I); i++)
     459      {
     460         sI[size(sI)+1] = sM[i][1];
     461      }
     462      matrix T = submat(sM,2..nrows(sM),1..ncols(sM)-size(I));
     463   }
     464
     465   setring R;
     466   return(imap(newR,sI),imap(newR,T));
     467}
     468example
     469{ "EXAMPLE:"; echo = 2;
     470   ring R = 0,(x,y,z),dp;
     471   poly f = x3+y7+z2+xyz;
     472   ideal i = jacob(f);
     473   matrix T;
     474   ideal sm = liftstd(i,T);
     475   sm;
     476   print(T);
     477   matrix(sm) - matrix(i)*T;
     478
     479
     480   ring S = 32003, x(1..6), lp;
     481   ideal I = cyclic(6);
     482   ideal sI;
     483   matrix T;
     484   sI,T = liftstd1(I);
     485   matrix(sI) - matrix(I)*T;
     486}
     487
     488////////////////////////////////////////////////////////////////////////////////
     489
     490proc modpStd(ideal I, int p, int variant, list #)
     491"USAGE:  modpStd(I,p,variant,#); I ideal, p integer, variant integer
     492ASSUME:  If size(#) > 0, then #[1] is an intvec describing the Hilbert series.
     493RETURN:  ideal - a standard basis of I mod p, integer - p (, matrix - if variant
     494         = 5, the transformation matrix obtained from liftstd)
     495NOTE:    The procedure computes a standard basis of the ideal I modulo p and
     496         fetches the result to the basering. If size(#) > 0 the Hilbert driven
     497         standard basis computation std(.,#[1]) is used instead of groebner.
     498         The standard basis computation modulo p does also vary depending on the
     499         integer variant, namely
     500@*       - variant = 1: std(.,#[1]) resp. groebner,
     501@*       - variant = 2: groebner,
     502@*       - variant = 3: homogenize - std(.,#[1]) resp. groebner - dehomogenize,
     503@*       - variant = 4: std(.,#[1]) resp. groebner,
     504@*       - variant = 5: liftstd.
     505EXAMPLE: example modpStd; shows an example
     506"
     507{
     508   def R0 = basering;
     509   list rl = ringlist(R0);
     510   rl[1] = p;
     511   def @r = ring(rl);
     512   setring @r;
     513   ideal i = fetch(R0,I);
     514
     515   option(redSB);
     516
     517   int t = timer;
     518   if((variant == 1) || (variant == 4))
     519   {
     520      if(size(#) > 0)
     521      {
     522         i = std(i, #[1]);
     523      }
     524      else
     525      {
     526         i = groebner(i);
     527      }
     528   }
     529   if(variant == 2)
     530   {
     531      i = groebner(i);
     532   }
     533
     534   if(variant == 3)
     535   {
     536      list rl = ringlist(@r);
     537      int nvar@r = nvars(@r);
     538
     539      int k;
     540      intvec w;
     541      for(k = 1; k <= nvar@r; k++)
     542      {
     543         w[k] = deg(var(k));
     544      }
     545      w[nvar@r + 1] = 1;
     546
     547      rl[2][nvar@r + 1] = "homvar";
     548      rl[3][2][2] = w;
     549
     550      def HomR = ring(rl);
     551      setring HomR;
     552      ideal i = imap(@r, i);
     553      i = homog(i, homvar);
     554
     555      if(size(#) > 0)
     556      {
     557         if(w == 1)
     558         {
     559            i = std(i, #[1]);
     560         }
     561         else
     562         {
     563            i = std(i, #[1], w);
     564         }
     565      }
     566      else
     567      {
     568         i = groebner(i);
     569      }
     570
     571      t = timer;
     572      i = subst(i, homvar, 1);
     573      i = simplify(i, 34);
     574
     575      setring @r;
     576      i = imap(HomR, i);
     577      i = interred(i);
     578      kill HomR;
     579   }
     580
     581   if(variant == 5)
     582   {
     583      matrix trans;
     584      i,trans = liftstd1(i);
     585      setring R0;
     586      return(list(fetch(@r,i),p,fetch(@r,trans)));
     587   }
     588
     589   setring R0;
     590   return(list(fetch(@r,i),p));
     591}
     592example
     593{ "EXAMPLE:"; echo = 2;
     594   ring r = 0, x(1..4), dp;
     595   ideal I = cyclic(4);
     596   int p = 181;
     597   list P = modpStd(I,p,5);
     598   P;
     599   matrix(P[1])-matrix(I)*P[3];
     600
     601   int q = 32003;
     602   list Q = modpStd(I,q,2);
     603   Q;
     604}
     605
     606////////////////////////////// main procedures /////////////////////////////////
     607
     608proc modStd(ideal I, list #)
     609"USAGE:  modStd(I); I ideal
     610ASSUME:  If size(#) > 0, then # contains either 1, 2 or 4 integers such that
     611@*       - #[1] is the number of available processors for the computation,
     612@*       - #[2] is an optional parameter for the exactness of the computation,
     613                if #[2] = 1, the procedure computes a standard basis for sure,
     614@*       - #[3] is the number of primes until the first lifting,
     615@*       - #[4] is the constant number of primes between two liftings until
     616           the computation stops.
    129617RETURN:  a standard basis of I if no warning appears;
    130 NOTE:    the procedure computes a standard basis of I (over the
    131          rational numbers) by using  modular methods. If a
    132          warning appears then the result is a standard basis
    133          containing I and with high probability a standard basis of I.
    134          For further experiments see procedure modS.
     618NOTE:    The procedure computes a standard basis of I (over the rational
     619         numbers) by using  modular methods. If a warning appears then the
     620         result is a standard basis containing I and with high probability
     621         a standard basis of I.
     622         By default the procedure computes a standard basis of I with high
     623         probability, but if the optional parameter #[2] = 1, it is exact.
     624         The procedure distinguishes between different variants for the standard
     625         basis computation in positive characteristic depending on the ordering
     626         of the basering, the parameter #[2] and if the ideal I is homogeneous.
     627@*       - variant = 1, if I is homogeneous and exactness = 0,
     628@*       - variant = 2, if I is not homogeneous, 1-block-ordering and
     629                        exactness = 0,
     630@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
     631                        > 1 block) and exactness = 0,
     632@*       - variant = 4, if I is homogeneous and exactness = 1,
     633@*       - variant = 5, if I is not homogeneous and exactness = 1.
    135634EXAMPLE: example modStd; shows an example
    136635"
    137636{
    138   def R0=basering;
    139   list rl=ringlist(R0);
    140   if((npars(R0)>0)||(rl[1]>0))
    141   {
    142      ERROR("characteristic of basering should be zero, basering should have
    143 no parameters");
    144   }
    145 
    146   int k,c;
    147   int pd=printlevel-voice+2;
    148   int j=1;
    149   int h=homog(I);
    150   int en=2134567879;
    151   int an=1000000000;
    152 
    153   intvec opt = option(get);     // Save current options
    154   intvec L=primeList(10);
    155   L[5]=prime(random(an,en));
    156   list T,TT,TH,LL;
    157 
    158 
    159   ideal J,K;
    160 
    161   option(redSB);
    162 
    163   while(1)
    164   {
    165      while(j<=size(L))
    166      {
    167         if(pd>2){c++;c;}
    168         rl[1]=L[j];
    169         def @r=ring(rl);
    170         setring @r;
    171         ideal i=fetch(R0,I);
    172         i=groebner(i);
    173         setring R0;
    174         T[size(T)+1]=fetch(@r,i);
    175         kill @r;
    176         j++;
    177      }
    178      if(pd>2){"lifting";}
    179   //================= delete unlucky primes ====================
    180   // unlucky if and only if the leading ideal is wrong
    181      LL=deleteUnluckyPrimes(T,L);
    182      TH=LL[1];
    183      L=LL[2];
    184   //============ now all leading ideals are the same ============
    185      for(j=1;j<=ncols(TH[1]);j++)
    186      {
    187          for(k=1;k<=size(L);k++)
    188          {
    189              TT[k]=TH[k][j];
    190          }
    191          J[j]=liftPoly(TT,L);
    192      }
    193      if(pd>2){"list of primes:";L;"pTest";}
    194      if(pTestSB(I,J,L))
    195      {
    196         attrib(J,"isSB",1);
    197         K=std(J);
    198         if(size(reduce(I,J))==0)
    199         {
    200            if(size(reduce(K,J))==0)
    201            {
    202                if(!((h)||(ord_test(R0)==-1)))
    203               {
    204 
    205 "==================================================================";
    206                  "    The input is not homogeneous and the ordering is not
    207 local.   ";
    208                  "WARNING: ideal generated by output may be greater then
    209 input ideal";
    210 
    211 "==================================================================";
    212               }
    213               option(set, opt);
    214               return(J);
    215            }
    216            if(pd>2){"pTest o.k. but result wrong";}
     637   int TT = timer;
     638   int RT = rtimer;
     639
     640   def R0 = basering;
     641   list rl = ringlist(R0);
     642   if((npars(R0) > 0) || (rl[1] > 0))
     643   {
     644      ERROR("characteristic of basering should be zero, basering should have no parameters");
     645   }
     646
     647   int index = 1;
     648   int i,k,c;
     649   int pd = printlevel-voice+2;
     650   int j = 1;
     651   int pTest;
     652   int en = 2134567879;
     653   int an = 1000000000;
     654   bigint N;
     655
     656//--------------------  Initialize optional parameters  ------------------------
     657   if(size(#) > 0)
     658   {
     659      if(size(#) == 1)
     660      {
     661         int n1 = #[1];
     662         if((n1 > 1) && (1 - system("with","MP")))
     663         {
     664            "========================================================================";
     665            "There is no MP available on your system. Since this is necessary to     ";
     666            "parallelize the algorithm, the computation will be done without forking.";
     667            "========================================================================";
     668            n1 = 1;
     669         }
     670         int exactness = 0;
     671         int n2 = 10;
     672         int n3 = 10;
     673      }
     674      if(size(#) == 2)
     675      {
     676         int n1 = #[1];
     677         if((n1 > 1) && (1 - system("with","MP")))
     678         {
     679            "========================================================================";
     680            "There is no MP available on your system. Since this is necessary to     ";
     681            "parallelize the algorithm, the computation will be done without forking.";
     682            "========================================================================";
     683            n1 = 1;
     684         }
     685         int exactness = #[2];
     686         int n2 = 10;
     687         int n3 = 10;
     688      }
     689      if(size(#) == 4)
     690      {
     691         int n1 = #[1];
     692         if((n1 > 1) && (1 - system("with","MP")))
     693         {
     694            "========================================================================";
     695            "There is no MP available on your system. Since this is necessary to     ";
     696            "parallelize the algorithm, the computation will be done without forking.";
     697            "========================================================================";
     698            n1 = 1;
     699         }
     700         int exactness = #[2];
     701         int n2 = #[3];
     702         int n3 = #[4];
     703      }
     704   }
     705   else
     706   {
     707      int n1 = 1;
     708      int exactness = 0;
     709      int n2 = 10;
     710      int n3 = 10;
     711   }
     712
     713//-------------------------  Save current options  -----------------------------
     714   intvec opt = option(get);
     715
     716   option(redSB);
     717
     718//--------------------  Initialize the list of primes  -------------------------
     719   intvec L = primeList(I,n2);
     720   L[5] = prime(random(an,en));
     721
     722//---------------------  Decide which variant to take  -------------------------
     723   int variant;
     724   int h = homog(I);
     725
     726   int tt = timer;
     727   int rt = rtimer;
     728
     729   if(h)
     730   {
     731      if(exactness == 0) { variant = 1; if(printlevel >= 10) { "variant = 1"; } }
     732      if(exactness == 1) { variant = 4; if(printlevel >= 10) { "variant = 4"; } }
     733      rl[1] = L[5];
     734      def @r = ring(rl);
     735      setring @r;
     736      def @s = changeord("dp");
     737      setring @s;
     738      ideal I = std(fetch(R0,I));
     739      intvec hi = hilb(I,1);
     740      setring R0;
     741      kill @r,@s;
     742   }
     743   else
     744   {
     745      if(exactness == 0)
     746      {
     747         string ordstr_R0 = ordstr(R0);
     748         int neg = 1 - attrib(R0,"global");
     749
     750         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
     751         {
     752            variant = 2;
     753            if(printlevel >= 10) { "variant = 2"; }
     754         }
     755         else
     756         {
     757            string order;
     758            if(system("nblocks") <= 2)
     759            {
     760               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") + find(ordstr_R0, "rp") <= 0)
     761               {
     762                  order = "simple";
     763               }
     764            }
     765
     766            if((order == "simple") || (size(rl) > 4))
     767            {
     768               variant = 2;
     769               if(printlevel >= 10) { "variant = 2"; }
     770            }
     771            else
     772            {
     773               variant = 3;
     774               if(printlevel >= 10) { "variant = 3"; }
     775               rl[1] = L[5];
     776               def @r = ring(rl);
     777               setring @r;
     778               int nvar@r = nvars(@r);
     779               intvec w;
     780               for(i = 1; i <= nvar@r; i++)
     781               {
     782                  w[i] = deg(var(i));
     783               }
     784               w[nvar@r + 1] = 1;
     785
     786               list hiRi = hilbRing(fetch(R0,I),w);
     787               intvec W = hiRi[2];
     788               def @s = hiRi[1];
     789               setring @s;
     790
     791               Id(1) = std(Id(1));
     792               intvec hi = hilb(Id(1), 1, W);
     793
     794               setring R0;
     795               kill @r,@s;
     796            }
     797         }
     798      }
     799
     800      if(exactness == 1)
     801      {
     802         variant = 5;
     803         if(printlevel >= 10) { "variant = 2"; }
     804         matrix trans;
     805      }
     806   }
     807
     808   list P,T1,T2,T3,LL;
     809
     810   ideal J,K,H;
     811
     812//-----  If there is more than one processor available, we parallelize the  ----
     813//-----  main standard basis computations in positive characteristic        ----
     814
     815   if(n1 > 1)
     816   {
     817      ideal I_for_fork = I;
     818      export(I_for_fork);           // I available for each link
     819
     820//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
     821//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
     822
     823      for(i = 1; i <= n1; i++)
     824      {
     825         link l(i) = "MPtcp:fork";
     826         open(l(i));
     827         if((variant == 1) || (variant == 3) || (variant == 4))
     828         {
     829            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), eval(variant), eval(hi))));
     830         }
     831         if((variant == 2) || (variant == 5))
     832         {
     833            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), eval(variant))));
     834         }
     835      }
     836
     837      int t = timer;
     838      if((variant == 1) || (variant == 3) || (variant == 4))
     839      {
     840         P = modpStd(I_for_fork, L[1], variant, hi);
     841      }
     842      if((variant == 2) || (variant == 5))
     843      {
     844         P = modpStd(I_for_fork, L[1], variant);
     845      }
     846      t = timer - t;
     847      if(t > 60) { t = 60; }
     848      int i_sleep = system("sh", "sleep "+string(t));
     849      T1[1] = P[1];
     850      T2[1] = bigint(P[2]);
     851      if(variant == 5) { T3[1] = P[3]; }
     852      index++;
     853
     854      j = j + n1 + 1;
     855   }
     856
     857//--------------  Main standard basis computations in positive  ----------------
     858//----------------------  characteristic start here  ---------------------------
     859
     860   while(1)
     861   {
     862      tt = timer; rt = rtimer;
     863
     864      if(n1 > 1)
     865      {
     866         if(printlevel >= 10) { "size(L) = "+string(size(L)); }
     867         while(j <= size(L) + 1)
     868         {
     869            for(i = 1; i <= n1; i++)
     870            {
     871               if(status(l(i), "read", "ready"))                         // ask if link l(i) is ready otherwise sleep for t seconds
     872               {
     873                  P = read(l(i));                                        // read the result from l(i)
     874                  T1[index] = P[1];
     875                  T2[index] = bigint(P[2]);
     876                  if(variant == 5) { T3[index] = P[3]; }
     877                  index++;
     878
     879                  if(j <= size(L))
     880                  {
     881                     if((variant == 1) || (variant == 3) || (variant == 4))
     882                     {
     883                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]), eval(variant), eval(hi))));
     884                        j++;
     885                     }
     886                     if((variant == 2) || (variant == 5))
     887                     {
     888                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]), eval(variant))));
     889                        j++;
     890                     }
     891                  }
     892                  else
     893                  {
     894                     k++;
     895                     close(l(i));
     896                  }
     897               }
     898            }
     899            if(k == n1)                                                   // k describes the number of closed links
     900            {
     901               j++;
     902            }
     903            i_sleep = system("sh", "sleep "+string(t));
     904         }
     905      }
     906      else
     907      {
     908         while(j <= size(L))
     909         {
     910            if((variant == 1) || (variant == 3) || (variant == 4))
     911            {
     912               P = modpStd(I, L[j], variant, hi);
     913            }
     914            if((variant == 2) || (variant == 5))
     915            {
     916               P = modpStd(I, L[j], variant);
     917            }
     918
     919            T1[index] = P[1];
     920            T2[index] = bigint(P[2]);
     921            if(variant == 5) { T3[index] = P[3]; }
     922            index++;
     923            j++;
     924         }
     925      }
     926
     927      if(printlevel >= 10)
     928      {
     929         "CPU-time for computing list is "+string(timer - tt)+" seconds.";
     930         "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
     931      }
     932
     933      if(pd>2){"lifting";}
     934
     935//------------------------  Delete unlucky primes  -----------------------------
     936//-------------  unlucky if and only if the leading ideal is wrong  ------------
     937
     938      if(variant == 5) { LL = deleteUnluckyPrimes(T1,T2,h,T3); T3 = LL[3]; }
     939      else { LL = deleteUnluckyPrimes(T1,T2,h); }
     940      T1 = LL[1];
     941      T2 = LL[2];
     942
     943//-------------------  Now all leading ideals are the same  --------------------
     944//-------------------  Lift results to basering via farey  ---------------------
     945
     946      N = T2[1];
     947      for(i = 2; i <= size(T2); i++){N = N*T2[i];}
     948      H = chinrem(T1,T2);
     949      J = farey(H,N);
     950      if(variant == 5) { trans = farey(chinrem(T3,T2), N); }
     951
     952//----------------  Test if we already have a standard basis of I --------------
     953
     954      tt = timer; rt = rtimer;
     955      if(pd > 2){ "list of primes:"; L; "pTest" ;}
     956      if((variant == 1) || (variant == 3) || (variant == 4)) { pTest = pTestSB(I,J,L,variant,hi); }
     957      if((variant == 2) || (variant == 5)) { pTest = pTestSB(I,J,L,variant); }
     958
     959      if(printlevel >= 10)
     960      {
     961         "CPU-time for pTest is "+string(timer - tt)+" seconds.";
     962         "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
     963      }
     964
     965      if(pTest)
     966      {
     967         if(printlevel >= 10)
     968         {
     969            "CPU-time for computation without final tests is "+string(timer - TT)+" seconds.";
     970            "Real-time for computation without final tests is "+string(rtimer - RT)+" seconds.";
     971         }
     972
     973         attrib(J,"isSB",1);
     974         tt = timer; rt = rtimer;
     975         int sizeTest = 1 - isIncluded(I,J,n1);
     976
     977         if(printlevel >= 10)
     978         {
     979            "CPU-time for checking if I subset <G> is "+string(timer - tt)+" seconds.";
     980            "Real-time for checking if I subset <G> is "+string(rtimer - rt)+" seconds.";
     981         }
     982
     983         if(sizeTest == 0)
     984         {
     985            if(variant == 1)
     986            {
     987               "===================================================================";
     988               "WARNING: Ideal generated by output may be greater than input ideal.";
     989               "===================================================================";
     990               option(set, opt);
     991               if(n1 > 1) { kill I_for_fork; }
     992               return(J);
     993            }
     994            if((variant == 2) || (variant == 3))
     995            {
     996               tt = timer; rt = rtimer;
     997               K = std(J);
     998
     999               if(printlevel >= 10)
     1000               {
     1001                  "CPU-time for last std-computation is "+string(timer - tt)+" seconds.";
     1002                  "Real-time for last std-computation is "+string(rtimer - rt)+" seconds.";
     1003               }
     1004
     1005               if(size(reduce(K,J)) == 0)
     1006               {
     1007                  "===================================================================";
     1008                  "WARNING: Ideal generated by output may be greater than input ideal.";
     1009                  "===================================================================";
     1010                  option(set, opt);
     1011                  if(n1 > 1) { kill I_for_fork; }
     1012                  return(J);
     1013               }
     1014            }
     1015            if(variant == 4)
     1016            {
     1017               tt = timer; rt = rtimer;
     1018               K = std(J);
     1019
     1020               if(printlevel >= 10)
     1021               {
     1022                  "CPU-time for last std-computation is "+string(timer - tt)+" seconds.";
     1023                  "Real-time for last std-computation is "+string(rtimer - rt)+" seconds.";
     1024               }
     1025
     1026               if(size(reduce(K,J)) == 0)
     1027               {
     1028                  option(set, opt);
     1029                  if(n1 > 1) { kill I_for_fork; }
     1030                  return(J);
     1031               }
     1032            }
     1033            if(variant == 5)
     1034            {
     1035               tt = timer; rt = rtimer;
     1036               K = std(J);
     1037
     1038               if(printlevel >= 10)
     1039               {
     1040                  "CPU-time for last std-computation is "+string(timer - tt)+" seconds.";
     1041                  "Real-time for last std-computation is "+string(rtimer - rt)+" seconds.";
     1042               }
     1043
     1044               if(size(reduce(K,J)) == 0)
     1045               {
     1046                  if(matrix(J) == matrix(I)*trans)
     1047                  {
     1048                     option(set, opt);
     1049                     if(n1 > 1) { kill I_for_fork; }
     1050                     return(J);
     1051                  }
     1052               }
     1053            }
     1054            if(pd>2){"pTest o.k. but result wrong";}
    2171055         }
    2181056         if(pd>2){"pTest o.k. but result wrong";}
    2191057      }
    220       j=size(L)+1;
    221       L=primeList(10,L);
     1058
     1059//--------------  We do not already have a standard basis of I  ----------------
     1060//-----------  Therefore do the main computation for more primes  --------------
     1061
     1062      T1 = H;
     1063      T2 = N;
     1064      index = 2;
     1065
     1066      j = size(L) + 1;
     1067      L = primeList(I,n3,L);
     1068
     1069      if(n1 > 1)
     1070      {
     1071         for(i = 1; i <= n1; i++)
     1072         {
     1073            open(l(i));
     1074            if((variant == 1) || (variant == 3) || (variant == 4))
     1075            {
     1076               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]), eval(variant), eval(hi))));
     1077            }
     1078            if((variant == 2) || (variant == 5))
     1079            {
     1080               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]), eval(variant))));
     1081            }
     1082         }
     1083         j = j + n1;
     1084         k = 0;
     1085      }
    2221086   }
    2231087}
    2241088example
    2251089{ "EXAMPLE:"; echo = 2;
    226    ring r=0,(x,y,z,t),dp;
    227    ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
    228    ideal J=modStd(I);
     1090   ring r = 0,(x,y,z,t),dp;
     1091   ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     1092   ideal J = modStd(I);
    2291093   J;
    230    I=homog(I,t);
    231    J=modStd(I);
     1094   I = homog(I,t);
     1095   J = modStd(I);
    2321096   J;
    233    ring s=0,(x,y,z),ds;
    234    ideal I=jacob(x5+y6+z7+xyz);
    235    ideal J=modStd(I);
    236    J;
    237 }
    238 ///////////////////////////////////////////////////////////////////////////////
    239 proc modS(ideal I, intvec L, list #)
     1097
     1098   ring s = 0,(x,y,z),ds;
     1099   ideal I = jacob(x5+y6+z7+xyz);
     1100   ideal J1 = modStd(I,1,1);
     1101   J1;
     1102   ideal J2 = modStd(I,3);
     1103   J2;
     1104   size(reduce(J1,J2));
     1105   size(reduce(J2,J1));
     1106
     1107   ring rr = 0,x(1..4),lp;
     1108   ideal I = cyclic(4);
     1109   ideal J1 = modStd(I,2);
     1110   ideal J2 = modStd(I,2,1);
     1111   size(reduce(J1,J2));
     1112   size(reduce(J2,J1));
     1113}
     1114
     1115////////////////////////////////////////////////////////////////////////////////
     1116
     1117proc modS(ideal I, list L, list #)
    2401118"USAGE:  modS(I,L); I ideal, L intvec of primes
    2411119         if size(#)>0 std is used instead of groebner
    … …  
    2471125"
    2481126{
    249   int j,k;
    250   list T,TT;
    251   def R0=basering;
    252   ideal J,cT,lT,K;
    253   ideal I0=I;
    254   list rl=ringlist(R0);
    255   if((npars(R0)>0)||(rl[1]>0))
    256   {
    257      ERROR("characteristic of basering should be zero");
    258   }
    259   for (j=1;j<=size(L);j++)
    260   {
    261     rl[1]=L[j];
    262     def @r=ring(rl);
    263     setring @r;
    264     ideal i=fetch(R0,I);
    265     option(redSB);
    266     if(size(#)>0)
    267     {
    268        i=std(i);
    269     }
    270     else
    271     {
    272        i=groebner(i);
    273     }
    274     setring R0;
    275     T[j]=fetch(@r,i);
    276     kill @r;
    277   }
    278   //================= delete unlucky primes ====================
    279   // unlucky if and only if the leading ideal is wrong
    280   list LL=deleteUnluckyPrimes(T,L);
    281   T=LL[1];
    282   L=LL[2];
    283   //============ now all leading ideals are the same ============
    284   for(j=1;j<=ncols(T[1]);j++)
    285   {
    286     for(k=1;k<=size(L);k++)
    287     {
    288       TT[k]=T[k][j];
    289     }
    290     J[j]=liftPoly(TT,L);
    291   }
    292   attrib(J,"isSB",1);
    293   return(J);
     1127   int j;
     1128   bigint N = 1;
     1129   def R0 = basering;
     1130   ideal J;
     1131   list T;
     1132   list rl = ringlist(R0);
     1133   if((npars(R0)>0) || (rl[1]>0)) { ERROR("characteristic of basering should be zero"); }
     1134   for(j = 1; j <= size(L); j++)
     1135   {
     1136      N = N*L[j];
     1137      rl[1] = L[j];
     1138      def @r = ring(rl);
     1139      setring @r;
     1140      ideal I = fetch(R0,I);
     1141      if(size(#) > 0)
     1142      {
     1143         I = std(I);
     1144      }
     1145      else
     1146      {
     1147         I = groebner(I);
     1148      }
     1149      setring R0;
     1150      T[j] = fetch(@r,I);
     1151      kill @r;
     1152   }
     1153   L = deleteUnluckyPrimes(T,L,homog(I)); // unlucky if and only if the leading ideal is wrong
     1154   J = farey(chinrem(L[1],L[2]),N);
     1155   attrib(J,"isSB",1);
     1156   return(J);
    2941157}
    2951158example
    2961159{ "EXAMPLE:"; echo = 2;
    297    intvec L=3,5,11,13,181;
    298    ring r=0,(x,y,z),dp;
    299    ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
    300    ideal J=modS(I,L);
     1160   list L = 3,5,11,13,181,32003;
     1161   ring r = 0,(x,y,z,t),dp;
     1162   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
     1163   I = homog(I,t);
     1164   ideal J = modS(I,L);
    3011165   J;
    3021166}
    303 ///////////////////////////////////////////////////////////////////////////////
    304 proc deleteUnluckyPrimes(list T,intvec L)
    305 "USAGE:  deleteUnluckyPrimes(T,L);T list of polys, L intvec of primes
    306 RETURN:  list L,T with T list of polys, L intvec of primes
    307 EXAMPLE: example deleteUnluckyPrimes; shows an example
    308 NOTE:    works only for homogeneous ideals with global orderings or
    309          for ideals with local orderings
     1167
     1168////////////////////////////////////////////////////////////////////////////////
     1169
     1170proc modHenselStd(ideal I, list #)
     1171"USAGE:  modHenselStd(I);
     1172RETURN:  a standard basis of I;
     1173NOTE:    The procedure computes a standard basis of I (over the rational
     1174         numbers) by using  modular computations and Hensellifting.
     1175         For further experiments see procedure modS.
     1176EXAMPLE: example modHenselStd; shows an example
    3101177"
    3111178{
    312   int j,k;
    313   intvec hl,hc;
    314   ideal cT,lT;
    315 
    316   lT=lead(T[size(T)]);
    317   attrib(lT,"isSB",1);
    318   hl=hilb(lT,1);
    319   for (j=1;j<size(T);j++)
    320   {
    321      cT=lead(T[j]);
    322      attrib(cT,"isSB",1);
    323      hc=hilb(cT,1);
    324      if(hl==hc)
    325      {
    326         for(k=1;k<=size(lT);k++)
    327         {
    328            if(lT[k]<cT[k]){lT=cT;break;}
    329            if(lT[k]>cT[k]){break;}
    330         }
    331      }
    332      else
    333      {
    334         if(hc<hl){lT=cT;hl=hilb(lT,1);}
    335      }
    336   }
    337   j=1;
    338   attrib(lT,"isSB",1);
    339   while(j<=size(T))
    340   {
    341      cT=lead(T[j]);
    342      attrib(cT,"isSB",1);
    343      if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
    344      {
    345         T=delete(T,j);
    346         if(j==1) { L=L[2..size(L)]; }
    347         else
    348         {
    349           if (j==size(L)) { L=L[1..size(L)-1]; }
    350           else { L=L[1..j-1],L[j+1..size(L)]; }
    351         }
    352         j--;
    353      }
    354      j++;
    355   }
    356   return(list(T,L,lT));
    357 }
    358 example
    359 { "EXAMPLE:"; echo = 2;
    360    intvec L=2,3,5,7,11;
    361    ring r=0,(y,x),Dp;
    362    ideal I1=y2x,y6;
    363    ideal I2=yx2,y3x,x5,y6;
    364    ideal I3=y2x,x3y,x5,y6;
    365    ideal I4=y2x,x3y,x5;
    366    ideal I5=y2x,yx3,x5,y6;
    367    list T=I1,I2,I3,I4,I5;
    368    list TT=deleteUnluckyPrimes(T,L);
    369    TT;
    370 }
    371 ///////////////////////////////////////////////////////////////////////////////
    372 proc liftPoly(list T, intvec L)
    373 "USAGE:  liftPoly(T,L); T list of polys, L intvec of primes
    374 RETURN:  poly p in Q[x] such that p mod L[i]=T[i]
    375 EXAMPLE: example liftPoly; shows an example
    376 "
    377 {
    378    int i;
    379    list TT;
    380    for(i=size(T);i>0;i--)
    381    { TT[i]=ideal(T[i]); }
    382    T=TT;
    383    ideal hh=chinrem(T,L);
    384    poly h=hh[1];
    385    poly p=lead(h);
    386    poly result;
    387    number n;
    388    bigint N=L[1];
    389    for(i=size(L);i>1;i--)
    390    {
    391       N=N*L[i];
    392    }
    393    while(h!=0)
    394    {
    395      n=Farey(N,bigint(leadcoef(h)));
    396      result=result+n*p;
    397      h=h-lead(h);
    398      p=leadmonom(h);
    399    }
    400    return(result);
    401 }
    402 example
    403 { "EXAMPLE:"; echo = 2;
    404    ring R = 0,(x,y),dp;
    405    intvec L=32003,181,241,499;
    406    list T=ideal(x2+7000x+13000),ideal(x2+100x+147y+40),ideal(x2+120x+191y+10),ideal(x2+x+67y+100);
    407    liftPoly(T,L);
    408 }
    409 ///////////////////////////////////////////////////////////////////////////
    410 proc liftPoly1(list T, intvec L)
    411 "USAGE:  liftPoly1(T,L); T list of polys, L intvec of primes
    412 RETURN:  poly p in Q[x] such that p mod L[i]=T[i]
    413 EXAMPLE: example liftPoly1; shows an example
    414 "
    415 {
    416    poly result;
    417    int i;
    418    poly p;
    419    list TT;
    420    number n;
    421 
    422    bigint N=L[1];
    423    for(i=2;i<=size(L);i++)
    424    {
    425       N=N*L[i];
    426    }
     1179   int i,j;
     1180
     1181   bigint p = 2134567879;
     1182   if(size(#)!=0) { p=#[1]; }
     1183   while(!primeTest(I,p))
     1184   {
     1185      p = prime(random(2000000000,2134567879));
     1186   }
     1187
     1188   def R = basering;
     1189   module F,PrevG,PrevZ,Z2;
     1190   ideal testG,testG1,G1,G2,G3,Gp;
     1191   list L = p;
     1192   list rl = ringlist(R);
     1193   rl[1] = int(p);
     1194
     1195   def S = ring(rl);
     1196   setring S;
     1197   option(redSB);
     1198   module Z,M,Z2;
     1199   ideal I = imap(R,I);
     1200   ideal Gp,G1,G2,G3;
     1201   Gp,Z = liftstd1(I);
     1202   attrib(Gp,"isSB",1);
     1203   module ZZ = syz(I);
     1204   attrib(ZZ,"isSB",1);
     1205   Z = reduce(Z,ZZ);
     1206
     1207   setring R;
     1208   Gp = imap(S,Gp);
     1209   PrevZ = imap(S,Z);
     1210   PrevG = module(Gp);
     1211   F = module(I);
     1212   testG = farey(Gp,p);
     1213   attrib(testG,"isSB",1);
    4271214   while(1)
    4281215   {
    429      p=leadmonom(T[1]);
    430      for(i=2;i<=size(T);i++)
    431      {
    432         if(leadmonom(T[i])>p)
    433         {
    434           p=leadmonom(T[i]);
    435         }
    436      }
    437      if (p==0) {return(result);}
    438      for(i=1;i<=size(T);i++)
    439      {
    440        if(p==leadmonom(T[i]))
    441        {
    442           TT[i]=leadcoef(T[i]);
    443           T[i]=T[i]-lead(T[i]);
    444        }
    445        else
    446        {
    447           TT[i]=0;
    448        }
    449      }
    450      n=chineseR(TT,L,N);
    451      n=Farey(N,bigint(n));
    452      result=result+n*p;
    453    }
    454 }
    455 example
    456 { "EXAMPLE:"; echo = 2;
    457    ring R = 0,(x,y),dp;
    458    intvec L=32003,181,241,499;
    459    list T=x2+7000x+13000,x2+100x+147y+40,x2+120x+191y+10,x2+x+67y+100;
    460    liftPoly1(T,L);
    461 }
    462  ///////////////////////////////////////////////////////////////////////////////
    463 proc fareyIdeal(ideal I,intvec L)
    464 {
    465    poly result,p;
    466    int i,j;
    467    number n;
    468    bigint N=L[1];
    469    for(i=2;i<=size(L);i++)
    470    {
    471       N=N*L[i];
    472    }
    473 
    474    for(i=1;i<=size(I);i++)
    475    {
    476      p=I[i];
    477      result=lead(p);
    478      while(1)
    479      {
    480         if (p==0) {break;}
    481         p=p-lead(p);
    482         n=Farey(N,bigint(leadcoef(p)));
    483         result=result+n*leadmonom(p);
    484      }
    485      I[i]=result;
    486    }
    487    return(I);
    488 }
    489 ///////////////////////////////////////////////////////////////////////////////
    490 proc Farey (bigint P, bigint N)
    491 "USAGE:  Farey (P,N); P, N number;
    492 RETURN:  a rational number a/b such that a/b=N mod P
    493          and |a|,|b|<(P/2)^{1/2}
    494 "
    495 {
    496    if (P<0){P=-P;}
    497    if (N<0){N=N+P;}
    498    bigint A,B,C,D,E;
    499    E=P;
    500    B=1;
    501    while (N!=0)
    502    {
    503         if (2*N^2<P)
    504         {
    505            return(number(N)/number(B));
    506         }
    507         D=E mod N;
    508         C=A-(E div N)*B;
    509         E=N;
    510         N=D;
    511         A=B;
    512         B=C;
    513    }
    514    return(0);
    515 }
    516 example
    517 { "EXAMPLE:"; echo = 2;
    518    ring R = 0,x,dp;
    519    Farey(32003,12345);
    520 }
    521  ///////////////////////////////////////////////////////////////////////////////
    522 proc chineseR(list T,intvec L,number N)
    523 "USAGE:  chineseR(T,L,N);
    524 RETURN: x such that x = T[i] mod L[i], N=product(L[i])
    525 NOTE:   chinese remainder theorem
    526 EXAMPLE:example chineseR; shows an example
    527 "
    528 {
    529    number x;
    530    if(size(L)==1)
    531    {
    532       x=T[1] mod L[1];
    533       return(x);
    534    }
    535    int i;
    536    int n=size(L);
    537    list M;
    538    for(i=1;i<=n;i++)
    539    {
    540       M[i]=N/L[i];
    541    }
    542    list S=eexgcdN(M);
    543    for(i=1;i<=n;i++)
    544    {
    545       x=x+S[i]*M[i]*T[i];
    546    }
    547    x=x mod N;
    548    return(x);
    549 }
    550 example
    551 { "EXAMPLE:"; echo = 2;
    552    ring R = 0,x,dp;
    553    chineseR(list(24,15,7),intvec(2,3,5),30);
    554 }
    555 
    556 ///////////////////////////////////////////////////////////////////////////////
    557 proc pStd(int p,ideal i)
    558 "USAGE:  pStd(p,i);p integer, i ideal;
    559 RETURN:  an ideal G which is the groebner base for i
    560 EXAMPLE: example pStd; shows an example
    561 "
    562 {
    563   def r=basering;
    564   list rl=ringlist(r);
    565   rl[1]=p;
    566   def r1=ring(rl);
    567   setring r1;
    568   option(redSB);
    569   ideal j=fetch(r,i);
    570   ideal GP=groebner(j);
    571   setring r;
    572   ideal G=fetch(r1,GP);
    573   attrib(G,"isSB",1);
    574   matrix Z=transmat(p,i,G);
    575   matrix G1=gstrich1(p,Z,i,G);
    576   ideal g1=G1;
    577   ideal g22=reduce(g1,G);
    578   matrix G22=transpose(matrix(g22));
    579   matrix M=redmat(G,G1,G22);
    580   matrix Z2=-M*Z;
    581   kill r1;
    582   number c=p;
    583   bigint cb=p;
    584   matrix G0=transpose(matrix(G));
    585   G0= MmodN(G0+ (c)* G22,cb^2);
    586   matrix GF=fareyMatrix(G0,cb^2);
    587   Z=MmodN(Z+(c)*Z2,cb^2);
    588   matrix C=transpose(G);
    589   int n=3;
    590   while(GF<>C)
    591   {
    592     C=GF;
    593     G1= gstrich2(c,Z,i,G0,n);
    594     g1=G1;
    595     g22=reduce(g1,G);
    596     G22=transpose(matrix(g22));
    597     M=redmat(G,G1,G22);
    598     Z2=-M*Z;
    599     Z=MmodN(Z+(c^(n-1))*Z2,cb^n);
    600     G0= MmodN(G0+ (c^(n-1))* G22,cb^n);
    601     GF=fareyMatrix(G0,cb^n);
    602     n++;
    603   }
    604   return(ideal(GF));
    605 }
    606 example
    607 { "EXAMPLE:"; echo = 2;
    608    ring r=0,(x,y,z),dp;
    609    ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
    610    ideal J=pStd(32003,I);
    611    J;
    612 }
    613 ///////////////////////////////////////////////////////////////////////////
    614 proc transmat(int p,ideal i,ideal G)
    615 "USAGE:  transmat(p,I,G); p integer, I,G ideal;
    616 RETURN:  the transformationmatrix Z for the ideal i mod p and the groebner base for i mod p
    617 EXAMPLE: example transmat; shows an example
    618 "
    619 {
    620   def r=basering;
    621   int n=nvars(r);
    622   list rl=ringlist(r);
    623   rl[1]=p;
    624   def r1=ring(rl);
    625   setring r1;
    626   ideal i=fetch(r,i);
    627   ideal G=fetch(r,G);
    628   attrib(G,"isSB",1);
    629   ring rhelp=p,x(1..n),dp;
    630   list lhelp=ringlist(rhelp);
    631   list l=lhelp[3];
    632   setring r;
    633   rl[3]=l;
    634   def r2=ring(rl);
    635   setring r2;
    636   ideal i=fetch(r,i);
    637   option(redSB);
    638   ideal j=std(i);
    639   matrix T=lift(i,j);
    640   setring r1;
    641   matrix T=fetch(r2,T);
    642   ideal j=fetch(r2,j);
    643   matrix M=lift(j,G);
    644   matrix Z=transpose(T*M);
    645   setring r;
    646   matrix Z=fetch(r1,Z);
    647   return(Z);
    648 }
    649 example
    650 { "EXAMPLE:"; echo = 2;
    651    ring r=0,(x,y,z),dp;
    652    ideal i=3x3+x2+1,11y5+y3+2,5z4+z2+4;
    653    ideal g=x3-60x2-60, z4-36z2+37, y5+33y3+66;
    654    int p=181;
    655    matrix Z=transmat(p,i,g);
    656    Z;
    657 }
    658 
    659 ///////////////////////////////////////////////////////////////////////////
    660 proc gstrich1(int p, matrix Z, ideal i, ideal gp)
    661 "USAGE:  gstrich1 (p,Z,i,gp); p integer, Z matrix, i,gp ideals;
    662 RETURN:  a matrix G such that (Z*F-GP)/p, where F and GP are the matrices of the ideals i and gp
    663 "
    664 {
    665   matrix F=transpose(matrix(i));
    666   matrix GP=transpose(matrix(gp));
    667   matrix G=(Z*F-GP)/p;
    668   return(G);
    669 }
    670 ///////////////////////////////////////////////////////////////////////////
    671 proc gstrich2(number p, matrix Z, ideal i, ideal gp, int n)
    672 "USAGE:  gstrich2 (p,Z,i,gp,n); p,n integer, Z matrix, i,gp ideals;
    673 RETURN:  a matrix G such that (Z*F-GP)/(p^(n-1)), where F and GP are the matrices of the ideals i and gp
    674 "
    675 {
    676   matrix F=transpose(matrix(i));
    677   matrix GP=transpose(matrix(gp));
    678   matrix G=(Z*F-GP)/(p^(n-1));
    679   return(G);
    680 }
    681 ///////////////////////////////////////////////////////////////////////////
    682 proc redmat(ideal i, matrix h, matrix g)
    683 "USAGE:  redmat(i,h,g); i ideal , h,g matrices;
    684 RETURN:  a matrix M such that i=M*h+g
    685 "
    686 {
    687   matrix c=h-g;
    688   ideal f=transpose(c);
    689   matrix N=lift(i,f);
    690   matrix M=transpose(N);
    691   return(M);
    692 }
    693 ///////////////////////////////////////////////////////////////////////////
    694 proc fareyMatrix(matrix m,bigint N)
    695 "USAGE:  fareyMatrix(m,y); m matrix, y integer;
    696 RETURN:  a matrix k of the matrix m with Farey rational numbers a/b as coefficients
    697 EXAMPLE: example fareyMatrix; shows an example
    698 "
    699 {
    700   ideal I=m;
    701   poly result,p;
    702   int i,j;
    703   number n;
    704   for(i=1;i<=size(I);i++)
    705   {
    706     p=I[i];
    707     result=lead(p);
    708     while(1)
    709     {
    710       if (p==0) {break;}
    711       p=p-lead(p);
    712       n=Farey(N,bigint(leadcoef(p)));
    713       result=result+n*leadmonom(p);
    714     }
    715     I[i]=result;
    716   }
    717   matrix k=transpose(I);
    718   return(k);
    719 }
    720 example
    721 {"EXAMPLE:"; echo = 2;
    722    ring r=0,(x,y,z),dp;
    723    matrix m[3][1]=x3+682794673x2+682794673,z4+204838402z2+819353608,    y5+186216729y3+372433458;
    724    int p=32003;
    725    matrix b=fareyMatrix(m,p^2);
    726    b;
    727 }
    728 ///////////////////////////////////////////////////////////////////////////
    729 proc MmodN(matrix Z,bigint N)
    730 "USAGE:  MmodN(Z,N);Z matrix, N bigint;
    731 RETURN:  the matrix Z mod N
    732 EXAMPLE: example MmodN;
    733 "
    734 {
    735   int i,j,k;
    736   poly m,p;
    737   number c;
    738   for(i=1;i<=nrows(Z);i++)
    739   {
    740     for(j=1;j<=ncols(Z);j++)
    741     {
    742       for(k=1;k<=size(Z[i,j]);k++)
    743       {
    744         m=leadmonom(Z[i,j][k]);
    745         c=bigint(leadcoef(Z[i,j][k])) mod N;
    746         p=p+c*m;
    747       }
    748       Z[i,j]=p;
    749       p=0;
    750     }
    751   }
    752   return(Z);
     1216      i++;
     1217      G1 = ideal(1/(p^i) * sum(F*PrevZ,(-1)*PrevG));
     1218      setring S;
     1219      G1 = imap(R,G1);
     1220      G2 = reduce(G1,Gp);
     1221      G3 = sum(G1,(-1)*G2);
     1222      M = lift(Gp,G3);
     1223      Z2 = (-1)*Z*M;
     1224
     1225      setring R;
     1226      G2 = imap(S,G2);
     1227      Z2 = imap(S,Z2);
     1228      PrevG = sum(PrevG, module(p^i*G2));
     1229      PrevZ = sum(PrevZ, multiply(poly(p^i),Z2));
     1230      testG1 = farey(ideal(PrevG),p^(i+1));
     1231      attrib(testG1,"isSB",1);
     1232      if(size(reduce(testG1,testG)) == 0)
     1233      {
     1234         if(size(reduce(I,testG1)) == 0)       // I is in testG1
     1235         {
     1236            if(pTestSB(I,testG1,L,2))
     1237            {
     1238               G3 = std(testG1);               // testG1 is SB
     1239               if(size(reduce(G3,testG1)) == 0)
     1240               {
     1241                  return(G3);
     1242               }
     1243            }
     1244         }
     1245      }
     1246      testG = testG1;
     1247      attrib(testG,"isSB",1);
     1248   }
    7531249}
    7541250example
    7551251{ "EXAMPLE:"; echo = 2;
    7561252   ring r = 0,(x,y,z),dp;
    757    matrix m[3][1]= x3+10668x2+10668, z4-12801z2+12802, y5-8728y3+14547;
    758    bigint p=32003;
    759    matrix b=MmodN(m,p^2);
    760    b;
    761 }
    762 ///////////////////////////////////////////////////////////////////////////////
     1253   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
     1254   ideal J = modHenselStd(I);
     1255   J;
     1256}
     1257
     1258////////////////////////////////////////////////////////////////////////////////
     1259
     1260static proc sum(list #)
     1261{
     1262   if(typeof(#[1])=="ideal")
     1263   {
     1264      ideal M;
     1265   }
     1266   else
     1267   {
     1268      module M;
     1269   }
     1270
     1271   int i;
     1272   for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; }
     1273   return(M);
     1274}
     1275
     1276////////////////////////////////////////////////////////////////////////////////
     1277
     1278static proc multiply(poly p, list #)
     1279{
     1280   if(typeof(#[1])=="ideal")
     1281   {
     1282      ideal M;
     1283   }
     1284   else
     1285   {
     1286      module M;
     1287   }
     1288
     1289   int i;
     1290   for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; }
     1291   return(M);
     1292}
     1293
     1294
     1295////////////////////////////// further examples ////////////////////////////////
     1296
    7631297/*
    764 ring r=0,(x,y,z),lp;
     1298ring r = 0, (x,y,z), lp;
    7651299poly s1 = 5x3y2z+3y3x2z+7xy2z2;
    7661300poly s2 = 3xy2z2+x5+11y2z2;
    … …  
    7691303ideal i =  s1, s2, s3, s4;
    7701304
    771 ring r=0,(x,y,z),lp;
     1305ring r = 0, (x,y,z), lp;
    7721306poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
    7731307poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
    … …  
    7761310ideal i =  s1, s2, s3, s4;
    7771311
    778 ring r=0,(x,y,z),lp;
     1312ring r = 0, (x,y,z), lp;
    7791313poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
    7801314poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
    7811315poly s3 = 8x3 + 12y3 + xz2 + 3;
    7821316poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
    783 ideal i =  s1, s2, s3, s4;
     1317ideal i = s1, s2, s3, s4;
    7841318
    7851319int n = 6;
    7861320ring r = 0,(x(1..n)),lp;
    7871321ideal i = cyclic(n);
    788 ring s=0,(x(1..n),t),lp;
    789 ideal i=imap(r,i);
    790 i=homog(i,t);
    791 
    792 ring r=0,(x(1..4),s),(dp(4),dp);
    793 poly s1 =1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
     1322ring s = 0, (x(1..n),t), lp;
     1323ideal i = imap(r,i);
     1324i = homog(i,t);
     1325
     1326ring r = 0, (x(1..4),s), (dp(4),dp);
     1327poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
    7941328poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
    7951329poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
    … …  
    7981332ideal i =  s1, s2, s3, s4, s5;
    7991333
    800 ring r=0,(x,y,z),ds;
    801 int a =16;
    802 int b =15;
    803 int c =4;
    804 int t =1;
    805 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;
    806 ideal i= jacob(f);
    807 
    808 ring r=0,(x,y,z),ds;
    809 int a =25;
    810 int b =25;
    811 int c =5;
    812 int t =1;
    813 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;
    814 ideal i= jacob(f),f;
    815 
    816 ring r=0,(x,y,z),ds;
    817 int a=10;
    818 poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
    819 ideal i= jacob(f);
    820 
    821 ring r=0,(x,y,z),ds;
    822 int a =6;
    823 int b =8;
    824 int c =10;
    825 int alpha =5;
    826 int beta= 5;
    827 int t= 1;
    828 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;
    829 ideal i= jacob(f);
    830 
     1334ring r = 0, (x,y,z), ds;
     1335int a = 16;
     1336int b = 15;
     1337int c = 4;
     1338int t = 1;
     1339poly 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;
     1340ideal i = jacob(f);
     1341
     1342ring r = 0, (x,y,z), ds;
     1343int a = 25;
     1344int b = 25;
     1345int c = 5;
     1346int t = 1;
     1347poly 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;
     1348ideal i = jacob(f),f;
     1349
     1350ring r = 0, (x,y,z), ds;
     1351int a = 10;
     1352poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
     1353ideal i = jacob(f);
     1354
     1355ring r = 0, (x,y,z), ds;
     1356int a = 6;
     1357int b = 8;
     1358int c = 10;
     1359int alpha = 5;
     1360int beta = 5;
     1361int t = 1;
     1362poly 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;
     1363ideal i = jacob(f);
    8311364*/
    8321365
    833 /*
    834 ring r=0,(x,y,z),lp;
    835 poly s1 = 5x3y2z+3y3x2z+7xy2z2;
    836 poly s2 = 3xy2z2+x5+11y2z2;
    837 poly s3 = 4xyz+7x3+12y3+1;
    838 poly s4 = 3x3-4y3+yz2;
    839 ideal i =  s1, s2, s3, s4;
    840 
    841 ring r=0,(x,y,z),lp;
    842 poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
    843 poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
    844 poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
    845 poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
    846 ideal i =  s1, s2, s3, s4;
    847 
    848 ring r=0,(x,y,z),lp;
    849 poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
    850 poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
    851 poly s3 = 8x3 + 12y3 + xz2 + 3;
    852 poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
    853 ideal i =  s1, s2, s3, s4;
    854 
    855 int n = 6;
    856 ring r = 0,(x(1..n)),lp;
    857 ideal i = cyclic(n);
    858 ring s=0,(x(1..n),t),lp;
    859 ideal i=imap(r,i);
    860 i=homog(i,t);
    861 
    862 ring r=0,(x(1..4),s),(dp(4),dp);
    863 poly s1 =1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
    864 poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
    865 poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
    866 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);
    867 poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
    868 ideal i =  s1, s2, s3, s4, s5;
    869 
    870 ring r=0,(x,y,z),ds;
    871 int a =16;
    872 int b =15;
    873 int c =4;
    874 int t =1;
    875 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;
    876 ideal i= jacob(f);
    877 
    878 ring r=0,(x,y,z),ds;
    879 int a =25;
    880 int b =25;
    881 int c =5;
    882 int t =1;
    883 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;
    884 ideal i= jacob(f),f;
    885 
    886 ring r=0,(x,y,z),ds;
    887 int a=10;
    888 poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
    889 ideal i= jacob(f);
    890 
    891 ring r=0,(x,y,z),ds;
    892 int a =6;
    893 int b =8;
    894 int c =10;
    895 int alpha =5;
    896 int beta= 5;
    897 int t= 1;
    898 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;
    899 ideal i= jacob(f);
    900 
    901 */
    902 
Note: See TracChangeset for help on using the changeset viewer.