Changeset 2ea186 in git for Singular


Ignore:
Timestamp:
Jul 25, 2013, 12:49:24 PM (11 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
3631a1348b10ecae02a13616277e1ba0458f46e7
Parents:
87977c77f3d3666b9c0306a8fd7e7e93a76b25e609830b6a4f7a61fb96ffe9bd1517c21777ec1096
Message:
Merge pull request #345 from gerhardpfister/spielwiese

Procedures for primariy decomposition of modules added.
Location:
Singular/LIB
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/gitfan.lib

    r09830b6 r2ea186  
    2727
    2828LIB "parallel.lib"; // for parallelWaitAll
     29LIB "general.lib";
    2930
    3031////////////////////////////////////////////////////
    31 proc mod_init()
    32 {
    33   LIB"gfanlib.so";
    34 }
    35 
    36 static proc int2face(int n, int r)
     32
     33proc int2face(int n, int r)
    3734{
    3835  int k = r-1;
     
    4340    while(2^k > n0){
    4441      k--;
    45       //v[size(v)+1] = 0;
    4642    }
    4743
     
    5652/////////////////////////////////
    5753
    58 proc isAface(ideal a, intvec gam0)
     54proc isAface(ideal a, intvec gam0, int n)
    5955"USAGE:  isAface(a,gam0); a: ideal, gam0:intvec
    6056PURPOSE: Checks whether the face of the positive orthant,
     
    6662"
    6763{
    68   poly pz;
    69 
    7064  // special case: gam0 is the zero-cone:
    7165  if (size(gam0) == 1 and gam0[1] == 0){
    72     ideal G;
     66    poly pz; ideal G;
    7367
    7468    // is an a-face if and only if RL0(0,...,0) = const
     
    9387  }
    9488
    95 
    96   // ring is too big: Switch to KK[T_i; e_i\in gam0] instead:
     89  // ring is too big: Switch to KK[T_i | e_i\in gam0] instead:
     90  intvec w=ringlist(basering)[3][1][2];
     91  intvec w0;
    9792  string initNewRing = "ring Rgam0 = 0,(";
    9893
    9994  for (int i=1; i<size(gam0); i++){
    10095    initNewRing = initNewRing + string(var(gam0[i])) + ",";
    101   }
    102 
     96    w0[i]=w[gam0[i]];
     97  }
     98
     99  w0 = w0,w[gam0[i]],1;
     100  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + ",U),Wp("+string(w0)+");";
     101  def R = basering;
     102  execute(initNewRing);
     103
     104  ideal agam0 = imap(R,a);
     105
     106  for (i = 1; i<=size(agam0); i=i+1)
     107  {
     108    if (size(agam0[i]) == 1)
     109    { return(0,0); }
     110  }
     111
     112  poly p = var(1); // first entry of g; p = prod T_i with i element of g
     113  for ( i = 2; i <= nvars(basering); i++ ) {
     114      p = p * var(i);
     115    }
     116  // p is now the product over all T_i, with e_i in gam0
     117
     118  agam0 = agam0, p - 1; // rad-membership
     119  agam0 = timeStd(agam0,5);
     120  // "timestd finished after: "+string(timer-t);
     121  // int timeout = 0;
     122  if (attrib(agam0,"isSB") < 1)
     123  {
     124    kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
     125    return(0,1);
     126    // // "timestd failed in "+string(gam0)+", falling back to saturation!";
     127    // setring R; kill Rgam0;
     128
     129    // initNewRing = "ring Rgam0 = 0,(";
     130
     131    // w0 = 0;
     132    // for (i=1; i<size(gam0); i++){
     133    //   initNewRing = initNewRing + string(var(gam0[i])) + ",";
     134    //   w0[i]=w[gam0[i]];
     135    // }
     136
     137    // w0 = w0,w[gam0[i]];
     138    // initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),Wp("+string(w0)+");";
     139    // execute(initNewRing);
     140
     141    // ideal G = imap(R,a);
     142
     143    // timeout = 1;
     144    // int t=rtimer;
     145    // for(int k=nvars(basering); k>0; k--) { G=sat(G,var(k))[1]; }
     146    // t = (rtimer - t) div 1000;
     147    // string(n)+": saturation successful after "+string(t)+" with: "+string(G<>1);
     148  }
     149
     150  // does G contain 1?, i.e. is G = 1?
     151  if(agam0 <> 1) {
     152    kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
     153    return(1,0); // true
     154  }
     155
     156  kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R;
     157  return(0,0); // false
     158}
     159example
     160{
     161  echo = 2;
     162
     163  ring R = 0,(T(1..4)),dp;
     164  ideal I = T(1)*T(2)-T(4);
     165
     166  intvec w = 1,4;
     167  intvec v = 1,2,4;
     168
     169  isAface(I,w); // should be 0
     170  "-----------";
     171  isAface(I,v); // should be 1
     172}
     173
     174
     175proc isAfaceNonZero(ideal a, intvec gam0)
     176{
     177  string initNewRing = "ring Rgam0 = 0,(";
     178  for (int i=1; i<size(gam0); i++)
     179    { initNewRing = initNewRing + string(var(gam0[i])) + ","; }
    103180  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;";
    104181  def R = basering;
     
    107184  ideal agam0 = imap(R,a);
    108185
    109   poly p = var(1); // first entry of g; p = prod T_i with i element of g
    110   for ( i = 2; i <= nvars(basering); i++ ) {
    111       p = p * var(i);
    112     }
    113   // p is now the product over all T_i, with e_i in gam0
    114 
    115   agam0 = agam0, p - 1; // rad-membership
     186  for ( i = 1; i<=size(agam0); i=i+1)
     187    { if (size(agam0[i]) == 1) { return(0); } }
     188
     189  poly p = var(1);
     190  for ( i = 2; i <= nvars(basering); i++ )
     191    { p = p * var(i); }
     192
     193  agam0 = agam0, p - 1;
    116194  ideal G = std(agam0);
    117195
    118   // does G contain 1?, i.e. is G = 1?
    119   if(G <> 1) {
    120     return(1); // true
    121   }
    122 
    123   return(0); // false
    124 }
    125 example
    126 {
    127   echo = 2;
    128 
    129   ring R = 0,(T(1..4)),dp;
    130   ideal I = T(1)*T(2)-T(4);
    131 
    132   intvec w = 1,4;
    133   intvec v = 1,2,4;
    134 
    135   isAface(I,w); // should be 0
    136   "-----------";
    137   isAface(I,v); // should be 1
    138 }
    139 
     196  if(G <> 1)
     197    { return(1); }
     198
     199  return(0);
     200}
    140201////////////////////////////////////////////////////
    141202
     
    145206  list AF;
    146207
    147   for(i = start; i <= end; i++){
     208  for(i = start; i <= end; i=i+1){
    148209    if(i < 2^r){
     210      string(i)+"    "+string(size(AF));
    149211      gam0 = int2face(i,r);
    150212
     
    165227
    166228proc afaces(ideal a, list #)
    167 "USAGE:  afaces(a, b, c); a: ideal, d: int, c: int
     229"USAGE:  afaces(a, b, c); a: ideal, b: int, c: int
    168230PURPOSE: Returns a list of all a-faces (represented by intvecs).
    169231         Moreover, it is possible to specify a dimensional bound b,
     
    210272
    211273  // do remaining ones:
    212   for(i = ncores * step +1; i < 2^r; i++){
     274  for(i = ncores * step +1; i < 2^r; i=i+1){
    213275    "another one needed";
    214276    gam0 = int2face(i,r);
     
    223285  }
    224286
     287  int l;
    225288  // read out afaces of out into AF:
    226   for(i = 1; i <= size(out); i++){
    227     AF = AF + out[i];
     289  for(l = 1; l <= size(out); l++){
     290    AF = AF + out[l];
    228291  }
    229292
  • Singular/LIB/primdecint.lib

    r87977c7 r2ea186  
    1717
    1818PROCEDURES:
    19  primdecZ(I);       compute the primary decomposition of I
     19 primdecZ(I);       compute the primary decomposition of ideal I
     20 primdecZM(I);      compute the primary decomposition of module I
    2021 minAssZ(I);        compute the minimal associated primes of I
    2122 radicalZ(I);       compute the radical of I
     
    224225               J=simplify(J,2);
    225226               attrib(J,"isSB",1);
    226                IS=maxIndependSet(J);
     227               IS=Primdec::maxIndependSet(J);
    227228               setring R;
    228229               //=== computing the pseudo primary and extract it
     
    809810                  J=simplify(J,2);
    810811                  attrib(J,"isSB",1);
    811                   IS=maxIndependSet(J);
     812                  IS=Primdec::maxIndependSet(J);
    812813                  setring R;
    813814                  //=== computing the pseudo primary and extract it
     
    10991100         //=== this is our way to obtain the coefficients in Z[u] of the
    11001101         //=== leading terms of the Groebner basis above
    1101          string quotring=prepareQuotientring(nvars(basering)-L[1][3]);
     1102         string quotring=Primdec::prepareQuotientring(nvars(basering)-L[1][3]);
    11021103         execute(quotring);
    11031104         ideal I=imap(Shelp,I);
     
    11841185   }
    11851186   return(J);
    1186 }
    1187 
    1188 ////////////////////////////////////////////////////////////////////////////////
    1189 
    1190 static proc prepareQuotientring (int nnp)
    1191 {
    1192 //=== this is from primdec.lib, it is static there, should be imported later
    1193 //=== if it is no more static
    1194   ideal @ih,@jh;
    1195   int npar=npars(basering);
    1196   int @n;
    1197 
    1198   string quotring= "ring quring = ("+charstr(basering);
    1199   for(@n=nnp+1;@n<=nvars(basering);@n++)
    1200   {
    1201      quotring=quotring+",var("+string(@n)+")";
    1202      @ih=@ih+var(@n);
    1203   }
    1204 
    1205   quotring=quotring+"),(var(1)";
    1206   @jh=@jh+var(1);
    1207   for(@n=2;@n<=nnp;@n++)
    1208   {
    1209      quotring=quotring+",var("+string(@n)+")";
    1210      @jh=@jh+var(@n);
    1211   }
    1212   quotring=quotring+"),(C,lp);";
    1213 
    1214   return(quotring);
    1215 }
    1216 
    1217 ////////////////////////////////////////////////////////////////////////////////
    1218 
    1219 static proc maxIndependSet (ideal j)
    1220 {
    1221 //=== this is from primdec.lib, it is static there, should be imported later
    1222 //=== if it is no more static
    1223   int n,k,di;
    1224 
    1225   list resu,hilf;
    1226   if(size(j)==0)
    1227   {
    1228      resu[1]=varstr(basering);
    1229      resu[2]=ordstr(basering);
    1230      resu[3]=0;
    1231      return(list(resu));
    1232   }
    1233   string var1,var2;
    1234   list v=indepSet(j,0);
    1235 
    1236   for(n=1;n<=size(v);n++)
    1237   {
    1238     di=0;
    1239     var1="";
    1240     var2="";
    1241     for(k=1;k<=size(v[n]);k++)
    1242     {
    1243       if(v[n][k]!=0)
    1244       {
    1245         di++;
    1246         var2=var2+"var("+string(k)+"),";
    1247       }
    1248       else
    1249       {
    1250         var1=var1+"var("+string(k)+"),";
    1251       }
    1252     }
    1253     if(di>0)
    1254     {
    1255       var1=var1+var2;
    1256       var1=var1[1..size(var1)-1];
    1257       hilf[1]=var1;
    1258       hilf[2]="lp";
    1259       hilf[3]=di;
    1260       resu[n]=hilf;
    1261     }
    1262     else
    1263     {
    1264       resu[n]=varstr(basering),ordstr(basering),0;
    1265     }
    1266   }
    1267   return(resu);
    12681187}
    12691188
     
    14021321   return(0);
    14031322}
     1323
     1324////////////////////////////////////////////////////////////////////////////////
     1325
     1326static proc pseudo_primdecZM(module N)
     1327{
     1328   ideal I=quotient(N,freemodule(nrows(N)));
     1329   if(size(I)==0){return(list(list(N,I)));}
     1330   
     1331   list B=minAssZ(I);
     1332   list S,R,L;
     1333   ideal K;
     1334   if(size(B)==0){return(S);}
     1335   for(int i=1;i<=size(B);i++)
     1336   {
     1337      S[i]=separatorsZ(i,B);
     1338   }
     1339   for(i=1;i<=size(B);i++)
     1340   {
     1341      L=sat(N,S[i]);
     1342      K[i]=S[i]^L[2];
     1343      R[i]=list(L[1],B[i]);
     1344   }
     1345   L=pseudo_primdecZM(N+K*freemodule(nrows(N)));
     1346   for(i=1;i<=size(L);i++)
     1347   {
     1348      R[size(R)+1]=L[i];
     1349   }
     1350   return(R);
     1351}
     1352
     1353
     1354
     1355static proc prepare_extractZM(list L)
     1356{
     1357   def R=basering;
     1358   module N=L[1];
     1359   ideal I=quotient(N,freemodule(nrows(N)));
     1360   list B=primdecZ(I);
     1361   list M;
     1362   if(size(B)==1){return(M);}
     1363   I=std(I);
     1364   list rl=ringlist(R);
     1365   if(deg(I[1])==0)
     1366   {
     1367      execute("int p="+string(I[1])+";");
     1368      rl[1]=p;
     1369   }
     1370   else
     1371   {
     1372      rl[1]=0;
     1373   }
     1374   def Shelp =ring(rl);
     1375   setring Shelp;
     1376   ideal I=imap(R,I);
     1377   I=std(I);
     1378   M=Primdec::maxIndependSet(I);
     1379   setring R;
     1380   return(M);
     1381}
     1382
     1383
     1384static proc extractZM(list M, list L)
     1385{
     1386   //=== M is a list of a pseudo primary module and the corresponding prime
     1387   //=== L is a list of maximal independent sets for P
     1388   def R=basering;
     1389   ideal P=M[2];
     1390   module I=M[1];
     1391   poly h=1;
     1392
     1393   //=== size(L)=0 means P is maximal ideal and I is primary
     1394   if(size(L)>0)
     1395   {
     1396      if(L[1][3]!=0)
     1397      {
     1398         //=== if u in x is an independent set of L then we compute a Groebner
     1399         //=== Basis in Z[u][x-u]
     1400         execute("ring S=integer,("+L[1][1]+"),lp;");
     1401         module I=imap(R,I);
     1402         I=std(I);
     1403         list rl=ringlist(S);
     1404         rl[1]=0;
     1405         def Shelp =ring(rl);
     1406         setring Shelp;
     1407         module I=imap(S,I);
     1408         //=== this is our way to obtain the coefficients in Z[u] of the
     1409         //=== leading terms of the Groebner basis above
     1410         string quotring=Primdec::prepareQuotientring(nvars(basering)-L[1][3]);
     1411         execute(quotring);
     1412         module I=imap(Shelp,I);
     1413         list C;
     1414         int i;
     1415         for(i=1;i<=size(I);i++)
     1416         {
     1417            C[i]=leadcoef(I[i]);
     1418         }
     1419         setring Shelp;
     1420         list C=imap(quring,C);
     1421         setring R;
     1422         list C=imap(Shelp,C);
     1423      }
     1424      else
     1425      {
     1426         // this is the case that P=<p>, p prime
     1427         I=std(I);
     1428         ideal IC=simplify(flatten(lead(I)),2);
     1429         list C;
     1430         int i;
     1431         for(i=1;i<=size(IC);i++)
     1432         {
     1433            C[i]=I[i];
     1434         }
     1435         list rl=ringlist(R);
     1436         rl[1]=0;
     1437         def Shelp =ring(rl);
     1438      }
     1439      for(i=1;i<=size(C);i++)
     1440      {
     1441         if(deg(C[i])>0){h=h*C[i];}  // das muss noch besser gemacht werden,
     1442                                     // nicht ausmultiplizieren!
     1443      }
     1444      setring Shelp;
     1445      poly h=imap(R,h);
     1446      ideal fac=factorize(h,1);
     1447      setring R;
     1448      list II;
     1449      h=1;
     1450      ideal fac=imap(Shelp,fac);
     1451      for(i=1;i<=size(fac);i++)
     1452      {
     1453         II=sat(I,fac[i]);
     1454          I=II[1];
     1455          h=h*fac[i]^II[2];
     1456      }
     1457   }
     1458   I=std(I);
     1459   return(list(I,h));
     1460}
     1461
     1462
     1463proc primdecZM(module N)
     1464"USAGE:  primdecZM(N); N module
     1465RETURN:  a list pr of primary modules and their associated primes:
     1466@format
     1467   pr[i][1]   the i-th primary component,
     1468   pr[i][2]   the i-th prime component.
     1469@end format
     1470EXAMPLE: example primdecZM; shows an example
     1471"
     1472{
     1473  list P,K,S;
     1474  int i,j;
     1475  list L=pseudo_primdecZM(N);
     1476  list M,O;
     1477  for(i=1;i<=size(L);i++)
     1478  {
     1479     if(size(L[i][2])!=0)
     1480     {
     1481        M=prepare_extractZM(L[i]);
     1482        O=extractZM(L[i],M);
     1483        P[size(P)+1]=list(O[1],L[i][2]);
     1484        K[size(K)+1]=L[i][1]+O[2]*freemodule(nrows(L[i][1]));
     1485     }
     1486     else
     1487     {
     1488        P[size(P)+1]=L[i];
     1489     }
     1490  }
     1491  for(j=1;j<=size(K);j++)
     1492  {
     1493     S=primdecZM(K[j]);
     1494     for(i=1;i<=size(S);i++)
     1495     {
     1496        P[size(P)+1]=S[i];
     1497     }
     1498  }
     1499  return(P);
     1500}
     1501example
     1502{ "EXAMPLE:";  echo = 2;
     1503   ring R=integer,(x,y),(c,lp);
     1504   module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18x];
     1505   primdecZM(N);
     1506}
     1507
    14041508
    14051509////////////////////////////////////////////////////////////////////////////////
     
    16021706        p^3*P^2-1,
    16031707        F^2*f^3-1;
     1708
     1709ring R=integer,(x,y),(c,lp);
     1710module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18x];
     1711
     1712ring R=integer,(x,y),(c,lp);
     1713module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18];
     1714
     1715ring R=integer,(x,y),(c,lp);
     1716module N=[-y,7,0],[2y3-y2],[3x,y2],[2y-y2,x],[4,5x3];
     1717
     1718ring r=integer,(x,y),(c,lp);
     1719module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,xy-x],[x,0,-xy],[5x,0,0];
     1720
     1721ring R2=integer,(a(1),a(2),a(3),b(1),b(2),b(3)),(c,lp);
     1722module N=[a(1)*b(1),a(2)*b(1),a(3)*b(1)],[a(1)*b(2),a(2)*b(2),a(3)*b(2)],[a(1)*b(3),a(2)*b(3),a(3)*b(3)];
     1723
     1724ring R3=integer,(x,y,z),(c,lp);
     1725module N=[y2+z2,xy,xz],[xy,x2+z2,yz],[xz,yz,x2+y2];
     1726
     1727ring R4=integer,(x,y,z,a,b,c),(c,lp);
     1728module N=[x3y2z2c,x2y3z2c,x2y2z3c],[x3y2z2b,x2y3z2b,x2y2z3b],[x3y2z2a,x2y3z2a,x2y2z3a];
     1729
    16041730*/
Note: See TracChangeset for help on using the changeset viewer.