Changeset 09830b6 in git


Ignore:
Timestamp:
Jul 25, 2013, 12:17:25 PM (11 years ago)
Author:
Gerhard Pfister <pfister@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '4188d308699580d975efd0f6cca8dcb41c396f70')
Children:
2ea18609e14482a4ddbe0fd0d68b8d5885268a0c
Parents:
0ab2b929de2eb4a16c8561981bbcef044441f401
Message:
Procedures for primariy decomposition of modules added.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/primdecint.lib

    r0ab2b92 r09830b6  
    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.