Changeset 38301d in git
- Timestamp:
- Jul 11, 2015, 10:21:42 AM (8 years ago)
- Branches:
- (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
- Children:
- 47362dcbb6b4fbb3b9d6cded63720e333b4c97b3
- Parents:
- d52203afadb18a3c05d4c08a2bc18bc327690cd9e9478bb30a6f3a2d77328b3805644a035ac3e934
- Files:
-
- 6 added
- 1 deleted
- 42 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/crypto.lib
rd52203 r38301d 69 69 merkle_hellman_transformation(list knapsack, int key, int mod1 generates a hard knapsack for the Merkle-Hellman Kryptosystem for a given easy knapsack , a multiplicator key and a modulus mod1 70 70 merkle_hellman_encryption(list knapsack, list message) encrypts a message with the Merkle-Hellman Kryptosystem, using a hard knapsack and a message encoded as binary list 71 merkle_hellman_decryption(list knapsack, int key, int mod1, int message) decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the hard knapsack, the key, the modulus mod1 and the message encoded as integer 71 merkle_hellman_decryption(list knapsack, int key, int mod1, int message) decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the hard knapsack, the key, the modulus mod1 and the message encoded as integer 72 72 super_increasing_knapsack(int ksize) Creates the smallest super-increasing knapsack of given size ksize 73 73 h_increasing_knapsack(int ksize, int h) Creates the smallest h-increasing knapsack of given size ksize and h … … 1272 1272 if((k mod 2)==0) 1273 1273 { 1274 resu=ellipticMult(N,a,b,P,k /2);1274 resu=ellipticMult(N,a,b,P,k div 2); 1275 1275 return(ellipticAdd(N,a,b,resu,resu)); 1276 1276 } -
Singular/LIB/gradedModules.lib
rd52203 r38301d 31 31 grsum(M,N) direct sum of two graded modules coker(M) + coker(N) 32 32 grpower(M,p) direct p-th power of graded module coker(M) 33 grtranspose(M) un-ordered graded transpose of map M 33 grtranspose(M) un-ordered graded transpose of map M 34 34 grgens(M) try to compute submodule generators of coker(M) 35 35 grpres(F) presentation of submodule generated by columns of F … … 51 51 mappingcone(M,N) mapping cone? 52 52 grlifting3(A,B) RND! chain lifting? probably wrong one 53 mappingcone3(A,B) mapping cone3? 53 mappingcone3(A,B) mapping cone3? 54 54 grrange(M) get the row-weightings 55 55 grneg(A) graded object given by -A 56 56 matrixpres(a) matrix presentation of direct sum of Omega^{a[i]}(i) 57 57 "; 58 58 59 59 // grisequal(A,B) check whether A is exactly eqal to B? TODO: isomorphic! 60 60 … … 86 86 87 87 v = -v; // NOTE: due to Mathematical meanings of Singular data 88 88 89 89 90 90 int lst = v[1]; 91 int cnt = 1; 91 int cnt = 1; 92 92 93 93 string p = R; … … 95 95 96 96 int k, d; 97 for (k = 2; k <= n; k++ ) 97 for (k = 2; k <= n; k++ ) 98 98 { 99 99 d = v[k]; 100 100 if( d == lst ) { cnt = cnt + 1; } 101 else 101 else 102 102 { 103 103 if (cnt > 1){ p = p + "^" + string(cnt); } … … 120 120 def E = grtwist(2, 0); 121 121 def v = grrange(E); // grdeg(E); 122 grsumstr("", v ); 122 grsumstr("", v ); 123 123 } 124 124 … … 180 180 { 181 181 int n = size(N); ASSUME(0, n > 0); 182 182 183 183 string msg1 = ""; 184 184 if( size(R) >= 2 ) … … 186 186 msg1 = msg1 + "(let R:="+R+")"; 187 187 R = "R"; // !!! 188 } 188 } 189 189 msg1 = msg1 + ": " ; 190 190 191 191 192 192 193 193 list arr; arr[n] = list(); 194 194 int exact = (1==1); 195 195 196 196 int i = 1; 197 197 198 198 ASSUME(1, grtest(N[i])); 199 199 200 200 string dst = grsumstr(R, grrange(N[i])); 201 201 string src = grsumstr(R, grdeg(N[i])); 202 202 203 203 arr[i] = list(dst, src); 204 204 205 205 i = i + 1; 206 206 207 207 while( i <= n ) 208 208 { 209 209 ASSUME(1, grtest(N[i])); 210 210 211 211 dst = grsumstr(R, grrange(N[i])); 212 212 213 213 if( exact && (src != dst) ) 214 214 { … … 218 218 219 219 src = grsumstr(R, grdeg(N[i])); 220 220 221 221 arr[i] = list(dst, src); 222 222 … … 236 236 msg = msg + newline + arr[i][1] + " <-- "+o+"_" + string(i) + " --"; 237 237 }; 238 238 239 239 msg = msg + newline + arr[n][2]; 240 msg = msg + ", given by maps: "; 240 msg = msg + ", given by maps: "; 241 241 } else 242 242 { 243 243 // print(arr); 244 244 245 msg = msg + "-object collection"; 245 msg = msg + "-object collection"; 246 246 o = "o"; 247 247 248 248 // for( i = 1; i <= n; i++ ) 249 249 // { 250 250 // msg = msg + newline + arr[i][1] + " <-- "+o+"_" + string(i) + " -- " + arr[i][2]; 251 // }; 252 msg = msg + ", given by the following maps (named here as "+o+"_[1 .. "+string(n)+"]): "; 251 // }; 252 msg = msg + ", given by the following maps (named here as "+o+"_[1 .. "+string(n)+"]): "; 253 253 } 254 254 … … 257 257 258 258 for( i = 1; i <= size(N); i++ ) 259 { 259 { 260 260 print( o+"_" + string(i) + " :" ); 261 grview( N[i] ); 261 grview( N[i] ); 262 262 }; 263 263 … … 268 268 269 269 ASSUME(1, grtest(N) ); 270 270 271 271 msg = msg + " homomorphism"; 272 272 if( size(R) >= 2 ) … … 280 280 intvec gr = grrange(N); // grading weights? 281 281 string dst = grsumstr(R, gr); 282 282 283 283 intvec G = grdeg(N); 284 284 string src = grsumstr(R, G); 285 285 286 286 if( ncols(N) == 0 ) 287 287 { 288 288 src = "0"; 289 } 289 } 290 290 291 291 lst = msg; 292 292 293 if( (size(lst) + size(dst) + size(src) + 4) > 80 ) 294 { 293 if( (size(lst) + size(dst) + size(src) + 4) > 80 ) 294 { 295 295 if( (size(lst) + size(dst)) > 80 ) { msg = msg + newline; lst = ""; } 296 296 … … 312 312 // lst = lst + ", given by "; 313 313 314 314 315 315 int nc = ncols(N); int nr = nrows(N); 316 316 … … 323 323 { 324 324 ASSUME(0, nc > 0); 325 325 326 326 matrix M = module(N); 327 327 328 328 int r,c; 329 329 int d = 1; // number of extra cols/rows for extra info around the central degree(N) block in D … … 361 361 } 362 362 363 } else 364 { 363 } else 364 { 365 365 msg = msg + "a matrix"; 366 366 } 367 367 368 print(msg + ", with degrees: " ); 369 draw(D, d); // print it nicely! 368 print(msg + ", with degrees: " ); 369 draw(D, d); // print it nicely! 370 370 } 371 371 } … … 462 462 " 463 463 { 464 ASSUME(1, grtest(M) ); 465 464 ASSUME(1, grtest(M) ); 465 466 466 if ( typeof(attrib(M, "degHomog")) == "intvec" ) 467 467 { … … 469 469 470 470 if( size(t) == 0 ){ return (t); } // ZERO! 471 471 472 472 ASSUME(2, ncols(M) == size(t) ); 473 473 return (t); … … 629 629 "Input degrees: "; grview(I); 630 630 631 def RR = grres(I, 0, 1); list L = RR; 631 def RR = grres(I, 0, 1); list L = RR; 632 632 633 633 " = Non-minimal betti numbers: "; print(betti(L, 0), "betti"); … … 830 830 831 831 ring r=32003,(x,y,z),dp; 832 832 833 833 grview( grtwists ( intvec(-4, 1, 6 )) ); 834 834 … … 843 843 " 844 844 { 845 ASSUME(0, a > 0); 845 ASSUME(0, a > 0); 846 846 def Z = grtwists( intvec(d:a) ); // will set the rank as well 847 847 // ASSUME(2, grisequal(Z, grpower( grshift(grzero(), d), a ) )); // optional check … … 931 931 932 932 def SS = grobj(S, c, dc); 933 933 934 934 ASSUME(0, size( grrange(SS) ) == (size(a) + size(b)) ); 935 935 ASSUME(0, size( grdeg(SS) ) == (size(da) + size(db)) ); 936 936 ASSUME(0, ncols( SS ) == size( grdeg(SS) ) ); 937 937 ASSUME(0, nrows( SS ) == size( grrange(SS) ) ); 938 938 939 939 return(SS); 940 940 } … … 984 984 { 985 985 ASSUME(1, grtest(M) ); 986 987 986 987 988 988 989 989 intvec a = grrange(M); … … 993 993 { 994 994 "!! Warning: shifting '0 <- 0' leaves it as it unchanged!"; 995 return (M); 996 } 997 995 return (M); 996 } 997 998 998 a = a - intvec(d:size(a)); 999 999 t = t - intvec(d:size(t)); … … 1056 1056 { 1057 1057 ASSUME(0, size(w) >= nrows(A) ); 1058 1058 1059 1059 module M = module(A); 1060 1060 … … 1063 1063 1064 1064 intvec @ww = 0:0; 1065 1065 1066 1066 if( size(#) > 0 ) 1067 1067 { 1068 1068 ASSUME(0, typeof(#[1]) == "intvec" ); 1069 1069 1070 1070 @ww = intvec( #[1] ); 1071 1071 … … 1077 1077 } 1078 1078 } 1079 1079 1080 1080 ASSUME(0, size(@ww) == ncols(M) ); 1081 1081 } … … 1091 1091 M = freemodule(0); 1092 1092 } 1093 1093 1094 1094 attrib( M, "rank", size(w) ); 1095 1095 attrib( M, "isHomog", w ); 1096 1096 1097 1097 // ASSUME(0, /* PROBLEM WITH ZERO COLUMNS / THEIR DEGREES! */ (ncols(M) == 0) ); 1098 1098 } 1099 1099 } 1100 1100 1101 1101 // type(@ww); type(M); 1102 1102 1103 1103 ASSUME(0, size(@ww) == ncols(M) ); // ?! 1104 1104 1105 1105 attrib(M, "degHomog", @ww); 1106 1106 … … 1132 1132 grview( grobj( freemodule(0), intvec(1,2,3) ) ); 1133 1133 1134 matrix z1[0][3]; grview( grobj( z1, 0:0, intvec(1,2,3) ) ); 1134 matrix z1[0][3]; grview( grobj( z1, 0:0, intvec(1,2,3) ) ); 1135 1135 grview( grobj( freemodule(0), 0:0, intvec(1,2,3) ) ); 1136 1136 1137 1137 matrix z0[0][0]; grview( grobj( z0, 0:0 ) ); 1138 1138 grview( grobj( freemodule(0), 0:0 ) ); 1139 1140 1139 1140 1141 1141 1142 1142 } … … 1145 1145 "USAGE: grtest(M[,b]), anyting M, optionally int b 1146 1146 RETURN: 1 if M is a valid graded object, 0 otherwise 1147 PURPOSE: validate a graded object. Print an invalid object message if b is not given 1147 PURPOSE: validate a graded object. Print an invalid object message if b is not given 1148 1148 NOTE: M should be an ideal or module or matrix, with weighting attribute 1149 1149 'isHomog' and optionally total graded degrees attribute 'degHomog'. … … 1169 1169 if ( nrows(N) != size(gr) ) 1170 1170 { 1171 if(b) { " ? grtest: Input has wrong number of rows!"; }; 1171 if(b) { " ? grtest: Input has wrong number of rows!"; }; 1172 1172 return (0); 1173 1173 }; … … 1177 1177 return(1); 1178 1178 } 1179 1179 1180 1180 // if( attrib(N, "rank") != size(gr) ){ return (0); } // wrong rank :( 1181 1181 … … 1189 1189 return (1); 1190 1190 } 1191 1191 1192 1192 if ( ncols(N) != size(T) ) 1193 1193 { 1194 if(b) { " ? grtest: Input has wrong number of cols!"; }; 1194 if(b) { " ? grtest: Input has wrong number of cols!"; }; 1195 1195 return (0); 1196 1196 }; 1197 1197 1198 1198 int k = nvars(basering) + 1; // index of mod. column in the leadexp 1199 1199 … … 1259 1259 { 1260 1260 ASSUME(0, d >= 0 ); 1261 1261 1262 1262 if( d == 0 ) { return (A); } 1263 1263 1264 1264 if( ncols(A) == 0 ) 1265 1265 { 1266 1266 matrix B[nrows(A) + d][0]; 1267 return (B); 1268 } 1269 1267 return (B); 1268 } 1269 1270 1270 module T; T[d] = 0; 1271 1271 T = T, module(transpose(A)); … … 1281 1281 matrix m[0][2]; 1282 1282 type( align(m, 3) ); 1283 } 1283 } 1284 1284 1285 1285 … … 1302 1302 module A = grobj( module([x+y, x, 0, 0], [0, x+y, y, 0]), intvec(0,0,0,1) ); 1303 1303 grview(A); 1304 1304 1305 1305 module B = grgroebner(A); 1306 1306 grview(B); … … 1324 1324 ring r=32003,(x,y,z),dp; 1325 1325 1326 module A = grgroebner( grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ) ); 1326 module A = grgroebner( grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ) ); 1327 1327 grview(A); 1328 1328 1329 1329 grview(grsyz(A)); 1330 1331 module X = grgroebner( grobj( module([x]), intvec(2) ) ); 1330 1331 module X = grgroebner( grobj( module([x]), intvec(2) ) ); 1332 1332 grview(X); 1333 1333 1334 1334 // syzygy module should be zero! 1335 1335 grview(grsyz(X)); 1336 1337 1336 1337 1338 1338 } 1339 1339 … … 1351 1351 intvec a = grdeg(A); 1352 1352 intvec b = grrange(B); 1353 1353 1354 1354 ASSUME(0, (size(a) == size(b)) && (a == b)); // == for intvec :( 1355 1355 … … 1361 1361 ring r=32003,(x,y,z),dp; 1362 1362 1363 module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 1363 module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 1364 1364 grview(A); 1365 1365 1366 1366 A = grgroebner(A); 1367 1367 grview(A); 1368 1368 1369 1369 module B = grsyz(A); 1370 1370 grview(B); 1371 1371 print(B); 1372 1372 1373 1373 module D = grprod( A, B ); 1374 1374 grview(D); … … 1384 1384 "USAGE: grres(M, l[, b]), graded object M, int l, int b 1385 1385 RETURN: graded resolution = list of graded objects 1386 PURPOSE: compute graded resolution of M (of length l) and minimise it if b was given 1386 PURPOSE: compute graded resolution of M (of length l) and minimise it if b was given 1387 1387 EXAMPLE: example grres; shows an example 1388 1388 " … … 1392 1392 1393 1393 intvec v = grrange(A); 1394 1394 1395 1395 int b = (size(#) > 0); 1396 1396 if(b) { list r = res(A, l, #[1]); } else { list r = res(A, l); } 1397 1397 1398 1398 l = size(r); 1399 1399 1400 1400 int i; 1401 1401 1402 1402 for ( i = 1; i <= l; i++ ) 1403 1403 { … … 1411 1411 r[i] = grobj(r[i], v); v = grdeg(r[i]); 1412 1412 } 1413 1413 1414 1414 i = i-1; 1415 1415 … … 1421 1421 ring r=32003,(x,y,z),dp; 1422 1422 1423 module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 1423 module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 1424 1424 grview(A); 1425 1425 1426 1426 module B = grgroebner(A); 1427 1427 grview(B); … … 1438 1438 RETURN: graded object 1439 1439 PURPOSE: graded transpose of M 1440 NOTE: no reordering is performend by this procedure 1440 NOTE: no reordering is performend by this procedure 1441 1441 EXAMPLE: example grtranspose; shows an example 1442 1442 " … … 1457 1457 1458 1458 module K = grsyz( L ); grview(K); 1459 1459 1460 1460 1461 1461 // Corner cases: 0 <- 0! … … 1463 1463 grview( grtranspose( Z ) ); 1464 1464 1465 1465 1466 1466 // Corner cases: * <- 0 1467 1467 matrix M1[3][0]; 1468 1468 1469 1469 module Z1 = grobj( M1, intvec(-1, 0, 1) ); grview(Z1); 1470 1470 grview( grtranspose( Z1 ) ); 1471 1472 1471 1472 1473 1473 // Corner cases: 0 <- * 1474 1474 matrix M2[0][3]; … … 1476 1476 module Z2 = grobj( M2, 0:0, intvec(-1, 0, 1) ); grview(Z2); 1477 1477 grview( grtranspose( Z2 ) ); 1478 1478 1479 1479 } 1480 1480 … … 1482 1482 proc grgens(def M) 1483 1483 "USAGE: grgens(M), graded object M (map) 1484 RETURN: graded object 1484 RETURN: graded object 1485 1485 PURPOSE: try compute graded generators of coker(M) and return them as columns 1486 1486 of a graded map. … … 1492 1492 1493 1493 module N = grtranspose( grsyz( grtranspose(M) ) ); 1494 1494 1495 1495 // ASSUME(3, grisequal( grgroebner(M), grgroebner( grpres( N ) ) ) ); // FIXME: not always true!? 1496 1496 1497 1497 return ( N ); 1498 1498 } … … 1505 1505 1506 1506 module N = grgens(M); 1507 1507 1508 1508 grview( N ); print(N); // fine == M 1509 1509 1510 1510 1511 module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 1511 module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 1512 1512 1513 1513 A = grgroebner(A); grview(A); … … 1518 1518 1519 1519 grview( grgens( grzero() ) ); 1520 1520 1521 1521 } 1522 1522 … … 1534 1534 1535 1535 // ASSUME(3, grisequal( M, grgens( N ) ) ); 1536 1536 1537 1537 return ( N ); 1538 1538 } … … 1545 1545 grview(A); 1546 1546 1547 "graded transpose: "; def B = grtranspose(A); grview( B ); print(B); 1547 "graded transpose: "; def B = grtranspose(A); grview( B ); print(B); 1548 1548 1549 1549 "... syzygy: "; def C = grsyz(B); grview(C); … … 1553 1553 "... and back to presentation: "; def E = grsyz( D ); grview(E); print(E); 1554 1554 1555 def F = grgens( E ); grview(F); print(F); 1555 def F = grgens( E ); grview(F); print(F); 1556 1556 def G = grpres( F ); grview(G); print(G); 1557 1557 1558 1558 1559 1559 def M = grtwists( intvec(-2, 0, 4, 4) ); grview(M); 1560 1560 1561 1561 def N = grgens(M); grview( N ); print(N); 1562 1563 def L = grpres( N ); grview( L ); print(L); 1562 1563 def L = grpres( N ); grview( L ); print(L); 1564 1564 } 1565 1565 … … 1652 1652 module N; 1653 1653 1654 if(i>n) 1654 if(i>n) 1655 1655 { // no middle part 1656 1656 if(a[1]>0) … … 1659 1659 1660 1660 if(a[n+1]>0) 1661 { N=grsum(N,grtwist(a[n+1],-1));} 1662 } 1661 { N=grsum(N,grtwist(a[n+1],-1));} 1662 } 1663 1663 else 1664 1664 { N=grtwist(a[n+1],-1);} 1665 1666 return (N); // grorder(N)); 1665 1666 return (N); // grorder(N)); 1667 1667 } 1668 else // i <= n: middle part is present, a_i != 0 1668 else // i <= n: middle part is present, a_i != 0 1669 1669 { // a = a1 ... | i:2, a_2 ..... i: n, a_n | .... i: n+1a_(n+1) 1670 1670 j = i - 1; … … 1697 1697 1698 1698 1699 return ((N)); // return (grorder(N)); 1699 return ((N)); // return (grorder(N)); 1700 1700 } 1701 1701 } … … 1704 1704 ring r = 32003,(x(0..4)),dp; 1705 1705 1706 def N1 = KeneshlouMatrixPresentation(intvec(2,0,0,0,0)); 1706 def N1 = KeneshlouMatrixPresentation(intvec(2,0,0,0,0)); 1707 1707 grview(N1); 1708 1708 … … 1738 1738 ASSUME(1, grtest(B)); 1739 1739 ASSUME(0, grrange(A)==grrange(B)); 1740 1740 1741 1741 intvec v = grrange(A); 1742 1742 intvec w=grdeg(A),grdeg(B); … … 1746 1746 { "EXAMPLE:"; echo = 2; 1747 1747 ring r; 1748 1748 1749 1749 module R=grobj(module([x,y,z]),intvec(0:3)); 1750 1750 grview(R); … … 1771 1771 // matrix U; 1772 1772 matrix L =lift(A,B/*,U*/); // module(B*U) = module(matrix(A)*L) 1773 1773 1774 1774 return(grobj(L, grdeg(A), grdeg(B))); 1775 1775 } … … 1782 1782 1783 1783 1784 module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 1784 module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 1785 1785 grview(D); 1786 1786 … … 1788 1788 grview(G); 1789 1789 1790 ASSUME(0, grisequal( grprod(D, G), P) ); 1790 ASSUME(0, grisequal( grprod(D, G), P) ); 1791 1791 } 1792 1792 … … 1807 1807 1808 1808 module Z = grobj(freemodule(0),intvec(0:0),intvec(0:0)); 1809 1809 1810 1810 grrange(Z); 1811 1811 grdeg(Z); 1812 1812 1813 1813 grview(Z); 1814 1814 … … 1838 1838 grtranspose( grprod( N, alpha1 ) ) ) 1839 1839 ) ); // alpha0! 1840 1840 1841 1841 } 1842 1842 example … … 1899 1899 1900 1900 ASSUME(0, t >= 2); 1901 1901 1902 1902 list P; 1903 1903 1904 1904 "t: ", t; 1905 1905 1906 1906 P[1]= grrndmap( rM[1], rN[1] ); // alpha1 1907 1907 1908 if(t==2){return(P[1]);} 1909 1908 if(t==2){return(P[1]);} 1909 1910 1910 for(k=2; k<=t; k++) 1911 1911 { 1912 P[k] = grlift( grprod(P[k-1],rM[k]), rN[k] ); 1913 grview(P[k]); 1914 1915 } 1916 1912 P[k] = grlift( grprod(P[k-1],rM[k]), rN[k] ); 1913 grview(P[k]); 1914 1915 } 1916 1917 1917 return(P); 1918 1918 1919 1919 } 1920 1920 example … … 1923 1923 ring r=32003,(x,y,z),dp; 1924 1924 1925 module P=grobj(module([xy,0,xz]),intvec(0,1,0)); 1925 module P=grobj(module([xy,0,xz]),intvec(0,1,0)); 1926 1926 grview(P); 1927 1927 1928 module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 1928 module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 1929 1929 grview(D); 1930 1930 … … 1937 1937 module D=grobj(module([y,0,z],[x2+y2,z,0], [z3, xy, xy2]),intvec(0,1,0)); 1938 1938 D = grgroebner(D); 1939 grview( grres(D, 0)); 1939 grview( grres(D, 0)); 1940 1940 1941 1941 def G=grlifting(D, D); … … 1951 1951 def M = KeneshlouMatrixPresentation(intvec(0,0,1,0)); 1952 1952 grview( grres(M, 0) ); 1953 1953 1954 1954 // module N = grshift(kos[3], 1); // psi, Syz_2(K(1)) 1955 1955 def N = KeneshlouMatrixPresentation(intvec(0,1,0,0)); … … 1966 1966 proc mappingcone(M,N) 1967 1967 "USAGE: mappingcone(M,N), M,N graded objects 1968 RETURN: chain complex (as a list) 1968 RETURN: chain complex (as a list) 1969 1969 PURPOSE: construct a free resolution of the cokernel of a random map between Img(M), and Img(N). 1970 1970 EXAMPLE: example mappingcone; shows an example … … 2026 2026 2027 2027 // correct 2028 proc grrndmap(def S, def D, list #) 2028 proc grrndmap(def S, def D, list #) 2029 2029 "USAGE: grrndmap(S,D), graded objects S and D 2030 2030 RETURN: graded object … … 2092 2092 // 0<---M<----F0<----F1<----F2<----F3<---- 2093 2093 // |p1 |p2 2094 // 2094 // 2095 2095 // 0<---N<----G0<----G1<----G2<----G3<---- 2096 2096 // B(g1) g2 g3 2097 // 2097 // 2098 2098 proc grlifting2(A,B) 2099 2099 "USAGE: grlifting2(A,B), graded objects A and B (matrices defining maps) 2100 RETURN: map of chain complexes (as a list) 2100 RETURN: map of chain complexes (as a list) 2101 2101 PURPOSE: construct a map of chain complexes between free resolution of 2102 2102 M=coker(A) and N=coker(B). … … 2127 2127 P[1]=grrndmap2(A,B); 2128 2128 2129 // A(or B)=0 2129 // A(or B)=0 2130 2130 if(t==1){return(P[1])}; 2131 2131 2132 2132 for(k=2;k<=t;k++) 2133 2133 { … … 2161 2161 def T=grlifting2(DD,PP); T; 2162 2162 2163 // def Z=grlifting2(P,D); Z; // WRONG!!! 2164 2163 // def Z=grlifting2(P,D); Z; // WRONG!!! 2164 2165 2165 } 2166 2166 … … 2171 2171 proc mappingcone2(A,B) 2172 2172 "USAGE: mappingcone2(A,B), graded objects A and B (matrices defining maps) 2173 RETURN: chain complex (as a list) 2173 RETURN: chain complex (as a list) 2174 2174 PURPOSE: construct the free resolution of a cokernel of a random map between M=coker(A), and N=coker(B) 2175 2175 EXAMPLE: example mappingcone2; … … 2234 2234 */ 2235 2235 2236 proc grlifting3(A,B) 2236 proc grlifting3(A,B) 2237 2237 "TODO: grlifting4 was newer and had more documentation than this proc, but was removed... Please verify and update! 2238 2238 " … … 2247 2247 list rN = grres(B,0,1); 2248 2248 print( betti(rN), "betti"); 2249 2249 2250 2250 int i,j,k; 2251 2251 … … 2260 2260 } 2261 2261 int t=min(i,j); 2262 2262 2263 2263 list P; 2264 2264 2265 2265 "t: ", t; 2266 2266 // grview(rM[t]); grview(rN[t]); 2267 2267 2268 2268 P[t]= grrndmap2(rM[t],rN[t]); 2269 2269 grview(P[t]); … … 2312 2312 //def I=KeneshlouMatrixPresentation(intvec(2,3,0,6,2)); 2313 2313 //def J=KeneshlouMatrixPresentation(intvec(4,0,1,2,1)); 2314 //def N=grlifting3(I,J); grview(N); 2314 //def N=grlifting3(I,J); grview(N); 2315 2315 } 2316 2316 … … 2318 2318 "USAGE: grneg(A), graded object A 2319 2319 RETURN: graded object 2320 PURPOSE: graded map defined by -A 2320 PURPOSE: graded map defined by -A 2321 2321 EXAMPLE: example grneg; shows an example 2322 2322 " … … 2341 2341 proc mappingcone3(A,B) 2342 2342 "USAGE: mappingcone3(A,B), graded objects A and B (matrices defining maps) 2343 RETURN: chain complex (as a list) 2343 RETURN: chain complex (as a list) 2344 2344 PURPOSE: construct a free resolution of the cokernel of a random map between M=coker(A), and N=coker(B) 2345 2345 EXAMPLE: example mappingcone3; shows an example … … 2376 2376 2377 2377 T[i]=grtranspose(D); 2378 2378 2379 2379 kill A, B, C, D, v, w, r, s, zero; 2380 2380 } … … 2394 2394 def F=grlifting3(A,T); grview(F); 2395 2395 2396 // BUG in the proc 2396 // BUG in the proc 2397 2397 def G=mappingcone3(A,T); grview(G); 2398 2398 … … 2403 2403 ideal P=groebner(flatten(U[2])); 2404 2404 resolution L=mres(P,0); 2405 print(betti(L),"betti"); 2405 print(betti(L),"betti"); 2406 2406 */ 2407 2407 … … 2421 2421 def I=KeneshlouMatrixPresentation(intvec(2,3,0,6,2)); 2422 2422 def J=KeneshlouMatrixPresentation(intvec(4,0,1,2,1)); 2423 // def N=grlifting3(I,J); 2423 // def N=grlifting3(I,J); 2424 2424 // 2nd module does not lie in the first: 2425 2425 // def NN=mappingcone3(I,J); // ???????? … … 2458 2458 module N; 2459 2459 2460 if(i>n) 2460 if(i>n) 2461 2461 { // no middle part 2462 2462 if(a[1]>0) … … 2465 2465 2466 2466 if(a[n+1]>0) 2467 { N=grsum(N,grtwist(a[n+1],0));} 2468 } 2467 { N=grsum(N,grtwist(a[n+1],0));} 2468 } 2469 2469 else 2470 2470 { N=grtwist(a[n+1],0);} 2471 2472 return (N); // grorder(N)); 2471 2472 return (N); // grorder(N)); 2473 2473 } 2474 2474 2475 else // i <= n: middle part is present, a_i != 0 2475 else // i <= n: middle part is present, a_i != 0 2476 2476 { // a = a1 ... | i:2, a_2 ..... i: n, a_n | .... i: n+1a_(n+1) 2477 2477 module I = maxideal(1); … … 2505 2505 2506 2506 2507 return ((N)); // return (grorder(N)); 2507 return ((N)); // return (grorder(N)); 2508 2508 } 2509 2509 } … … 2518 2518 grview(S); 2519 2519 2520 def N1 = matrixpres(intvec(2,0,0,0,0)); 2520 def N1 = matrixpres(intvec(2,0,0,0,0)); 2521 2521 grview(N1); 2522 2522 -
Singular/LIB/grwalk.lib
re9478b r38301d 255 255 } 256 256 257 proc gwalk(ideal Go, list #)258 "SYNTAX: gwalk(ideal i );259 gwalk(ideal i, int vec v, intvec w);257 proc gwalk(ideal Go, int reduction,int printout, list #) 258 "SYNTAX: gwalk(ideal i, int reduction, int printout); 259 gwalk(ideal i, int reduction, int printout, intvec v, intvec w); 260 260 TYPE: ideal 261 261 PURPOSE: compute the standard basis of the ideal, calculated via … … 274 274 string ord_str = OSCTW[2]; 275 275 intvec curr_weight = OSCTW[3]; /* original weight vector */ 276 intvec target_weight = OSCTW[4]; /* t erget weight vector */276 intvec target_weight = OSCTW[4]; /* target weight vector */ 277 277 kill OSCTW; 278 278 option(redSB); … … 284 284 //print("//** help ring = " + string(basering)); 285 285 ideal G = fetch(xR, Go); 286 G = system("Mwalk", G, curr_weight, target_weight,basering );286 G = system("Mwalk", G, curr_weight, target_weight,basering,reduction,printout); 287 287 288 288 setring xR; … … 299 299 //** compute a Groebner basis of I w.r.t. lp. 300 300 ring r = 32003,(z,y,x), lp; 301 ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z; 302 gwalk(I); 301 ideal I = zy2+yx2+yx+3, 302 z3x+y3+zyx-yx2-yx-3, 303 z2yx3-y5+z2yx2+y3x2+y2x3+y3x+y2x2+3z2x+3y2+3yx, 304 zyx5+y6-y4x2-y3x3+2zyx4-y4x-y3x2+zyx3-3z2yx+3zx3-3y3-3y2x+3zx2, 305 yx7-y7+y5x2+y4x3+3yx6+y5x+y4x2+3yx5-6zyx3+yx4+3x5+3y4+3y3x-6zyx2+6x4+3x3-9zx; 306 gwalk(I,0,1); 303 307 } 304 308 … … 346 350 } 347 351 348 proc fwalk(ideal Go, list #)349 "SYNTAX: fwalk(ideal i );350 fwalk(ideal i, int vec v, intvec w);352 proc fwalk(ideal Go, int reduction, int printout, list #) 353 "SYNTAX: fwalk(ideal i,int reductioin); 354 fwalk(ideal i, int reduction intvec v, intvec w); 351 355 TYPE: ideal 352 356 PURPOSE: compute the standard basis of the ideal w.r.t. the … … 372 376 373 377 ideal G = fetch(xR, Go); 374 G = system("Mfwalk", G, curr_weight, target_weight );378 G = system("Mfwalk", G, curr_weight, target_weight, reduction, printout); 375 379 376 380 setring xR; … … 387 391 ring r = 32003,(z,y,x), lp; 388 392 ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z; 389 fwalk(I); 393 int reduction = 1; 394 int printout = 1; 395 fwalk(I,reduction,printout); 390 396 } 391 397 … … 437 443 } 438 444 439 proc pwalk(ideal Go, int n1, int n2, list #)440 "SYNTAX: pwalk(int d, ideal i, int n1, int n2 );441 pwalk(int d, ideal i, int n1, int n2, int vec v, intvec w);445 proc pwalk(ideal Go, int n1, int n2, int reduction, int printout, list #) 446 "SYNTAX: pwalk(int d, ideal i, int n1, int n2, int reduction, int printout); 447 pwalk(int d, ideal i, int n1, int n2, int reduction, int printout, intvec v, intvec w); 442 448 TYPE: ideal 443 449 PURPOSE: compute the standard basis of the ideal, calculated via … … 477 483 ideal G = fetch(xR, Go); 478 484 479 G = system("Mpwalk", G, n1, n2, curr_weight, target_weight,nP);480 485 G = system("Mpwalk",G,n1,n2,curr_weight,target_weight,nP,reduction,printout); 486 481 487 setring xR; 482 //kill Go; 488 //kill Go; //unused 483 489 484 490 keepring basering; … … 492 498 ring r = 32003,(z,y,x), lp; 493 499 ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z; 494 //I = std(I); 495 //ring rr = 32003,(z,y,x),lp; 496 //ideal I = fetch(r,I); 497 pwalk(I,2,2); 500 int reduction = 1; 501 int printout = 2; 502 pwalk(I,2,2,reduction,printout); 498 503 } 499 504 -
Singular/LIB/modwalk.lib
-
Property
mode
changed from
100644
to100755
re9478b r38301d 16 16 17 17 PROCEDURES: 18 modWalk(ideal,int,int[,int,int,int,int]); standard basis conversion of I using modular methods (chinese remainder) 18 19 modWalk(I,#); standard basis conversion of I by Groebner Walk using modular methods 20 modrWalk(I,radius,pertdeg,#); standard basis conversion of I by Random Walk using modular methods 21 modfWalk(I,#); standard basis conversion of I by Fractal Walk using modular methods 22 modfrWalk(I,radius,#); standard basis conversion of I by Random Fractal Walk using modular methods 19 23 "; 20 24 21 LIB "poly.lib";22 LIB "ring.lib";23 LIB "parallel.lib";24 25 LIB "rwalk.lib"; 25 26 LIB "grwalk.lib"; 26 LIB "modstd.lib"; 27 28 29 //////////////////////////////////////////////////////////////////////////////// 30 31 static proc modpWalk(def II, int p, int variant, list #) 32 "USAGE: modpWalk(I,p,#); I ideal, p integer, variant integer 33 ASSUME: If size(#) > 0, then 34 #[1] is an intvec describing the current weight vector 35 #[2] is an intvec describing the target weight vector 36 RETURN: ideal - a standard basis of I mod p, integer - p 37 NOTE: The procedure computes a standard basis of the ideal I modulo p and 38 fetches the result to the basering. 39 EXAMPLE: example modpWalk; shows an example 40 " 41 { 42 option(redSB); 43 int k,nvar@r; 44 def R0 = basering; 45 string ordstr_R0 = ordstr(R0); 46 list rl = ringlist(R0); 47 int sizerl = size(rl); 48 int neg = 1 - attrib(R0,"global"); 49 if(typeof(II) == "ideal") 50 { 51 ideal I = II; 27 LIB "modular.lib"; 28 29 proc modWalk(ideal I, list #) 30 "USAGE: modWalk(I, [, v, w]); I ideal, v intvec, w intvec 31 RETURN: a standard basis of I 32 NOTE: The procedure computes a standard basis of I (over the rational 33 numbers) by using modular methods. 34 SEE ALSO: modular 35 EXAMPLE: example modWalk; shows an example" 36 { 37 /* read optional parameter */ 38 if (size(#) > 0) { 39 if (size(#) == 1) { 40 intvec w = #[1]; 41 } 42 if (size(#) == 2) { 43 intvec v = #[1]; 44 intvec w = #[2]; 45 } 46 if (size(#) > 2 || typeof(#[1]) != "intvec") { 47 ERROR("wrong optional parameter"); 48 } 49 } 50 51 /* save options */ 52 intvec opt = option(get); 53 option(redSB); 54 55 /* set additional parameters */ 56 int reduction = 1; 57 int printout = 0; 58 59 /* call modular() */ 60 if (size(#) > 0) { 61 I = modular("gwalk", list(I,reduction,printout,#)); 62 } 63 else { 64 I = modular("gwalk", list(I,reduction,printout)); 65 } 66 67 /* return the result */ 68 attrib(I, "isSB", 1); 69 option(set, opt); 70 return(I); 71 } 72 example 73 { 74 "EXAMPLE:"; 75 echo = 2; 76 ring R1 = 0, (x,y,z,t), dp; 77 ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; 78 I = std(I); 79 ring R2 = 0, (x,y,z,t), lp; 80 ideal I = fetch(R1, I); 81 ideal J = modWalk(I); 82 J; 83 } 84 85 proc modrWalk(ideal I, int radius, int pertdeg, list #) 86 "USAGE: modrWalk(I, radius, pertdeg[, v, w]); 87 I ideal, radius int, pertdeg int, v intvec, w intvec 88 RETURN: a standard basis of I 89 NOTE: The procedure computes a standard basis of I (over the rational 90 numbers) by using modular methods. 91 SEE ALSO: modular 92 EXAMPLE: example modrWalk; shows an example" 93 { 94 /* read optional parameter */ 95 if (size(#) > 0) { 96 if (size(#) == 1) { 97 intvec w = #[1]; 98 } 99 if (size(#) == 2) { 100 intvec v = #[1]; 101 intvec w = #[2]; 102 } 103 if (size(#) > 2 || typeof(#[1]) != "intvec") { 104 ERROR("wrong optional parameter"); 105 } 106 } 107 108 /* save options */ 109 intvec opt = option(get); 110 option(redSB); 111 112 /* set additional parameters */ 113 int reduction = 1; 114 int printout = 0; 115 116 /* call modular() */ 117 if (size(#) > 0) { 118 I = modular("rwalk", list(I,radius,pertdeg,reduction,printout,#)); 119 } 120 else { 121 I = modular("rwalk", list(I,radius,pertdeg,reduction,printout)); 122 } 123 124 /* return the result */ 125 attrib(I, "isSB", 1); 126 option(set, opt); 127 return(I); 128 } 129 example 130 { 131 "EXAMPLE:"; 132 echo = 2; 133 ring R1 = 0, (x,y,z,t), dp; 134 ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; 135 I = std(I); 136 ring R2 = 0, (x,y,z,t), lp; 137 ideal I = fetch(R1, I); 52 138 int radius = 2; 53 int pert_deg = 2; 54 } 55 if(typeof(II) == "list" && typeof(II[1]) == "ideal") 56 { 57 ideal I = II[1]; 58 if(size(II) == 2) 59 { 60 int radius = II[2]; 61 int pert_deg = 2; 62 } 63 if(size(II) == 3) 64 { 65 int radius = II[2]; 66 int pert_deg = II[3]; 67 } 68 } 69 rl[1] = p; 70 int h = homog(I); 71 def @r = ring(rl); 72 setring @r; 73 ideal i = fetch(R0,I); 74 string order; 75 if(system("nblocks") <= 2) 76 { 77 if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") + find(ordstr_R0, "rp") <= 0) 78 { 79 order = "simple"; 80 } 81 } 82 83 //------------------------- make i homogeneous ----------------------------- 84 if(!mixedTest() && !h) 85 { 86 if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)) 87 { 88 if(!((order == "simple") || (sizerl > 4))) 89 { 90 list rl@r = ringlist(@r); 91 nvar@r = nvars(@r); 92 intvec w; 93 for(k = 1; k <= nvar@r; k++) 94 { 95 w[k] = deg(var(k)); 96 } 97 w[nvar@r + 1] = 1; 98 rl@r[2][nvar@r + 1] = "homvar"; 99 rl@r[3][2][2] = w; 100 def HomR = ring(rl@r); 101 setring HomR; 102 ideal i = imap(@r, i); 103 i = homog(i, homvar); 104 } 105 } 106 } 107 108 //------------------------- compute a standard basis mod p ----------------------------- 109 110 if(variant == 1) 111 { 112 if(size(#)>0) 113 { 114 i = rwalk(i,radius,pert_deg,#); 115 // rwalk(i,radius,pert_deg,#); std(i); 116 } 117 else 118 { 119 i = rwalk(i,radius,pert_deg); 120 } 121 } 122 if(variant == 2) 123 { 124 if(size(#) == 2) 125 { 126 i = gwalk(i,#); 127 } 128 else 129 { 130 i = gwalk(i); 131 } 132 } 133 if(variant == 3) 134 { 135 if(size(#) == 2) 136 { 137 i = frandwalk(i,radius,#); 138 } 139 else 140 { 141 i = frandwalk(i,radius); 142 } 143 } 144 if(variant == 4) 145 { 146 if(size(#) == 2) 147 { 148 i=fwalk(i,#); 149 } 150 else 151 { 152 i=fwalk(i); 153 } 154 } 155 if(variant == 5) 156 { 157 if(size(#) == 2) 158 { 159 i=prwalk(i,radius,pert_deg,pert_deg,#); 160 } 161 else 162 { 163 i=prwalk(i,radius,pert_deg,pert_deg); 164 } 165 } 166 if(variant == 6) 167 { 168 if(size(#) == 2) 169 { 170 i=pwalk(i,pert_deg,pert_deg,#); 171 } 172 else 173 { 174 i=pwalk(i,pert_deg,pert_deg); 175 } 176 } 177 178 if(!mixedTest() && !h) 179 { 180 if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)) 181 { 182 if(!((order == "simple") || (sizerl > 4))) 183 { 184 i = subst(i, homvar, 1); 185 i = simplify(i, 34); 186 setring @r; 187 i = imap(HomR, i); 188 i = interred(i); 189 kill HomR; 190 } 191 } 192 } 193 setring R0; 194 return(list(fetch(@r,i),p)); 195 } 196 example 197 { 198 "EXAMPLE:"; echo = 2; 199 option(redSB); 200 201 int p = 181; 202 intvec a = 2,1,3,4; 203 intvec b = 1,9,1,1; 204 ring ra = 0,x(1..4),(a(a),lp); 205 ideal I = std(cyclic(4)); 206 ring rb = 0,x(1..4),(a(b),lp); 207 ideal I = imap(ra,I); 208 modpWalk(I,p,1,a,b); 209 std(I); 210 } 211 212 //////////////////////////////////////////////////////////////////////////////// 213 214 proc modWalk(def II, int variant, list #) 215 "USAGE: modWalk(II); II ideal or list(ideal,int) 216 ASSUME: If variant = 217 @* - 1 the Random Walk algorithm with radius II[2] is applied 218 to II[1] if II = list(ideal, int). It is applied to II with radius 2 219 if II is an ideal. 220 @* - 2, the Groebner Walk algorithm is applied to II[1] or to II, respectively. 221 @* - 3, the Fractal Walk algorithm with random element is applied to II[1] or II, 222 respectively. 223 @* - 4, the Fractal Walk algorithm is applied to II[1] or II, respectively. 224 @* - 5, the Perturbation Walk algorithm with random element is applied to II[1] 225 or II, respectively, with radius II[3] and perturbation degree II[2]. 226 @* - 6, the Perturbation Walk algorithm is applied to II[1] or II, respectively, 227 with perturbation degree II[3]. 228 If size(#) > 0, then # contains either 1, 2 or 4 integers such that 229 @* - #[1] is the number of available processors for the computation, 230 @* - #[2] is an optional parameter for the exactness of the computation, 231 if #[2] = 1, the procedure computes a standard basis for sure, 232 @* - #[3] is the number of primes until the first lifting, 233 @* - #[4] is the constant number of primes between two liftings until 234 the computation stops. 235 RETURN: a standard basis of I if no warning appears. 236 NOTE: The procedure converts a standard basis of I (over the rational 237 numbers) from the ordering \"a(v),lp\", "dp\" or \"Dp\" to the ordering 238 \"(a(w),lp\" or \"a(1,0,...,0),lp\" by using modular methods. 239 By default the procedure computes a standard basis of I for sure, but 240 if the optional parameter #[2] = 0, it computes a standard basis of I 241 with high probability. 242 EXAMPLE: example modWalk; shows an example 243 " 244 { 245 int TT = timer; 246 int RT = rtimer; 247 int i,j,pTest,sizeTest,weighted,n1; 248 bigint N; 249 250 def R0 = basering; 251 list rl = ringlist(R0); 252 if((npars(R0) > 0) || (rl[1] > 0)) 253 { 254 ERROR("Characteristic of basering should be zero, basering should have no parameters."); 255 } 256 257 if(typeof(II) == "ideal") 258 { 259 ideal I = II; 260 kill II; 261 list II; 262 II[1] = I; 263 II[2] = 2; 264 II[3] = 2; 265 } 266 else 267 { 268 if(typeof(II) == "list" && typeof(II[1]) == "ideal") 269 { 270 ideal I = II[1]; 271 if(size(II) == 1) 272 { 273 II[2] = 2; 274 II[3] = 2; 275 } 276 if(size(II) == 2) 277 { 278 II[3] = 2; 279 } 280 281 } 282 else 283 { 284 ERROR("Unexpected type of input."); 285 } 286 } 287 288 //-------------------- Initialize optional parameters ------------------------ 289 n1 = system("--cpus"); 290 if(size(#) == 0) 291 { 292 int exactness = 1; 293 int n2 = 10; 294 int n3 = 10; 295 } 296 else 297 { 298 if(size(#) == 1) 299 { 300 if(typeof(#[1]) == "int") 301 { 302 if(#[1] < n1) 303 { 304 n1 = #[1]; 305 } 306 int exactness = 1; 307 if(n1 >= 10) 308 { 309 int n2 = n1 + 1; 310 int n3 = n1; 311 } 312 else 313 { 314 int n2 = 10; 315 int n3 = 10; 316 } 317 } 318 else 319 { 320 ERROR("Unexpected type of input."); 321 } 322 } 323 if(size(#) == 2) 324 { 325 if(typeof(#[1]) == "int" && typeof(#[2]) == "int") 326 { 327 if(#[1] < n1) 328 { 329 n1 = #[1]; 330 } 331 int exactness = #[2]; 332 if(n1 >= 10) 333 { 334 int n2 = n1 + 1; 335 int n3 = n1; 336 } 337 else 338 { 339 int n2 = 10; 340 int n3 = 10; 341 } 342 } 343 else 344 { 345 if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec") 346 { 347 intvec curr_weight = #[1]; 348 intvec target_weight = #[2]; 349 weighted = 1; 350 int n2 = 10; 351 int n3 = 10; 352 int exactness = 1; 353 } 354 else 355 { 356 ERROR("Unexpected type of input."); 357 } 358 } 359 } 360 if(size(#) == 3) 361 { 362 if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int") 363 { 364 intvec curr_weight = #[1]; 365 intvec target_weight = #[2]; 366 weighted = 1; 367 n1 = #[3]; 368 int n2 = 10; 369 int n3 = 10; 370 int exactness = 1; 371 } 372 else 373 { 374 ERROR("Unexpected type of input."); 375 } 376 } 377 if(size(#) == 4) 378 { 379 if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int" 380 && typeof(#[4]) == "int") 381 { 382 intvec curr_weight = #[1]; 383 intvec target_weight = #[2]; 384 weighted = 1; 385 if(#[1] < n1) 386 { 387 n1 = #[3]; 388 } 389 int exactness = #[4]; 390 if(n1 >= 10) 391 { 392 int n2 = n1 + 1; 393 int n3 = n1; 394 } 395 else 396 { 397 int n2 = 10; 398 int n3 = 10; 399 } 400 } 401 else 402 { 403 if(typeof(#[1]) == "int" && typeof(#[2]) == "int" && typeof(#[3]) == "int" && typeof(#[4]) == "int") 404 { 405 if(#[1] < n1) 406 { 407 n1 = #[1]; 408 } 409 int exactness = #[2]; 410 if(n1 >= #[3]) 411 { 412 int n2 = n1 + 1; 413 } 414 else 415 { 416 int n2 = #[3]; 417 } 418 if(n1 >= #[4]) 419 { 420 int n3 = n1; 421 } 422 else 423 { 424 int n3 = #[4]; 425 } 426 } 427 else 428 { 429 ERROR("Unexpected type of input."); 430 } 431 } 432 } 433 if(size(#) == 6) 434 { 435 if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int" && typeof(#[4]) == "int" && typeof(#[5]) == "int" && typeof(#[6]) == "int") 436 { 437 intvec curr_weight = #[1]; 438 intvec target_weight = #[2]; 439 weighted = 1; 440 if(#[3] < n1) 441 { 442 n1 = #[3]; 443 } 444 int exactness = #[4]; 445 if(n1 >= #[5]) 446 { 447 int n2 = n1 + 1; 448 } 449 else 450 { 451 int n2 = #[5]; 452 } 453 if(n1 >= #[6]) 454 { 455 int n3 = n1; 456 } 457 else 458 { 459 int n3 = #[6]; 460 } 461 } 462 else 463 { 464 ERROR("Expected list(intvec,intvec,int,int,int,int) as optional parameter list."); 465 } 466 } 467 if(size(#) == 1 || size(#) == 5 || size(#) > 6) 468 { 469 ERROR("Expected 0,2,3,4 or 5 optional arguments."); 470 } 471 } 472 if(printlevel >= 10) 473 { 474 "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)+", exactness = "+string(exactness); 475 } 476 477 //------------------------- Save current options ----------------------------- 478 intvec opt = option(get); 479 //option(redSB); 480 481 //-------------------- Initialize the list of primes ------------------------- 482 int tt = timer; 483 int rt = rtimer; 484 int en = 2134567879; 485 int an = 1000000000; 486 intvec L = primeList(I,n2); 487 if(n2 > 4) 488 { 489 // L[5] = prime(random(an,en)); 490 } 491 if(printlevel >= 10) 492 { 493 "CPU-time for primeList: "+string(timer-tt)+" seconds."; 494 "Real-time for primeList: "+string(rtimer-rt)+" seconds."; 495 } 496 int h = homog(I); 497 list P,T1,T2,LL,Arguments,PP; 498 ideal J,K,H; 499 500 //------------------- parallelized Groebner Walk in positive characteristic -------------------- 501 502 if(weighted) 503 { 504 for(i=1; i<=size(L); i++) 505 { 506 Arguments[i] = list(II,L[i],variant,list(curr_weight,target_weight)); 507 } 508 } 509 else 510 { 511 for(i=1; i<=size(L); i++) 512 { 513 Arguments[i] = list(II,L[i],variant); 514 } 515 } 516 P = parallelWaitAll("modpWalk",Arguments); 517 for(i=1; i<=size(P); i++) 518 { 519 T1[i] = P[i][1]; 520 T2[i] = bigint(P[i][2]); 521 } 522 523 while(1) 524 { 525 LL = deleteUnluckyPrimes(T1,T2,h); 526 T1 = LL[1]; 527 T2 = LL[2]; 528 //------------------- Now all leading ideals are the same -------------------- 529 //------------------- Lift results to basering via farey --------------------- 530 531 tt = timer; rt = rtimer; 532 N = T2[1]; 533 for(i=2; i<=size(T2); i++) 534 { 535 N = N*T2[i]; 536 } 537 H = chinrem(T1,T2); 538 J = parallelFarey(H,N,n1); 539 //J=farey(H,N); 540 if(printlevel >= 10) 541 { 542 "CPU-time for lifting-process is "+string(timer - tt)+" seconds."; 543 "Real-time for lifting-process is "+string(rtimer - rt)+" seconds."; 544 } 545 546 //---------------- Test if we already have a standard basis of I -------------- 547 548 tt = timer; rt = rtimer; 549 pTest = pTestSB(I,J,L,variant); 550 //pTest = primeTestSB(I,J,L,variant); 551 if(printlevel >= 10) 552 { 553 "CPU-time for pTest is "+string(timer - tt)+" seconds."; 554 "Real-time for pTest is "+string(rtimer - rt)+" seconds."; 555 } 556 if(pTest) 557 { 558 if(printlevel >= 10) 559 { 560 "CPU-time for computation without final tests is "+string(timer - TT)+" seconds."; 561 "Real-time for computation without final tests is "+string(rtimer - RT)+" seconds."; 562 } 563 attrib(J,"isSB",1); 564 if(exactness == 0) 565 { 566 option(set, opt); 567 return(J); 568 } 569 else 570 { 571 tt = timer; 572 rt = rtimer; 573 sizeTest = 1 - isIdealIncluded(I,J,n1); 574 if(printlevel >= 10) 575 { 576 "CPU-time for checking if I subset <G> is "+string(timer - tt)+" seconds."; 577 "Real-time for checking if I subset <G> is "+string(rtimer - rt)+" seconds."; 578 } 579 if(sizeTest == 0) 580 { 581 tt = timer; 582 rt = rtimer; 583 K = std(J); 584 if(printlevel >= 10) 585 { 586 "CPU-time for last std-computation is "+string(timer - tt)+" seconds."; 587 "Real-time for last std-computation is "+string(rtimer - rt)+" seconds."; 588 } 589 if(size(reduce(K,J)) == 0) 590 { 591 option(set, opt); 592 return(J); 593 } 594 } 595 } 596 } 597 //-------------- We do not already have a standard basis of I, therefore do the main computation for more primes -------------- 598 599 T1 = H; 600 T2 = N; 601 j = size(L)+1; 602 tt = timer; rt = rtimer; 603 L = primeList(I,n3,L,n1); 604 if(printlevel >= 10) 605 { 606 "CPU-time for primeList: "+string(timer-tt)+" seconds."; 607 "Real-time for primeList: "+string(rtimer-rt)+" seconds."; 608 } 609 Arguments = list(); 610 PP = list(); 611 if(weighted) 612 { 613 for(i=j; i<=size(L); i++) 614 { 615 //Arguments[i-j+1] = list(II,L[i],variant,list(curr_weight,target_weight)); 616 Arguments[size(Arguments)+1] = list(II,L[i],variant,list(curr_weight,target_weight)); 617 } 618 } 619 else 620 { 621 for(i=j; i<=size(L); i++) 622 { 623 //Arguments[i-j+1] = list(II,L[i],variant); 624 Arguments[size(Arguments)+1] = list(II,L[i],variant); 625 } 626 } 627 PP = parallelWaitAll("modpWalk",Arguments); 628 if(printlevel >= 10) 629 { 630 "parallel modpWalk"; 631 } 632 for(i=1; i<=size(PP); i++) 633 { 634 //P[size(P) + 1] = PP[i]; 635 T1[size(T1) + 1] = PP[i][1]; 636 T2[size(T2) + 1] = bigint(PP[i][2]); 637 } 638 } 639 if(printlevel >= 10) 640 { 641 "CPU-time for computation with final tests is "+string(timer - TT)+" seconds."; 642 "Real-time for computation with final tests is "+string(rtimer - RT)+" seconds."; 643 } 644 } 645 646 example 647 { 648 "EXAMPLE:"; 649 echo = 2; 650 ring R=0,(x,y,z),lp; 651 ideal I=-x+y2z-z,xz+1,x2+y2-1; 652 // I is a standard basis in dp 653 ideal J = modWalk(I,1); 654 J; 655 } 656 657 //////////////////////////////////////////////////////////////////////////////// 658 static proc isIdealIncluded(ideal I, ideal J, int n1) 659 "USAGE: isIdealIncluded(I,J,int n1); I ideal, J ideal, n1 integer 660 " 661 { 662 if(n1 > 1) 663 { 664 int k; 665 list args,results; 666 for(k=1; k<=size(I); k++) 667 { 668 args[k] = list(ideal(I[k]),J,1); 669 } 670 results = parallelWaitAll("reduce",args); 671 for(k=1; k<=size(results); k++) 672 { 673 if(results[k] == 0) 674 { 675 return(1); 676 } 677 } 678 return(0); 679 } 680 else 681 { 682 if(reduce(I,J,1) == 0) 683 { 684 return(1); 685 } 686 else 687 { 688 return(0); 689 } 690 } 691 } 692 693 //////////////////////////////////////////////////////////////////////////////// 694 static proc parallelChinrem(list T1, list T2, int n1) 695 "USAGE: parallelChinrem(T1,T2); T1 list of ideals, T2 list of primes, n1 integer" 696 { 697 int i,j,k; 698 699 ideal H,J; 700 701 list arguments_chinrem,results_chinrem; 702 for(i=1; i<=size(T1); i++) 703 { 704 J = ideal(T1[i]); 705 attrib(J,"isSB",1); 706 arguments_chinrem[size(arguments_chinrem)+1] = list(list(J),T2); 707 } 708 results_chinrem = parallelWaitAll("chinrem",arguments_chinrem); 709 for(j=1; j <= size(results_chinrem); j++) 710 { 711 J = results_chinrem[j]; 712 attrib(J,"isSB",1); 713 if(isIdealIncluded(J,H,n1) == 0) 714 { 715 if(H == 0) 716 { 717 H = J; 718 } 719 else 720 { 721 H = H,J; 722 } 723 } 724 } 725 return(H); 726 } 727 728 //////////////////////////////////////////////////////////////////////////////// 729 static proc parallelFarey(ideal H, bigint N, int n1) 730 "USAGE: parallelFarey(H,N,n1); H ideal, N bigint, n1 integer 731 " 732 { 733 int i,j; 734 int ii = 1; 735 list arguments_farey,results_farey; 736 for(i=1; i<=size(H); i++) 737 { 738 for(j=1; j<=size(H[i]); j++) 739 { 740 arguments_farey[size(arguments_farey)+1] = list(H[i][j],N); 741 } 742 } 743 results_farey = parallelWaitAll("farey",arguments_farey); 744 ideal J,K; 745 poly f_farey; 746 while(ii<=size(results_farey)) 747 { 748 for(i=1; i<=size(H); i++) 749 { 750 f_farey = 0; 751 for(j=1; j<=size(H[i]); j++) 752 { 753 f_farey = f_farey + results_farey[ii][1]; 754 ii++; 755 } 756 K = ideal(f_farey); 757 attrib(K,"isSB",1); 758 attrib(J,"isSB",1); 759 if(isIdealIncluded(K,J,n1) == 0) 760 { 761 if(J == 0) 762 { 763 J = K; 764 } 765 else 766 { 767 J = J,K; 768 } 769 } 770 } 771 } 772 return(J); 773 } 774 ////////////////////////////////////////////////////////////////////////////////// 775 static proc primeTestSB(def II, ideal J, list L, int variant, list #) 776 "USAGE: primeTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int 777 RETURN: 1 (resp. 0) if for a randomly chosen prime p that is not in L 778 J mod p is (resp. is not) a standard basis of I mod p 779 EXAMPLE: example primeTestSB; shows an example 780 " 781 { 782 if(typeof(II) == "ideal") 783 { 784 ideal I = II; 785 int radius = 2; 786 } 787 if(typeof(II) == "list") 788 { 789 ideal I = II[1]; 790 int radius = II[2]; 791 } 792 793 int i,j,k,p; 794 def R = basering; 795 list r = ringlist(R); 796 797 while(!j) 798 { 799 j = 1; 800 p = prime(random(1000000000,2134567879)); 801 for(i = 1; i <= size(L); i++) 802 { 803 if(p == L[i]) 804 { 805 j = 0; 806 break; 807 } 808 } 809 if(j) 810 { 811 for(i = 1; i <= ncols(I); i++) 812 { 813 for(k = 2; k <= size(I[i]); k++) 814 { 815 if((denominator(leadcoef(I[i][k])) mod p) == 0) 816 { 817 j = 0; 818 break; 819 } 820 } 821 if(!j) 822 { 823 break; 824 } 825 } 826 } 827 if(j) 828 { 829 if(!primeTest(I,p)) 830 { 831 j = 0; 832 } 833 } 834 } 835 r[1] = p; 836 def @R = ring(r); 837 setring @R; 838 ideal I = imap(R,I); 839 ideal J = imap(R,J); 840 attrib(J,"isSB",1); 841 842 int t = timer; 843 j = 1; 844 if(isIncluded(I,J) == 0) 845 { 846 j = 0; 847 } 848 if(printlevel >= 11) 849 { 850 "isIncluded(I,J) takes "+string(timer - t)+" seconds"; 851 "j = "+string(j); 852 } 853 t = timer; 854 if(j) 855 { 856 if(size(#) > 0) 857 { 858 ideal K = modpWalk(I,p,variant,#)[1]; 859 } 860 else 861 { 862 ideal K = modpWalk(I,p,variant)[1]; 863 } 864 t = timer; 865 if(isIncluded(J,K) == 0) 866 { 867 j = 0; 868 } 869 if(printlevel >= 11) 870 { 871 "isIncluded(K,J) takes "+string(timer - t)+" seconds"; 872 "j = "+string(j); 873 } 874 } 875 setring R; 876 877 return(j); 878 } 879 example 880 { "EXAMPLE:"; echo = 2; 881 intvec L = 2,3,5; 882 ring r = 0,(x,y,z),lp; 883 ideal I = x+1,x+y+1; 884 ideal J = x+1,y; 885 primeTestSB(I,I,L,1); 886 primeTestSB(I,J,L,1); 887 } 888 889 //////////////////////////////////////////////////////////////////////////////// 890 static proc pTestSB(ideal I, ideal J, list L, int variant, list #) 891 "USAGE: pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int 892 RETURN: 1 (resp. 0) if for a randomly chosen prime p that is not in L 893 J mod p is (resp. is not) a standard basis of I mod p 894 EXAMPLE: example pTestSB; shows an example 895 " 896 { 897 int i,j,k,p; 898 def R = basering; 899 list r = ringlist(R); 900 901 while(!j) 902 { 903 j = 1; 904 p = prime(random(1000000000,2134567879)); 905 for(i = 1; i <= size(L); i++) 906 { 907 if(p == L[i]) { j = 0; break; } 908 } 909 if(j) 910 { 911 for(i = 1; i <= ncols(I); i++) 912 { 913 for(k = 2; k <= size(I[i]); k++) 914 { 915 if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; } 916 } 917 if(!j){ break; } 918 } 919 } 920 if(j) 921 { 922 if(!primeTest(I,p)) { j = 0; } 923 } 924 } 925 r[1] = p; 926 def @R = ring(r); 927 setring @R; 928 ideal I = imap(R,I); 929 ideal J = imap(R,J); 930 attrib(J,"isSB",1); 931 932 int t = timer; 933 j = 1; 934 if(isIncluded(I,J) == 0) { j = 0; } 935 936 if(printlevel >= 11) 937 { 938 "isIncluded(I,J) takes "+string(timer - t)+" seconds"; 939 "j = "+string(j); 940 } 941 942 t = timer; 943 if(j) 944 { 945 if(size(#) > 0) 946 { 947 ideal K = modpStd(I,p,variant,#[1])[1]; 948 } 949 else 950 { 951 ideal K = groebner(I); 952 } 953 t = timer; 954 if(isIncluded(J,K) == 0) { j = 0; } 955 956 if(printlevel >= 11) 957 { 958 "isIncluded(J,K) takes "+string(timer - t)+" seconds"; 959 "j = "+string(j); 960 } 961 } 962 setring R; 963 return(j); 964 } 965 example 966 { "EXAMPLE:"; echo = 2; 967 intvec L = 2,3,5; 968 ring r = 0,(x,y,z),dp; 969 ideal I = x+1,x+y+1; 970 ideal J = x+1,y; 971 pTestSB(I,I,L,2); 972 pTestSB(I,J,L,2); 973 } 974 //////////////////////////////////////////////////////////////////////////////// 975 static proc mixedTest() 976 "USAGE: mixedTest(); 977 RETURN: 1 if ordering of basering is mixed, 0 else 978 EXAMPLE: example mixedTest(); shows an example 979 " 980 { 981 int i,p,m; 982 for(i = 1; i <= nvars(basering); i++) 983 { 984 if(var(i) > 1) 985 { 986 p++; 987 } 988 else 989 { 990 m++; 991 } 992 } 993 if((p > 0) && (m > 0)) { return(1); } 994 return(0); 995 } 996 example 997 { "EXAMPLE:"; echo = 2; 998 ring R1 = 0,(x,y,z),dp; 999 mixedTest(); 1000 ring R2 = 31,(x(1..4),y(1..3)),(ds(4),lp(3)); 1001 mixedTest(); 1002 ring R3 = 181,x(1..9),(dp(5),lp(4)); 1003 mixedTest(); 1004 } 139 int pertdeg = 3; 140 ideal J = modrWalk(I,radius,pertdeg); 141 J; 142 } 143 144 proc modfWalk(ideal I, list #) 145 "USAGE: modfWalk(I, [, v, w]); I ideal, v intvec, w intvec 146 RETURN: a standard basis of I 147 NOTE: The procedure computes a standard basis of I (over the rational 148 numbers) by using modular methods. 149 SEE ALSO: modular 150 EXAMPLE: example modfWalk; shows an example" 151 { 152 /* read optional parameter */ 153 if (size(#) > 0) { 154 if (size(#) == 1) { 155 intvec w = #[1]; 156 } 157 if (size(#) == 2) { 158 intvec v = #[1]; 159 intvec w = #[2]; 160 } 161 if (size(#) > 2 || typeof(#[1]) != "intvec") { 162 ERROR("wrong optional parameter"); 163 } 164 } 165 166 /* save options */ 167 intvec opt = option(get); 168 option(redSB); 169 170 /* set additional parameters */ 171 int reduction = 1; 172 int printout = 0; 173 174 /* call modular() */ 175 if (size(#) > 0) { 176 I = modular("fwalk", list(I,reduction,printout,#)); 177 } 178 else { 179 I = modular("fwalk", list(I,reduction,printout)); 180 } 181 182 /* return the result */ 183 attrib(I, "isSB", 1); 184 option(set, opt); 185 return(I); 186 } 187 example 188 { 189 "EXAMPLE:"; 190 echo = 2; 191 ring R1 = 0, (x,y,z,t), dp; 192 ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; 193 I = std(I); 194 ring R2 = 0, (x,y,z,t), lp; 195 ideal I = fetch(R1, I); 196 ideal J = modfWalk(I); 197 J; 198 } 199 200 proc modfrWalk(ideal I, int radius, list #) 201 "USAGE: modfrWalk(I, radius [, v, w]); I ideal, radius int, v intvec, w intvec 202 RETURN: a standard basis of I 203 NOTE: The procedure computes a standard basis of I (over the rational 204 numbers) by using modular methods. 205 SEE ALSO: modular 206 EXAMPLE: example modfrWalk; shows an example" 207 { 208 /* read optional parameter */ 209 if (size(#) > 0) { 210 if (size(#) == 1) { 211 intvec w = #[1]; 212 } 213 if (size(#) == 2) { 214 intvec v = #[1]; 215 intvec w = #[2]; 216 } 217 if (size(#) > 2 || typeof(#[1]) != "intvec") { 218 ERROR("wrong optional parameter"); 219 } 220 } 221 222 /* save options */ 223 intvec opt = option(get); 224 option(redSB); 225 226 /* set additional parameters */ 227 int reduction = 1; 228 int printout = 0; 229 230 /* call modular() */ 231 if (size(#) > 0) { 232 I = modular("frandwalk", list(I,radius,reduction,printout,#)); 233 } 234 else { 235 I = modular("frandwalk", list(I,radius,reduction,printout)); 236 } 237 238 /* return the result */ 239 attrib(I, "isSB", 1); 240 option(set, opt); 241 return(I); 242 } 243 example 244 { 245 "EXAMPLE:"; 246 echo = 2; 247 ring R1 = 0, (x,y,z,t), dp; 248 ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; 249 I = std(I); 250 ring R2 = 0, (x,y,z,t), lp; 251 ideal I = fetch(R1, I); 252 int radius = 2; 253 ideal J = modfrWalk(I,radius); 254 J; 255 } -
Property
mode
changed from
-
Singular/LIB/nfmodstd.lib
rd52203 r38301d 524 524 I; 525 525 } 526 527 ////////////////////////////////////////////////////////////////////////////// 528 529 /* 530 Benchmark Problems from 531 532 Boku, Decker, Fieker, Steenpass: Groebner Bases over Algebraic Number 533 Fields. 534 535 // 1 536 ring R = (0,a), (x,y,z), dp; 537 minpoly = (a^2+1); 538 poly f1 = (a+8)*x^2*y^2+5*x*y^3+(-a+3)*x^3*z 539 +x^2*y*z; 540 poly f2 = x^5+2*y^3*z^2+13*y^2*z^3+5*y*z^4; 541 poly f3 = 8*x^3+(a+12)*y^3+x*z^2+3; 542 poly f4 = (-a+7)*x^2*y^4+y^3*z^3+18*y^3*z^2; 543 ideal I1 = f1,f2,f3,f4; 544 545 // 2 546 ring R = (0,a), (x,y,z), dp; 547 minpoly = (a^5+a^2+2); 548 poly f1 = 2*x*y^4*z^2+(a-1)*x^2*y^3*z 549 +(2*a)*x*y*z^2+7*y^3+(7*a+1); 550 poly f2 = 2*x^2*y^4*z+(a)*x^2*y*z^2-x*y^2*z^2 551 +(2*a+3)*x^2*y*z-12*x+(12*a)*y; 552 poly f3 = (2*a)*y^5*z+x^2*y^2*z-x*y^3*z 553 +(-a)*x*y^3+y^4+2*y^2*z; 554 poly f4 = (3*a)*x*y^4*z^3+(a+1)*x^2*y^2*z 555 -x*y^3*z+4*y^3*z^2+(3*a)*x*y*z^3 556 +4*z^2-x+(a)*y; 557 ideal I2 = f1,f2,f3,f4; 558 559 // 3a 560 ring R = (0,a), (v,w,x,y,z), dp; 561 minpoly = (a^7-7*a+3); 562 poly f1 = (a)*v+(a-1)*w+x+(a+2)*y+z; 563 poly f2 = v*w+(a-1)*w*x+(a+2)*v*y+x*y+(a)*y*z; 564 poly f3 = (a)*v*w*x+(a+5)*w*x*y+(a)*v*w*z 565 +(a+2)*v*y*z+(a)*x*y*z; 566 poly f4 = (a-11)*v*w*x*y+(a+5)*v*w*x*z 567 +(a)*v*w*y*z+(a)*v*x*y*z 568 +(a)*w*x*y*z; 569 poly f5 = (a+3)*v*w*x*y*z+(a+23); 570 ideal I3a = f1,f2,f3,f4,f5; 571 572 // 3b 573 ring R = (0,a), (u,v,w,x,y,z), dp; 574 minpoly = (a^7-7*a+3); 575 poly f1 = (a)*u+(a+2)*v+w+x+y+z; 576 poly f2 = u*v+v*w+w*x+x*y+(a+3)*u*z+y*z; 577 poly f3 = u*v*w+v*w*x+(a+1)*w*x*y+u*v*z+u*y*z 578 +x*y*z; 579 poly f4 = (a-1)*u*v*w*x+v*w*x*y+u*v*w*z 580 +u*v*y*z+u*x*y*z+w*x*y*z; 581 poly f5 = u*v*w*x*y+(a+1)*u*v*w*x*z+u*v*w*y*z 582 +u*v*x*y*z+u*w*x*y*z+v*w*x*y*z; 583 poly f6 = u*v*w*x*y*z+(-a+2); 584 ideal I3b = f1,f2,f3,f4,f5,f6; 585 586 // 4 587 ring R = (0,a), (w,x,y,z), dp; 588 minpoly = (a^6+a^5+a^4+a^3+a^2+a+1); 589 poly f1 = (a+5)*w^3*x^2*y+(a-3)*w^2*x^3*y 590 +(a+7)*w*x^2*y^2; 591 poly f2 = (a)*w^5+(a+3)*w*x^2*y^2 592 +(a^2+11)*x^2*y^2*z; 593 poly f3 = (a+7)*w^3+12*x^3+4*w*x*y+(a)*z^3; 594 poly f4 = 3*w^3+(a-4)*x^3+x*y^2; 595 ideal I4 = f1,f2,f3,f4; 596 597 // 5 598 ring R = (0,a), (w,x,y,z), dp; 599 minpoly = (a^12-5*a^11+24*a^10-115*a^9+551*a^8 600 -2640*a^7+12649*a^6-2640*a^5+551*a^4 601 -115*a^3+24*a^2-5*a+1); 602 poly f1 = (2*a+3)*w*x^4*y^2+(a+1)*w^2*x^3*y*z 603 +2*w*x*y^2*z^3+(7*a-1)*x^3*z^4; 604 poly f2 = 2*w^2*x^4*y+w^2*x*y^2*z^2 605 +(-a)*w*x^2*y^2*z^2 606 +(a+11)*w^2*x*y*z^3-12*w*z^6 607 +12*x*z^6; 608 poly f3 = 2*x^5*y+w^2*x^2*y*z-w*x^3*y*z 609 -w*x^3*z^2+(a)*x^4*z^2+2*x^2*y*z^3; 610 poly f4 = 3*w*x^4*y^3+w^2*x^2*y*z^3 611 -w*x^3*y*z^3+(a+4)*x^3*y^2*z^3 612 +3*w*x*y^3*z^3+(4*a)*y^2*z^6-w*z^7 613 +x*z^7; 614 ideal I5 = f1,f2,f3,f4; 615 616 // 6 617 ring R = (0,a), (u,v,w,x,y,z), dp; 618 minpoly = (a^2+5*a+1); 619 poly f1 = u+v+w+x+y+z+(a); 620 poly f2 = u*v+v*w+w*x+x*y+y*z+(a)*u+(a)*z; 621 poly f3 = u*v*w+v*w*x+w*x*y+x*y*z+(a)*u*v 622 +(a)*u*z+(a)*y*z; 623 poly f4 = u*v*w*x+v*w*x*y+w*x*y*z+(a)*u*v*w 624 +(a)*u*v*z+(a)*u*y*z+(a)*x*y*z; 625 poly f5 = u*v*w*x*y+v*w*x*y*z+(a)*u*v*w*x 626 +(a)*u*v*w*z+(a)*u*v*y*z+(a)*u*x*y*z 627 +(a)*w*x*y*z; 628 poly f6 = u*v*w*x*y*z+(a)*u*v*w*x*y 629 +(a)*u*v*w*x*z+(a)*u*v*w*y*z 630 +(a)*u*v*x*y*z+(a)*u*w*x*y*z 631 +(a)*v*w*x*y*z; 632 poly f7 = (a)*u*v*w*x*y*z-1; 633 ideal I6 = f1,f2,f3,f4,f5,f6,f7; 634 635 // 7 636 ring R = (0,a), (w,x,y,z), dp; 637 minpoly = (a^8-16*a^7+19*a^6-a^5-5*a^4+13*a^3 638 -9*a^2+13*a+17); 639 poly f1 = (-a^2-1)*x^2*y+2*w*x*z-2*w 640 +(a^2+1)*y; 641 poly f2 = (a^3-a-3)*w^3*y+4*w*x^2*y+4*w^2*x*z 642 +2*x^3*z+(a)*w^2-10*x^2+4*w*y-10*x*z 643 +(2*a^2+a); 644 poly f3 = (a^2+a+11)*x*y*z+w*z^2-w-2*y; 645 poly f4 = -w*y^3+4*x*y^2*z+4*w*y*z^2+2*x*z^3 646 +(2*a^3+a^2)*w*y+4*y^2-10*x*z-10*z^2 647 +(3*a^2+5); 648 ideal I7 = f1,f2,f3,f4; 649 650 // 8 651 ring R = (0,a), (t,u,v,w,x,y,z), dp; 652 minpoly = (a^7+10*a^5+5*a^3+10*a+1); 653 poly f1 = v*x+w*y-x*z-w-y; 654 poly f2 = v*w-u*x+x*y-w*z+v+x+z; 655 poly f3 = t*w-w^2+x^2-t; 656 poly f4 = (-a)*v^2-u*y+y^2-v*z-z^2+u; 657 poly f5 = t*v+v*w+(-a^2-a-5)*x*y-t*z+w*z+v+x+z 658 +(a+1); 659 poly f6 = t*u+u*w+(-a-11)*v*x-t*y+w*y-x*z-t-u 660 +w+y; 661 poly f7 = w^2*y^3-w*x*y^3+x^2*y^3+w^2*y^2*z 662 -w*x*y^2*z+x^2*y^2*z+w^2*y*z^2 663 -w*x*y*z^2+x^2*y*z^2+w^2*z^3-w*x*z^3 664 +x^2*z^3; 665 poly f8 = t^2*u^3+t^2*u^2*v+t^2*u*v^2+t^2*v^3 666 -t*u^3*x-t*u^2*v*x-t*u*v^2*x-t*v^3*x 667 +u^3*x^2+u^2*v*x^2+u*v^2*x^2 668 +v^3*x^2; 669 ideal I8 = f1,f2,f3,f4,f5,f6,f7,f8; 670 */ -
Singular/LIB/primdec.lib
rd52203 r38301d 807 807 if((size(sact)==1)&&(m==2)) 808 808 { 809 l[2*i]=l[2*i-1]; 810 attrib(l[2*i],"isSB",1); 809 l[2*i]=std(l[2*i-1],sact[1]); 811 810 } 812 811 if((size(sact)==1)&&(m>2)) … … 875 874 int uytrewq; 876 875 int nva = nvars(basering); 877 int @k,@s,@n,@k1, zz;876 int @k,@s,@n,@k1,@zz; 878 877 list primary,lres0,lres1,act,@lh,@wh; 879 878 map phi,psi,phi1,psi1; … … 1056 1055 { 1057 1056 poly randp; 1058 for( zz=1;zz<nvars(basering);zz++)1057 for(@zz=1;@zz<nvars(basering);@zz++) 1059 1058 { 1060 1059 randp=randp 1061 +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var( zz);1060 +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(@zz); 1062 1061 } 1063 1062 randp=randp+var(nvars(basering)); … … 1069 1068 if (size(primary[2*@k])==0) 1070 1069 { 1071 for( zz=1;zz<size(primary[2*@k-1])-1;zz++)1070 for(@zz=1;@zz<size(primary[2*@k-1])-1;@zz++) 1072 1071 { 1073 1072 attrib(primary[2*@k-1],"isSB",1); 1074 if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][ zz]))1073 if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][@zz])) 1075 1074 { 1076 1075 primary[2*@k]=primary[2*@k-1]; … … 1100 1099 if(deg(lead(primary[2*@k-1][@n]))==1) 1101 1100 { 1102 for( zz=1;zz<=nva;zz++)1101 for(@zz=1;@zz<=nva;@zz++) 1103 1102 { 1104 if(lead(primary[2*@k-1][@n])/var( zz)!=0)1103 if(lead(primary[2*@k-1][@n])/var(@zz)!=0) 1105 1104 { 1106 jmap1[ zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]1105 jmap1[@zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n] 1107 1106 +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]); 1108 jmap2[ zz]=primary[2*@k-1][@n];1109 @qht[@n]=var( zz);1107 jmap2[@zz]=primary[2*@k-1][@n]; 1108 @qht[@n]=var(@zz); 1110 1109 } 1111 1110 } … … 9024 9023 int uytrewq; 9025 9024 int nva = nvars(basering); 9026 int @k,@s,@n,@k1, zz;9025 int @k,@s,@n,@k1,@zz; 9027 9026 list primary,lres0,lres1,act,@lh,@wh; 9028 9027 map phi,psi,phi1,psi1; … … 9207 9206 { 9208 9207 poly randp; 9209 for( zz=1;zz<nvars(basering);zz++)9208 for(@zz=1;@zz<nvars(basering);@zz++) 9210 9209 { 9211 9210 randp=randp 9212 +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var( zz);9211 +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(@zz); 9213 9212 } 9214 9213 randp=randp+var(nvars(basering)); … … 9220 9219 if (size(primary[2*@k])==0) 9221 9220 { 9222 for( zz=1;zz<size(primary[2*@k-1])-1;zz++)9221 for(@zz=1;@zz<size(primary[2*@k-1])-1;@zz++) 9223 9222 { 9224 9223 attrib(primary[2*@k-1],"isSB",1); 9225 if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][ zz]))9224 if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][@zz])) 9226 9225 { 9227 9226 primary[2*@k]=primary[2*@k-1]; … … 9251 9250 if(deg(lead(primary[2*@k-1][@n]))==1) 9252 9251 { 9253 for( zz=1;zz<=nva;zz++)9252 for(@zz=1;@zz<=nva;@zz++) 9254 9253 { 9255 if(lead(primary[2*@k-1][@n])/var( zz)!=0)9254 if(lead(primary[2*@k-1][@n])/var(@zz)!=0) 9256 9255 { 9257 jmap1[ zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]9256 jmap1[@zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n] 9258 9257 +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]); 9259 jmap2[ zz]=primary[2*@k-1][@n];9260 @qht[@n]=var( zz);9258 jmap2[@zz]=primary[2*@k-1][@n]; 9259 @qht[@n]=var(@zz); 9261 9260 } 9262 9261 } -
Singular/LIB/rwalk.lib
-
Property
mode
changed from
100644
to100755
re9478b r38301d 10 10 rwalk(ideal,int,int[,intvec,intvec]); standard basis of ideal via Random Walk algorithm 11 11 rwalk(ideal,int[,intvec,intvec]); standard basis of ideal via Random Perturbation Walk algorithm 12 fr walk(ideal,int[,intvec,intvec]); standard basis of ideal via Random Fractal Walk algorithm12 frandwalk(ideal,int[,intvec,intvec]); standard basis of ideal via Random Fractal Walk algorithm 13 13 "; 14 14 … … 141 141 * Random Walk * 142 142 ****************/ 143 proc rwalk(ideal Go, int radius, int pert_deg, list #)143 proc rwalk(ideal Go, int radius, int pert_deg, int reduction, int printout, list #) 144 144 "SYNTAX: rwalk(ideal i, int radius); 145 145 if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w); 146 intermediate Groebner bases are not reduced if reduction = 0 146 147 TYPE: ideal 147 148 PURPOSE: compute the standard basis of the ideal, calculated via … … 178 179 179 180 ideal G = fetch(xR, Go); 180 G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, basering);181 G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, reduction, printout); 181 182 182 183 setring xR; … … 196 197 int radius = 1; 197 198 int perturb_deg = 2; 198 rwalk(I,radius,perturb_deg); 199 int reduction = 0; 200 int printout = 1; 201 rwalk(I,radius,perturb_deg,reduction,printout); 199 202 } 200 203 … … 202 205 * Perturbation Walk with random element * 203 206 *****************************************/ 204 proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, list #)207 proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, int reduction, int printout, list #) 205 208 "SYNTAX: rwalk(ideal i, int radius); 206 209 if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w); … … 227 230 OSCTW= OrderStringalp_NP("al", #); 228 231 } 232 int nP = OSCTW[1]; 229 233 string ord_str = OSCTW[2]; 230 234 intvec curr_weight = OSCTW[3]; // original weight vector … … 238 242 239 243 ideal G = fetch(xR, Go); 240 G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg, basering); 244 G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg, 245 nP, reduction, printout); 241 246 242 247 setring xR; … … 257 262 int o_perturb_deg = 2; 258 263 int t_perturb_deg = 2; 259 prwalk(I,radius,o_perturb_deg,t_perturb_deg); 264 int reduction = 0; 265 int printout = 2; 266 prwalk(I,radius,o_perturb_deg,t_perturb_deg,reduction,printout); 260 267 } 261 268 … … 263 270 * Fractal Walk with random element * 264 271 ************************************/ 265 proc frandwalk(ideal Go, int radius, list #)266 "SYNTAX: frwalk(ideal i, int radius );267 frwalk(ideal i, int radius, int vec v, intvec w);272 proc frandwalk(ideal Go, int radius, int reduction, int printout, list #) 273 "SYNTAX: frwalk(ideal i, int radius, int reduction, int printout); 274 frwalk(ideal i, int radius, int reduction, int printout, intvec v, intvec w); 268 275 TYPE: ideal 269 276 PURPOSE: compute the standard basis of the ideal, calculated via … … 299 306 ideal G = fetch(xR, Go); 300 307 int pert_deg = 2; 301 G = system("Mfrwalk", G, curr_weight, target_weight, radius); 308 309 G = system("Mfrwalk", G, curr_weight, target_weight, radius, reduction, printout); 302 310 303 311 setring xR; … … 314 322 ring r = 0,(z,y,x), lp; 315 323 ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z; 316 frandwalk(I,2); 317 } 324 int reduction = 0; 325 frandwalk(I,2,0,1); 326 } -
Property
mode
changed from
-
Singular/LIB/surf.lib
rd52203 r38301d 274 274 275 275 string surf_call; i = 0; 276 276 277 277 if (isWindows()) 278 278 { … … 303 303 else 304 304 { 305 surf_call = "surfer"; 305 surf_call = "surfer"; 306 306 surf_call = surf_call + " " + l + " >/dev/null 2>&1"; 307 307 "Close window to exit from `surfer`."; 308 308 i = system("sh", surf_call); 309 310 if ( (i != 0) && isMacOSX() ) 309 310 if ( (i != 0) && isMacOSX() ) 311 311 { 312 312 "*!* Sorry: calling `surfer` failed ['"+surf_call+"']" + newline 313 313 + " (The shell returned the error code " + string(i) + "." + newline 314 + "But since we are on Mac OS X, let us try to open SURFER.app instead..." + newline 314 + "But since we are on Mac OS X, let us try to open SURFER.app instead..." + newline 315 315 + "Appropriate SURFER.app is available for instance at http://www.mathematik.uni-kl.de/~motsak/files/SURFER.dmg"; 316 316 317 317 // fallback, will only work if SURFER.app is available (e.g. in /Applications) 318 318 // get SURFER.app e.g. from http://www.mathematik.uni-kl.de/~motsak/files/SURFER.dmg 319 319 // note that the newer (Java-based) variant of Surfer may not support command line usage yet :( 320 320 321 321 surf_call = "open -a SURFER -W --args -t -s"; 322 322 surf_call = surf_call + " " + l + " >/dev/null 2>&1"; … … 324 324 i = system("sh", surf_call); 325 325 } 326 327 328 326 327 328 329 329 } 330 330 system("sh", "/bin/rm " + l); … … 381 381 { 382 382 string s = system("uname"); 383 383 384 384 for (int i = 1; i <= size(s)-2; i = i + 1) 385 385 { -
Singular/LIB/swalk.lib
-
Property
mode
changed from
100644
to100755
-
Property
mode
changed from
-
Singular/Makefile.am
rd52203 r38301d 98 98 mmalloc.h \ 99 99 mod_lib.h \ 100 omSingularConfig.h \101 100 links/ndbm.h \ 102 101 newstruct.h \ -
Singular/dyn_modules/singmathic/singmathic.cc
rd52203 r38301d 468 468 result->rtyp=NONE; 469 469 470 if (arg == NULL || arg->next != NULL || 470 if (arg == NULL || arg->next != NULL || 471 471 ((arg->Typ() != IDEAL_CMD) &&(arg->Typ() != MODUL_CMD))){ 472 472 WerrorS("Syntax: mathicgb(<ideal>/<module>)"); -
Singular/dyn_modules/syzextra/mod_main.cc
rd52203 r38301d 82 82 for (int j=0; j<l; j++) 83 83 if (id->m[j] != NULL && p_GetComp(id->m[j], r) > 0) 84 return TRUE; 84 return TRUE; 85 85 86 86 return FALSE; // rank: 1, only zero or no entries? can be an ideal OR module... BUT in the use-case should better be an ideal! … … 91 91 92 92 93 93 94 94 95 95 static inline void NoReturn(leftv& res) … … 1550 1550 const int iLimit = r->typ[pos].data.is.limit; 1551 1551 const ideal F = r->typ[pos].data.is.F; 1552 1552 1553 1553 ideal FF = id_Copy(F, r); 1554 1554 -
Singular/extra.cc
rd52203 r38301d 3794 3794 } 3795 3795 else 3796 #endif 3796 #endif 3797 3797 /*==================== Error =================*/ 3798 3798 Werror( "(extended) system(\"%s\",...) %s", sys_cmd, feNotImplemented ); -
Singular/maps_ip.cc
rd52203 r38301d 73 73 if (P!=0) 74 74 { 75 // WerrorS("Sorry 'napPermNumber' was lost in the refactoring process (due to Frank): needs to be fixed");76 // return TRUE;77 #if 178 75 // poly n_PermNumber(const number z, const int *par_perm, const int OldPar, const ring src, const ring dst); 79 76 res->data= (void *) n_PermNumber((number)data, par_perm, P, preimage_r, currRing); 80 #endif81 77 res->rtyp=POLY_CMD; 82 78 if (nCoeff_is_algExt(currRing->cf)) … … 87 83 { 88 84 assume( nMap != NULL ); 89 90 85 number a = nMap((number)data, preimage_r->cf, currRing->cf); 91 86 if (nCoeff_is_Extension(currRing->cf)) … … 150 145 } 151 146 } 152 else 153 if ( (what==IMAP_CMD) || /*(*/ (what==FETCH_CMD) /*)*/) /* && (nMap!=nCopy)*/ 147 else if ((what==IMAP_CMD) || (what==FETCH_CMD)) 154 148 { 155 149 for (i=R*C-1;i>=0;i--) … … 160 154 } 161 155 } 162 else /* if(what==MAP_CMD) */ 163 { 156 else /* (what==MAP_CMD) */ 157 { 158 assume(what==MAP_CMD); 164 159 matrix s=mpNew(N,maMaxDeg_Ma((ideal)data,preimage_r)); 165 160 for (i=R*C-1;i>=0;i--) … … 170 165 idDelete((ideal *)&s); 171 166 } 172 if (nCoeff_is_ Extension(currRing->cf))167 if (nCoeff_is_algExt(currRing->cf)) 173 168 { 174 169 for (i=R*C-1;i>=0;i--) -
Singular/misc_ip.cc
rd52203 r38301d 716 716 } while (v!=NULL); 717 717 718 #ifdef OM_SINGULAR_CONFIG_H719 718 // set global variable to show memory usage 720 719 extern int om_sing_opt_show_mem; 721 720 if (BVERBOSE(V_SHOW_MEM)) om_sing_opt_show_mem = 1; 722 721 else om_sing_opt_show_mem = 0; 723 #endif724 722 725 723 return FALSE; … … 1095 1093 #endif // HAVE_SIMPLEIPC 1096 1094 fe_reset_input_mode(); 1095 monitor(NULL,0); 1097 1096 #ifdef PAGE_TEST 1098 1097 mmEndStat(); -
Singular/test.cc
rd52203 r38301d 194 194 #include <Singular/links/ndbm.h> 195 195 #include <Singular/newstruct.h> 196 #include <Singular/omSingularConfig.h>197 196 #include <Singular/pcv.h> 198 197 #include <Singular/links/pipeLink.h> -
Singular/walk.cc
-
Property
mode
changed from
100644
to100755
re9478b r38301d 27 27 28 28 #define FIRST_STEP_FRACTAL // to define the first step of the fractal 29 //#define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if29 #define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if 30 30 // tau doesn't stay in the correct cone 31 31 … … 42 42 #include <Singular/ipshell.h> 43 43 #include <Singular/ipconv.h> 44 #include <coeffs/ffields.h> 44 45 #include <coeffs/coeffs.h> 45 46 #include <Singular/subexpr.h> 47 #include <polys/templates/p_Procs.h> 46 48 47 49 #include <polys/monomials/maps.h> … … 429 431 #endif 430 432 431 #ifdef CHECK_IDEAL_MWALK433 //#ifdef CHECK_IDEAL_MWALK 432 434 static void idString(ideal L, const char* st) 433 435 { … … 441 443 Print(" %s;", pString(L->m[nL-1])); 442 444 } 443 #endif445 //#endif 444 446 445 447 #if defined(CHECK_IDEAL_MWALK) || defined(ENDWALKS) … … 511 513 512 514 //unused 513 #if 0515 //#if 0 514 516 static void MivString(intvec* iva, intvec* ivb, intvec* ivc) 515 517 { … … 534 536 Print("%d)", (*ivc)[nV]); 535 537 } 536 #endif538 //#endif 537 539 538 540 /******************************************************************** … … 558 560 } 559 561 return p0; 562 } 563 564 /***************************************************************************** 565 * compute the gcd of the entries of the vectors curr_weight and diff_weight * 566 *****************************************************************************/ 567 static int simplify_gcd(intvec* curr_weight, intvec* diff_weight) 568 { 569 int j; 570 int nRing = currRing->N; 571 int gcd_tmp = (*curr_weight)[0]; 572 for (j=1; j<nRing; j++) 573 { 574 gcd_tmp = gcd(gcd_tmp, (*curr_weight)[j]); 575 if(gcd_tmp == 1) 576 { 577 break; 578 } 579 } 580 if(gcd_tmp != 1) 581 { 582 for (j=0; j<nRing; j++) 583 { 584 gcd_tmp = gcd(gcd_tmp, (*diff_weight)[j]); 585 if(gcd_tmp == 1) 586 { 587 break; 588 } 589 } 590 } 591 return gcd_tmp; 560 592 } 561 593 … … 774 806 for(i=nG-1; i>=0; i--) 775 807 { 776 mi = MpolyInitialForm(G->m[i], iv); 777 gi = G->m[i]; 778 808 mi = pHead(MpolyInitialForm(G->m[i], iv)); 809 //Print("\n **// test_w_in_ConeCC: lm(initial)= %s \n",pString(mi)); 810 gi = pHead(G->m[i]); 811 //Print("\n **// test_w_in_ConeCC: lm(ideal)= %s \n",pString(gi)); 779 812 if(mi == NULL) 780 813 { … … 953 986 } 954 987 955 /***************************************************************************** 956 * create a weight matrix order as intvec of an extra weight vector (a(iv), lp)*957 ****************************************************************************** /988 /********************************************************************************* 989 * create a weight matrix order as intvec of an extra weight vector (a(iv),M(iw)) * 990 **********************************************************************************/ 958 991 intvec* MivMatrixOrderRefine(intvec* iv, intvec* iw) 959 992 { 960 assume( iv->length() == iw->length());961 int i, nR = iv->length();962 993 assume((iv->length())*(iv->length()) == iw->length()); 994 int i,j, nR = iv->length(); 995 963 996 intvec* ivm = new intvec(nR*nR); 964 997 … … 966 999 { 967 1000 (*ivm)[i] = (*iv)[i]; 968 (*ivm)[i+nR] = (*iw)[i]; 969 } 970 for(i=2; i<nR; i++) 971 { 972 (*ivm)[i*nR+i-2] = 1; 1001 } 1002 for(i=1; i<nR; i++) 1003 { 1004 for(j=0; j<nR; j++) 1005 { 1006 (*ivm)[j+i*nR] = (*iw)[j+i*nR]; 1007 } 973 1008 } 974 1009 return ivm; … … 1612 1647 (* result)[i] = mpz_get_si(pert_vector[i]); 1613 1648 } 1614 1649 /* 1615 1650 j = 0; 1616 for(i=0; i<n V; i++)1651 for(i=0; i<niv; i++) 1617 1652 { 1618 1653 (* result1)[i] = mpz_get_si(pert_vector[i]); … … 1624 1659 } 1625 1660 } 1626 if(j > n V- 1)1661 if(j > niv - 1) 1627 1662 { 1628 1663 // Print("\n// MfPertwalk: geaenderter vector gleich Null! \n"); … … 1630 1665 goto CHECK_OVERFLOW; 1631 1666 } 1632 1667 /* 1633 1668 // check that the perturbed weight vector lies in the Groebner cone 1669 Print("\n========================================\n//** MfPertvector: test in Cone.\n"); 1634 1670 if(test_w_in_ConeCC(G,result1) != 0) 1635 1671 { … … 1647 1683 // Print("\n// Mfpertwalk: geaenderter vector liegt nicht in Groebnerkegel! \n"); 1648 1684 } 1649 1685 Print("\n ========================================\n");*/ 1650 1686 CHECK_OVERFLOW: 1651 1687 … … 1859 1895 } 1860 1896 1897 1898 /************************************************************** 1899 * Look for the position of the smallest absolut value in vec * 1900 **************************************************************/ 1901 static int MivAbsMaxArg(intvec* vec) 1902 { 1903 int k = MivAbsMax(vec); 1904 int i=0; 1905 while(1) 1906 { 1907 if((*vec)[i] == k || (*vec)[i] == -k) 1908 { 1909 break; 1910 } 1911 i++; 1912 } 1913 return i; 1914 } 1915 1916 1861 1917 /********************************************************************** 1862 1918 * Compute a next weight vector between curr_weight and target_weight * … … 1873 1929 1874 1930 int nRing = currRing->N; 1875 int checkRed, j, kkk,nG = IDELEMS(G);1931 int checkRed, j, nG = IDELEMS(G); 1876 1932 intvec* ivtemp; 1877 1933 … … 1911 1967 mpz_init(dcw); 1912 1968 1913 //int tn0, tn1, tz1, ncmp, gcd_tmp, ntmp;1914 1969 int gcd_tmp; 1915 1970 intvec* diff_weight = MivSub(target_weight, curr_weight); … … 1917 1972 intvec* diff_weight1 = MivSub(target_weight, curr_weight); 1918 1973 poly g; 1919 //poly g, gw; 1974 1920 1975 for (j=0; j<nG; j++) 1921 1976 { … … 1934 1989 mpz_set_si(MwWd, MivDotProduct(ivtemp, curr_weight)); 1935 1990 mpz_sub(s_zaehler, deg_w0_p1, MwWd); 1936 1937 1991 if(mpz_cmp(s_zaehler, t_null) != 0) 1938 1992 { 1939 1993 mpz_set_si(MwWd, MivDotProduct(ivtemp, diff_weight)); 1940 1994 mpz_sub(s_nenner, MwWd, deg_d0_p1); 1941 1942 1995 // check for 0 < s <= 1 1943 1996 if( (mpz_cmp(s_zaehler,t_null) > 0 && … … 1979 2032 } 1980 2033 } 1981 //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));2034 //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t)); 1982 2035 mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t)); 1983 2036 1984 2037 1985 // there is no 0<t<1 and define the next weight vector that is equal to the current weight vector 2038 // there is no 0<t<1 and define the next weight vector that is equal 2039 // to the current weight vector 1986 2040 if(mpz_cmp(t_nenner, t_null) == 0) 1987 2041 { … … 2054 2108 #endif 2055 2109 2056 // BOOLEAN isdwpos; 2057 2058 // construct a new weight vector 2110 // construct a new weight vector and check whether vec[j] is overflow, 2111 // i.e. vec[j] > 2^31. 2112 // If vec[j] doesn't overflow, define a weight vector. Otherwise, 2113 // report that overflow appears. In the second case, test whether the 2114 // the correctness of the new vector plays an important role 2115 2059 2116 for (j=0; j<nRing; j++) 2060 2117 { … … 2100 2157 } 2101 2158 } 2102 2159 // reduce the vector with the gcd 2160 if(mpz_cmp_si(ggt,1) != 0) 2161 { 2162 for (j=0; j<nRing; j++) 2163 { 2164 mpz_divexact(vec[j], vec[j], ggt); 2165 } 2166 } 2103 2167 #ifdef NEXT_VECTORS_CC 2104 2168 PrintS("\n// gcd of elements of the vector: "); … … 2106 2170 #endif 2107 2171 2108 /**********************************************************************2109 * construct a new weight vector and check whether vec[j] is overflow, *2110 * i.e. vec[j] > 2^31. *2111 * If vec[j] doesn't overflow, define a weight vector. Otherwise, *2112 * report that overflow appears. In the second case, test whether the *2113 * the correctness of the new vector plays an important role *2114 **********************************************************************/2115 kkk=0;2116 2172 for(j=0; j<nRing; j++) 2117 2173 { … … 2129 2185 2130 2186 REDUCTION: 2187 checkRed = 1; 2131 2188 for (j=0; j<nRing; j++) 2132 2189 { 2133 (*diff_weight)[j] = mpz_get_si(vec[j]); 2134 } 2135 while(MivAbsMax(diff_weight) >= 5) 2136 { 2137 for (j=0; j<nRing; j++) 2138 { 2139 if(mpz_cmp_si(ggt,1)==0) 2140 { 2141 (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5); 2142 // Print("\n// vector[%d] = %d \n",j+1, (*diff_weight1)[j]); 2143 } 2144 else 2145 { 2146 mpz_divexact(vec[j], vec[j], ggt); 2147 (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5); 2148 // Print("\n// vector[%d] = %d \n",j+1, (*diff_weight1)[j]); 2149 } 2150 /* 2151 if((*diff_weight1)[j] == 0) 2152 { 2153 kkk = kkk + 1; 2154 } 2155 */ 2156 } 2157 2158 2159 /* 2160 if(kkk > nRing - 1) 2161 { 2162 // diff_weight was reduced to zero 2163 // Print("\n // MwalkNextWeightCC: geaenderter Vector gleich Null! \n"); 2164 goto TEST_OVERFLOW; 2165 } 2166 */ 2167 2168 if(test_w_in_ConeCC(G,diff_weight1) != 0) 2169 { 2170 Print("\n// MwalkNextWeightCC: geaenderter vector liegt in Groebnerkegel! \n"); 2171 for (j=0; j<nRing; j++) 2172 { 2173 (*diff_weight)[j] = (*diff_weight1)[j]; 2174 } 2175 if(MivAbsMax(diff_weight) < 5) 2176 { 2177 checkRed = 1; 2178 goto SIMPLIFY_GCD; 2179 } 2180 } 2181 else 2182 { 2183 // Print("\n// MwalkNextWeightCC: geaenderter vector liegt nicht in Groebnerkegel! \n"); 2184 break; 2185 } 2190 (*diff_weight1)[j] = mpz_get_si(vec[j]); 2191 } 2192 while(test_w_in_ConeCC(G,diff_weight1)) 2193 { 2194 for(j=0; j<nRing; j++) 2195 { 2196 (*diff_weight)[j] = (*diff_weight1)[j]; 2197 mpz_set_si(vec[j], (*diff_weight)[j]); 2198 } 2199 for(j=0; j<nRing; j++) 2200 { 2201 (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5); 2202 } 2203 } 2204 if(MivAbsMax(diff_weight)>10000) 2205 { 2206 for(j=0; j<nRing; j++) 2207 { 2208 (*diff_weight1)[j] = (*diff_weight)[j]; 2209 } 2210 j = 0; 2211 while(test_w_in_ConeCC(G,diff_weight1)) 2212 { 2213 (*diff_weight)[j] = (*diff_weight1)[j]; 2214 mpz_set_si(vec[j], (*diff_weight)[j]); 2215 j = MivAbsMaxArg(diff_weight1); 2216 (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5); 2217 } 2218 goto SIMPLIFY_GCD; 2186 2219 } 2187 2220 … … 2222 2255 mpz_clear(t_null); 2223 2256 2224 2225 2226 2257 if(Overflow_Error == FALSE) 2227 2258 { 2228 2259 Overflow_Error = nError; 2229 2260 } 2230 rComplete(currRing);2231 for( kkk=0; kkk<IDELEMS(G);kkk++)2232 { 2233 poly p=G->m[ kkk];2261 rComplete(currRing); 2262 for(j=0; j<IDELEMS(G); j++) 2263 { 2264 poly p=G->m[j]; 2234 2265 while(p!=NULL) 2235 2266 { … … 2271 2302 } 2272 2303 2273 /************************************************************** 2304 /******************************************************************** 2274 2305 * define and execute a new ring which order is (a(vb),a(va),lp,C) * 2275 * ************************************************************/ 2276 #if 0 2277 // unused 2306 * ******************************************************************/ 2278 2307 static void VMrHomogeneous(intvec* va, intvec* vb) 2279 2308 { … … 2353 2382 rChangeCurrRing(r); 2354 2383 } 2355 #endif 2384 2356 2385 2357 2386 /************************************************************** … … 2428 2457 } 2429 2458 2430 static ring VMrDefault1(intvec* va)2431 {2432 2433 if ((currRing->ppNoether)!=NULL)2434 {2435 pDelete(&(currRing->ppNoether));2436 }2437 if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||2438 ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))2439 {2440 sLastPrinted.CleanUp();2441 }2442 2443 ring r = (ring) omAlloc0Bin(sip_sring_bin);2444 int i, nv = currRing->N;2445 2446 r->cf = currRing->cf;2447 r->N = currRing->N;2448 2449 int nb = 4;2450 2451 //names2452 char* Q; // In order to avoid the corrupted memory, do not change.2453 r->names = (char **) omAlloc0(nv * sizeof(char_ptr));2454 for(i=0; i<nv; i++)2455 {2456 Q = currRing->names[i];2457 r->names[i] = omStrDup(Q);2458 }2459 2460 /*weights: entries for 3 blocks: NULL Made:???*/2461 r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));2462 r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));2463 for(i=0; i<nv; i++)2464 r->wvhdl[0][i] = (*va)[i];2465 2466 /* order: a,lp,C,0 */2467 r->order = (int *) omAlloc(nb * sizeof(int *));2468 r->block0 = (int *)omAlloc0(nb * sizeof(int *));2469 r->block1 = (int *)omAlloc0(nb * sizeof(int *));2470 2471 // ringorder a for the first block: var 1..nv2472 r->order[0] = ringorder_a;2473 r->block0[0] = 1;2474 r->block1[0] = nv;2475 2476 // ringorder lp for the second block: var 1..nv2477 r->order[1] = ringorder_lp;2478 r->block0[1] = 1;2479 r->block1[1] = nv;2480 2481 // ringorder C for the third block2482 // it is very important within "idLift",2483 // especially, by ring syz_ring=rCurrRingAssure_SyzComp();2484 // therefore, nb must be (nBlocks(currRing) + 1)2485 r->order[2] = ringorder_C;2486 2487 // the last block: everything is 02488 r->order[3] = 0;2489 2490 // polynomial ring2491 r->OrdSgn = 1;2492 2493 // complete ring intializations2494 2495 rComplete(r);2496 2497 //rChangeCurrRing(r);2498 return r;2499 }2500 2501 2459 /**************************************************************** 2502 2460 * define and execute a new ring with ordering (a(va),Wp(vb),C) * 2503 2461 * **************************************************************/ 2504 2505 2462 static ring VMrRefine(intvec* va, intvec* vb) 2506 2463 { … … 2576 2533 2577 2534 // complete ring intializations 2578 2535 2579 2536 rComplete(r); 2580 2537 … … 2806 2763 } 2807 2764 2808 2809 /* define a ring with parameters und change to it */ 2810 /* DefRingPar and DefRingParlp corrupt still memory */ 2765 /*************************************************** 2766 * define a ring with parameters und change to it * 2767 * DefRingPar and DefRingParlp corrupt still memory * 2768 ****************************************************/ 2811 2769 static void DefRingPar(intvec* va) 2812 2770 { … … 2956 2914 } 2957 2915 2958 //unused 2959 /************************************************************** 2960 * check wheather one or more components of a vector are zero * 2961 **************************************************************/ 2962 #if 0 2916 /************************************************************* 2917 * check whether one or more components of a vector are zero * 2918 *************************************************************/ 2963 2919 static int isNolVector(intvec* hilb) 2964 2920 { … … 2973 2929 return 0; 2974 2930 } 2975 #endif 2931 2932 /************************************************************* 2933 * check whether one or more components of a vector are <= 0 * 2934 *************************************************************/ 2935 static int isNegNolVector(intvec* hilb) 2936 { 2937 int i; 2938 for(i=hilb->length()-1; i>=0; i--) 2939 { 2940 if((* hilb)[i]<=0) 2941 { 2942 return 1; 2943 } 2944 } 2945 return 0; 2946 } 2947 2948 /************************************************************************** 2949 * Gomega is the initial ideal of G w. r. t. the current weight vector * 2950 * curr_weight. Check whether curr_weight lies on a border of the Groebner * 2951 * cone, i. e. check whether a monomial is divisible by a leading monomial * 2952 ***************************************************************************/ 2953 static ideal middleOfCone(ideal G, ideal Gomega) 2954 { 2955 BOOLEAN middle = FALSE; 2956 int i,j,N = IDELEMS(Gomega); 2957 poly p,lm,factor1,factor2; 2958 //PrintS("\n//** idCopy\n"); 2959 ideal Go = idCopy(G); 2960 2961 //PrintS("\n//** jetzt for-Loop!\n"); 2962 2963 // check whether leading monomials of G and Gomega coincide 2964 // and return NULL if not 2965 for(i=0; i<N; i++) 2966 { 2967 p = pCopy(Gomega->m[i]); 2968 lm = pCopy(pHead(G->m[i])); 2969 if(!pIsConstant(pSub(p,lm))) 2970 { 2971 //pDelete(&p); 2972 //pDelete(&lm); 2973 idDelete(&Go); 2974 return NULL; 2975 } 2976 //pDelete(&p); 2977 //pDelete(&lm); 2978 } 2979 for(i=0; i<N; i++) 2980 { 2981 for(j=0; j<N; j++) 2982 { 2983 if(i!=j) 2984 { 2985 p = pCopy(Gomega->m[i]); 2986 lm = pCopy(Gomega->m[j]); 2987 pIter(p); 2988 while(p!=NULL) 2989 { 2990 if(pDivisibleBy(lm,p)) 2991 { 2992 if(middle == FALSE) 2993 { 2994 middle = TRUE; 2995 } 2996 factor1 = singclap_pdivide(pHead(p),lm,currRing); 2997 factor2 = pMult(pCopy(factor1),pCopy(Go->m[j])); 2998 pDelete(&factor1); 2999 Go->m[i] = pAdd((Go->m[i]),pNeg(pCopy(factor2))); 3000 pDelete(&factor2); 3001 } 3002 pIter(p); 3003 } 3004 pDelete(&lm); 3005 pDelete(&p); 3006 } 3007 } 3008 } 3009 3010 //PrintS("\n//** jetzt Delete!\n"); 3011 //pDelete(&p); 3012 //pDelete(&factor); 3013 //pDelete(&lm); 3014 if(middle == TRUE) 3015 { 3016 //PrintS("\n//** middle TRUE!\n"); 3017 return Go; 3018 } 3019 //PrintS("\n//** middle FALSE!\n"); 3020 idDelete(&Go); 3021 return NULL; 3022 } 2976 3023 2977 3024 /****************************** Februar 2002 **************************** … … 3128 3175 else 3129 3176 { 3130 rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung3177 rChangeCurrRing(VMrDefault(curr_weight)); 3131 3178 } 3132 3179 newRing = currRing; … … 3280 3327 } 3281 3328 return 0; 3329 } 3330 3331 /***************************************** 3332 * return maximal polynomial length of G * 3333 *****************************************/ 3334 static int maxlengthpoly(ideal G) 3335 { 3336 int i,k,length=0; 3337 for(i=IDELEMS(G)-1; i>=0; i--) 3338 { 3339 k = pLength(G->m[i]); 3340 if(k>length) 3341 { 3342 length = k; 3343 } 3344 } 3345 return length; 3282 3346 } 3283 3347 … … 3884 3948 else 3885 3949 { 3886 rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung3950 rChangeCurrRing(VMrDefault(curr_weight)); 3887 3951 } 3888 3952 newRing = currRing; … … 4145 4209 else 4146 4210 { 4147 rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung4211 rChangeCurrRing(VMrDefault(curr_weight)); 4148 4212 } 4149 4213 newRing = currRing; … … 4287 4351 intvec* Xivlp; 4288 4352 4289 #if 0 4353 4290 4354 /******************************** 4291 4355 * compute a next weight vector * 4292 4356 ********************************/ 4293 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg) 4357 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, 4358 intvec* target_weight, int weight_rad, int pert_deg) 4294 4359 { 4295 int i, weight_norm; 4296 int nV = currRing->N; 4360 assume(currRing != NULL && curr_weight != NULL && 4361 target_weight != NULL && G->m[0] != NULL); 4362 4363 int i,weight_norm,nV = currRing->N; 4297 4364 intvec* next_weight2; 4298 4365 intvec* next_weight22 = new intvec(nV); 4299 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G); 4300 if(MivComp(next_weight, target_weight) == 1) 4301 { 4302 return(next_weight); 4303 } 4304 else 4305 { 4306 //compute a perturbed next weight vector "next_weight1" 4307 intvec* next_weight1 = MkInterRedNextWeight(MPertVectors(G, MivMatrixOrder(curr_weight), pert_deg), target_weight, G); 4308 //Print("\n // size of next_weight1 = %d", sizeof((*next_weight1))); 4309 4310 //compute a random next weight vector "next_weight2" 4311 while(1) 4312 { 4313 weight_norm = 0; 4314 while(weight_norm == 0) 4366 intvec* result = new intvec(nV); 4367 4368 intvec* next_weight1 =MkInterRedNextWeight(curr_weight,target_weight,G); 4369 //compute a random next weight vector "next_weight2" 4370 while(1) 4371 { 4372 weight_norm = 0; 4373 while(weight_norm == 0) 4374 { 4375 4376 for(i=0; i<nV; i++) 4377 { 4378 (*next_weight22)[i] = rand() % 60000 - 30000; 4379 weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i]; 4380 } 4381 weight_norm = 1 + floor(sqrt(weight_norm)); 4382 } 4383 4384 for(i=0; i<nV; i++) 4385 { 4386 if((*next_weight22)[i] < 0) 4387 { 4388 (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 4389 } 4390 else 4391 { 4392 (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 4393 } 4394 } 4395 4396 if(test_w_in_ConeCC(G, next_weight22) == 1) 4397 { 4398 next_weight2 = MkInterRedNextWeight(next_weight22,target_weight,G); 4399 if(MivAbsMax(next_weight2)>1147483647) 4315 4400 { 4316 4401 for(i=0; i<nV; i++) 4317 4402 { 4318 //Print("\n// next_weight[%d] = %d", i, (*next_weight)[i]); 4319 (*next_weight22)[i] = rand() % 60000 - 30000; 4320 weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i]; 4403 (*next_weight22)[i] = (*next_weight2)[i]; 4321 4404 } 4322 weight_norm = 1 + floor(sqrt(weight_norm)); 4323 } 4324 4325 for(i=nV-1; i>=0; i--) 4326 { 4327 if((*next_weight22)[i] < 0) 4405 i = 0; 4406 /* reduce the size of the maximal entry of the vector*/ 4407 while(test_w_in_ConeCC(G,next_weight22)) 4328 4408 { 4329 (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 4409 (*next_weight2)[i] = (*next_weight22)[i]; 4410 i = MivAbsMaxArg(next_weight22); 4411 (*next_weight22)[i] = floor(0.1*(*next_weight22)[i] + 0.5); 4330 4412 } 4331 else 4332 { 4333 (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 4334 } 4335 //Print("\n// next_weight22[%d] = %d", i, (*next_weight22)[i]); 4336 } 4337 4338 if(test_w_in_ConeCC(G, next_weight22) == 1) 4339 { 4340 //Print("\n//MWalkRandomNextWeight: next_weight2 im Kegel\n"); 4341 next_weight2 = MkInterRedNextWeight(next_weight22, target_weight, G); 4342 delete next_weight22; 4343 break; 4344 } 4345 } 4346 intvec* result = new intvec(nV); 4347 ideal G_test = MwalkInitialForm(G, next_weight); 4413 } 4414 delete next_weight22; 4415 break; 4416 } 4417 } 4418 4419 // compute "usual" next weight vector 4420 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G); 4421 ideal G_test = MwalkInitialForm(G, next_weight); 4422 ideal G_test2 = MwalkInitialForm(G, next_weight2); 4423 4424 // compare next weights 4425 if(Overflow_Error == FALSE) 4426 { 4348 4427 ideal G_test1 = MwalkInitialForm(G, next_weight1); 4349 ideal G_test2 = MwalkInitialForm(G, next_weight2); 4350 4351 // compare next_weights 4352 if(IDELEMS(G_test1) < IDELEMS(G_test)) 4353 { 4354 if(IDELEMS(G_test2) <= IDELEMS(G_test1)) // |G_test2| <= |G_test1| < |G_test| 4428 if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test))//if(IDELEMS(G_test1) < IDELEMS(G_test)) 4429 { 4430 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test1))//if(IDELEMS(G_test2) < IDELEMS(G_test1)) 4355 4431 { 4356 4432 for(i=0; i<nV; i++) … … 4359 4435 } 4360 4436 } 4361 else // |G_test1| < |G_test|, |G_test1| < |G_test2|4437 else 4362 4438 { 4363 4439 for(i=0; i<nV; i++) … … 4365 4441 (*result)[i] = (*next_weight1)[i]; 4366 4442 } 4367 } 4443 } 4368 4444 } 4369 4445 else 4370 4446 { 4371 if( IDELEMS(G_test2) <= IDELEMS(G_test)) // |G_test2| <=|G_test| <= |G_test1|4447 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))//if(IDELEMS(G_test2) < IDELEMS(G_test)) // |G_test2| < |G_test| <= |G_test1| 4372 4448 { 4373 4449 for(i=0; i<nV; i++) … … 4376 4452 } 4377 4453 } 4378 else // |G_test| <= |G_test1|, |G_test| < |G_test2|4454 else 4379 4455 { 4380 4456 for(i=0; i<nV; i++) … … 4384 4460 } 4385 4461 } 4386 delete next_weight;4387 delete next_weight1;4388 idDelete(&G_test);4389 4462 idDelete(&G_test1); 4390 idDelete(&G_test2); 4391 if(test_w_in_ConeCC(G, result) == 1) 4392 { 4393 delete next_weight2; 4394 return result; 4463 } 4464 else 4465 { 4466 Overflow_Error = FALSE; 4467 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))//if(IDELEMS(G_test2) < IDELEMS(G_test)) 4468 { 4469 for(i=1; i<nV; i++) 4470 { 4471 (*result)[i] = (*next_weight2)[i]; 4472 } 4395 4473 } 4396 4474 else 4397 4475 { 4398 delete result;4399 return next_weight2;4400 }4401 }4402 }4403 #endif4404 4405 /********************************4406 * compute a next weight vector *4407 ********************************/4408 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg)4409 {4410 int i, weight_norm;4411 //int randCount=0;4412 int nV = currRing->N;4413 intvec* next_weight2;4414 intvec* next_weight22 = new intvec(nV);4415 intvec* result = new intvec(nV);4416 4417 //compute a perturbed next weight vector "next_weight1"4418 //intvec* next_weight1 = MkInterRedNextWeight(MPertVectors(G,MivMatrixOrderRefine(curr_weight,target_weight),pert_deg),target_weight,G);4419 intvec* next_weight1 =MkInterRedNextWeight(curr_weight,target_weight,G);4420 //compute a random next weight vector "next_weight2"4421 while(1)4422 {4423 weight_norm = 0;4424 while(weight_norm == 0)4425 {4426 4476 for(i=0; i<nV; i++) 4427 4477 { 4428 (*next_weight22)[i] = rand() % 60000 - 30000;4429 weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i];4430 }4431 weight_norm = 1 + floor(sqrt(weight_norm));4432 }4433 for(i=0; i<nV; i++)4434 {4435 if((*next_weight22)[i] < 0)4436 {4437 (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);4438 }4439 else4440 {4441 (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);4442 }4443 }4444 if(test_w_in_ConeCC(G, next_weight22) == 1)4445 {4446 next_weight2 = MkInterRedNextWeight(next_weight22,target_weight,G);4447 delete next_weight22;4448 break;4449 }4450 }4451 // compute "usual" next weight vector4452 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);4453 ideal G_test = MwalkInitialForm(G, next_weight);4454 ideal G_test2 = MwalkInitialForm(G, next_weight2);4455 4456 // compare next weights4457 if(Overflow_Error == FALSE)4458 {4459 ideal G_test1 = MwalkInitialForm(G, next_weight1);4460 if(IDELEMS(G_test1) < IDELEMS(G_test))4461 {4462 if(IDELEMS(G_test2) < IDELEMS(G_test1))4463 {4464 // |G_test2| < |G_test1| < |G_test|4465 for(i=0; i<nV; i++)4466 {4467 (*result)[i] = (*next_weight2)[i];4468 }4469 }4470 else4471 {4472 // |G_test1| < |G_test|, |G_test1| <= |G_test2|4473 for(i=0; i<nV; i++)4474 {4475 (*result)[i] = (*next_weight1)[i];4476 }4477 }4478 }4479 else4480 {4481 if(IDELEMS(G_test2) < IDELEMS(G_test)) // |G_test2| < |G_test| <= |G_test1|4482 {4483 for(i=0; i<nV; i++)4484 {4485 (*result)[i] = (*next_weight2)[i];4486 }4487 }4488 else4489 {4490 // |G_test| < |G_test1|, |G_test| <= |G_test2|4491 for(i=0; i<nV; i++)4492 {4493 (*result)[i] = (*next_weight)[i];4494 }4495 }4496 }4497 idDelete(&G_test1);4498 }4499 else4500 {4501 Overflow_Error = FALSE;4502 if(IDELEMS(G_test2) < IDELEMS(G_test))4503 {4504 for(i=1; i<nV; i++)4505 {4506 (*result)[i] = (*next_weight2)[i];4507 }4508 }4509 else4510 {4511 for(i=0; i<nV; i++)4512 {4513 4478 (*result)[i] = (*next_weight)[i]; 4514 4479 } 4515 4480 } 4516 4481 } 4482 //PrintS("\n MWalkRandomNextWeight: Ende ok!\n"); 4517 4483 idDelete(&G_test); 4518 4484 idDelete(&G_test2); … … 4575 4541 else 4576 4542 { 4577 rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung4543 rChangeCurrRing(VMrDefault(orig_target_weight)); 4578 4544 } 4579 4545 TargetRing = currRing; … … 4646 4612 else 4647 4613 { 4648 rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung4614 rChangeCurrRing(VMrDefault(curr_weight)); 4649 4615 } 4650 4616 newRing = currRing; … … 4755 4721 else 4756 4722 { 4757 rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung4723 rChangeCurrRing(VMrDefault(orig_target_weight)); 4758 4724 } 4759 4725 F1 = idrMoveR(G, newRing,currRing); … … 4786 4752 else 4787 4753 { 4788 rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung4754 rChangeCurrRing(VMrDefault(orig_target_weight)); 4789 4755 } 4790 4756 KSTD_Finish: … … 4884 4850 tim = clock(); 4885 4851 /* 4886 Print("\n// **** Gr ᅵbnerwalk took %d steps and ", nwalk);4852 Print("\n// **** Groebnerwalk took %d steps and ", nwalk); 4887 4853 PrintS("\n// **** call the rec. Pert. Walk to compute a red GB of:"); 4888 4854 idElements(Gomega, "G_omega"); … … 4914 4880 oldRing = currRing; 4915 4881 4916 / * create a new ring newRing */4882 // create a new ring newRing 4917 4883 if (rParameter(currRing) != NULL) 4918 4884 { … … 4921 4887 else 4922 4888 { 4923 rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung4889 rChangeCurrRing(VMrDefault(curr_weight)); 4924 4890 } 4925 4891 newRing = currRing; … … 4947 4913 else 4948 4914 { 4949 rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung4915 rChangeCurrRing(VMrDefault(curr_weight)); 4950 4916 } 4951 4917 newRing = currRing; … … 4959 4925 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight); 4960 4926 delete hilb_func; 4961 #endif // BUCHBERGER_ALG4927 #endif 4962 4928 tstd = tstd + clock() - to; 4963 4929 … … 4968 4934 4969 4935 to = clock(); 4970 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega). Gomega is a reduced Groebner basis w.r.t. the current ring. 4936 // compute a representation of the generators of submod (M) with respect 4937 // to those of mod (Gomega). 4938 // Gomega is a reduced Groebner basis w.r.t. the current ring. 4971 4939 F = MLifttwoIdeal(Gomega2, M1, G); 4972 4940 tlift = tlift + clock() - to; … … 5018 4986 else 5019 4987 { 5020 rChangeCurrRing(VMrDefault(target_weight)); // Aenderung4988 rChangeCurrRing(VMrDefault(target_weight)); 5021 4989 } 5022 4990 F1 = idrMoveR(G, newRing,currRing); … … 5065 5033 * THE GROEBNER WALK ALGORITHM * 5066 5034 *******************************/ 5067 ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, ring baseRing) 5035 ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, 5036 ring baseRing, int reduction, int printout) 5068 5037 { 5069 BITSET save1 = si_opt_1; // save current options 5070 //si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 5071 //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5072 //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB)); 5038 // save current options 5039 BITSET save1 = si_opt_1; 5040 if(reduction == 0) 5041 { 5042 si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 5043 si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5044 } 5073 5045 Set_Error(FALSE); 5074 5046 Overflow_Error = FALSE; 5047 //BOOLEAN endwalks = FALSE; 5075 5048 #ifdef TIME_TEST 5076 5049 clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0; … … 5080 5053 #endif 5081 5054 nstep=0; 5082 int i,nwalk ,endwalks = 0;5055 int i,nwalk; 5083 5056 int nV = baseRing->N; 5084 5057 5085 ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;5058 ideal Gomega, M, F, FF, Gomega1, Gomega2, M1; 5086 5059 ring newRing; 5087 5060 ring XXRing = baseRing; 5061 ring targetRing; 5088 5062 intvec* ivNull = new intvec(nV); 5089 5063 intvec* curr_weight = new intvec(nV); 5090 5064 intvec* target_weight = new intvec(nV); 5091 5065 intvec* exivlp = Mivlp(nV); 5066 /* 5092 5067 intvec* tmp_weight = new intvec(nV); 5093 5068 for(i=0; i<nV; i++) 5094 5069 { 5095 (*tmp_weight)[i] = (*target_M)[i]; 5096 } 5070 (*tmp_weight)[i] = (*orig_M)[i]; 5071 } 5072 */ 5097 5073 for(i=0; i<nV; i++) 5098 5074 { … … 5111 5087 #endif 5112 5088 rComplete(currRing); 5113 #ifdef CHECK_IDEAL_MWALK 5114 idString(Go,"Go"); 5115 #endif 5089 //#ifdef CHECK_IDEAL_MWALK 5090 if(printout > 2) 5091 { 5092 idString(Go,"//** Mwalk: Go"); 5093 } 5094 //#endif 5095 5096 if(target_M->length() == nV) 5097 { 5098 // define the target ring 5099 targetRing = VMrDefault(target_weight); 5100 } 5101 else 5102 { 5103 targetRing = VMatrDefault(target_M); 5104 } 5105 if(orig_M->length() == nV) 5106 { 5107 // define a new ring with ordering "(a(curr_weight),lp) 5108 newRing = VMrDefault(curr_weight); 5109 } 5110 else 5111 { 5112 newRing = VMatrDefault(orig_M); 5113 } 5114 rChangeCurrRing(newRing); 5115 5116 5116 #ifdef TIME_TEST 5117 5117 to = clock(); 5118 5118 #endif 5119 if(orig_M->length() == nV) 5120 { 5121 newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp) 5122 } 5123 else 5124 { 5125 newRing = VMatrDefault(orig_M); 5126 } 5127 rChangeCurrRing(newRing); 5119 5128 5120 ideal G = MstdCC(idrMoveR(Go,baseRing,currRing)); 5129 baseRing = currRing; 5121 5130 5122 #ifdef TIME_TEST 5131 5123 tostd = clock()-to; 5132 5124 #endif 5133 5125 5126 baseRing = currRing; 5134 5127 nwalk = 0; 5128 5135 5129 while(1) 5136 5130 { 5137 5131 nwalk ++; 5138 5132 nstep ++; 5133 5139 5134 #ifdef TIME_TEST 5140 5135 to = clock(); 5141 5136 #endif 5142 #ifdef CHECK_IDEAL_MWALK 5143 idString(G,"G"); 5144 #endif 5145 Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector" 5137 // compute an initial form ideal of <G> w.r.t. "curr_vector" 5138 Gomega = MwalkInitialForm(G, curr_weight); 5146 5139 #ifdef TIME_TEST 5147 tif = tif + clock()-to; //time for computing initial form ideal 5148 #endif 5149 #ifdef CHECK_IDEAL_MWALK 5150 idString(Gomega,"Gomega"); 5151 #endif 5140 tif = tif + clock()-to; 5141 #endif 5142 5143 //#ifdef CHECK_IDEAL_MWALK 5144 if(printout > 1) 5145 { 5146 idString(Gomega,"//** Mwalk: Gomega"); 5147 } 5148 //#endif 5149 5150 if(reduction == 0) 5151 { 5152 FF = middleOfCone(G,Gomega); 5153 if( FF != NULL) 5154 { 5155 idDelete(&G); 5156 G = idCopy(FF); 5157 idDelete(&FF); 5158 goto NEXT_VECTOR; 5159 } 5160 } 5161 5152 5162 #ifndef BUCHBERGER_ALG 5153 5163 if(isNolVector(curr_weight) == 0) 5154 5164 { 5155 5165 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 5156 } 5166 } 5157 5167 else 5158 5168 { … … 5160 5170 } 5161 5171 #endif 5172 5162 5173 if(nwalk == 1) 5163 5174 { 5164 5175 if(orig_M->length() == nV) 5165 5176 { 5166 newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp) 5177 // define a new ring with ordering "(a(curr_weight),lp) 5178 newRing = VMrDefault(curr_weight); 5167 5179 } 5168 5180 else … … 5175 5187 if(target_M->length() == nV) 5176 5188 { 5177 newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))" 5189 //define a new ring with ordering "(a(curr_weight),lp)" 5190 newRing = VMrDefault(curr_weight); 5178 5191 } 5179 5192 else 5180 5193 { 5194 //define a new ring with matrix ordering 5181 5195 newRing = VMatrRefine(target_M,curr_weight); 5182 5196 } … … 5199 5213 #endif 5200 5214 idSkipZeroes(M); 5201 #ifdef CHECK_IDEAL_MWALK 5202 PrintS("\n//** Mwalk: computed M.\n"); 5203 idString(M, "M"); 5204 #endif 5215 //#ifdef CHECK_IDEAL_MWALK 5216 if(printout > 2) 5217 { 5218 idString(M, "//** Mwalk: M"); 5219 } 5220 //#endif 5205 5221 //change the ring to baseRing 5206 5222 rChangeCurrRing(baseRing); … … 5212 5228 to = clock(); 5213 5229 #endif 5214 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega), where Gomega is a reduced Groebner basis w.r.t. the current ring 5230 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega), 5231 // where Gomega is a reduced Groebner basis w.r.t. the current ring 5215 5232 F = MLifttwoIdeal(Gomega2, M1, G); 5233 5216 5234 #ifdef TIME_TEST 5217 5235 tlift = tlift + clock() - to; 5218 5236 #endif 5219 #ifdef CHECK_IDEAL_MWALK 5220 idString(F, "F"); 5221 #endif 5237 //#ifdef CHECK_IDEAL_MWALK 5238 if(printout > 2) 5239 { 5240 idString(F, "//** Mwalk: F"); 5241 } 5242 //#endif 5222 5243 idDelete(&Gomega2); 5223 5244 idDelete(&M1); 5245 5224 5246 rChangeCurrRing(newRing); // change the ring to newRing 5225 5247 G = idrMoveR(F,baseRing,currRing); 5226 5248 idDelete(&F); 5249 idSkipZeroes(G); 5250 5251 //#ifdef CHECK_IDEAL_MWALK 5252 if(printout > 2) 5253 { 5254 idString(G, "//** Mwalk: G"); 5255 } 5256 //#endif 5257 5258 rChangeCurrRing(targetRing); 5259 G = idrMoveR(G,newRing,currRing); 5260 // test whether target cone is reached 5261 if(reduction !=0 && test_w_in_ConeCC(G,curr_weight) == 1) 5262 { 5263 baseRing = currRing; 5264 break; 5265 //endwalks = TRUE; 5266 } 5267 5268 rChangeCurrRing(newRing); 5269 G = idrMoveR(G,targetRing,currRing); 5227 5270 baseRing = currRing; 5271 5272 /* 5228 5273 #ifdef TIME_TEST 5229 5274 to = clock(); 5230 5275 #endif 5231 //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL); 5276 5232 5277 #ifdef TIME_TEST 5233 5278 tstd = tstd + clock() - to; 5234 5279 #endif 5235 idSkipZeroes(G); 5236 #ifdef CHECK_IDEAL_MWALK 5237 idString(G, "G"); 5238 #endif 5280 */ 5281 5282 5239 5283 #ifdef TIME_TEST 5240 5284 to = clock(); 5241 5285 #endif 5286 NEXT_VECTOR: 5242 5287 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 5243 5288 #ifdef TIME_TEST 5244 5289 tnw = tnw + clock() - to; 5245 5290 #endif 5246 #ifdef PRINT_VECTORS 5247 MivString(curr_weight, target_weight, next_weight); 5248 #endif 5249 if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1)// || test_w_in_ConeCC(G, target_weight) == 1 || MivComp(next_weight,curr_weight) == 1) 5250 { 5251 #ifdef CHECK_IDEAL_MWALK 5252 PrintS("\n//** Mwalk: entering last cone.\n"); 5253 #endif 5291 //#ifdef PRINT_VECTORS 5292 if(printout > 0) 5293 { 5294 MivString(curr_weight, target_weight, next_weight); 5295 } 5296 //#endif 5297 if(MivComp(target_weight,curr_weight) == 1)// || endwalks == TRUE) 5298 {/* 5299 //#ifdef CHECK_IDEAL_MWALK 5300 if(printout > 0) 5301 { 5302 PrintS("\n//** Mwalk: entering last cone.\n"); 5303 } 5304 //#endif 5305 5254 5306 Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector" 5255 5307 if(target_M->length() == nV) … … 5264 5316 Gomega1 = idrMoveR(Gomega, baseRing,currRing); 5265 5317 idDelete(&Gomega); 5266 #ifdef CHECK_IDEAL_MWALK 5267 idString(Gomega1, "Gomega"); 5268 #endif 5318 //#ifdef CHECK_IDEAL_MWALK 5319 if(printout > 1) 5320 { 5321 idString(Gomega1, "//** Mwalk: Gomega"); 5322 } 5323 PrintS("\n //** Mwalk: kStd(Gomega)"); 5324 //#endif 5269 5325 M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL); 5270 #ifdef CHECK_IDEAL_MWALK 5271 idString(M,"M"); 5272 #endif 5326 //#ifdef CHECK_IDEAL_MWALK 5327 if(printout > 1) 5328 { 5329 idString(M,"//** Mwalk: M"); 5330 } 5331 //#endif 5273 5332 rChangeCurrRing(baseRing); 5274 5333 M1 = idrMoveR(M, newRing,currRing); … … 5276 5335 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 5277 5336 idDelete(&Gomega1); 5337 //PrintS("\n //** Mwalk: MLifttwoIdeal"); 5278 5338 F = MLifttwoIdeal(Gomega2, M1, G); 5279 #ifdef CHECK_IDEAL_MWALK 5280 idString(F,"F"); 5281 #endif 5339 //#ifdef CHECK_IDEAL_MWALK 5340 if(printout > 2) 5341 { 5342 idString(F,"//** Mwalk: F"); 5343 } 5344 //#endif 5282 5345 idDelete(&Gomega2); 5283 5346 idDelete(&M1); … … 5291 5354 to = clock(); 5292 5355 #endif 5293 // if(si_opt_1 == (Sy_bit(OPT_REDSB))) 5294 // { 5295 G = kInterRedCC(G,NULL); //reduce the Groebner basis <G> w.r.t. currRing, if option(redSB) is set 5296 // } 5356 //PrintS("\n //**Mwalk: Interreduce"); 5357 //interreduce the Groebner basis <G> w.r.t. currRing 5358 //G = kInterRedCC(G,NULL); 5297 5359 #ifdef TIME_TEST 5298 5360 tred = tred + clock() - to; 5299 5361 #endif 5300 5362 idSkipZeroes(G); 5301 delete next_weight; 5363 delete next_weight; */ 5302 5364 break; 5303 #ifdef CHECK_IDEAL_MWALK 5304 PrintS("\n//** Mwalk: last cone.\n"); 5305 #endif 5306 } 5307 #ifdef CHECK_IDEAL_MWALK 5308 PrintS("\n//** Mwalk: update weight vectors.\n"); 5309 #endif 5365 } 5366 5310 5367 for(i=nV-1; i>=0; i--) 5311 5368 { 5312 (*tmp_weight)[i] = (*curr_weight)[i];5369 //(*tmp_weight)[i] = (*curr_weight)[i]; 5313 5370 (*curr_weight)[i] = (*next_weight)[i]; 5314 5371 } … … 5317 5374 rChangeCurrRing(XXRing); 5318 5375 ideal result = idrMoveR(G,baseRing,currRing); 5376 idDelete(&Go); 5319 5377 idDelete(&G); 5320 /*#ifdef CHECK_IDEAL_MWALK 5321 pDelete(&p); 5322 #endif*/ 5323 delete tmp_weight; 5378 //delete tmp_weight; 5324 5379 delete ivNull; 5325 5380 delete exivlp; … … 5328 5383 #endif 5329 5384 #ifdef TIME_TEST 5330 Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);5331 5385 TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep); 5332 Print("\n//** Mwalk: Ergebnis.\n");5333 5386 //Print("\n// pSetm_Error = (%d)", ErrorCheck()); 5334 5387 //Print("\n// Overflow_Error? (%d)\n", Overflow_Error); 5335 5388 #endif 5389 Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep); 5336 5390 return(result); 5337 5391 } 5338 5392 5339 // 07.11.20125340 // THE RANDOM WALK ALGORITHM ideal Go, intvec* orig_M, intvec* target_M, ring baseRing 5341 ideal Mrwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int pert_deg, ring baseRing)5393 // THE RANDOM WALK ALGORITHM 5394 ideal Mrwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int pert_deg, 5395 int reduction, int printout) 5342 5396 { 5343 5397 BITSET save1 = si_opt_1; // save current options 5344 //si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 5345 //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5346 //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB)); 5398 if(reduction == 0) 5399 { 5400 si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 5401 si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5402 } 5347 5403 Set_Error(FALSE); 5348 5404 Overflow_Error = FALSE; 5405 //BOOLEAN endwalks = FALSE; 5349 5406 #ifdef TIME_TEST 5350 5407 clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0; … … 5354 5411 #endif 5355 5412 nstep=0; 5356 int i, nwalk,endwalks = 0;5357 int nV = baseRing->N;5358 5359 ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;5413 int i,polylength,nwalk; 5414 int nV = currRing->N; 5415 5416 ideal Gomega, M, F,FF, Gomega1, Gomega2, M1; 5360 5417 ring newRing; 5361 ring XXRing = baseRing; 5418 ring targetRing; 5419 ring baseRing = currRing; 5420 ring XXRing = currRing; 5362 5421 intvec* ivNull = new intvec(nV); 5363 5422 intvec* curr_weight = new intvec(nV); 5364 5423 intvec* target_weight = new intvec(nV); 5365 5424 intvec* exivlp = Mivlp(nV); 5425 intvec* next_weight= new intvec(nV); 5426 /* 5366 5427 intvec* tmp_weight = new intvec(nV); 5367 5428 for(i=0; i<nV; i++) … … 5369 5430 (*tmp_weight)[i] = (*target_M)[i]; 5370 5431 } 5432 */ 5371 5433 for(i=0; i<nV; i++) 5372 5434 { … … 5385 5447 #endif 5386 5448 rComplete(currRing); 5387 #ifdef CHECK_IDEAL_MWALK5388 idString(Go,"Go");5389 #endif5390 5449 #ifdef TIME_TEST 5391 5450 to = clock(); 5392 5451 #endif 5393 if(orig_M->length() == nV) 5394 { 5395 newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp) 5396 } 5397 else 5398 { 5399 newRing = VMatrDefault(orig_M); 5400 } 5452 5453 if(target_M->length() == nV) 5454 { 5455 // define the target ring 5456 targetRing = VMrDefault(target_weight); 5457 } 5458 else 5459 { 5460 targetRing = VMatrDefault(target_M); 5461 } 5462 if(orig_M->length() == nV) 5463 { 5464 newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp) 5465 } 5466 else 5467 { 5468 newRing = VMatrDefault(orig_M); 5469 } 5401 5470 rChangeCurrRing(newRing); 5402 5471 ideal G = MstdCC(idrMoveR(Go,baseRing,currRing)); … … 5414 5483 to = clock(); 5415 5484 #endif 5416 #ifdef CHECK_IDEAL_MWALK 5417 idString(G,"G"); 5418 #endif 5485 5419 5486 Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector" 5487 //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise 5488 polylength = lengthpoly(Gomega); 5420 5489 #ifdef TIME_TEST 5421 5490 tif = tif + clock()-to; //time for computing initial form ideal 5422 5491 #endif 5423 #ifdef CHECK_IDEAL_MWALK 5424 idString(Gomega,"Gomega"); 5425 #endif 5492 //#ifdef CHECK_IDEAL_MWALK 5493 if(printout > 1) 5494 { 5495 idString(Gomega,"//** Mrwalk: Gomega"); 5496 } 5497 //#endif 5498 // test whether target cone is reached 5499 /* if(test_w_in_ConeCC(G,target_weight) == 1) 5500 { 5501 endwalks = TRUE; 5502 }*/ 5503 if(reduction == 0) 5504 { 5505 5506 FF = middleOfCone(G,Gomega); 5507 5508 if(FF != NULL) 5509 { 5510 idDelete(&G); 5511 G = idCopy(FF); 5512 idDelete(&FF); 5513 5514 goto NEXT_VECTOR; 5515 } 5516 } 5426 5517 #ifndef BUCHBERGER_ALG 5427 5518 if(isNolVector(curr_weight) == 0) 5428 5519 { 5429 5520 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 5430 } 5521 } 5431 5522 else 5432 5523 { … … 5449 5540 if(target_M->length() == nV) 5450 5541 { 5451 newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))" 5542 newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp) 5543 //newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))" 5452 5544 } 5453 5545 else … … 5473 5565 #endif 5474 5566 idSkipZeroes(M); 5475 #ifdef CHECK_IDEAL_MWALK 5476 PrintS("\n//** Mwalk: computed M.\n"); 5477 idString(M, "M"); 5478 #endif 5567 //#ifdef CHECK_IDEAL_MWALK 5568 if(printout > 2) 5569 { 5570 idString(M, "//** Mrwalk: M"); 5571 } 5572 //#endif 5479 5573 //change the ring to baseRing 5480 5574 rChangeCurrRing(baseRing); … … 5486 5580 to = clock(); 5487 5581 #endif 5488 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega), where Gomega is a reduced Groebner basis w.r.t. the current ring 5582 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega), 5583 // where Gomega is a reduced Groebner basis w.r.t. the current ring 5489 5584 F = MLifttwoIdeal(Gomega2, M1, G); 5490 5585 #ifdef TIME_TEST 5491 5586 tlift = tlift + clock() - to; 5492 5587 #endif 5493 #ifdef CHECK_IDEAL_MWALK 5494 idString(F, "F"); 5495 #endif 5588 //#ifdef CHECK_IDEAL_MWALK 5589 if(printout > 2) 5590 { 5591 idString(F, "//** Mrwalk: F"); 5592 } 5593 //#endif 5496 5594 idDelete(&Gomega2); 5497 5595 idDelete(&M1); … … 5502 5600 #ifdef TIME_TEST 5503 5601 to = clock(); 5504 #endif5505 //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL);5506 #ifdef TIME_TEST5507 5602 tstd = tstd + clock() - to; 5508 5603 #endif 5509 5604 idSkipZeroes(G); 5510 #ifdef CHECK_IDEAL_MWALK 5511 idString(G, "G"); 5512 #endif 5605 //#ifdef CHECK_IDEAL_MWALK 5606 if(printout > 2) 5607 { 5608 idString(G, "//** Mrwalk: G"); 5609 } 5610 //#endif 5611 5612 rChangeCurrRing(targetRing); 5613 G = idrMoveR(G,newRing,currRing); 5614 // test whether target cone is reached 5615 if(reduction !=0 && test_w_in_ConeCC(G,curr_weight) == 1) 5616 { 5617 baseRing = currRing; 5618 break; 5619 } 5620 5621 rChangeCurrRing(newRing); 5622 G = idrMoveR(G,targetRing,currRing); 5623 baseRing = currRing; 5624 5625 5626 NEXT_VECTOR: 5513 5627 #ifdef TIME_TEST 5514 5628 to = clock(); 5515 5629 #endif 5516 intvec* next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg);//next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 5630 next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 5631 if(polylength > 0) 5632 { 5633 //there is a polynomial in Gomega with at least 3 monomials, 5634 //low-dimensional facet of the cone 5635 delete next_weight; 5636 next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg); 5637 } 5517 5638 #ifdef TIME_TEST 5518 5639 tnw = tnw + clock() - to; 5519 5640 #endif 5520 #ifdef PRINT_VECTORS 5521 MivString(curr_weight, target_weight, next_weight); 5522 #endif 5523 if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1)// || test_w_in_ConeCC(G, target_weight) == 1 || MivComp(next_weight,curr_weight) == 1) 5524 { 5525 #ifdef CHECK_IDEAL_MWALK 5526 PrintS("\n//** Mwalk: entering last cone.\n"); 5527 #endif 5641 //#ifdef PRINT_VECTORS 5642 if(printout > 0) 5643 { 5644 MivString(curr_weight, target_weight, next_weight); 5645 } 5646 //#endif 5647 if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1)// || endwalks == TRUE) 5648 {/* 5649 //#ifdef CHECK_IDEAL_MWALK 5650 if(printout > 0) 5651 { 5652 PrintS("\n//** Mrwalk: entering last cone.\n"); 5653 } 5654 //#endif 5655 5528 5656 Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector" 5529 5657 if(target_M->length() == nV) … … 5538 5666 Gomega1 = idrMoveR(Gomega, baseRing,currRing); 5539 5667 idDelete(&Gomega); 5540 #ifdef CHECK_IDEAL_MWALK 5541 idString(Gomega1, "Gomega"); 5542 #endif 5668 //#ifdef CHECK_IDEAL_MWALK 5669 if(printout > 1) 5670 { 5671 idString(Gomega1, "//** Mrwalk: Gomega"); 5672 } 5673 PrintS("\n //** Mrwalk: kStd(Gomega)"); 5674 //#endif 5543 5675 M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL); 5544 #ifdef CHECK_IDEAL_MWALK 5545 idString(M,"M"); 5546 #endif 5676 //#ifdef CHECK_IDEAL_MWALK 5677 if(printout > 1) 5678 { 5679 idString(M,"//** Mrwalk: M"); 5680 } 5681 //#endif 5547 5682 rChangeCurrRing(baseRing); 5548 5683 M1 = idrMoveR(M, newRing,currRing); … … 5550 5685 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 5551 5686 idDelete(&Gomega1); 5687 //PrintS("\n //** Mrwalk: MLifttwoIdeal"); 5552 5688 F = MLifttwoIdeal(Gomega2, M1, G); 5553 #ifdef CHECK_IDEAL_MWALK 5554 idString(F,"F"); 5555 #endif 5689 //#ifdef CHECK_IDEAL_MWALK 5690 if(printout > 2) 5691 { 5692 idString(F,"//** Mrwalk: F"); 5693 } 5694 //#endif 5556 5695 idDelete(&Gomega2); 5557 5696 idDelete(&M1); … … 5565 5704 to = clock(); 5566 5705 #endif 5567 // if(si_opt_1 == (Sy_bit(OPT_REDSB))) 5568 // { 5569 //G = kInterRedCC(G,NULL); //reduce the Groebner basis <G> w.r.t. currRing, if option(redSB) is set 5570 // } 5706 //PrintS("\n //**Mrwalk: Interreduce"); 5707 //interreduce the Groebner basis <G> w.r.t. currRing 5708 //G = kInterRedCC(G,NULL); 5571 5709 #ifdef TIME_TEST 5572 5710 tred = tred + clock() - to; 5573 5711 #endif 5574 5712 idSkipZeroes(G); 5575 delete next_weight; 5713 delete next_weight;*/ 5576 5714 break; 5577 #ifdef CHECK_IDEAL_MWALK 5578 PrintS("\n//** Mwalk: last cone.\n"); 5579 #endif 5580 } 5581 #ifdef CHECK_IDEAL_MWALK 5582 PrintS("\n//** Mwalk: update weight vectors.\n"); 5583 #endif 5715 } 5716 5584 5717 for(i=nV-1; i>=0; i--) 5585 5718 { 5586 (*tmp_weight)[i] = (*curr_weight)[i];5719 //(*tmp_weight)[i] = (*curr_weight)[i]; 5587 5720 (*curr_weight)[i] = (*next_weight)[i]; 5588 5721 } 5589 5722 delete next_weight; 5590 5723 } 5724 baseRing = currRing; 5591 5725 rChangeCurrRing(XXRing); 5592 5726 ideal result = idrMoveR(G,baseRing,currRing); 5593 5727 idDelete(&G); 5594 /*#ifdef CHECK_IDEAL_MWALK 5595 pDelete(&p); 5596 #endif*/ 5597 delete tmp_weight; 5728 si_opt_1 = save1; //set original options, e. g. option(RedSB) 5729 //delete tmp_weight; 5598 5730 delete ivNull; 5599 5731 delete exivlp; … … 5601 5733 delete last_omega; 5602 5734 #endif 5735 Print("\n//** Mrwalk: Groebner Walk took %d steps.\n", nstep); 5603 5736 #ifdef TIME_TEST 5604 Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);5605 5737 TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep); 5606 Print("\n//** Mwalk: Ergebnis.\n");5607 5738 //Print("\n// pSetm_Error = (%d)", ErrorCheck()); 5608 5739 //Print("\n// Overflow_Error? (%d)\n", Overflow_Error); … … 5751 5882 // use kStd, if nP = 0, else call LastGB 5752 5883 ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight, 5753 intvec* target_weight, int nP )5884 intvec* target_weight, int nP, int reduction, int printout) 5754 5885 { 5886 BITSET save1 = si_opt_1; // save current options 5887 if(reduction == 0) 5888 { 5889 si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 5890 //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5891 //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB)); 5892 } 5755 5893 Set_Error(FALSE ); 5756 5894 Overflow_Error = FALSE; … … 5766 5904 nstep = 0; 5767 5905 int i, ntwC=1, ntestw=1, nV = currRing->N; 5768 int endwalks=0;5769 5770 ideal Gomega, M, F, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;5906 BOOLEAN endwalks = FALSE; 5907 5908 ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG; 5771 5909 ring newRing, oldRing, TargetRing; 5772 5910 intvec* iv_M_dp; … … 5790 5928 ring XXRing = currRing; 5791 5929 5792 5793 5930 to = clock(); 5794 / * perturbs the original vector */5931 // perturbs the original vector 5795 5932 if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp" 5796 5933 { … … 5809 5946 DefRingPar(curr_weight); 5810 5947 else 5811 rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 15948 rChangeCurrRing(VMrDefault(curr_weight)); 5812 5949 5813 5950 G = idrMoveR(Go, XXRing,currRing); … … 5824 5961 ring HelpRing = currRing; 5825 5962 5826 / * perturbs the target weight vector */5963 // perturbs the target weight vector 5827 5964 if(tp_deg > 1 && tp_deg <= nV) 5828 5965 { … … 5830 5967 DefRingPar(target_weight); 5831 5968 else 5832 rChangeCurrRing(VMrDefault(target_weight)); // Aenderung 25969 rChangeCurrRing(VMrDefault(target_weight)); 5833 5970 5834 5971 TargetRing = currRing; … … 5837 5974 { 5838 5975 iv_M_lp = MivMatrixOrderlp(nV); 5839 //ivString(iv_M_lp, "iv_M_lp");5840 //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);5841 5976 target_weight = MPertVectors(ssG, iv_M_lp, tp_deg); 5842 5977 } … … 5844 5979 { 5845 5980 iv_M_lp = MivMatrixOrder(target_weight); 5846 //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);5847 5981 target_weight = MPertVectors(ssG, iv_M_lp, tp_deg); 5848 5982 } … … 5852 5986 G = idrMoveR(ssG, TargetRing,currRing); 5853 5987 } 5854 /* 5855 Print("\n// Perturbationwalkalg. vom Gradpaar (%d,%d):",op_deg,tp_deg); 5856 ivString(curr_weight, "new sigma"); 5857 ivString(target_weight, "new tau"); 5858 */ 5988 if(printout > 0) 5989 { 5990 Print("\n//** Mpwalk: Perturbation Walk of degree (%d,%d):",op_deg,tp_deg); 5991 ivString(curr_weight, "//** Mpwalk: new current weight"); 5992 ivString(target_weight, "//** Mpwalk: new target weight"); 5993 } 5859 5994 while(1) 5860 5995 { 5861 5996 nstep ++; 5862 5997 to = clock(); 5863 / *compute an initial form ideal of <G> w.r.t. the weight vector5864 "curr_weight" */5998 // compute an initial form ideal of <G> w.r.t. the weight vector 5999 // "curr_weight" 5865 6000 Gomega = MwalkInitialForm(G, curr_weight); 5866 6001 //#ifdef CHECK_IDEAL_MWALK 6002 if(printout > 1) 6003 { 6004 idString(Gomega,"//** Mpwalk: Gomega"); 6005 } 6006 //#endif 6007 if(reduction == 0 && nstep > 1) 6008 { 6009 // check whether weight vector is in the interior of the cone 6010 while(1) 6011 { 6012 FF = middleOfCone(G,Gomega); 6013 if(FF != NULL) 6014 { 6015 idDelete(&G); 6016 G = idCopy(FF); 6017 idDelete(&FF); 6018 next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 6019 if(printout > 0) 6020 { 6021 MivString(curr_weight, target_weight, next_weight); 6022 } 6023 } 6024 else 6025 { 6026 break; 6027 } 6028 for(i=nV-1; i>=0; i--) 6029 { 6030 (*curr_weight)[i] = (*next_weight)[i]; 6031 } 6032 Gomega = MwalkInitialForm(G, curr_weight); 6033 if(printout > 1) 6034 { 6035 idString(Gomega,"//** Mpwalk: Gomega"); 6036 } 6037 } 6038 } 5867 6039 5868 6040 #ifdef ENDWALKS 5869 if(endwalks == 1){ 6041 if(endwalks == TRUE) 6042 { 5870 6043 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 5871 6044 idElements(G, "G"); 5872 // idElements(Gomega, "Gw");5873 6045 headidString(G, "G"); 5874 //headidString(Gomega, "Gw");5875 6046 } 5876 6047 #endif … … 5891 6062 DefRingPar(curr_weight); 5892 6063 else 5893 rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 36064 rChangeCurrRing(VMrDefault(curr_weight)); 5894 6065 5895 6066 newRing = currRing; … … 5897 6068 5898 6069 #ifdef ENDWALKS 5899 if(endwalks== 1)6070 if(endwalks==TRUE) 5900 6071 { 5901 6072 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); … … 5912 6083 tim = clock(); 5913 6084 to = clock(); 5914 / * compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */6085 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 5915 6086 #ifdef BUCHBERGER_ALG 5916 6087 M = MstdhomCC(Gomega1); … … 5918 6089 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight); 5919 6090 delete hilb_func; 5920 #endif // BUCHBERGER_ALG 5921 5922 if(endwalks == 1){ 6091 #endif 6092 //#ifdef CHECK_IDEAL_MWALK 6093 if(printout > 2) 6094 { 6095 idString(M,"//** Mpwalk: M"); 6096 } 6097 //#endif 6098 6099 if(endwalks == TRUE){ 5923 6100 xtstd = xtstd+clock()-to; 5924 6101 #ifdef ENDWALKS … … 5930 6107 tstd=tstd+clock()-to; 5931 6108 5932 / * change the ring to oldRing */6109 // change the ring to oldRing 5933 6110 rChangeCurrRing(oldRing); 5934 6111 M1 = idrMoveR(M, newRing,currRing); 5935 6112 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 5936 5937 //if(endwalks==1) PrintS("\n// Lifting is working:..");5938 6113 5939 6114 to=clock(); … … 5942 6117 Gomega is a reduced Groebner basis w.r.t. the current ring */ 5943 6118 F = MLifttwoIdeal(Gomega2, M1, G); 5944 if(endwalks != 1)6119 if(endwalks == FALSE) 5945 6120 tlift = tlift+clock()-to; 5946 6121 else 5947 6122 xtlift=clock()-to; 5948 6123 6124 //#ifdef CHECK_IDEAL_MWALK 6125 if(printout > 2) 6126 { 6127 idString(F,"//** Mpwalk: F"); 6128 } 6129 //#endif 6130 5949 6131 idDelete(&M1); 5950 6132 idDelete(&Gomega2); 5951 6133 idDelete(&G); 5952 6134 5953 / * change the ring to newRing */6135 // change the ring to newRing 5954 6136 rChangeCurrRing(newRing); 5955 F1 = idrMoveR(F, oldRing,currRing); 5956 5957 //if(endwalks==1)PrintS("\n// InterRed is working now:"); 6137 if(reduction == 0) 6138 { 6139 G = idrMoveR(F,oldRing,currRing); 6140 } 6141 else 6142 { 6143 F1 = idrMoveR(F, oldRing,currRing); 6144 if(printout > 2) 6145 { 6146 PrintS("\n //** Mpwalk: reduce the Groebner basis.\n"); 6147 } 6148 to=clock(); 6149 G = kInterRedCC(F1, NULL); 6150 if(endwalks == FALSE) 6151 tred = tred+clock()-to; 6152 else 6153 xtred=clock()-to; 6154 idDelete(&F1); 6155 } 6156 if(endwalks == TRUE) 6157 break; 5958 6158 5959 6159 to=clock(); 5960 /* reduce the Groebner basis <G> w.r.t. new ring */ 5961 G = kInterRedCC(F1, NULL); 5962 if(endwalks != 1) 5963 tred = tred+clock()-to; 5964 else 5965 xtred=clock()-to; 5966 5967 idDelete(&F1); 5968 5969 if(endwalks == 1) 5970 break; 5971 5972 to=clock(); 5973 /* compute a next weight vector */ 6160 // compute a next weight vector 5974 6161 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 5975 6162 tnw=tnw+clock()-to; 5976 #ifdef PRINT_VECTORS 5977 MivString(curr_weight, target_weight, next_weight); 5978 #endif 6163 //#ifdef PRINT_VECTORS 6164 if(printout > 0) 6165 { 6166 MivString(curr_weight, target_weight, next_weight); 6167 } 6168 //#endif 5979 6169 5980 6170 if(Overflow_Error == TRUE) … … 5994 6184 } 5995 6185 if(MivComp(next_weight, target_weight) == 1) 5996 endwalks = 1;6186 endwalks = TRUE; 5997 6187 5998 6188 for(i=nV-1; i>=0; i--) … … 6000 6190 6001 6191 delete next_weight; 6002 }// while6192 }//end of while-loop 6003 6193 6004 6194 if(tp_deg != 1) … … 6014 6204 DefRingPar(orig_target); 6015 6205 else 6016 rChangeCurrRing(VMrDefault(orig_target)); //Aenderung6206 rChangeCurrRing(VMrDefault(orig_target)); 6017 6207 6018 6208 TargetRing=currRing; … … 6030 6220 if( ntestw != 1 || ntwC == 0) 6031 6221 { 6032 /*6033 if(ntestw != 1){6222 if(ntestw != 1 && printout >2) 6223 { 6034 6224 ivString(pert_target_vector, "tau"); 6035 6225 PrintS("\n// ** perturbed target vector doesn't stay in cone!!"); 6036 6226 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6037 idElements(F1, "G"); 6038 } 6039 */ 6227 //idElements(F1, "G"); 6228 } 6040 6229 // LastGB is "better" than the kStd subroutine 6041 6230 to=clock(); … … 6068 6257 Eresult = idrMoveR(G, newRing,currRing); 6069 6258 } 6259 si_opt_1 = save1; //set original options, e. g. option(RedSB) 6070 6260 delete ivNull; 6071 6261 if(tp_deg != 1) … … 6082 6272 tnw+xtnw); 6083 6273 6084 Print("\n// pSetm_Error = (%d)", ErrorCheck()); 6085 Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep, Overflow_Error); 6086 #endif 6274 //Print("\n// pSetm_Error = (%d)", ErrorCheck()); 6275 //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep, Overflow_Error); 6276 #endif 6277 Print("\n//** Mpwalk: Perturbation Walk took %d steps.\n", nstep); 6278 return(Eresult); 6279 } 6280 6281 /******************************************************* 6282 * THE PERTURBATION WALK ALGORITHM WITH RANDOM ELEMENT * 6283 *******************************************************/ 6284 ideal Mprwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, 6285 int op_deg, int tp_deg, int nP, int reduction, int printout) 6286 { 6287 BITSET save1 = si_opt_1; // save current options 6288 if(reduction == 0) 6289 { 6290 si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 6291 //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 6292 //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB)); 6293 } 6294 Set_Error(FALSE); 6295 Overflow_Error = FALSE; 6296 //Print("// pSetm_Error = (%d)", ErrorCheck()); 6297 6298 clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0; 6299 xtextra=0; 6300 xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; 6301 tinput = clock(); 6302 6303 clock_t tim; 6304 6305 nstep = 0; 6306 int i, ntwC=1, ntestw=1, polylength, nV = currRing->N; 6307 BOOLEAN endwalks = FALSE; 6308 6309 ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG; 6310 ring newRing, oldRing, TargetRing; 6311 intvec* iv_M_dp; 6312 intvec* iv_M_lp; 6313 intvec* exivlp = Mivlp(nV); 6314 intvec* curr_weight = new intvec(nV); 6315 intvec* target_weight = new intvec(nV); 6316 for(i=0; i<nV; i++) 6317 { 6318 (*curr_weight)[i] = (*orig_M)[i]; 6319 (*target_weight)[i] = (*target_M)[i]; 6320 } 6321 intvec* orig_target = target_weight; 6322 intvec* pert_target_vector = target_weight; 6323 intvec* ivNull = new intvec(nV); 6324 intvec* iv_dp = MivUnit(nV);// define (1,1,...,1) 6325 #ifndef BUCHBERGER_ALG 6326 intvec* hilb_func; 6327 #endif 6328 intvec* next_weight; 6329 6330 // to avoid (1,0,...,0) as the target vector 6331 intvec* last_omega = new intvec(nV); 6332 for(i=nV-1; i>0; i--) 6333 (*last_omega)[i] = 1; 6334 (*last_omega)[0] = 10000; 6335 6336 ring XXRing = currRing; 6337 6338 to = clock(); 6339 // perturbs the original vector 6340 if(orig_M->length() == nV) 6341 { 6342 if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp" 6343 { 6344 G = MstdCC(Go); 6345 tostd = clock()-to; 6346 if(op_deg != 1) 6347 { 6348 iv_M_dp = MivMatrixOrderdp(nV); 6349 curr_weight = MPertVectors(G, iv_M_dp, op_deg); 6350 } 6351 } 6352 else 6353 { 6354 //define ring order := (a(curr_weight),lp); 6355 if (rParameter(currRing) != NULL) 6356 DefRingPar(curr_weight); 6357 else 6358 rChangeCurrRing(VMrDefault(curr_weight)); 6359 6360 G = idrMoveR(Go, XXRing,currRing); 6361 G = MstdCC(G); 6362 tostd = clock()-to; 6363 if(op_deg != 1) 6364 { 6365 iv_M_dp = MivMatrixOrder(curr_weight); 6366 curr_weight = MPertVectors(G, iv_M_dp, op_deg); 6367 } 6368 } 6369 } 6370 else 6371 { 6372 rChangeCurrRing(VMatrDefault(orig_M)); 6373 G = idrMoveR(Go, XXRing,currRing); 6374 G = MstdCC(G); 6375 tostd = clock()-to; 6376 if(op_deg != 1) 6377 { 6378 curr_weight = MPertVectors(G, orig_M, op_deg); 6379 } 6380 } 6381 6382 delete iv_dp; 6383 if(op_deg != 1) delete iv_M_dp; 6384 6385 ring HelpRing = currRing; 6386 6387 // perturbs the target weight vector 6388 if(target_M->length() == nV) 6389 { 6390 if(tp_deg > 1 && tp_deg <= nV) 6391 { 6392 if (rParameter(currRing) != NULL) 6393 DefRingPar(target_weight); 6394 else 6395 rChangeCurrRing(VMrDefault(target_weight)); 6396 6397 TargetRing = currRing; 6398 ssG = idrMoveR(G,HelpRing,currRing); 6399 if(MivSame(target_weight, exivlp) == 1) 6400 { 6401 iv_M_lp = MivMatrixOrderlp(nV); 6402 target_weight = MPertVectors(ssG, iv_M_lp, tp_deg); 6403 } 6404 else 6405 { 6406 iv_M_lp = MivMatrixOrder(target_weight); 6407 target_weight = MPertVectors(ssG, iv_M_lp, tp_deg); 6408 } 6409 delete iv_M_lp; 6410 pert_target_vector = target_weight; 6411 rChangeCurrRing(HelpRing); 6412 G = idrMoveR(ssG, TargetRing,currRing); 6413 } 6414 } 6415 else 6416 { 6417 if(tp_deg > 1 && tp_deg <= nV) 6418 { 6419 rChangeCurrRing(VMatrDefault(target_M)); 6420 TargetRing = currRing; 6421 ssG = idrMoveR(G,HelpRing,currRing); 6422 target_weight = MPertVectors(ssG, target_M, tp_deg); 6423 } 6424 } 6425 if(printout > 0) 6426 { 6427 Print("\n//** Mprwalk: Random Perturbation Walk of degree (%d,%d):",op_deg,tp_deg); 6428 ivString(curr_weight, "//** Mprwalk: new current weight"); 6429 ivString(target_weight, "//** Mprwalk: new target weight"); 6430 } 6431 while(1) 6432 { 6433 nstep ++; 6434 to = clock(); 6435 // compute an initial form ideal of <G> w.r.t. the weight vector 6436 // "curr_weight" 6437 Gomega = MwalkInitialForm(G, curr_weight); 6438 polylength = lengthpoly(Gomega); 6439 //#ifdef CHECK_IDEAL_MWALK 6440 if(printout > 1) 6441 { 6442 idString(Gomega,"//** Mprwalk: Gomega"); 6443 } 6444 //#endif 6445 6446 if(reduction == 0 && nstep > 1) 6447 { 6448 // check whether weight vector is in the interior of the cone 6449 while(1) 6450 { 6451 FF = middleOfCone(G,Gomega); 6452 if(FF != NULL) 6453 { 6454 idDelete(&G); 6455 G = idCopy(FF); 6456 idDelete(&FF); 6457 next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 6458 if(printout > 0) 6459 { 6460 MivString(curr_weight, target_weight, next_weight); 6461 } 6462 } 6463 else 6464 { 6465 break; 6466 } 6467 for(i=nV-1; i>=0; i--) 6468 { 6469 (*curr_weight)[i] = (*next_weight)[i]; 6470 } 6471 Gomega = MwalkInitialForm(G, curr_weight); 6472 if(printout > 1) 6473 { 6474 idString(Gomega,"//** Mprwalk: Gomega"); 6475 } 6476 } 6477 } 6478 6479 #ifdef ENDWALKS 6480 if(endwalks == TRUE) 6481 { 6482 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6483 idElements(G, "G"); 6484 headidString(G, "G"); 6485 } 6486 #endif 6487 6488 tif = tif + clock()-to; 6489 6490 #ifndef BUCHBERGER_ALG 6491 if(isNolVector(curr_weight) == 0) 6492 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 6493 else 6494 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 6495 #endif // BUCHBERGER_ALG 6496 6497 oldRing = currRing; 6498 6499 if(target_M->length() == nV) 6500 { 6501 // define a new ring with ordering "(a(curr_weight),lp) 6502 if (rParameter(currRing) != NULL) 6503 DefRingPar(curr_weight); 6504 else 6505 rChangeCurrRing(VMrDefault(curr_weight)); 6506 } 6507 else 6508 { 6509 rChangeCurrRing(VMatrRefine(target_M,curr_weight)); 6510 } 6511 newRing = currRing; 6512 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 6513 #ifdef ENDWALKS 6514 if(endwalks == TRUE) 6515 { 6516 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6517 idElements(Gomega1, "Gw"); 6518 headidString(Gomega1, "headGw"); 6519 PrintS("\n// compute a rGB of Gw:\n"); 6520 6521 #ifndef BUCHBERGER_ALG 6522 ivString(hilb_func, "w"); 6523 #endif 6524 } 6525 #endif 6526 6527 tim = clock(); 6528 to = clock(); 6529 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 6530 #ifdef BUCHBERGER_ALG 6531 M = MstdhomCC(Gomega1); 6532 #else 6533 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight); 6534 delete hilb_func; 6535 #endif 6536 //#ifdef CHECK_IDEAL_MWALK 6537 if(printout > 2) 6538 { 6539 idString(M,"//** Mprwalk: M"); 6540 } 6541 //#endif 6542 6543 if(endwalks == TRUE) 6544 { 6545 xtstd = xtstd+clock()-to; 6546 #ifdef ENDWALKS 6547 Print("\n// time for the last std(Gw) = %.2f sec\n", 6548 ((double) clock())/1000000 -((double)tim) /1000000); 6549 #endif 6550 } 6551 else 6552 tstd=tstd+clock()-to; 6553 6554 /* change the ring to oldRing */ 6555 rChangeCurrRing(oldRing); 6556 M1 = idrMoveR(M, newRing,currRing); 6557 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 6558 6559 to=clock(); 6560 /* compute a representation of the generators of submod (M) 6561 with respect to those of mod (Gomega). 6562 Gomega is a reduced Groebner basis w.r.t. the current ring */ 6563 F = MLifttwoIdeal(Gomega2, M1, G); 6564 if(endwalks == FALSE) 6565 tlift = tlift+clock()-to; 6566 else 6567 xtlift=clock()-to; 6568 6569 //#ifdef CHECK_IDEAL_MWALK 6570 if(printout > 2) 6571 { 6572 idString(F,"//** Mprwalk: F"); 6573 } 6574 //#endif 6575 6576 idDelete(&M1); 6577 idDelete(&Gomega2); 6578 idDelete(&G); 6579 6580 // change the ring to newRing 6581 rChangeCurrRing(newRing); 6582 if(reduction == 0) 6583 { 6584 G = idrMoveR(F,oldRing,currRing); 6585 } 6586 else 6587 { 6588 F1 = idrMoveR(F, oldRing,currRing); 6589 if(printout > 2) 6590 { 6591 PrintS("\n //** Mprwalk: reduce the Groebner basis.\n"); 6592 } 6593 to=clock(); 6594 G = kInterRedCC(F1, NULL); 6595 if(endwalks == FALSE) 6596 tred = tred+clock()-to; 6597 else 6598 xtred=clock()-to; 6599 idDelete(&F1); 6600 } 6601 6602 if(endwalks == TRUE) 6603 break; 6604 6605 to=clock(); 6606 6607 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 6608 if(polylength > 0) 6609 { 6610 if(printout > 2) 6611 { 6612 Print("\n //**Mprwalk: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone.\n"); 6613 } 6614 delete next_weight; 6615 next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg); 6616 } 6617 tnw=tnw+clock()-to; 6618 //#ifdef PRINT_VECTORS 6619 if(printout > 0) 6620 { 6621 MivString(curr_weight, target_weight, next_weight); 6622 } 6623 //#endif 6624 6625 if(Overflow_Error == TRUE) 6626 { 6627 ntwC = 0; 6628 //ntestomega = 1; 6629 //Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6630 //idElements(G, "G"); 6631 delete next_weight; 6632 goto FINISH_160302; 6633 } 6634 if(MivComp(next_weight, ivNull) == 1){ 6635 newRing = currRing; 6636 delete next_weight; 6637 //Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6638 break; 6639 } 6640 if(MivComp(next_weight, target_weight) == 1) 6641 endwalks = TRUE; 6642 6643 for(i=nV-1; i>=0; i--) 6644 (*curr_weight)[i] = (*next_weight)[i]; 6645 6646 delete next_weight; 6647 }//while 6648 6649 if(tp_deg != 1) 6650 { 6651 FINISH_160302: 6652 if(target_M->length() == nV) 6653 { 6654 if(MivSame(orig_target, exivlp) == 1) 6655 if (rParameter(currRing) != NULL) 6656 DefRingParlp(); 6657 else 6658 VMrDefaultlp(); 6659 else 6660 if (rParameter(currRing) != NULL) 6661 DefRingPar(orig_target); 6662 else 6663 rChangeCurrRing(VMrDefault(orig_target)); 6664 } 6665 else 6666 { 6667 rChangeCurrRing(VMatrDefault(target_M)); 6668 } 6669 TargetRing=currRing; 6670 F1 = idrMoveR(G, newRing,currRing); 6671 #ifdef CHECK_IDEAL 6672 headidString(G, "G"); 6673 #endif 6674 6675 // check whether the pertubed target vector stays in the correct cone 6676 if(ntwC != 0) 6677 { 6678 ntestw = test_w_in_ConeCC(F1, pert_target_vector); 6679 } 6680 if(ntestw != 1 || ntwC == 0) 6681 { 6682 if(ntestw != 1 && printout > 2) 6683 { 6684 ivString(pert_target_vector, "tau"); 6685 PrintS("\n// **Mprwalk: perturbed target vector doesn't stay in cone."); 6686 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6687 //idElements(F1, "G"); 6688 } 6689 // LastGB is "better" than the kStd subroutine 6690 to=clock(); 6691 ideal eF1; 6692 if(nP == 0 || tp_deg == 1 || MivSame(orig_target, exivlp) != 1 || target_M->length() != nV) 6693 { 6694 if(printout > 2) 6695 { 6696 PrintS("\n// ** Mprwalk: Call \"std\" to compute a Groebner basis.\n"); 6697 } 6698 eF1 = MstdCC(F1); 6699 idDelete(&F1); 6700 } 6701 else 6702 { 6703 if(printout > 2) 6704 { 6705 PrintS("\n// **Mprwalk: Call \"LastGB\" to compute a Groebner basis.\n"); 6706 } 6707 rChangeCurrRing(newRing); 6708 ideal F2 = idrMoveR(F1, TargetRing,currRing); 6709 eF1 = LastGB(F2, curr_weight, tp_deg-1); 6710 F2=NULL; 6711 } 6712 xtextra=clock()-to; 6713 ring exTargetRing = currRing; 6714 6715 rChangeCurrRing(XXRing); 6716 Eresult = idrMoveR(eF1, exTargetRing,currRing); 6717 } 6718 else{ 6719 rChangeCurrRing(XXRing); 6720 Eresult = idrMoveR(F1, TargetRing,currRing); 6721 } 6722 } 6723 else { 6724 rChangeCurrRing(XXRing); 6725 Eresult = idrMoveR(G, newRing,currRing); 6726 } 6727 si_opt_1 = save1; //set original options, e. g. option(RedSB) 6728 delete ivNull; 6729 if(tp_deg != 1) 6730 delete target_weight; 6731 6732 if(op_deg != 1 ) 6733 delete curr_weight; 6734 6735 delete exivlp; 6736 delete last_omega; 6737 6738 #ifdef TIME_TEST 6739 TimeStringFractal(tinput, tostd, tif+xtif, tstd+xtstd,0, tlift+xtlift, tred+xtred, 6740 tnw+xtnw); 6741 6742 //Print("\n// pSetm_Error = (%d)", ErrorCheck()); 6743 //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep, Overflow_Error); 6744 #endif 6745 Print("\n//** Mprwalk: Perturbation Walk took %d steps.\n", nstep); 6087 6746 return(Eresult); 6088 6747 } … … 6110 6769 * Perturb the start weight vector at the top level, i.e. nlev = 1 * 6111 6770 ***********************************************************************/ 6112 static ideal rec_fractal_call(ideal G, int nlev, intvec* omtmp) 6771 static ideal rec_fractal_call(ideal G, int nlev, intvec* ivtarget, 6772 int reduction, int printout) 6113 6773 { 6114 6774 Overflow_Error = FALSE; 6115 //Print("\n\n// Entering the %d-th recursion:", nlev); 6116 6775 if(printout >0) 6776 { 6777 Print("\n\n// Entering the %d-th recursion:", nlev); 6778 } 6117 6779 int i, nV = currRing->N; 6118 6780 ring new_ring, testring; 6119 6781 //ring extoRing; 6120 ideal Gomega, Gomega1, Gomega2, F , F1, Gresult, Gresult1, G1, Gt;6782 ideal Gomega, Gomega1, Gomega2, FF, F, F1, Gresult, Gresult1, G1, Gt; 6121 6783 int nwalks = 0; 6122 6784 intvec* Mwlp; … … 6127 6789 intvec* next_vect; 6128 6790 intvec* omega2 = new intvec(nV); 6129 intvec* altomega = new intvec(nV); 6130 6791 intvec* omtmp = new intvec(nV); 6792 //intvec* altomega = new intvec(nV); 6793 6794 for(i = nV -1; i>0; i--) 6795 { 6796 (*omtmp)[i] = (*ivtarget)[i]; 6797 } 6131 6798 //BOOLEAN isnewtarget = FALSE; 6132 6799 … … 6169 6836 NEXT_VECTOR_FRACTAL: 6170 6837 to=clock(); 6171 / * determine the next border */6838 // determine the next border 6172 6839 next_vect = MkInterRedNextWeight(omega,omega2,G); 6173 6840 xtnw=xtnw+clock()-to; 6174 #ifdef PRINT_VECTORS 6175 MivString(omega, omega2, next_vect); 6176 #endif 6841 6177 6842 oRing = currRing; 6178 6843 6179 / * We only perturb the current target vector at the recursion level 1 */6844 // We only perturb the current target vector at the recursion level 1 6180 6845 if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3 6181 if (MivComp(next_vect, omega2) != 1) 6182 { 6183 /* to dispense with taking initial (and lifting/interreducing 6184 after the call of recursion */ 6185 //Print("\n\n// ** Perturb the both vectors with degree %d with",nlev); 6186 //idElements(G, "G"); 6846 if (MivComp(next_vect, omega2) == 1) 6847 { 6848 // to dispense with taking initial (and lifting/interreducing 6849 // after the call of recursion 6850 if(printout > 0) 6851 { 6852 Print("\n//** rec_fractal_call: Perturb the both vectors with degree %d.",nlev); 6853 //idElements(G, "G"); 6854 } 6187 6855 6188 6856 Xngleich = 1; 6189 6857 nlev +=1; 6190 6858 6191 if (rParameter(currRing) != NULL) 6192 DefRingPar(omtmp); 6859 if(ivtarget->length() == nV) 6860 { 6861 if (rParameter(currRing) != NULL) 6862 DefRingPar(omtmp); 6863 else 6864 rChangeCurrRing(VMrDefault(omtmp)); 6865 } 6193 6866 else 6194 rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung3 6195 6867 { 6868 rChangeCurrRing(VMatrDefault(ivtarget)); 6869 } 6196 6870 testring = currRing; 6197 6871 Gt = idrMoveR(G, oRing,currRing); 6198 6872 6199 /* perturb the original target vector w.r.t. the current GB */ 6200 delete Xtau; 6201 Xtau = NewVectorlp(Gt); 6873 // perturb the original target vector w.r.t. the current GB 6874 if(ivtarget->length() == nV) 6875 { 6876 delete Xtau; 6877 Xtau = NewVectorlp(Gt); 6878 } 6879 else 6880 { 6881 delete Xtau; 6882 Xtau = Mfpertvector(Gt,ivtarget); 6883 } 6202 6884 6203 6885 rChangeCurrRing(oRing); 6204 6886 G = idrMoveR(Gt, testring,currRing); 6205 6887 6206 / * perturb the current vector w.r.t. the current GB */6888 // perturb the current vector w.r.t. the current GB 6207 6889 Mwlp = MivWeightOrderlp(omega); 6208 6890 Xsigma = Mfpertvector(G, Mwlp); … … 6217 6899 to=clock(); 6218 6900 6219 / * to avoid the value of Overflow_Error that occur in Mfpertvector*/6901 // to avoid the value of Overflow_Error that occur in Mfpertvector 6220 6902 Overflow_Error = FALSE; 6221 6222 6903 next_vect = MkInterRedNextWeight(omega,omega2,G); 6223 6904 xtnw=xtnw+clock()-to; 6224 6225 #ifdef PRINT_VECTORS 6905 }// end of (if MivComp(next_vect, omega2) == 1) 6906 6907 //#ifdef PRINT_VECTORS 6908 if(printout > 0) 6909 { 6226 6910 MivString(omega, omega2, next_vect); 6227 #endif 6228 } 6229 6230 6231 /* check whether the the computed vector is in the correct cone */ 6232 /* If no, the reduced GB of an omega-homogeneous ideal will be 6233 computed by Buchberger algorithm and stop this recursion step*/ 6234 //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6 6235 if(Overflow_Error == TRUE) 6911 } 6912 //#endif 6913 6914 // check whether the the computed vector is in the correct cone. 6915 // If no, compute the reduced Groebner basis of an omega-homogeneous 6916 // ideal with Buchberger's algorithm and stop this recursion step 6917 if(Overflow_Error == TRUE || test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6 6236 6918 { 6237 6919 delete next_vect; 6238 if (rParameter(currRing) != NULL) 6239 { 6240 DefRingPar(omtmp); 6920 if(ivtarget->length() == nV) 6921 { 6922 if (rParameter(currRing) != NULL) 6923 DefRingPar(omtmp); 6924 else 6925 rChangeCurrRing(VMrDefault(omtmp)); 6241 6926 } 6242 6927 else 6243 6928 { 6244 rChangeCurrRing(VM rDefault1(omtmp)); // Aenderung46929 rChangeCurrRing(VMatrDefault(ivtarget)); 6245 6930 } 6246 6931 #ifdef TEST_OVERFLOW … … 6248 6933 Gt = NULL; return(Gt); 6249 6934 #endif 6250 6251 //Print("\n\n// apply BB's alg. in ring r = %s;", rString(currRing)); 6935 if(printout > 0) 6936 { 6937 Print("\n//** rec_fractal_call: Applying Buchberger's algorithm in ring r = %s;", 6938 rString(currRing)); 6939 } 6252 6940 to=clock(); 6253 6941 Gt = idrMoveR(G, oRing,currRing); … … 6257 6945 6258 6946 delete omega2; 6259 delete altomega; 6260 6261 //Print("\n// Leaving the %d-th recursion with %d steps", nlev, nwalks); 6262 //Print(" ** Overflow_Error? (%d)", Overflow_Error); 6947 //delete altomega; 6948 if(printout > 0) 6949 { 6950 Print("\n//** rec_fractal_call: Overflow. Leaving the %d-th recursion with %d steps.\n", 6951 nlev, nwalks); 6952 //Print(" ** Overflow_Error? (%d)", Overflow_Error); 6953 } 6954 6263 6955 nnflow ++; 6264 6265 6956 Overflow_Error = FALSE; 6266 6957 return (G1); … … 6277 6968 if (MivComp(next_vect, XivNull) == 1) 6278 6969 { 6279 if (rParameter(currRing) != NULL) 6280 DefRingPar(omtmp); 6970 if(ivtarget->length() == nV) 6971 { 6972 if (rParameter(currRing) != NULL) 6973 DefRingPar(omtmp); 6974 else 6975 rChangeCurrRing(VMrDefault(omtmp)); 6976 } 6281 6977 else 6282 rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung5 6978 { 6979 rChangeCurrRing(VMatrDefault(ivtarget)); 6980 } 6283 6981 6284 6982 testring = currRing; 6285 6983 Gt = idrMoveR(G, oRing,currRing); 6286 6984 6287 if(test_w_in_ConeCC(Gt, omega2) == 1) { 6985 if(test_w_in_ConeCC(Gt, omega2) == 1) 6986 { 6987 delete omega2; 6988 delete next_vect; 6989 //delete altomega; 6990 if(printout > 0) 6991 { 6992 Print("\n//** rec_fractal_call: Correct cone. Leaving the %d-th recursion with %d steps.\n", 6993 nlev, nwalks); 6994 } 6995 return (Gt); 6996 } 6997 else 6998 { 6999 if(printout > 0) 7000 { 7001 Print("\n//** rec_fractal_call: Wrong cone. Tau doesn't stay in the correct cone.\n"); 7002 } 7003 7004 #ifndef MSTDCC_FRACTAL 7005 intvec* Xtautmp; 7006 if(ivtarget->length() == nV) 7007 { 7008 Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp)); 7009 } 7010 else 7011 { 7012 Xtautmp = Mfpertvector(Gt, ivtarget); 7013 } 7014 #ifdef TEST_OVERFLOW 7015 if(Overflow_Error == TRUE) 7016 Gt = NULL; return(Gt); 7017 #endif 7018 7019 if(MivSame(Xtau, Xtautmp) == 1) 7020 { 7021 if(printout > 0) 7022 { 7023 Print("\n//** rec_fractal_call: Updated vectors are equal to the old vectors.\n"); 7024 } 7025 delete Xtautmp; 7026 goto FRACTAL_MSTDCC; 7027 } 7028 7029 Xtau = Xtautmp; 7030 Xtautmp = NULL; 7031 7032 for(i=nV-1; i>=0; i--) 7033 (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i]; 7034 7035 rChangeCurrRing(oRing); 7036 G = idrMoveR(Gt, testring,currRing); 7037 7038 goto NEXT_VECTOR_FRACTAL; 7039 #endif 7040 7041 FRACTAL_MSTDCC: 7042 if(printout > 0) 7043 { 7044 Print("\n//** rec_fractal_call: Wrong cone. Applying Buchberger's algorithm in ring = %s.\n", 7045 rString(currRing)); 7046 } 7047 to=clock(); 7048 G = MstdCC(Gt); 7049 xtextra=xtextra+clock()-to; 7050 7051 oRing = currRing; 7052 7053 // update the original target vector w.r.t. the current GB 7054 if(ivtarget->length() == nV) 7055 { 7056 if(MivSame(Xivinput, Xivlp) == 1) 7057 if (rParameter(currRing) != NULL) 7058 DefRingParlp(); 7059 else 7060 VMrDefaultlp(); 7061 else 7062 if (rParameter(currRing) != NULL) 7063 DefRingPar(Xivinput); 7064 else 7065 rChangeCurrRing(VMrDefault(Xivinput)); 7066 } 7067 else 7068 { 7069 rChangeCurrRing(VMatrRefine(ivtarget,Xivinput)); 7070 } 7071 testring = currRing; 7072 Gt = idrMoveR(G, oRing,currRing); 7073 7074 // perturb the original target vector w.r.t. the current GB 7075 if(ivtarget->length() == nV) 7076 { 7077 delete Xtau; 7078 Xtau = NewVectorlp(Gt); 7079 } 7080 else 7081 { 7082 delete Xtau; 7083 Xtau = Mfpertvector(Gt,ivtarget); 7084 } 7085 7086 rChangeCurrRing(oRing); 7087 G = idrMoveR(Gt, testring,currRing); 7088 7089 delete omega2; 7090 delete next_vect; 7091 //delete altomega; 7092 if(printout > 0) 7093 { 7094 Print("\n//** rec_fractal_call: Vectors updated. Leaving the %d-th recursion with %d steps.\n", 7095 nlev, nwalks); 7096 //Print(" ** Overflow_Error? (%d)", Overflow_Error); 7097 } 7098 if(Overflow_Error == TRUE) 7099 nnflow ++; 7100 7101 Overflow_Error = FALSE; 7102 return(G); 7103 } 7104 }// end of (if next_vect==nullvector) 7105 7106 for(i=nV-1; i>=0; i--) { 7107 //(*altomega)[i] = (*omega)[i]; 7108 (*omega)[i] = (*next_vect)[i]; 7109 } 7110 delete next_vect; 7111 7112 to=clock(); 7113 // Take the initial form of <G> w.r.t. omega 7114 Gomega = MwalkInitialForm(G, omega); 7115 xtif=xtif+clock()-to; 7116 if(printout > 1) 7117 { 7118 idString(Gomega,"//** rec_fractal_call: Gomega"); 7119 } 7120 7121 if(reduction == 0) 7122 { 7123 // Check whether the intermediate weight vector lies in the interior of the cone. 7124 // If so, only perform reductions. Otherwise apply Buchberger's algorithm. 7125 FF = middleOfCone(G,Gomega); 7126 if( FF != NULL) 7127 { 7128 idDelete(&G); 7129 G = idCopy(FF); 7130 idDelete(&FF); 7131 // Compue next vector. 7132 goto NEXT_VECTOR_FRACTAL; 7133 } 7134 } 7135 7136 #ifndef BUCHBERGER_ALG 7137 if(isNolVector(omega) == 0) 7138 hilb_func = hFirstSeries(Gomega,NULL,NULL,omega,currRing); 7139 else 7140 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 7141 #endif 7142 7143 if(ivtarget->length() == nV) 7144 { 7145 if (rParameter(currRing) != NULL) 7146 DefRingPar(omega); 7147 else 7148 rChangeCurrRing(VMrDefault(omega)); 7149 } 7150 else 7151 { 7152 rChangeCurrRing(VMatrRefine(ivtarget,omega)); 7153 } 7154 Gomega1 = idrMoveR(Gomega, oRing,currRing); 7155 7156 // Maximal recursion depth, to compute a red. GB 7157 // Fractal walk with the alternative recursion 7158 // alternative recursion 7159 if(nlev == Xnlev || lengthpoly(Gomega1) == 0) 7160 { 7161 if(printout > 1) 7162 { 7163 Print("\n//** rec_fractal_call: Maximal recursion depth.\n"); 7164 } 7165 7166 to=clock(); 7167 #ifdef BUCHBERGER_ALG 7168 Gresult = MstdhomCC(Gomega1); 7169 #else 7170 Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega); 7171 delete hilb_func; 7172 #endif 7173 xtstd=xtstd+clock()-to; 7174 } 7175 else 7176 { 7177 rChangeCurrRing(oRing); 7178 Gomega1 = idrMoveR(Gomega1, oRing,currRing); 7179 Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega,reduction,printout); 7180 } 7181 if(printout > 2) 7182 { 7183 idString(Gresult,"//** rec_fractal_call: M"); 7184 } 7185 //convert a Groebner basis from a ring to another ring 7186 new_ring = currRing; 7187 7188 rChangeCurrRing(oRing); 7189 Gresult1 = idrMoveR(Gresult, new_ring,currRing); 7190 Gomega2 = idrMoveR(Gomega1, new_ring,currRing); 7191 7192 to=clock(); 7193 // Lifting process 7194 F = MLifttwoIdeal(Gomega2, Gresult1, G); 7195 xtlift=xtlift+clock()-to; 7196 if(printout > 2) 7197 { 7198 idString(F,"//** rec_fractal_call: F"); 7199 } 7200 idDelete(&Gresult1); 7201 idDelete(&Gomega2); 7202 idDelete(&G); -
Property
mode
changed from