Changeset c10f46 in git for Singular/LIB/modstd.lib


Ignore:
Timestamp:
Aug 3, 2012, 1:06:22 PM (12 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
b0732eb4254fbe76540d72c337f3d7f985adf948
Parents:
b8610796eaca5ba7db8dc75f8bc88d698e077b365abb0e7fdca82a2bad9778e410148c74d7b49a0b
Message:
Merge pull request #162 from steenpass/libs

changes in LIBs from master
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/modstd.lib

    rb86107 rc10f46  
    11////////////////////////////////////////////////////////////////////////////////
    2 version="$Id$";
     2version="$Id: modstd.lib 14375 2011-08-23 09:29:47Z steidel $";
    33category = "Commutative Algebra";
    44info="
    … …  
    88@*        G. Pfister      pfister@mathematik.uni-kl.de
    99@*        H. Schoenemann  hannes@mathematik.uni-kl.de
     10@*        A. Steenpass    steenpass@mathematik.uni-kl.de
    1011@*        S. Steidel      steidel@mathematik.uni-kl.de
    1112
    … …  
    2627LIB "poly.lib";
    2728LIB "ring.lib";
     29LIB "parallel.lib";
    2830
    2931////////////////////////////////////////////////////////////////////////////////
    … …  
    243245      if(printlevel >= 11)
    244246      {
    245          "isIncluded(K,J) takes "+string(timer - t)+" seconds";
     247         "isIncluded(J,K) takes "+string(timer - t)+" seconds";
    246248         "j = "+string(j);
    247249      }
    … …  
    373375////////////////////////////////////////////////////////////////////////////////
    374376
    375 proc primeTest(ideal I, bigint p)
     377proc primeTest(def II, bigint p)
    376378{
     379   if(typeof(II) == "string")
     380   {
     381      execute("ideal I = "+II+";");
     382   }
     383   else
     384   {
     385      ideal I = II;
     386   }
     387
     388   I = simplify(I, 2);   // erase zero generators
     389
    377390   int i,j;
     391   poly f;
     392   number cnt;
    378393   for(i = 1; i <= size(I); i++)
    379394   {
    380       for(j = 1; j <= size(I[i]); j++)
    381       {
    382          if((numerator(leadcoef(I[i][j])) mod p) == 0) { return(0); }
    383          if((denominator(leadcoef(I[i][j])) mod p) == 0) { return(0); }
     395      f = cleardenom(I[i]);
     396      if(f == 0) { return(0); }
     397      cnt = leadcoef(I[i])/leadcoef(f);
     398      if((numerator(cnt) mod p) == 0) { return(0); }
     399      if((denominator(cnt) mod p) == 0) { return(0); }
     400      for(j = size(f); j > 0; j--)
     401      {
     402         if((leadcoef(f[j]) mod p) == 0) { return(0); }
    384403      }
    385404   }
    … …  
    390409
    391410proc primeList(ideal I, int n, list #)
    392 "USAGE:  primeList(I,n); ( resp. primeList(I,n,L); ) I ideal, n integer
    393 RETURN:  the intvec of n greatest primes  <= 2147483647 (resp. n greatest primes
     411"USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
     412         n integer
     413RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
    394414         < L[size(L)] union with L) such that none of these primes divides any
    395415         coefficient occuring in I
     416NOTE:    The number of cores to use can be defined by ncores, default is 1.
    396417EXAMPLE: example primList; shows an example
    397418"
    … …  
    399420   intvec L;
    400421   int i,p;
     422   int ncores = 1;
     423   
     424//-----------------  Initialize optional parameter ncores  ---------------------
     425   if(size(#) > 0)
     426   {
     427      if(size(#) == 1)
     428      {
     429         if(typeof(#[1]) == "int")
     430         {
     431            ncores = #[1];
     432            # = list();
     433         }
     434      }
     435      else
     436      {
     437         ncores = #[2];
     438      }
     439   }
     440   
    401441   if(size(#) == 0)
    402442   {
    … …  
    421461   }
    422462   if(p == 2) { ERROR("no more primes"); }
    423    for(i = 2; i <= n; i++)
    424    {
    425       p = prime(p-1);
    426       while(!primeTest(I,p))
     463   if(ncores == 1)
     464   {
     465      for(i = 2; i <= n; i++)
    427466      {
    428467         p = prime(p-1);
    429          if(p == 2) { ERROR("no more primes"); }
    430       }
    431       L[size(L)+1] = p;
     468         while(!primeTest(I,p))
     469         {
     470            p = prime(p-1);
     471            if(p == 2) { ERROR("no more primes"); }
     472         }
     473         L[size(L)+1] = p;
     474      }
     475   }
     476   else
     477   {
     478      int neededSize = size(L)+n-1;;
     479      list parallelResults;
     480      list arguments;
     481      int neededPrimes = neededSize-size(L);
     482      while(neededPrimes > 0)
     483      {
     484         arguments = list();
     485         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
     486            *ncores; i > 0; i--)
     487         {
     488            p = prime(p-1);
     489            if(p == 2) { ERROR("no more primes"); }
     490            arguments[i] = list("I", p);
     491         }
     492         parallelResults = parallelWaitAll("primeTest", arguments,
     493            list(list(list(ncores))));
     494         for(i = size(arguments); i > 0; i--)
     495         {
     496            if(parallelResults[i])
     497            {
     498               L[size(L)+1] = arguments[i][2];
     499            }
     500         }
     501         neededPrimes = neededSize-size(L);
     502      }
     503      if(size(L) > neededSize)
     504      {
     505         L = L[1..neededSize];
     506      }
    432507   }
    433508   return(L);
    … …  
    470545   ideal I = imap(R,I);
    471546
     547   intvec opt = option(get);
    472548   option(none);
    473549   option(prompt);
    … …  
    504580
    505581   setring R;
     582   option(set, opt);
    506583   return(imap(newR,sI),imap(newR,T));
    507584}
    … …  
    537614         The standard basis computation modulo p does also vary depending on the
    538615         integer variant, namely
    539 @*       - variant = 1/4: std(.,#[1]) resp. groebner,
    540 @*       - variant = 2/5: groebner,
    541 @*       - variant = 3/6: homog. - std(.,#[1]) resp. groebner - dehomog..
     616@*       - variant = 1: std(.,#[1]) resp. groebner,
     617@*       - variant = 2: groebner,
     618@*       - variant = 3: homog. - std(.,#[1]) resp. groebner - dehomog.,
     619@*       - variant = 4: fglm.
    542620EXAMPLE: example modpStd; shows an example
    543621"
    … …  
    552630   option(redSB);
    553631
    554    int t = timer;
    555    if((variant == 1) || (variant == 4))
     632   if(variant == 1)
    556633   {
    557634      if(size(#) > 0)
    … …  
    565642   }
    566643
    567    if((variant == 2) || (variant == 5))
     644   if(variant == 2)
    568645   {
    569646      i = groebner(i);
    570647   }
    571648
    572    if((variant == 3) || (variant == 6))
     649   if(variant == 3)
    573650   {
    574651      list rl = ringlist(@r);
    … …  
    607684      }
    608685
    609       t = timer;
    610686      i = subst(i, homvar, 1);
    611687      i = simplify(i, 34);
    … …  
    615691      i = interred(i);
    616692      kill HomR;
     693   }
     694   
     695   if(variant == 4)
     696   {
     697      def R1 = changeord("dp");
     698      setring R1;
     699      ideal i = fetch(@r,i);
     700      i = std(i);
     701      setring @r;
     702      i = fglm(R1,i);
    617703   }
    618704
    … …  
    646732RETURN:  a standard basis of I if no warning appears;
    647733NOTE:    The procedure computes a standard basis of I (over the rational
    648          numbers) by using  modular methods. If a warning appears then the
    649          result is a standard basis containing I and with high probability
    650          a standard basis of I.
     734         numbers) by using  modular methods.
    651735         By default the procedure computes a standard basis of I for sure, but
    652736         if the optional parameter #[2] = 0, it computes a standard basis of I
    653          with high probability and a warning appears at the end.
     737         with high probability.
    654738         The procedure distinguishes between different variants for the standard
    655739         basis computation in positive characteristic depending on the ordering
    656740         of the basering, the parameter #[2] and if the ideal I is homogeneous.
    657 @*       - variant = 1, if I is homogeneous and exactness = 0,
    658 @*       - variant = 2, if I is not homogeneous, 1-block-ordering and
    659                         exactness = 0,
     741@*       - variant = 1, if I is homogeneous,
     742@*       - variant = 2, if I is not homogeneous, 1-block-ordering,
    660743@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
    661                         > 1 block) and exactness = 0,
    662 @*       - variant = 4, if I is homogeneous and exactness = 1,
    663 @*       - variant = 5, if I is not homogeneous, 1-block-ordering and
    664                         exactness = 1,
    665 @*       - variant = 6, if I is not homogeneous, complicated ordering (lp or
    666                         > 1 block) and exactness = 1.
     744                        > 1 block),
     745@*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
    667746EXAMPLE: example modStd; shows an example
    668747"
    … …  
    752831   if(printlevel >= 10)
    753832   {
    754       "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3);
     833      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
     834       +", exactness = "+string(exactness);
    755835   }
    756836
    … …  
    761841
    762842//--------------------  Initialize the list of primes  -------------------------
    763    intvec L = primeList(I,n2);
     843   int tt = timer;
     844   int rt = rtimer;
     845   intvec L = primeList(I,n2,n1);
     846   if(printlevel >= 10)
     847   {
     848      "CPU-time for primeList: "+string(timer-tt)+" seconds.";
     849      "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
     850   }
    764851   L[5] = prime(random(an,en));
    765852
    … …  
    768855   int h = homog(I);
    769856
    770    int tt = timer;
    771    int rt = rtimer;
     857   tt = timer;
     858   rt = rtimer;
    772859
    773860   if(!mixedTest())
    … …  
    775862      if(h)
    776863      {
    777          if(exactness == 0)
    778          {
    779             variant = 1;
    780             if(printlevel >= 10) { "variant = 1"; }
    781          }
    782          if(exactness == 1)
    783          {
    784             variant = 4;
    785             if(printlevel >= 10) { "variant = 4"; }
    786          }
     864         variant = 1;
     865         if(printlevel >= 10) { "variant = 1"; }
     866         
    787867         rl[1] = L[5];
    788868         def @r = ring(rl);
    … …  
    802882         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
    803883         {
    804             if(exactness == 0)
    805             {
    806                variant = 2;
    807                if(printlevel >= 10) { "variant = 2"; }
    808             }
    809             if(exactness == 1)
    810             {
    811                variant = 5;
    812                if(printlevel >= 10) { "variant = 5"; }
    813             }
     884            variant = 2;
     885            if(printlevel >= 10) { "variant = 2"; }
    814886         }
    815887         else
    … …  
    827899            if((order == "simple") || (size(rl) > 4))
    828900            {
    829                if(exactness == 0)
     901               variant = 2;
     902               if(printlevel >= 10) { "variant = 2"; }
     903            }
     904            else
     905            {
     906               rl[1] = L[5];
     907               def @r = ring(rl);
     908               setring @r;
     909               
     910               def @s = changeord("dp");
     911               setring @s;
     912               ideal I = std(fetch(R0,I));
     913               if(dim(I) == 0)
    830914               {
    831                   variant = 2;
    832                   if(printlevel >= 10) { "variant = 2"; }
     915                  variant = 4;
     916                  if(printlevel >= 10) { "variant = 4"; }
    833917               }
    834                if(exactness == 1)
    835                {
    836                   variant = 5;
    837                   if(printlevel >= 10) { "variant = 5"; }
    838                }
    839             }
    840             else
    841             {
    842                if(exactness == 0)
     918               else
    843919               {
    844920                  variant = 3;
    845921                  if(printlevel >= 10) { "variant = 3"; }
     922                 
     923                  int nvar@r = nvars(@r);
     924                  intvec w;
     925                  for(i = 1; i <= nvar@r; i++)
     926                  {
     927                     w[i] = deg(var(i));
     928                  }
     929                  w[nvar@r + 1] = 1;
     930
     931                  list hiRi = hilbRing(fetch(R0,I),w);
     932                  intvec W = hiRi[2];
     933                  @s = hiRi[1];
     934                  setring @s;
     935
     936                  Id(1) = std(Id(1));
     937                  intvec hi = hilb(Id(1), 1, W);
    846938               }
    847                if(exactness == 1)
    848                {
    849                   variant = 6;
    850                   if(printlevel >= 10) { "variant = 6"; }
    851                }
    852 
    853                rl[1] = L[5];
    854                def @r = ring(rl);
    855                setring @r;
    856                int nvar@r = nvars(@r);
    857                intvec w;
    858                for(i = 1; i <= nvar@r; i++)
    859                {
    860                   w[i] = deg(var(i));
    861                }
    862                w[nvar@r + 1] = 1;
    863 
    864                list hiRi = hilbRing(fetch(R0,I),w);
    865                intvec W = hiRi[2];
    866                def @s = hiRi[1];
    867                setring @s;
    868 
    869                Id(1) = std(Id(1));
    870                intvec hi = hilb(Id(1), 1, W);
    871 
     939               
    872940               setring R0;
    873941               kill @r,@s;
    … …  
    9711039         link l(i) = "ssi:fork";
    9721040         open(l(i));
    973          if((variant == 1) || (variant == 3) ||
    974             (variant == 4) || (variant == 6))
     1041         if((variant == 1) || (variant == 3))
    9751042         {
    9761043            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
    9771044                                                  eval(variant), eval(hi))));
    9781045         }
    979          if((variant == 2) || (variant == 5))
     1046         if((variant == 2) || (variant == 4))
    9801047         {
    9811048            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
    … …  
    9851052
    9861053      int t = timer;
    987       if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
     1054      if((variant == 1) || (variant == 3))
    9881055      {
    9891056         P = modpStd(I_for_fork, L[1], variant, hi);
    9901057      }
    991       if((variant == 2) || (variant == 5))
     1058      if((variant == 2) || (variant == 4))
    9921059      {
    9931060         P = modpStd(I_for_fork, L[1], variant);
    … …  
    10051072//--------------  Main standard basis computations in positive  ----------------
    10061073//----------------------  characteristic start here  ---------------------------
     1074
     1075   list arguments_farey, results_farey;
    10071076
    10081077   while(1)
    … …  
    10291098                  if(j <= size(L))
    10301099                  {
    1031                      if((variant == 1) || (variant == 3) ||
    1032                         (variant == 4) || (variant == 6))
     1100                     if((variant == 1) || (variant == 3))
    10331101                     {
    10341102                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
    … …  
    10361104                        j++;
    10371105                     }
    1038                      if((variant == 2) || (variant == 5))
     1106                     if((variant == 2) || (variant == 4))
    10391107                     {
    10401108                        write(l(i), quote(modpStd(I_for_fork,
    … …  
    10621130         while(j <= size(L))
    10631131         {
    1064             if((variant == 1) || (variant == 3) ||
    1065                (variant == 4) || (variant == 6))
     1132            if((variant == 1) || (variant == 3))
    10661133            {
    10671134               P = modpStd(I, L[j], variant, hi);
    10681135            }
    1069             if((variant == 2) || (variant == 5))
     1136            if((variant == 2) || (variant == 4))
    10701137            {
    10711138               P = modpStd(I, L[j], variant);
    … …  
    10951162//-------------------  Lift results to basering via farey  ---------------------
    10961163
    1097       tt = timer;
     1164      tt = timer; rt = rtimer;
    10981165      N = T2[1];
    10991166      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
    11001167      H = chinrem(T1,T2);
    1101       J = farey(H,N);
    1102       if(printlevel >= 10) { "Lifting-process takes "+string(timer - tt)
    1103                              +" seconds"; }
     1168      if(n1 == 1)
     1169      {
     1170         J = farey(H,N);
     1171      }
     1172      else
     1173      {
     1174         for(i = ncols(H); i > 0; i--)
     1175         {
     1176            arguments_farey[i] = list(ideal(H[i]), N);
     1177         }
     1178         results_farey = parallelWaitAll("farey", arguments_farey,
     1179                                         list(list(list(n1))));
     1180         for(i = ncols(H); i > 0; i--)
     1181         {
     1182            J[i] = results_farey[i][1];
     1183         }
     1184      }
     1185      if(printlevel >= 10)
     1186      {
     1187         "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
     1188         "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
     1189      }
    11041190
    11051191//----------------  Test if we already have a standard basis of I --------------
    11061192
    11071193      tt = timer; rt = rtimer;
    1108       if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
     1194      if((variant == 1) || (variant == 3))
    11091195      {
    11101196         pTest = pTestSB(I,J,L,variant,hi);
    11111197      }
    1112       if((variant == 2) || (variant == 5))
     1198      if((variant == 2) || (variant == 4))
    11131199      {
    11141200         pTest = pTestSB(I,J,L,variant);
    … …  
    11321218
    11331219         attrib(J,"isSB",1);
    1134          tt = timer; rt = rtimer;
    1135          sizeTest = 1 - isIncluded(I,J,n1);
    1136 
    1137          if(printlevel >= 10)
    1138          {
    1139             "CPU-time for checking if I subset <G> is "
    1140             +string(timer - tt)+" seconds.";
    1141             "Real-time for checking if I subset <G> is "
    1142             +string(rtimer - rt)+" seconds.";
    1143          }
    1144 
    1145          if(sizeTest == 0)
    1146          {
    1147             if((variant == 1) || (variant == 2) || (variant == 3))
    1148             {
    1149                "===========================================================";
    1150                "WARNING: Output might not be a standard basis of the input.";
    1151                "===========================================================";
    1152                option(set, opt);
    1153                if(n1 > 1) { kill I_for_fork; }
    1154                return(J);
    1155             }
    1156             if((variant == 4) || (variant == 5) || (variant == 6))
     1220         
     1221         if(exactness == 0)
     1222         {
     1223            option(set, opt);
     1224            if(n1 > 1) { kill I_for_fork; }
     1225            return(J);
     1226         }
     1227         
     1228         if(exactness == 1)
     1229         {
     1230            tt = timer; rt = rtimer;
     1231            sizeTest = 1 - isIncluded(I,J,n1);
     1232
     1233            if(printlevel >= 10)
     1234            {
     1235               "CPU-time for checking if I subset <G> is "
     1236               +string(timer - tt)+" seconds.";
     1237               "Real-time for checking if I subset <G> is "
     1238               +string(rtimer - rt)+" seconds.";
     1239            }
     1240
     1241            if(sizeTest == 0)
    11571242            {
    11581243               tt = timer; rt = rtimer;
    … …  
    11851270
    11861271      j = size(L) + 1;
    1187       L = primeList(I,n3,L);
     1272      tt = timer; rt = rtimer;
     1273      L = primeList(I,n3,L,n1);
     1274      if(printlevel >= 10)
     1275      {
     1276         "CPU-time for primeList: "+string(timer-tt)+" seconds.";
     1277         "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
     1278      }
    11881279
    11891280      if(n1 > 1)
    … …  
    11921283         {
    11931284            open(l(i));
    1194             if((variant == 1) || (variant == 3) ||
    1195                (variant == 4) || (variant == 6))
     1285            if((variant == 1) || (variant == 3))
    11961286            {
    11971287               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
    11981288                                                     eval(variant), eval(hi))));
    11991289            }
    1200             if((variant == 2) || (variant == 5))
     1290            if((variant == 2) || (variant == 4))
    12011291            {
    12021292               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
    … …  
    13181408   def S = ring(rl);
    13191409   setring S;
     1410   intvec opt = option(get);
    13201411   option(redSB);
    13211412   module Z,M,Z2;
    … …  
    13621453               if(size(reduce(G3,testG1)) == 0)
    13631454               {
     1455                  option(set, opt);
    13641456                  return(G3);
    13651457               }
Note: See TracChangeset for help on using the changeset viewer.