Changeset 36b81ac in git for Singular/LIB/modwalk.lib
- Timestamp:
- Sep 1, 2015, 3:13:14 PM (9 years ago)
- Branches:
- (u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
- Children:
- 510257fc6c7ff9c0c28e81f62f6fe06e63f024f8b430ca2b1337f52e1932f689e451e9abb002ec4a
- Parents:
- 74fe3587102ed73b3858258975abbe6d63139abfcd38fbd954495900afc188a0e871f8586db00535
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/modwalk.lib
-
Property
mode
changed from
100644
to100755
rcd38fbd r36b81ac 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
Note: See TracChangeset
for help on using the changeset viewer.