Changeset c10f46 in git for Singular/LIB/modstd.lib
 Timestamp:
 Aug 3, 2012, 1:06:22 PM (11 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0d7ffeded36f5660db047d76670b184bc97cecd3')
 Children:
 b0732eb4254fbe76540d72c337f3d7f985adf948
 Parents:
 b8610796eaca5ba7db8dc75f8bc88d698e077b365abb0e7fdca82a2bad9778e410148c74d7b49a0b
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/modstd.lib
rb86107 rc10f46 1 1 //////////////////////////////////////////////////////////////////////////////// 2 version="$Id $";2 version="$Id: modstd.lib 14375 20110823 09:29:47Z steidel $"; 3 3 category = "Commutative Algebra"; 4 4 info=" … … 8 8 @* G. Pfister pfister@mathematik.unikl.de 9 9 @* H. Schoenemann hannes@mathematik.unikl.de 10 @* A. Steenpass steenpass@mathematik.unikl.de 10 11 @* S. Steidel steidel@mathematik.unikl.de 11 12 … … 26 27 LIB "poly.lib"; 27 28 LIB "ring.lib"; 29 LIB "parallel.lib"; 28 30 29 31 //////////////////////////////////////////////////////////////////////////////// … … 243 245 if(printlevel >= 11) 244 246 { 245 "isIncluded( K,J) takes "+string(timer  t)+" seconds";247 "isIncluded(J,K) takes "+string(timer  t)+" seconds"; 246 248 "j = "+string(j); 247 249 } … … 373 375 //////////////////////////////////////////////////////////////////////////////// 374 376 375 proc primeTest( idealI, bigint p)377 proc primeTest(def II, bigint p) 376 378 { 379 if(typeof(II) == "string") 380 { 381 execute("ideal I = "+II+";"); 382 } 383 else 384 { 385 ideal I = II; 386 } 387 388 I = simplify(I, 2); // erase zero generators 389 377 390 int i,j; 391 poly f; 392 number cnt; 378 393 for(i = 1; i <= size(I); i++) 379 394 { 380 for(j = 1; j <= size(I[i]); j++) 381 { 382 if((numerator(leadcoef(I[i][j])) mod p) == 0) { return(0); } 383 if((denominator(leadcoef(I[i][j])) mod p) == 0) { return(0); } 395 f = cleardenom(I[i]); 396 if(f == 0) { return(0); } 397 cnt = leadcoef(I[i])/leadcoef(f); 398 if((numerator(cnt) mod p) == 0) { return(0); } 399 if((denominator(cnt) mod p) == 0) { return(0); } 400 for(j = size(f); j > 0; j) 401 { 402 if((leadcoef(f[j]) mod p) == 0) { return(0); } 384 403 } 385 404 } … … 390 409 391 410 proc primeList(ideal I, int n, list #) 392 "USAGE: primeList(I,n); ( resp. primeList(I,n,L); ) I ideal, n integer 393 RETURN: the intvec of n greatest primes <= 2147483647 (resp. n greatest primes 411 "USAGE: primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal, 412 n integer 413 RETURN: the intvec of n greatest primes <= 2147483647 (resp. n greatest primes 394 414 < L[size(L)] union with L) such that none of these primes divides any 395 415 coefficient occuring in I 416 NOTE: The number of cores to use can be defined by ncores, default is 1. 396 417 EXAMPLE: example primList; shows an example 397 418 " … … 399 420 intvec L; 400 421 int i,p; 422 int ncores = 1; 423 424 // Initialize optional parameter ncores  425 if(size(#) > 0) 426 { 427 if(size(#) == 1) 428 { 429 if(typeof(#[1]) == "int") 430 { 431 ncores = #[1]; 432 # = list(); 433 } 434 } 435 else 436 { 437 ncores = #[2]; 438 } 439 } 440 401 441 if(size(#) == 0) 402 442 { … … 421 461 } 422 462 if(p == 2) { ERROR("no more primes"); } 423 for(i = 2; i <= n; i++) 424 { 425 p = prime(p1); 426 while(!primeTest(I,p)) 463 if(ncores == 1) 464 { 465 for(i = 2; i <= n; i++) 427 466 { 428 467 p = prime(p1); 429 if(p == 2) { ERROR("no more primes"); } 430 } 431 L[size(L)+1] = p; 468 while(!primeTest(I,p)) 469 { 470 p = prime(p1); 471 if(p == 2) { ERROR("no more primes"); } 472 } 473 L[size(L)+1] = p; 474 } 475 } 476 else 477 { 478 int neededSize = size(L)+n1;; 479 list parallelResults; 480 list arguments; 481 int neededPrimes = neededSizesize(L); 482 while(neededPrimes > 0) 483 { 484 arguments = list(); 485 for(i = ((neededPrimes div ncores)+1(neededPrimes%ncores == 0)) 486 *ncores; i > 0; i) 487 { 488 p = prime(p1); 489 if(p == 2) { ERROR("no more primes"); } 490 arguments[i] = list("I", p); 491 } 492 parallelResults = parallelWaitAll("primeTest", arguments, 493 list(list(list(ncores)))); 494 for(i = size(arguments); i > 0; i) 495 { 496 if(parallelResults[i]) 497 { 498 L[size(L)+1] = arguments[i][2]; 499 } 500 } 501 neededPrimes = neededSizesize(L); 502 } 503 if(size(L) > neededSize) 504 { 505 L = L[1..neededSize]; 506 } 432 507 } 433 508 return(L); … … 470 545 ideal I = imap(R,I); 471 546 547 intvec opt = option(get); 472 548 option(none); 473 549 option(prompt); … … 504 580 505 581 setring R; 582 option(set, opt); 506 583 return(imap(newR,sI),imap(newR,T)); 507 584 } … … 537 614 The standard basis computation modulo p does also vary depending on the 538 615 integer variant, namely 539 @*  variant = 1/4: std(.,#[1]) resp. groebner, 540 @*  variant = 2/5: groebner, 541 @*  variant = 3/6: homog.  std(.,#[1]) resp. groebner  dehomog.. 616 @*  variant = 1: std(.,#[1]) resp. groebner, 617 @*  variant = 2: groebner, 618 @*  variant = 3: homog.  std(.,#[1]) resp. groebner  dehomog., 619 @*  variant = 4: fglm. 542 620 EXAMPLE: example modpStd; shows an example 543 621 " … … 552 630 option(redSB); 553 631 554 int t = timer; 555 if((variant == 1)  (variant == 4)) 632 if(variant == 1) 556 633 { 557 634 if(size(#) > 0) … … 565 642 } 566 643 567 if( (variant == 2)  (variant == 5))644 if(variant == 2) 568 645 { 569 646 i = groebner(i); 570 647 } 571 648 572 if( (variant == 3)  (variant == 6))649 if(variant == 3) 573 650 { 574 651 list rl = ringlist(@r); … … 607 684 } 608 685 609 t = timer;610 686 i = subst(i, homvar, 1); 611 687 i = simplify(i, 34); … … 615 691 i = interred(i); 616 692 kill HomR; 693 } 694 695 if(variant == 4) 696 { 697 def R1 = changeord("dp"); 698 setring R1; 699 ideal i = fetch(@r,i); 700 i = std(i); 701 setring @r; 702 i = fglm(R1,i); 617 703 } 618 704 … … 646 732 RETURN: a standard basis of I if no warning appears; 647 733 NOTE: The procedure computes a standard basis of I (over the rational 648 numbers) by using modular methods. If a warning appears then the 649 result is a standard basis containing I and with high probability 650 a standard basis of I. 734 numbers) by using modular methods. 651 735 By default the procedure computes a standard basis of I for sure, but 652 736 if the optional parameter #[2] = 0, it computes a standard basis of I 653 with high probability and a warning appears at the end.737 with high probability. 654 738 The procedure distinguishes between different variants for the standard 655 739 basis computation in positive characteristic depending on the ordering 656 740 of the basering, the parameter #[2] and if the ideal I is homogeneous. 657 @*  variant = 1, if I is homogeneous and exactness = 0, 658 @*  variant = 2, if I is not homogeneous, 1blockordering and 659 exactness = 0, 741 @*  variant = 1, if I is homogeneous, 742 @*  variant = 2, if I is not homogeneous, 1blockordering, 660 743 @*  variant = 3, if I is not homogeneous, complicated ordering (lp or 661 > 1 block) and exactness = 0, 662 @*  variant = 4, if I is homogeneous and exactness = 1, 663 @*  variant = 5, if I is not homogeneous, 1blockordering and 664 exactness = 1, 665 @*  variant = 6, if I is not homogeneous, complicated ordering (lp or 666 > 1 block) and exactness = 1. 744 > 1 block), 745 @*  variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0. 667 746 EXAMPLE: example modStd; shows an example 668 747 " … … 752 831 if(printlevel >= 10) 753 832 { 754 "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3); 833 "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3) 834 +", exactness = "+string(exactness); 755 835 } 756 836 … … 761 841 762 842 // Initialize the list of primes  763 intvec L = primeList(I,n2); 843 int tt = timer; 844 int rt = rtimer; 845 intvec L = primeList(I,n2,n1); 846 if(printlevel >= 10) 847 { 848 "CPUtime for primeList: "+string(timertt)+" seconds."; 849 "Realtime for primeList: "+string(rtimerrt)+" seconds."; 850 } 764 851 L[5] = prime(random(an,en)); 765 852 … … 768 855 int h = homog(I); 769 856 770 inttt = timer;771 intrt = rtimer;857 tt = timer; 858 rt = rtimer; 772 859 773 860 if(!mixedTest()) … … 775 862 if(h) 776 863 { 777 if(exactness == 0) 778 { 779 variant = 1; 780 if(printlevel >= 10) { "variant = 1"; } 781 } 782 if(exactness == 1) 783 { 784 variant = 4; 785 if(printlevel >= 10) { "variant = 4"; } 786 } 864 variant = 1; 865 if(printlevel >= 10) { "variant = 1"; } 866 787 867 rl[1] = L[5]; 788 868 def @r = ring(rl); … … 802 882 if((find(ordstr_R0, "M") > 0)  (find(ordstr_R0, "a") > 0)  neg) 803 883 { 804 if(exactness == 0) 805 { 806 variant = 2; 807 if(printlevel >= 10) { "variant = 2"; } 808 } 809 if(exactness == 1) 810 { 811 variant = 5; 812 if(printlevel >= 10) { "variant = 5"; } 813 } 884 variant = 2; 885 if(printlevel >= 10) { "variant = 2"; } 814 886 } 815 887 else … … 827 899 if((order == "simple")  (size(rl) > 4)) 828 900 { 829 if(exactness == 0) 901 variant = 2; 902 if(printlevel >= 10) { "variant = 2"; } 903 } 904 else 905 { 906 rl[1] = L[5]; 907 def @r = ring(rl); 908 setring @r; 909 910 def @s = changeord("dp"); 911 setring @s; 912 ideal I = std(fetch(R0,I)); 913 if(dim(I) == 0) 830 914 { 831 variant = 2;832 if(printlevel >= 10) { "variant = 2"; }915 variant = 4; 916 if(printlevel >= 10) { "variant = 4"; } 833 917 } 834 if(exactness == 1) 835 { 836 variant = 5; 837 if(printlevel >= 10) { "variant = 5"; } 838 } 839 } 840 else 841 { 842 if(exactness == 0) 918 else 843 919 { 844 920 variant = 3; 845 921 if(printlevel >= 10) { "variant = 3"; } 922 923 int nvar@r = nvars(@r); 924 intvec w; 925 for(i = 1; i <= nvar@r; i++) 926 { 927 w[i] = deg(var(i)); 928 } 929 w[nvar@r + 1] = 1; 930 931 list hiRi = hilbRing(fetch(R0,I),w); 932 intvec W = hiRi[2]; 933 @s = hiRi[1]; 934 setring @s; 935 936 Id(1) = std(Id(1)); 937 intvec hi = hilb(Id(1), 1, W); 846 938 } 847 if(exactness == 1) 848 { 849 variant = 6; 850 if(printlevel >= 10) { "variant = 6"; } 851 } 852 853 rl[1] = L[5]; 854 def @r = ring(rl); 855 setring @r; 856 int nvar@r = nvars(@r); 857 intvec w; 858 for(i = 1; i <= nvar@r; i++) 859 { 860 w[i] = deg(var(i)); 861 } 862 w[nvar@r + 1] = 1; 863 864 list hiRi = hilbRing(fetch(R0,I),w); 865 intvec W = hiRi[2]; 866 def @s = hiRi[1]; 867 setring @s; 868 869 Id(1) = std(Id(1)); 870 intvec hi = hilb(Id(1), 1, W); 871 939 872 940 setring R0; 873 941 kill @r,@s; … … 971 1039 link l(i) = "ssi:fork"; 972 1040 open(l(i)); 973 if((variant == 1)  (variant == 3)  974 (variant == 4)  (variant == 6)) 1041 if((variant == 1)  (variant == 3)) 975 1042 { 976 1043 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), 977 1044 eval(variant), eval(hi)))); 978 1045 } 979 if((variant == 2)  (variant == 5))1046 if((variant == 2)  (variant == 4)) 980 1047 { 981 1048 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), … … 985 1052 986 1053 int t = timer; 987 if((variant == 1)  (variant == 3)  (variant == 4)  (variant == 6))1054 if((variant == 1)  (variant == 3)) 988 1055 { 989 1056 P = modpStd(I_for_fork, L[1], variant, hi); 990 1057 } 991 if((variant == 2)  (variant == 5))1058 if((variant == 2)  (variant == 4)) 992 1059 { 993 1060 P = modpStd(I_for_fork, L[1], variant); … … 1005 1072 // Main standard basis computations in positive  1006 1073 // characteristic start here  1074 1075 list arguments_farey, results_farey; 1007 1076 1008 1077 while(1) … … 1029 1098 if(j <= size(L)) 1030 1099 { 1031 if((variant == 1)  (variant == 3)  1032 (variant == 4)  (variant == 6)) 1100 if((variant == 1)  (variant == 3)) 1033 1101 { 1034 1102 write(l(i), quote(modpStd(I_for_fork, eval(L[j]), … … 1036 1104 j++; 1037 1105 } 1038 if((variant == 2)  (variant == 5))1106 if((variant == 2)  (variant == 4)) 1039 1107 { 1040 1108 write(l(i), quote(modpStd(I_for_fork, … … 1062 1130 while(j <= size(L)) 1063 1131 { 1064 if((variant == 1)  (variant == 3)  1065 (variant == 4)  (variant == 6)) 1132 if((variant == 1)  (variant == 3)) 1066 1133 { 1067 1134 P = modpStd(I, L[j], variant, hi); 1068 1135 } 1069 if((variant == 2)  (variant == 5))1136 if((variant == 2)  (variant == 4)) 1070 1137 { 1071 1138 P = modpStd(I, L[j], variant); … … 1095 1162 // Lift results to basering via farey  1096 1163 1097 tt = timer; 1164 tt = timer; rt = rtimer; 1098 1165 N = T2[1]; 1099 1166 for(i = 2; i <= size(T2); i++) { N = N*T2[i]; } 1100 1167 H = chinrem(T1,T2); 1101 J = farey(H,N); 1102 if(printlevel >= 10) { "Liftingprocess takes "+string(timer  tt) 1103 +" seconds"; } 1168 if(n1 == 1) 1169 { 1170 J = farey(H,N); 1171 } 1172 else 1173 { 1174 for(i = ncols(H); i > 0; i) 1175 { 1176 arguments_farey[i] = list(ideal(H[i]), N); 1177 } 1178 results_farey = parallelWaitAll("farey", arguments_farey, 1179 list(list(list(n1)))); 1180 for(i = ncols(H); i > 0; i) 1181 { 1182 J[i] = results_farey[i][1]; 1183 } 1184 } 1185 if(printlevel >= 10) 1186 { 1187 "CPUtime for liftingprocess is "+string(timer  tt)+" seconds."; 1188 "Realtime for liftingprocess is "+string(rtimer  rt)+" seconds."; 1189 } 1104 1190 1105 1191 // Test if we already have a standard basis of I  1106 1192 1107 1193 tt = timer; rt = rtimer; 1108 if((variant == 1)  (variant == 3)  (variant == 4)  (variant == 6))1194 if((variant == 1)  (variant == 3)) 1109 1195 { 1110 1196 pTest = pTestSB(I,J,L,variant,hi); 1111 1197 } 1112 if((variant == 2)  (variant == 5))1198 if((variant == 2)  (variant == 4)) 1113 1199 { 1114 1200 pTest = pTestSB(I,J,L,variant); … … 1132 1218 1133 1219 attrib(J,"isSB",1); 1134 tt = timer; rt = rtimer; 1135 sizeTest = 1  isIncluded(I,J,n1); 1136 1137 if(printlevel >= 10) 1138 { 1139 "CPUtime for checking if I subset <G> is " 1140 +string(timer  tt)+" seconds."; 1141 "Realtime for checking if I subset <G> is " 1142 +string(rtimer  rt)+" seconds."; 1143 } 1144 1145 if(sizeTest == 0) 1146 { 1147 if((variant == 1)  (variant == 2)  (variant == 3)) 1148 { 1149 "==========================================================="; 1150 "WARNING: Output might not be a standard basis of the input."; 1151 "==========================================================="; 1152 option(set, opt); 1153 if(n1 > 1) { kill I_for_fork; } 1154 return(J); 1155 } 1156 if((variant == 4)  (variant == 5)  (variant == 6)) 1220 1221 if(exactness == 0) 1222 { 1223 option(set, opt); 1224 if(n1 > 1) { kill I_for_fork; } 1225 return(J); 1226 } 1227 1228 if(exactness == 1) 1229 { 1230 tt = timer; rt = rtimer; 1231 sizeTest = 1  isIncluded(I,J,n1); 1232 1233 if(printlevel >= 10) 1234 { 1235 "CPUtime for checking if I subset <G> is " 1236 +string(timer  tt)+" seconds."; 1237 "Realtime for checking if I subset <G> is " 1238 +string(rtimer  rt)+" seconds."; 1239 } 1240 1241 if(sizeTest == 0) 1157 1242 { 1158 1243 tt = timer; rt = rtimer; … … 1185 1270 1186 1271 j = size(L) + 1; 1187 L = primeList(I,n3,L); 1272 tt = timer; rt = rtimer; 1273 L = primeList(I,n3,L,n1); 1274 if(printlevel >= 10) 1275 { 1276 "CPUtime for primeList: "+string(timertt)+" seconds."; 1277 "Realtime for primeList: "+string(rtimerrt)+" seconds."; 1278 } 1188 1279 1189 1280 if(n1 > 1) … … 1192 1283 { 1193 1284 open(l(i)); 1194 if((variant == 1)  (variant == 3)  1195 (variant == 4)  (variant == 6)) 1285 if((variant == 1)  (variant == 3)) 1196 1286 { 1197 1287 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i1]), 1198 1288 eval(variant), eval(hi)))); 1199 1289 } 1200 if((variant == 2)  (variant == 5))1290 if((variant == 2)  (variant == 4)) 1201 1291 { 1202 1292 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i1]), … … 1318 1408 def S = ring(rl); 1319 1409 setring S; 1410 intvec opt = option(get); 1320 1411 option(redSB); 1321 1412 module Z,M,Z2; … … 1362 1453 if(size(reduce(G3,testG1)) == 0) 1363 1454 { 1455 option(set, opt); 1364 1456 return(G3); 1365 1457 }
Note: See TracChangeset
for help on using the changeset viewer.