Changeset 7f30e2 in git for Singular/LIB/modstd.lib
 Timestamp:
 Nov 25, 2011, 5:22:38 PM (11 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
 Children:
 f1cfef36a45967a6dc88672e6efc050bf9c592b1
 Parents:
 4093f96acfcc6419c56ee983c447099c3b87a36e
 gitauthor:
 Andreas Steenpass <steenpass@mathematik.unikl.de>20111125 17:22:38+01:00
 gitcommitter:
 Andreas Steenpass <steenpass@mathematik.unikl.de>20120802 18:26:02+02:00
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/modstd.lib
r4093f9 r7f30e2 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 377 388 int i,j; 389 poly f; 390 number cnt; 378 391 for(i = 1; i <= size(I); i++) 379 392 { 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); } 393 f = cleardenom(I[i]); 394 if(f == 0) { return(0); } 395 cnt = leadcoef(I[i])/leadcoef(f); 396 if((numerator(cnt) mod p) == 0) { return(0); } 397 if((denominator(cnt) mod p) == 0) { return(0); } 398 for(j = size(f); j > 0; j) 399 { 400 if((leadcoef(f[j]) mod p) == 0) { return(0); } 384 401 } 385 402 } … … 390 407 391 408 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 409 "USAGE: primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal, 410 n integer 411 RETURN: the intvec of n greatest primes <= 2147483647 (resp. n greatest primes 394 412 < L[size(L)] union with L) such that none of these primes divides any 395 413 coefficient occuring in I 414 NOTE: The number of cores to use can be defined by ncores, default is 1. 396 415 EXAMPLE: example primList; shows an example 397 416 " … … 399 418 intvec L; 400 419 int i,p; 420 int ncores = 1; 421 422 // Initialize optional parameter ncores  423 if(size(#) > 0) 424 { 425 if(size(#) == 1) 426 { 427 if(typeof(#[1]) == "int") 428 { 429 ncores = #[1]; 430 # = list(); 431 } 432 } 433 else 434 { 435 ncores = #[2]; 436 } 437 } 438 401 439 if(size(#) == 0) 402 440 { … … 421 459 } 422 460 if(p == 2) { ERROR("no more primes"); } 423 for(i = 2; i <= n; i++) 424 { 425 p = prime(p1); 426 while(!primeTest(I,p)) 461 if(ncores == 1) 462 { 463 for(i = 2; i <= n; i++) 427 464 { 428 465 p = prime(p1); 429 if(p == 2) { ERROR("no more primes"); } 430 } 431 L[size(L)+1] = p; 466 while(!primeTest(I,p)) 467 { 468 p = prime(p1); 469 if(p == 2) { ERROR("no more primes"); } 470 } 471 L[size(L)+1] = p; 472 } 473 } 474 else 475 { 476 int neededSize = size(L)+n1;; 477 list parallelResults; 478 list arguments; 479 int neededPrimes = neededSizesize(L); 480 while(neededPrimes > 0) 481 { 482 arguments = list(); 483 for(i = ((neededPrimes div ncores)+1(neededPrimes%ncores == 0)) 484 *ncores; i > 0; i) 485 { 486 p = prime(p1); 487 if(p == 2) { ERROR("no more primes"); } 488 arguments[i] = list("I", p); 489 } 490 parallelResults = parallelWaitAll("primeTest", arguments, 491 list(list(list(ncores)))); 492 for(i = size(arguments); i > 0; i) 493 { 494 if(parallelResults[i]) 495 { 496 L[size(L)+1] = arguments[i][2]; 497 } 498 } 499 neededPrimes = neededSizesize(L); 500 } 501 if(size(L) > neededSize) 502 { 503 L = L[1..neededSize]; 504 } 432 505 } 433 506 return(L); … … 470 543 ideal I = imap(R,I); 471 544 545 intvec opt = option(get); 472 546 option(none); 473 547 option(prompt); … … 504 578 505 579 setring R; 580 option(set, opt); 506 581 return(imap(newR,sI),imap(newR,T)); 507 582 } … … 537 612 The standard basis computation modulo p does also vary depending on the 538 613 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.. 614 @*  variant = 1: std(.,#[1]) resp. groebner, 615 @*  variant = 2: groebner, 616 @*  variant = 3: homog.  std(.,#[1]) resp. groebner  dehomog., 617 @*  variant = 4: fglm. 542 618 EXAMPLE: example modpStd; shows an example 543 619 " … … 552 628 option(redSB); 553 629 554 int t = timer; 555 if((variant == 1)  (variant == 4)) 630 if(variant == 1) 556 631 { 557 632 if(size(#) > 0) … … 565 640 } 566 641 567 if( (variant == 2)  (variant == 5))642 if(variant == 2) 568 643 { 569 644 i = groebner(i); 570 645 } 571 646 572 if( (variant == 3)  (variant == 6))647 if(variant == 3) 573 648 { 574 649 list rl = ringlist(@r); … … 607 682 } 608 683 609 t = timer;610 684 i = subst(i, homvar, 1); 611 685 i = simplify(i, 34); … … 615 689 i = interred(i); 616 690 kill HomR; 691 } 692 693 if(variant == 4) 694 { 695 def R1 = changeord("dp"); 696 setring R1; 697 ideal i = fetch(@r,i); 698 i = std(i); 699 setring @r; 700 i = fglm(R1,i); 617 701 } 618 702 … … 646 730 RETURN: a standard basis of I if no warning appears; 647 731 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. 732 numbers) by using modular methods. 651 733 By default the procedure computes a standard basis of I for sure, but 652 734 if the optional parameter #[2] = 0, it computes a standard basis of I 653 with high probability and a warning appears at the end.735 with high probability. 654 736 The procedure distinguishes between different variants for the standard 655 737 basis computation in positive characteristic depending on the ordering 656 738 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, 739 @*  variant = 1, if I is homogeneous, 740 @*  variant = 2, if I is not homogeneous, 1blockordering, 660 741 @*  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. 742 > 1 block), 743 @*  variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0. 667 744 EXAMPLE: example modStd; shows an example 668 745 " … … 752 829 if(printlevel >= 10) 753 830 { 754 "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3); 831 "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3) 832 +", exactness = "+string(exactness); 755 833 } 756 834 … … 761 839 762 840 // Initialize the list of primes  763 intvec L = primeList(I,n2); 841 int tt = timer; 842 int rt = rtimer; 843 intvec L = primeList(I,n2,n1); 844 if(printlevel >= 10) 845 { 846 "CPUtime for primeList: "+string(timertt)+" seconds."; 847 "Realtime for primeList: "+string(rtimerrt)+" seconds."; 848 } 764 849 L[5] = prime(random(an,en)); 765 850 … … 768 853 int h = homog(I); 769 854 770 inttt = timer;771 intrt = rtimer;855 tt = timer; 856 rt = rtimer; 772 857 773 858 if(!mixedTest()) … … 775 860 if(h) 776 861 { 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 } 862 variant = 1; 863 if(printlevel >= 10) { "variant = 1"; } 864 787 865 rl[1] = L[5]; 788 866 def @r = ring(rl); … … 802 880 if((find(ordstr_R0, "M") > 0)  (find(ordstr_R0, "a") > 0)  neg) 803 881 { 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 } 882 variant = 2; 883 if(printlevel >= 10) { "variant = 2"; } 814 884 } 815 885 else … … 827 897 if((order == "simple")  (size(rl) > 4)) 828 898 { 829 if(exactness == 0) 899 variant = 2; 900 if(printlevel >= 10) { "variant = 2"; } 901 } 902 else 903 { 904 rl[1] = L[5]; 905 def @r = ring(rl); 906 setring @r; 907 908 def @s = changeord("dp"); 909 setring @s; 910 ideal I = std(fetch(R0,I)); 911 if(dim(I) == 0) 830 912 { 831 variant = 2;832 if(printlevel >= 10) { "variant = 2"; }913 variant = 4; 914 if(printlevel >= 10) { "variant = 4"; } 833 915 } 834 if(exactness == 1) 835 { 836 variant = 5; 837 if(printlevel >= 10) { "variant = 5"; } 838 } 839 } 840 else 841 { 842 if(exactness == 0) 916 else 843 917 { 844 918 variant = 3; 845 919 if(printlevel >= 10) { "variant = 3"; } 920 921 int nvar@r = nvars(@r); 922 intvec w; 923 for(i = 1; i <= nvar@r; i++) 924 { 925 w[i] = deg(var(i)); 926 } 927 w[nvar@r + 1] = 1; 928 929 list hiRi = hilbRing(fetch(R0,I),w); 930 intvec W = hiRi[2]; 931 @s = hiRi[1]; 932 setring @s; 933 934 Id(1) = std(Id(1)); 935 intvec hi = hilb(Id(1), 1, W); 846 936 } 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 937 872 938 setring R0; 873 939 kill @r,@s; … … 971 1037 link l(i) = "ssi:fork"; 972 1038 open(l(i)); 973 if((variant == 1)  (variant == 3)  974 (variant == 4)  (variant == 6)) 1039 if((variant == 1)  (variant == 3)) 975 1040 { 976 1041 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), 977 1042 eval(variant), eval(hi)))); 978 1043 } 979 if((variant == 2)  (variant == 5))1044 if((variant == 2)  (variant == 4)) 980 1045 { 981 1046 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), … … 985 1050 986 1051 int t = timer; 987 if((variant == 1)  (variant == 3)  (variant == 4)  (variant == 6))1052 if((variant == 1)  (variant == 3)) 988 1053 { 989 1054 P = modpStd(I_for_fork, L[1], variant, hi); 990 1055 } 991 if((variant == 2)  (variant == 5))1056 if((variant == 2)  (variant == 4)) 992 1057 { 993 1058 P = modpStd(I_for_fork, L[1], variant); … … 1005 1070 // Main standard basis computations in positive  1006 1071 // characteristic start here  1072 1073 list arguments_farey, results_farey; 1007 1074 1008 1075 while(1) … … 1029 1096 if(j <= size(L)) 1030 1097 { 1031 if((variant == 1)  (variant == 3)  1032 (variant == 4)  (variant == 6)) 1098 if((variant == 1)  (variant == 3)) 1033 1099 { 1034 1100 write(l(i), quote(modpStd(I_for_fork, eval(L[j]), … … 1036 1102 j++; 1037 1103 } 1038 if((variant == 2)  (variant == 5))1104 if((variant == 2)  (variant == 4)) 1039 1105 { 1040 1106 write(l(i), quote(modpStd(I_for_fork, … … 1062 1128 while(j <= size(L)) 1063 1129 { 1064 if((variant == 1)  (variant == 3)  1065 (variant == 4)  (variant == 6)) 1130 if((variant == 1)  (variant == 3)) 1066 1131 { 1067 1132 P = modpStd(I, L[j], variant, hi); 1068 1133 } 1069 if((variant == 2)  (variant == 5))1134 if((variant == 2)  (variant == 4)) 1070 1135 { 1071 1136 P = modpStd(I, L[j], variant); … … 1095 1160 // Lift results to basering via farey  1096 1161 1097 tt = timer; 1162 tt = timer; rt = rtimer; 1098 1163 N = T2[1]; 1099 1164 for(i = 2; i <= size(T2); i++) { N = N*T2[i]; } 1100 1165 H = chinrem(T1,T2); 1101 J = farey(H,N); 1102 if(printlevel >= 10) { "Liftingprocess takes "+string(timer  tt) 1103 +" seconds"; } 1166 if(n1 == 1) 1167 { 1168 J = farey(H,N); 1169 } 1170 else 1171 { 1172 for(i = ncols(H); i > 0; i) 1173 { 1174 arguments_farey[i] = list(ideal(H[i]), N); 1175 } 1176 results_farey = parallelWaitAll("farey", arguments_farey, 1177 list(list(list(n1)))); 1178 for(i = ncols(H); i > 0; i) 1179 { 1180 J[i] = results_farey[i][1]; 1181 } 1182 } 1183 if(printlevel >= 10) 1184 { 1185 "CPUtime for liftingprocess is "+string(timer  tt)+" seconds."; 1186 "Realtime for liftingprocess is "+string(rtimer  rt)+" seconds."; 1187 } 1104 1188 1105 1189 // Test if we already have a standard basis of I  1106 1190 1107 1191 tt = timer; rt = rtimer; 1108 if((variant == 1)  (variant == 3)  (variant == 4)  (variant == 6))1192 if((variant == 1)  (variant == 3)) 1109 1193 { 1110 1194 pTest = pTestSB(I,J,L,variant,hi); 1111 1195 } 1112 if((variant == 2)  (variant == 5))1196 if((variant == 2)  (variant == 4)) 1113 1197 { 1114 1198 pTest = pTestSB(I,J,L,variant); … … 1132 1216 1133 1217 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)) 1218 1219 if(exactness == 0) 1220 { 1221 option(set, opt); 1222 if(n1 > 1) { kill I_for_fork; } 1223 return(J); 1224 } 1225 1226 if(exactness == 1) 1227 { 1228 tt = timer; rt = rtimer; 1229 sizeTest = 1  isIncluded(I,J,n1); 1230 1231 if(printlevel >= 10) 1232 { 1233 "CPUtime for checking if I subset <G> is " 1234 +string(timer  tt)+" seconds."; 1235 "Realtime for checking if I subset <G> is " 1236 +string(rtimer  rt)+" seconds."; 1237 } 1238 1239 if(sizeTest == 0) 1157 1240 { 1158 1241 tt = timer; rt = rtimer; … … 1185 1268 1186 1269 j = size(L) + 1; 1187 L = primeList(I,n3,L); 1270 tt = timer; rt = rtimer; 1271 L = primeList(I,n3,L,n1); 1272 if(printlevel >= 10) 1273 { 1274 "CPUtime for primeList: "+string(timertt)+" seconds."; 1275 "Realtime for primeList: "+string(rtimerrt)+" seconds."; 1276 } 1188 1277 1189 1278 if(n1 > 1) … … 1192 1281 { 1193 1282 open(l(i)); 1194 if((variant == 1)  (variant == 3)  1195 (variant == 4)  (variant == 6)) 1283 if((variant == 1)  (variant == 3)) 1196 1284 { 1197 1285 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i1]), 1198 1286 eval(variant), eval(hi)))); 1199 1287 } 1200 if((variant == 2)  (variant == 5))1288 if((variant == 2)  (variant == 4)) 1201 1289 { 1202 1290 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i1]), … … 1318 1406 def S = ring(rl); 1319 1407 setring S; 1408 intvec opt = option(get); 1320 1409 option(redSB); 1321 1410 module Z,M,Z2; … … 1362 1451 if(size(reduce(G3,testG1)) == 0) 1363 1452 { 1453 option(set, opt); 1364 1454 return(G3); 1365 1455 }
Note: See TracChangeset
for help on using the changeset viewer.