- Timestamp:
- Jul 25, 2013, 12:49:24 PM (11 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- 3631a1348b10ecae02a13616277e1ba0458f46e7
- Parents:
- 87977c77f3d3666b9c0306a8fd7e7e93a76b25e609830b6a4f7a61fb96ffe9bd1517c21777ec1096
- Location:
- Singular/LIB
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/gitfan.lib
r09830b6 r2ea186 27 27 28 28 LIB "parallel.lib"; // for parallelWaitAll 29 LIB "general.lib"; 29 30 30 31 //////////////////////////////////////////////////// 31 proc mod_init() 32 { 33 LIB"gfanlib.so"; 34 } 35 36 static proc int2face(int n, int r) 32 33 proc int2face(int n, int r) 37 34 { 38 35 int k = r-1; … … 43 40 while(2^k > n0){ 44 41 k--; 45 //v[size(v)+1] = 0;46 42 } 47 43 … … 56 52 ///////////////////////////////// 57 53 58 proc isAface(ideal a, intvec gam0 )54 proc isAface(ideal a, intvec gam0, int n) 59 55 "USAGE: isAface(a,gam0); a: ideal, gam0:intvec 60 56 PURPOSE: Checks whether the face of the positive orthant, … … 66 62 " 67 63 { 68 poly pz;69 70 64 // special case: gam0 is the zero-cone: 71 65 if (size(gam0) == 1 and gam0[1] == 0){ 72 ideal G;66 poly pz; ideal G; 73 67 74 68 // is an a-face if and only if RL0(0,...,0) = const … … 93 87 } 94 88 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; 97 92 string initNewRing = "ring Rgam0 = 0,("; 98 93 99 94 for (int i=1; i<size(gam0); i++){ 100 95 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 } 159 example 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 175 proc 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])) + ","; } 103 180 initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;"; 104 181 def R = basering; … … 107 184 ideal agam0 = imap(R,a); 108 185 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; 116 194 ideal G = std(agam0); 117 195 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 } 140 201 //////////////////////////////////////////////////// 141 202 … … 145 206 list AF; 146 207 147 for(i = start; i <= end; i ++){208 for(i = start; i <= end; i=i+1){ 148 209 if(i < 2^r){ 210 string(i)+" "+string(size(AF)); 149 211 gam0 = int2face(i,r); 150 212 … … 165 227 166 228 proc afaces(ideal a, list #) 167 "USAGE: afaces(a, b, c); a: ideal, d: int, c: int229 "USAGE: afaces(a, b, c); a: ideal, b: int, c: int 168 230 PURPOSE: Returns a list of all a-faces (represented by intvecs). 169 231 Moreover, it is possible to specify a dimensional bound b, … … 210 272 211 273 // 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){ 213 275 "another one needed"; 214 276 gam0 = int2face(i,r); … … 223 285 } 224 286 287 int l; 225 288 // 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]; 228 291 } 229 292 -
Singular/LIB/primdecint.lib
r87977c7 r2ea186 17 17 18 18 PROCEDURES: 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 20 21 minAssZ(I); compute the minimal associated primes of I 21 22 radicalZ(I); compute the radical of I … … 224 225 J=simplify(J,2); 225 226 attrib(J,"isSB",1); 226 IS= maxIndependSet(J);227 IS=Primdec::maxIndependSet(J); 227 228 setring R; 228 229 //=== computing the pseudo primary and extract it … … 809 810 J=simplify(J,2); 810 811 attrib(J,"isSB",1); 811 IS= maxIndependSet(J);812 IS=Primdec::maxIndependSet(J); 812 813 setring R; 813 814 //=== computing the pseudo primary and extract it … … 1099 1100 //=== this is our way to obtain the coefficients in Z[u] of the 1100 1101 //=== 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]); 1102 1103 execute(quotring); 1103 1104 ideal I=imap(Shelp,I); … … 1184 1185 } 1185 1186 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 later1193 //=== if it is no more static1194 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 later1222 //=== if it is no more static1223 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 else1249 {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 else1263 {1264 resu[n]=varstr(basering),ordstr(basering),0;1265 }1266 }1267 return(resu);1268 1187 } 1269 1188 … … 1402 1321 return(0); 1403 1322 } 1323 1324 //////////////////////////////////////////////////////////////////////////////// 1325 1326 static 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 1355 static 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 1384 static 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 1463 proc primdecZM(module N) 1464 "USAGE: primdecZM(N); N module 1465 RETURN: 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 1470 EXAMPLE: 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 } 1501 example 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 1404 1508 1405 1509 //////////////////////////////////////////////////////////////////////////////// … … 1602 1706 p^3*P^2-1, 1603 1707 F^2*f^3-1; 1708 1709 ring R=integer,(x,y),(c,lp); 1710 module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18x]; 1711 1712 ring R=integer,(x,y),(c,lp); 1713 module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18]; 1714 1715 ring R=integer,(x,y),(c,lp); 1716 module N=[-y,7,0],[2y3-y2],[3x,y2],[2y-y2,x],[4,5x3]; 1717 1718 ring r=integer,(x,y),(c,lp); 1719 module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,xy-x],[x,0,-xy],[5x,0,0]; 1720 1721 ring R2=integer,(a(1),a(2),a(3),b(1),b(2),b(3)),(c,lp); 1722 module 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 1724 ring R3=integer,(x,y,z),(c,lp); 1725 module N=[y2+z2,xy,xz],[xy,x2+z2,yz],[xz,yz,x2+y2]; 1726 1727 ring R4=integer,(x,y,z,a,b,c),(c,lp); 1728 module N=[x3y2z2c,x2y3z2c,x2y2z3c],[x3y2z2b,x2y3z2b,x2y2z3b],[x3y2z2a,x2y3z2a,x2y2z3a]; 1729 1604 1730 */
Note: See TracChangeset
for help on using the changeset viewer.