Changeset 604d8b in git
 Timestamp:
 Jun 10, 2010, 3:02:22 PM (13 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
 Children:
 b41ecb2f02b86cc0b679c2ecb470330dae622a51
 Parents:
 6e0750067834212a8910aad022e669ca646b5786
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/modstd.lib
r6e0750 r604d8b 1 1 //GP, last modified 23.10.06 2 /////////////////////////////////////////////////////////////////////////////// 2 //////////////////////////////////////////////////////////////////////////////// 3 3 version="$Id$"; 4 category ="Commutative Algebra";4 category = "Commutative Algebra"; 5 5 info=" 6 LIBRARY: modstd.lib Grobner basis of ideals7 AUTHORS: A. Hashemi, Amir.Hashemi@lip6.fr 8 @* G. Pfister pfister@mathematik.unikl.de 9 @* H. Schoenemann hannes@mathematik.unikl.de10 11 NOTE: 12 A library for computing the Grobner basis of an ideal in the polynomial 13 ring over the rational numbers using modular methods. The procedures are 14 inspired by the following paper: 15 Elizabeth A. Arnold:16 Modular Algorithms for Computing Groebner Bases , Journal of Symbolic17 Computation , April 2003, Volume 35, (4), p. 403419.18 19 6 LIBRARY: modstd.lib Groebner basis of ideals 7 8 AUTHORS: A. Hashemi Amir.Hashemi@lip6.fr 9 @* G. Pfister pfister@mathematik.unikl.de 10 @* H. Schoenemann hannes@mathematik.unikl.de 11 @* S. Steidel steidel@mathematik.unikl.de 12 13 OVERVIEW: 14 15 A library for computing the Groebner basis of an ideal in the polynomial 16 ring over the rational numbers using modular methods. The procedures are 17 inspired by the following paper: 18 Elizabeth A. Arnold: Modular algorithms for computing Groebner bases. 19 Journal of Symbolic Computation 35, 403419 (2003). 20 20 21 21 PROCEDURES: 22 modStd(I); compute a standard basis of I using modular methods 23 modS(I,L); liftings to Q of standard bases of I mod p for p in L 24 primeList(n); intvec of n primes <= 2134567879 in decreasing order 22 modStd(I); standard basis of I using modular methods (chinese remainder) 23 modHenselStd(I); standard basis of I using modular methods (hensel lifting) 24 modS(I,L); liftings to Q of standard bases of I mod p for p in L 25 25 "; 26 26 27 27 LIB "poly.lib"; 28 LIB "crypto.lib"; 29 /////////////////////////////////////////////////////////////////////////////// 30 proc pTestSB(ideal I, ideal J, list L) 31 "USAGE: pTestSB(I,J,L); I,J ideals, L intvec of primes 32 RETURN: 1 (resp. 0) if for a randomly chosen prime p not in L 33 J mod p is (resp. is not) a standard basis of I mod p 34 EXAMPLE:example primList; shows an example 28 LIB "ring.lib"; 29 30 //////////////////////////////////////////////////////////////////////////////// 31 32 static proc redFork(ideal I, ideal J, int n) 33 { 34 attrib(J,"isSB",1); 35 return(reduce(I,J,1)); 36 } 37 38 //////////////////////////////////////////////////////////////////////////////// 39 40 proc isIncluded(ideal I, ideal J, list #) 41 "USAGE: isIncluded(I,J); I,J ideals 42 RETURN: 1 if J includes I, 43 0 if there is an element f in I which does not reduce to 0 w.r.t. J. 44 EXAMPLE: example isIncluded; shows an example 35 45 " 36 46 { 37 int i,j,k,p; 38 def R=basering; 39 list r= ringlist(R); 40 41 while(!j) 42 { 43 j=1; 44 p=prime(random(1000000000,2134567879)); 45 for(i=1;i<=size(L);i++) 46 { 47 if(p==L[i]){j=0;break} 48 } 49 if(j) 50 { 51 for(i=1;i<=ncols(J);i++) 52 { 53 for(k=2;k<=size(J[i]);k++) 54 { 55 if((denominator(leadcoef(J[i][k])) mod 56 p)==0){j=0;break} 57 } 58 if(!j){break;} 59 } 60 } 61 } 62 r[1]=p; 63 def @R=ring(r); 64 setring @R; 65 ideal I=imap(R,I); 66 ideal J=imap(R,J); 67 attrib(J,"isSB",1); 68 j=1; 69 if(size(reduce(I,J))!=0){j=0;} 70 if(j) 71 { 72 ideal K=std(J); 73 if(size(reduce(K,J))!=0){j=0;} 74 } 75 setring R; 76 return(j); 47 attrib(J,"isSB",1); 48 int i,j,k; 49 50 if(size(#) > 0) 51 { 52 int n = #[1]; 53 if((n > 1) && (n < ncols(I)) && system("with","MP")) 54 { 55 for(i = 1; i <= n  1; i++) 56 { 57 link l(i) = "MPtcp:fork"; 58 open(l(i)); 59 write(l(i), quote(redFork(eval(I[ncols(I)i]), eval(J), 1))); 60 } 61 62 int t = timer; 63 if(reduce(I[ncols(I)], J, 1) != 0) 64 { 65 for(i = 1; i <= n  1; i++) 66 { 67 close(l(i)); 68 } 69 return(0); 70 } 71 t = timer  t; 72 if(t > 60) { t = 60; } 73 int i_sleep = system("sh", "sleep "+string(t)); 74 75 j = ncols(I)  n; 76 77 while(j >= 0) 78 { 79 for(i = 1; i <= n  1; i++) 80 { 81 if(status(l(i), "read", "ready")) 82 { 83 if(read(l(i)) != 0) 84 { 85 for(i = 1; i <= n  1; i++) 86 { 87 close(l(i)); 88 } 89 return(0); 90 } 91 else 92 { 93 if(j >= 1) 94 { 95 write(l(i), quote(redFork(eval(I[j]), eval(J), 1))); 96 j; 97 } 98 else 99 { 100 k++; 101 close(l(i)); 102 } 103 } 104 } 105 } 106 if(k == n  1) 107 { 108 j; 109 } 110 i_sleep = system("sh", "sleep "+string(t)); 111 } 112 return(1); 113 } 114 } 115 116 for(i = ncols(I); i >= 1; i) 117 { 118 if(reduce(I[i],J,1) != 0){ return(0); } 119 } 120 return(1); 77 121 } 78 122 example 79 123 { "EXAMPLE:"; echo = 2; 80 intvec L=2,3,5;81 124 ring r=0,(x,y,z),dp; 82 ideal i=x+1,x+y+1; 83 ideal j=x+1,y; 84 pTestSB(i,i,L); 85 pTestSB(i,j,L); 86 } 87 88 proc primeList(int n, list #) 89 "USAGE: primeList(n); (resp. primeList(n,L); ) 90 RETURN: the intvec of n greatest primes <= 2147483647 (resp. n greatest 91 primes < L[size(L)] union with L) 92 EXAMPLE:example primList; shows an example 125 ideal I = x+1,x+y+1; 126 ideal J = x+1,y; 127 isIncluded(I,J); 128 isIncluded(J,I); 129 isIncluded(I,J,4); 130 131 ring R = 0, x(1..5), dp; 132 ideal I1 = cyclic(4); 133 ideal I2 = I1,x(5)^2; 134 isIncluded(I1,I2,4); 135 } 136 137 //////////////////////////////////////////////////////////////////////////////// 138 139 proc pTestSB(ideal I, ideal J, list L, int variant, list #) 140 "USAGE: pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int 141 RETURN: 1 (resp. 0) if for a randomly chosen prime p that is not in L 142 J mod p is (resp. is not) a standard basis of I mod p 143 EXAMPLE: example pTestSB; shows an example 93 144 " 94 145 { 95 intvec L; 96 int i,p; 97 if(size(#)==0) 98 { 99 p=2147483647; 100 L[1]=p; 101 } 102 else 103 { 104 L=#[1]; 105 p=prime(L[size(L)]1); 106 L[size(L)+1]=p; 107 } 108 if(p==2){ERROR("no more primes");} 109 for(i=2;i<=n;i++) 110 { 111 p=prime(p1); 112 L[size(L)+1]=p; 113 } 114 return(L); 146 int i,j,k,p; 147 def R = basering; 148 list r = ringlist(R); 149 150 while(!j) 151 { 152 j = 1; 153 p = prime(random(1000000000,2134567879)); 154 for(i = 1; i <= size(L); i++) 155 { 156 if(p == L[i]) { j = 0; break; } 157 } 158 if(j) 159 { 160 for(i = 1; i <= ncols(I); i++) 161 { 162 for(k = 2; k <= size(I[i]); k++) 163 { 164 if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; } 165 } 166 if(!j){ break; } 167 } 168 } 169 if(j) 170 { 171 if(!primeTest(I,p)) { j = 0; } 172 } 173 } 174 r[1] = p; 175 def @R = ring(r); 176 setring @R; 177 ideal I = imap(R,I); 178 ideal J = imap(R,J); 179 attrib(J,"isSB",1); 180 181 int t = timer; 182 j = 1; 183 if(isIncluded(I,J) == 0) { j = 0; } 184 185 if(printlevel >= 10) 186 { 187 "isIncluded(I,J) takes "+string(timer  t)+" seconds"; 188 "j = "+string(j); 189 } 190 191 t = timer; 192 if(j) 193 { 194 if(size(#) > 0) 195 { 196 ideal K = modpStd(I,p,variant,#[1])[1]; 197 } 198 else 199 { 200 ideal K = groebner(I); 201 } 202 t = timer; 203 if(isIncluded(J,K) == 0) { j = 0; } 204 205 if(printlevel >= 10) 206 { 207 "isIncluded(K,J) takes "+string(timer  t)+" seconds"; 208 "j = "+string(j); 209 } 210 } 211 setring R; 212 return(j); 115 213 } 116 214 example 117 215 { "EXAMPLE:"; echo = 2; 118 intvec L=primeList(10); 216 intvec L = 2,3,5; 217 ring r = 0,(x,y,z),dp; 218 ideal I = x+1,x+y+1; 219 ideal J = x+1,y; 220 pTestSB(I,I,L,2); 221 pTestSB(I,J,L,2); 222 } 223 224 //////////////////////////////////////////////////////////////////////////////// 225 226 proc deleteUnluckyPrimes(list T, list L, int ho, list #) 227 "USAGE: deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer 228 RETURN: lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #), lT ideal 229 NOTE:  if ho = 1, the polynomials in T are homogeneous, else ho = 0, 230  lT is prevalent, i.e. the most appearing leading ideal in T 231 EXAMPLE: example deleteUnluckyPrimes; shows an example 232 " 233 { 234 ho = ((ho)(ord_test(basering) == 1)); 235 int j,k,c; 236 intvec hl,hc; 237 ideal cT,lT,cK; 238 lT = lead(T[size(T)]); 239 attrib(lT,"isSB",1); 240 if(!ho) 241 { 242 for(j = 1; j < size(T); j++) 243 { 244 cT = lead(T[j]); 245 attrib(cT,"isSB",1); 246 if((size(reduce(cT,lT))!=0)(size(reduce(lT,cT))!=0)) 247 { 248 cK = cT; 249 c++; 250 } 251 } 252 if(c > size(T)/2){ lT = cK; } 253 } 254 else 255 { 256 hl = hilb(lT,1); 257 for(j = 1; j < size(T); j++) 258 { 259 cT = lead(T[j]); 260 attrib(cT,"isSB",1); 261 hc = hilb(cT,1); 262 if(hl == hc) 263 { 264 for(k = 1; k <= size(lT); k++) 265 { 266 if(lT[k] < cT[k]) { lT = cT; c++; break; } 267 if(lT[k] > cT[k]) { c++; break; } 268 } 269 } 270 else 271 { 272 if(hc < hl){ lT = cT; hl = hilb(lT,1); c++ } 273 } 274 } 275 } 276 277 int addList; 278 if(size(#) > 0) { list M = #; addList = 1; } 279 j = 1; 280 attrib(lT,"isSB",1); 281 while((j <= size(T))&&(c > 0)) 282 { 283 cT = lead(T[j]); 284 attrib(cT,"isSB",1); 285 if((size(reduce(cT,lT)) != 0)(size(reduce(lT,cT)) != 0)) 286 { 287 T = delete(T,j); 288 if(j == 1) 289 { 290 L = L[2..size(L)]; 291 if(addList == 1) { M = M[2..size(M)]; } 292 } 293 else 294 { 295 if(j == size(L)) 296 { 297 L = L[1..size(L)1]; 298 if(addList == 1) { M = M[1..size(M)1]; } 299 } 300 else 301 { 302 L = L[1..j1],L[j+1..size(L)]; 303 if(addList == 1) { M = M[1..j1],M[j+1..size(M)]; } 304 } 305 } 306 j; 307 } 308 j++; 309 } 310 311 for(j = 1; j <= size(L); j++) 312 { 313 L[j] = bigint(L[j]); 314 } 315 316 if(addList == 0) { return(list(T,L,lT)); } 317 if(addList == 1) { return(list(T,L,M,lT)); } 318 } 319 example 320 { "EXAMPLE:"; echo = 2; 321 list L = 2,3,5,7,11; 322 ring r = 0,(y,x),Dp; 323 ideal I1 = 2y2x,y6; 324 ideal I2 = yx2,y3x,x5,y6; 325 ideal I3 = y2x,x3y,x5,y6; 326 ideal I4 = y2x,11x3y,x5; 327 ideal I5 = y2x,yx3,x5,7y6; 328 list T = I1,I2,I3,I4,I5; 329 deleteUnluckyPrimes(T,L,1); 330 list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5); 331 deleteUnluckyPrimes(T,L,1,P); 332 } 333 334 //////////////////////////////////////////////////////////////////////////////// 335 336 proc primeTest(ideal I, bigint p) 337 { 338 int i,j; 339 for(i = 1; i <= size(I); i++) 340 { 341 for(j = 1; j <= size(I[i]); j++) 342 { 343 if((leadcoef(I[i][j]) mod p) == 0) { return(0); } 344 } 345 } 346 return(1); 347 } 348 349 //////////////////////////////////////////////////////////////////////////////// 350 351 proc primeList(ideal I, int n, list #) 352 "USAGE: primeList(I,n); ( resp. primeList(I,n,L); ) I ideal, n integer 353 RETURN: the intvec of n greatest primes <= 2147483647 (resp. n greatest primes 354 < L[size(L)] union with L) such that none ot these primes divides any 355 coefficient occuring in I 356 EXAMPLE: example primList; shows an example 357 " 358 { 359 intvec L; 360 int i,p; 361 if(size(#) == 0) 362 { 363 p = 2147483647; 364 while(!primeTest(I,p)) 365 { 366 p = prime(p1); 367 if(p == 2) { ERROR("no more primes"); } 368 } 369 L[1] = p; 370 } 371 else 372 { 373 L = #[1]; 374 p = prime(L[size(L)]1); 375 while(!primeTest(I,p)) 376 { 377 p = prime(p1); 378 if(p == 2) { ERROR("no more primes"); } 379 } 380 L[size(L)+1] = p; 381 } 382 if(p == 2) { ERROR("no more primes"); } 383 for(i = 2; i <= n; i++) 384 { 385 p = prime(p1); 386 while(!primeTest(I,p)) 387 { 388 p = prime(p1); 389 if(p == 2) { ERROR("no more primes"); } 390 } 391 L[size(L)+1] = p; 392 } 393 return(L); 394 } 395 example 396 { "EXAMPLE:"; echo = 2; 397 ring r = 0,(x,y,z),dp; 398 ideal I = 2147483647x+y, z181; 399 intvec L = primeList(I,10); 400 size(L); 401 L[1]; 402 L[size(L)]; 403 L = primeList(I,5,L); 119 404 size(L); 120 405 L[size(L)]; 121 L=primeList(5,L); 122 size(L); 123 L[size(L)]; 124 } 125 126 127 proc modStd(ideal I) 128 "USAGE: modStd(I); 406 } 407 408 //////////////////////////////////////////////////////////////////////////////// 409 410 proc liftstd1(ideal I) 411 { 412 def R = basering; 413 list rl = ringlist(R); 414 list ordl = rl[3]; 415 416 int i; 417 for(i = 1; i <= size(ordl); i++) 418 { 419 if((ordl[i][1] == "C")  (ordl[i][1] == "c")) 420 { 421 ordl = delete(ordl, i); 422 break; 423 } 424 } 425 426 ordl = insert(ordl, list("c", 0)); 427 rl[3] = ordl; 428 def newR = ring(rl); 429 setring newR; 430 ideal I = imap(R,I); 431 432 option(none); 433 option(prompt); 434 435 module M; 436 for(i = 1; i <= size(I); i++) 437 { 438 M = M + module(I[i]*gen(1) + gen(i+1)); 439 M = M + module(gen(i+1)); 440 } 441 442 module sM = std(M); 443 444 ideal sI; 445 if(attrib(R,"global")) 446 { 447 for(i = size(I)+1; i <= size(sM); i++) 448 { 449 sI[size(sI)+1] = sM[i][1]; 450 } 451 matrix T = submat(sM,2..nrows(sM),size(I)+1..ncols(sM)); 452 } 453 else 454 { 455 "=========================================================="; 456 "WARNING: Algorithm is not applicable if ordering is mixed."; 457 "=========================================================="; 458 for(i = 1; i <= size(sM)size(I); i++) 459 { 460 sI[size(sI)+1] = sM[i][1]; 461 } 462 matrix T = submat(sM,2..nrows(sM),1..ncols(sM)size(I)); 463 } 464 465 setring R; 466 return(imap(newR,sI),imap(newR,T)); 467 } 468 example 469 { "EXAMPLE:"; echo = 2; 470 ring R = 0,(x,y,z),dp; 471 poly f = x3+y7+z2+xyz; 472 ideal i = jacob(f); 473 matrix T; 474 ideal sm = liftstd(i,T); 475 sm; 476 print(T); 477 matrix(sm)  matrix(i)*T; 478 479 480 ring S = 32003, x(1..6), lp; 481 ideal I = cyclic(6); 482 ideal sI; 483 matrix T; 484 sI,T = liftstd1(I); 485 matrix(sI)  matrix(I)*T; 486 } 487 488 //////////////////////////////////////////////////////////////////////////////// 489 490 proc modpStd(ideal I, int p, int variant, list #) 491 "USAGE: modpStd(I,p,variant,#); I ideal, p integer, variant integer 492 ASSUME: If size(#) > 0, then #[1] is an intvec describing the Hilbert series. 493 RETURN: ideal  a standard basis of I mod p, integer  p (, matrix  if variant 494 = 5, the transformation matrix obtained from liftstd) 495 NOTE: The procedure computes a standard basis of the ideal I modulo p and 496 fetches the result to the basering. If size(#) > 0 the Hilbert driven 497 standard basis computation std(.,#[1]) is used instead of groebner. 498 The standard basis computation modulo p does also vary depending on the 499 integer variant, namely 500 @*  variant = 1: std(.,#[1]) resp. groebner, 501 @*  variant = 2: groebner, 502 @*  variant = 3: homogenize  std(.,#[1]) resp. groebner  dehomogenize, 503 @*  variant = 4: std(.,#[1]) resp. groebner, 504 @*  variant = 5: liftstd. 505 EXAMPLE: example modpStd; shows an example 506 " 507 { 508 def R0 = basering; 509 list rl = ringlist(R0); 510 rl[1] = p; 511 def @r = ring(rl); 512 setring @r; 513 ideal i = fetch(R0,I); 514 515 option(redSB); 516 517 int t = timer; 518 if((variant == 1)  (variant == 4)) 519 { 520 if(size(#) > 0) 521 { 522 i = std(i, #[1]); 523 } 524 else 525 { 526 i = groebner(i); 527 } 528 } 529 if(variant == 2) 530 { 531 i = groebner(i); 532 } 533 534 if(variant == 3) 535 { 536 list rl = ringlist(@r); 537 int nvar@r = nvars(@r); 538 539 int k; 540 intvec w; 541 for(k = 1; k <= nvar@r; k++) 542 { 543 w[k] = deg(var(k)); 544 } 545 w[nvar@r + 1] = 1; 546 547 rl[2][nvar@r + 1] = "homvar"; 548 rl[3][2][2] = w; 549 550 def HomR = ring(rl); 551 setring HomR; 552 ideal i = imap(@r, i); 553 i = homog(i, homvar); 554 555 if(size(#) > 0) 556 { 557 if(w == 1) 558 { 559 i = std(i, #[1]); 560 } 561 else 562 { 563 i = std(i, #[1], w); 564 } 565 } 566 else 567 { 568 i = groebner(i); 569 } 570 571 t = timer; 572 i = subst(i, homvar, 1); 573 i = simplify(i, 34); 574 575 setring @r; 576 i = imap(HomR, i); 577 i = interred(i); 578 kill HomR; 579 } 580 581 if(variant == 5) 582 { 583 matrix trans; 584 i,trans = liftstd1(i); 585 setring R0; 586 return(list(fetch(@r,i),p,fetch(@r,trans))); 587 } 588 589 setring R0; 590 return(list(fetch(@r,i),p)); 591 } 592 example 593 { "EXAMPLE:"; echo = 2; 594 ring r = 0, x(1..4), dp; 595 ideal I = cyclic(4); 596 int p = 181; 597 list P = modpStd(I,p,5); 598 P; 599 matrix(P[1])matrix(I)*P[3]; 600 601 int q = 32003; 602 list Q = modpStd(I,q,2); 603 Q; 604 } 605 606 ////////////////////////////// main procedures ///////////////////////////////// 607 608 proc modStd(ideal I, list #) 609 "USAGE: modStd(I); I ideal 610 ASSUME: If size(#) > 0, then # contains either 1, 2 or 4 integers such that 611 @*  #[1] is the number of available processors for the computation, 612 @*  #[2] is an optional parameter for the exactness of the computation, 613 if #[2] = 1, the procedure computes a standard basis for sure, 614 @*  #[3] is the number of primes until the first lifting, 615 @*  #[4] is the constant number of primes between two liftings until 616 the computation stops. 129 617 RETURN: a standard basis of I if no warning appears; 130 NOTE: the procedure computes a standard basis of I (over the 131 rational numbers) by using modular methods. If a 132 warning appears then the result is a standard basis 133 containing I and with high probability a standard basis of I. 134 For further experiments see procedure modS. 618 NOTE: The procedure computes a standard basis of I (over the rational 619 numbers) by using modular methods. If a warning appears then the 620 result is a standard basis containing I and with high probability 621 a standard basis of I. 622 By default the procedure computes a standard basis of I with high 623 probability, but if the optional parameter #[2] = 1, it is exact. 624 The procedure distinguishes between different variants for the standard 625 basis computation in positive characteristic depending on the ordering 626 of the basering, the parameter #[2] and if the ideal I is homogeneous. 627 @*  variant = 1, if I is homogeneous and exactness = 0, 628 @*  variant = 2, if I is not homogeneous, 1blockordering and 629 exactness = 0, 630 @*  variant = 3, if I is not homogeneous, complicated ordering (lp or 631 > 1 block) and exactness = 0, 632 @*  variant = 4, if I is homogeneous and exactness = 1, 633 @*  variant = 5, if I is not homogeneous and exactness = 1. 135 634 EXAMPLE: example modStd; shows an example 136 635 " 137 636 { 138 def R0=basering; 139 list rl=ringlist(R0); 140 if((npars(R0)>0)(rl[1]>0)) 141 { 142 ERROR("characteristic of basering should be zero, basering should have 143 no parameters"); 144 } 145 146 int k,c; 147 int pd=printlevelvoice+2; 148 int j=1; 149 int h=homog(I); 150 int en=2134567879; 151 int an=1000000000; 152 153 intvec opt = option(get); // Save current options 154 intvec L=primeList(10); 155 L[5]=prime(random(an,en)); 156 list T,TT,TH,LL; 157 158 159 ideal J,K; 160 161 option(redSB); 162 163 while(1) 164 { 165 while(j<=size(L)) 166 { 167 if(pd>2){c++;c;} 168 rl[1]=L[j]; 169 def @r=ring(rl); 170 setring @r; 171 ideal i=fetch(R0,I); 172 i=groebner(i); 173 setring R0; 174 T[size(T)+1]=fetch(@r,i); 175 kill @r; 176 j++; 177 } 178 if(pd>2){"lifting";} 179 //================= delete unlucky primes ==================== 180 // unlucky if and only if the leading ideal is wrong 181 LL=deleteUnluckyPrimes(T,L); 182 TH=LL[1]; 183 L=LL[2]; 184 //============ now all leading ideals are the same ============ 185 for(j=1;j<=ncols(TH[1]);j++) 186 { 187 for(k=1;k<=size(L);k++) 188 { 189 TT[k]=TH[k][j]; 190 } 191 J[j]=liftPoly(TT,L); 192 } 193 if(pd>2){"list of primes:";L;"pTest";} 194 if(pTestSB(I,J,L)) 195 { 196 attrib(J,"isSB",1); 197 K=std(J); 198 if(size(reduce(I,J))==0) 199 { 200 if(size(reduce(K,J))==0) 201 { 202 if(!((h)(ord_test(R0)==1))) 203 { 204 205 "=================================================================="; 206 " The input is not homogeneous and the ordering is not 207 local. "; 208 "WARNING: ideal generated by output may be greater then 209 input ideal"; 210 211 "=================================================================="; 212 } 213 option(set, opt); 214 return(J); 215 } 216 if(pd>2){"pTest o.k. but result wrong";} 637 int TT = timer; 638 int RT = rtimer; 639 640 def R0 = basering; 641 list rl = ringlist(R0); 642 if((npars(R0) > 0)  (rl[1] > 0)) 643 { 644 ERROR("characteristic of basering should be zero, basering should have no parameters"); 645 } 646 647 int index = 1; 648 int i,k,c; 649 int pd = printlevelvoice+2; 650 int j = 1; 651 int pTest; 652 int en = 2134567879; 653 int an = 1000000000; 654 bigint N; 655 656 // Initialize optional parameters  657 if(size(#) > 0) 658 { 659 if(size(#) == 1) 660 { 661 int n1 = #[1]; 662 if((n1 > 1) && (1  system("with","MP"))) 663 { 664 "========================================================================"; 665 "There is no MP available on your system. Since this is necessary to "; 666 "parallelize the algorithm, the computation will be done without forking."; 667 "========================================================================"; 668 n1 = 1; 669 } 670 int exactness = 0; 671 int n2 = 10; 672 int n3 = 10; 673 } 674 if(size(#) == 2) 675 { 676 int n1 = #[1]; 677 if((n1 > 1) && (1  system("with","MP"))) 678 { 679 "========================================================================"; 680 "There is no MP available on your system. Since this is necessary to "; 681 "parallelize the algorithm, the computation will be done without forking."; 682 "========================================================================"; 683 n1 = 1; 684 } 685 int exactness = #[2]; 686 int n2 = 10; 687 int n3 = 10; 688 } 689 if(size(#) == 4) 690 { 691 int n1 = #[1]; 692 if((n1 > 1) && (1  system("with","MP"))) 693 { 694 "========================================================================"; 695 "There is no MP available on your system. Since this is necessary to "; 696 "parallelize the algorithm, the computation will be done without forking."; 697 "========================================================================"; 698 n1 = 1; 699 } 700 int exactness = #[2]; 701 int n2 = #[3]; 702 int n3 = #[4]; 703 } 704 } 705 else 706 { 707 int n1 = 1; 708 int exactness = 0; 709 int n2 = 10; 710 int n3 = 10; 711 } 712 713 // Save current options  714 intvec opt = option(get); 715 716 option(redSB); 717 718 // Initialize the list of primes  719 intvec L = primeList(I,n2); 720 L[5] = prime(random(an,en)); 721 722 // Decide which variant to take  723 int variant; 724 int h = homog(I); 725 726 int tt = timer; 727 int rt = rtimer; 728 729 if(h) 730 { 731 if(exactness == 0) { variant = 1; if(printlevel >= 10) { "variant = 1"; } } 732 if(exactness == 1) { variant = 4; if(printlevel >= 10) { "variant = 4"; } } 733 rl[1] = L[5]; 734 def @r = ring(rl); 735 setring @r; 736 def @s = changeord("dp"); 737 setring @s; 738 ideal I = std(fetch(R0,I)); 739 intvec hi = hilb(I,1); 740 setring R0; 741 kill @r,@s; 742 } 743 else 744 { 745 if(exactness == 0) 746 { 747 string ordstr_R0 = ordstr(R0); 748 int neg = 1  attrib(R0,"global"); 749 750 if((find(ordstr_R0, "M") > 0)  (find(ordstr_R0, "a") > 0)  neg) 751 { 752 variant = 2; 753 if(printlevel >= 10) { "variant = 2"; } 754 } 755 else 756 { 757 string order; 758 if(system("nblocks") <= 2) 759 { 760 if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") + find(ordstr_R0, "rp") <= 0) 761 { 762 order = "simple"; 763 } 764 } 765 766 if((order == "simple")  (size(rl) > 4)) 767 { 768 variant = 2; 769 if(printlevel >= 10) { "variant = 2"; } 770 } 771 else 772 { 773 variant = 3; 774 if(printlevel >= 10) { "variant = 3"; } 775 rl[1] = L[5]; 776 def @r = ring(rl); 777 setring @r; 778 int nvar@r = nvars(@r); 779 intvec w; 780 for(i = 1; i <= nvar@r; i++) 781 { 782 w[i] = deg(var(i)); 783 } 784 w[nvar@r + 1] = 1; 785 786 list hiRi = hilbRing(fetch(R0,I),w); 787 intvec W = hiRi[2]; 788 def @s = hiRi[1]; 789 setring @s; 790 791 Id(1) = std(Id(1)); 792 intvec hi = hilb(Id(1), 1, W); 793 794 setring R0; 795 kill @r,@s; 796 } 797 } 798 } 799 800 if(exactness == 1) 801 { 802 variant = 5; 803 if(printlevel >= 10) { "variant = 2"; } 804 matrix trans; 805 } 806 } 807 808 list P,T1,T2,T3,LL; 809 810 ideal J,K,H; 811 812 // If there is more than one processor available, we parallelize the  813 // main standard basis computations in positive characteristic  814 815 if(n1 > 1) 816 { 817 ideal I_for_fork = I; 818 export(I_for_fork); // I available for each link 819 820 // Create n1 links l(1),...,l(n1), open all of them and compute  821 // standard basis for the primes L[2],...,L[n1 + 1].  822 823 for(i = 1; i <= n1; i++) 824 { 825 link l(i) = "MPtcp:fork"; 826 open(l(i)); 827 if((variant == 1)  (variant == 3)  (variant == 4)) 828 { 829 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), eval(variant), eval(hi)))); 830 } 831 if((variant == 2)  (variant == 5)) 832 { 833 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), eval(variant)))); 834 } 835 } 836 837 int t = timer; 838 if((variant == 1)  (variant == 3)  (variant == 4)) 839 { 840 P = modpStd(I_for_fork, L[1], variant, hi); 841 } 842 if((variant == 2)  (variant == 5)) 843 { 844 P = modpStd(I_for_fork, L[1], variant); 845 } 846 t = timer  t; 847 if(t > 60) { t = 60; } 848 int i_sleep = system("sh", "sleep "+string(t)); 849 T1[1] = P[1]; 850 T2[1] = bigint(P[2]); 851 if(variant == 5) { T3[1] = P[3]; } 852 index++; 853 854 j = j + n1 + 1; 855 } 856 857 // Main standard basis computations in positive  858 // characteristic start here  859 860 while(1) 861 { 862 tt = timer; rt = rtimer; 863 864 if(n1 > 1) 865 { 866 if(printlevel >= 10) { "size(L) = "+string(size(L)); } 867 while(j <= size(L) + 1) 868 { 869 for(i = 1; i <= n1; i++) 870 { 871 if(status(l(i), "read", "ready")) // ask if link l(i) is ready otherwise sleep for t seconds 872 { 873 P = read(l(i)); // read the result from l(i) 874 T1[index] = P[1]; 875 T2[index] = bigint(P[2]); 876 if(variant == 5) { T3[index] = P[3]; } 877 index++; 878 879 if(j <= size(L)) 880 { 881 if((variant == 1)  (variant == 3)  (variant == 4)) 882 { 883 write(l(i), quote(modpStd(I_for_fork, eval(L[j]), eval(variant), eval(hi)))); 884 j++; 885 } 886 if((variant == 2)  (variant == 5)) 887 { 888 write(l(i), quote(modpStd(I_for_fork, eval(L[j]), eval(variant)))); 889 j++; 890 } 891 } 892 else 893 { 894 k++; 895 close(l(i)); 896 } 897 } 898 } 899 if(k == n1) // k describes the number of closed links 900 { 901 j++; 902 } 903 i_sleep = system("sh", "sleep "+string(t)); 904 } 905 } 906 else 907 { 908 while(j <= size(L)) 909 { 910 if((variant == 1)  (variant == 3)  (variant == 4)) 911 { 912 P = modpStd(I, L[j], variant, hi); 913 } 914 if((variant == 2)  (variant == 5)) 915 { 916 P = modpStd(I, L[j], variant); 917 } 918 919 T1[index] = P[1]; 920 T2[index] = bigint(P[2]); 921 if(variant == 5) { T3[index] = P[3]; } 922 index++; 923 j++; 924 } 925 } 926 927 if(printlevel >= 10) 928 { 929 "CPUtime for computing list is "+string(timer  tt)+" seconds."; 930 "Realtime for computing list is "+string(rtimer  rt)+" seconds."; 931 } 932 933 if(pd>2){"lifting";} 934 935 // Delete unlucky primes  936 // unlucky if and only if the leading ideal is wrong  937 938 if(variant == 5) { LL = deleteUnluckyPrimes(T1,T2,h,T3); T3 = LL[3]; } 939 else { LL = deleteUnluckyPrimes(T1,T2,h); } 940 T1 = LL[1]; 941 T2 = LL[2]; 942 943 // Now all leading ideals are the same  944 // Lift results to basering via farey  945 946 N = T2[1]; 947 for(i = 2; i <= size(T2); i++){N = N*T2[i];} 948 H = chinrem(T1,T2); 949 J = farey(H,N); 950 if(variant == 5) { trans = farey(chinrem(T3,T2), N); } 951 952 // Test if we already have a standard basis of I  953 954 tt = timer; rt = rtimer; 955 if(pd > 2){ "list of primes:"; L; "pTest" ;} 956 if((variant == 1)  (variant == 3)  (variant == 4)) { pTest = pTestSB(I,J,L,variant,hi); } 957 if((variant == 2)  (variant == 5)) { pTest = pTestSB(I,J,L,variant); } 958 959 if(printlevel >= 10) 960 { 961 "CPUtime for pTest is "+string(timer  tt)+" seconds."; 962 "Realtime for pTest is "+string(rtimer  rt)+" seconds."; 963 } 964 965 if(pTest) 966 { 967 if(printlevel >= 10) 968 { 969 "CPUtime for computation without final tests is "+string(timer  TT)+" seconds."; 970 "Realtime for computation without final tests is "+string(rtimer  RT)+" seconds."; 971 } 972 973 attrib(J,"isSB",1); 974 tt = timer; rt = rtimer; 975 int sizeTest = 1  isIncluded(I,J,n1); 976 977 if(printlevel >= 10) 978 { 979 "CPUtime for checking if I subset <G> is "+string(timer  tt)+" seconds."; 980 "Realtime for checking if I subset <G> is "+string(rtimer  rt)+" seconds."; 981 } 982 983 if(sizeTest == 0) 984 { 985 if(variant == 1) 986 { 987 "==================================================================="; 988 "WARNING: Ideal generated by output may be greater than input ideal."; 989 "==================================================================="; 990 option(set, opt); 991 if(n1 > 1) { kill I_for_fork; } 992 return(J); 993 } 994 if((variant == 2)  (variant == 3)) 995 { 996 tt = timer; rt = rtimer; 997 K = std(J); 998 999 if(printlevel >= 10) 1000 { 1001 "CPUtime for last stdcomputation is "+string(timer  tt)+" seconds."; 1002 "Realtime for last stdcomputation is "+string(rtimer  rt)+" seconds."; 1003 } 1004 1005 if(size(reduce(K,J)) == 0) 1006 { 1007 "==================================================================="; 1008 "WARNING: Ideal generated by output may be greater than input ideal."; 1009 "==================================================================="; 1010 option(set, opt); 1011 if(n1 > 1) { kill I_for_fork; } 1012 return(J); 1013 } 1014 } 1015 if(variant == 4) 1016 { 1017 tt = timer; rt = rtimer; 1018 K = std(J); 1019 1020 if(printlevel >= 10) 1021 { 1022 "CPUtime for last stdcomputation is "+string(timer  tt)+" seconds."; 1023 "Realtime for last stdcomputation is "+string(rtimer  rt)+" seconds."; 1024 } 1025 1026 if(size(reduce(K,J)) == 0) 1027 { 1028 option(set, opt); 1029 if(n1 > 1) { kill I_for_fork; } 1030 return(J); 1031 } 1032 } 1033 if(variant == 5) 1034 { 1035 tt = timer; rt = rtimer; 1036 K = std(J); 1037 1038 if(printlevel >= 10) 1039 { 1040 "CPUtime for last stdcomputation is "+string(timer  tt)+" seconds."; 1041 "Realtime for last stdcomputation is "+string(rtimer  rt)+" seconds."; 1042 } 1043 1044 if(size(reduce(K,J)) == 0) 1045 { 1046 if(matrix(J) == matrix(I)*trans) 1047 { 1048 option(set, opt); 1049 if(n1 > 1) { kill I_for_fork; } 1050 return(J); 1051 } 1052 } 1053 } 1054 if(pd>2){"pTest o.k. but result wrong";} 217 1055 } 218 1056 if(pd>2){"pTest o.k. but result wrong";} 219 1057 } 220 j=size(L)+1; 221 L=primeList(10,L); 1058 1059 // We do not already have a standard basis of I  1060 // Therefore do the main computation for more primes  1061 1062 T1 = H; 1063 T2 = N; 1064 index = 2; 1065 1066 j = size(L) + 1; 1067 L = primeList(I,n3,L); 1068 1069 if(n1 > 1) 1070 { 1071 for(i = 1; i <= n1; i++) 1072 { 1073 open(l(i)); 1074 if((variant == 1)  (variant == 3)  (variant == 4)) 1075 { 1076 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i1]), eval(variant), eval(hi)))); 1077 } 1078 if((variant == 2)  (variant == 5)) 1079 { 1080 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i1]), eval(variant)))); 1081 } 1082 } 1083 j = j + n1; 1084 k = 0; 1085 } 222 1086 } 223 1087 } 224 1088 example 225 1089 { "EXAMPLE:"; echo = 2; 226 ring r =0,(x,y,z,t),dp;227 ideal I =3x3+x2+1,11y5+y3+2,5z4+z2+4;228 ideal J =modStd(I);1090 ring r = 0,(x,y,z,t),dp; 1091 ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; 1092 ideal J = modStd(I); 229 1093 J; 230 I =homog(I,t);231 J =modStd(I);1094 I = homog(I,t); 1095 J = modStd(I); 232 1096 J; 233 ring s=0,(x,y,z),ds; 234 ideal I=jacob(x5+y6+z7+xyz); 235 ideal J=modStd(I); 236 J; 237 } 238 /////////////////////////////////////////////////////////////////////////////// 239 proc modS(ideal I, intvec L, list #) 1097 1098 ring s = 0,(x,y,z),ds; 1099 ideal I = jacob(x5+y6+z7+xyz); 1100 ideal J1 = modStd(I,1,1); 1101 J1; 1102 ideal J2 = modStd(I,3); 1103 J2; 1104 size(reduce(J1,J2)); 1105 size(reduce(J2,J1)); 1106 1107 ring rr = 0,x(1..4),lp; 1108 ideal I = cyclic(4); 1109 ideal J1 = modStd(I,2); 1110 ideal J2 = modStd(I,2,1); 1111 size(reduce(J1,J2)); 1112 size(reduce(J2,J1)); 1113 } 1114 1115 //////////////////////////////////////////////////////////////////////////////// 1116 1117 proc modS(ideal I, list L, list #) 240 1118 "USAGE: modS(I,L); I ideal, L intvec of primes 241 1119 if size(#)>0 std is used instead of groebner … … 247 1125 " 248 1126 { 249 int j,k; 250 list T,TT; 251 def R0=basering; 252 ideal J,cT,lT,K; 253 ideal I0=I; 254 list rl=ringlist(R0); 255 if((npars(R0)>0)(rl[1]>0)) 256 { 257 ERROR("characteristic of basering should be zero"); 258 } 259 for (j=1;j<=size(L);j++) 260 { 261 rl[1]=L[j]; 262 def @r=ring(rl); 263 setring @r; 264 ideal i=fetch(R0,I); 265 option(redSB); 266 if(size(#)>0) 267 { 268 i=std(i); 269 } 270 else 271 { 272 i=groebner(i); 273 } 274 setring R0; 275 T[j]=fetch(@r,i); 276 kill @r; 277 } 278 //================= delete unlucky primes ==================== 279 // unlucky if and only if the leading ideal is wrong 280 list LL=deleteUnluckyPrimes(T,L); 281 T=LL[1]; 282 L=LL[2]; 283 //============ now all leading ideals are the same ============ 284 for(j=1;j<=ncols(T[1]);j++) 285 { 286 for(k=1;k<=size(L);k++) 287 { 288 TT[k]=T[k][j]; 289 } 290 J[j]=liftPoly(TT,L); 291 } 292 attrib(J,"isSB",1); 293 return(J); 1127 int j; 1128 bigint N = 1; 1129 def R0 = basering; 1130 ideal J; 1131 list T; 1132 list rl = ringlist(R0); 1133 if((npars(R0)>0)  (rl[1]>0)) { ERROR("characteristic of basering should be zero"); } 1134 for(j = 1; j <= size(L); j++) 1135 { 1136 N = N*L[j]; 1137 rl[1] = L[j]; 1138 def @r = ring(rl); 1139 setring @r; 1140 ideal I = fetch(R0,I); 1141 if(size(#) > 0) 1142 { 1143 I = std(I); 1144 } 1145 else 1146 { 1147 I = groebner(I); 1148 } 1149 setring R0; 1150 T[j] = fetch(@r,I); 1151 kill @r; 1152 } 1153 L = deleteUnluckyPrimes(T,L,homog(I)); // unlucky if and only if the leading ideal is wrong 1154 J = farey(chinrem(L[1],L[2]),N); 1155 attrib(J,"isSB",1); 1156 return(J); 294 1157 } 295 1158 example 296 1159 { "EXAMPLE:"; echo = 2; 297 intvec L=3,5,11,13,181; 298 ring r=0,(x,y,z),dp; 299 ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4; 300 ideal J=modS(I,L); 1160 list L = 3,5,11,13,181,32003; 1161 ring r = 0,(x,y,z,t),dp; 1162 ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4; 1163 I = homog(I,t); 1164 ideal J = modS(I,L); 301 1165 J; 302 1166 } 303 /////////////////////////////////////////////////////////////////////////////// 304 proc deleteUnluckyPrimes(list T,intvec L) 305 "USAGE: deleteUnluckyPrimes(T,L);T list of polys, L intvec of primes 306 RETURN: list L,T with T list of polys, L intvec of primes 307 EXAMPLE: example deleteUnluckyPrimes; shows an example 308 NOTE: works only for homogeneous ideals with global orderings or 309 for ideals with local orderings 1167 1168 //////////////////////////////////////////////////////////////////////////////// 1169 1170 proc modHenselStd(ideal I, list #) 1171 "USAGE: modHenselStd(I); 1172 RETURN: a standard basis of I; 1173 NOTE: The procedure computes a standard basis of I (over the rational 1174 numbers) by using modular computations and Hensellifting. 1175 For further experiments see procedure modS. 1176 EXAMPLE: example modHenselStd; shows an example 310 1177 " 311 1178 { 312 int j,k; 313 intvec hl,hc; 314 ideal cT,lT; 315 316 lT=lead(T[size(T)]); 317 attrib(lT,"isSB",1); 318 hl=hilb(lT,1); 319 for (j=1;j<size(T);j++) 320 { 321 cT=lead(T[j]); 322 attrib(cT,"isSB",1); 323 hc=hilb(cT,1); 324 if(hl==hc) 325 { 326 for(k=1;k<=size(lT);k++) 327 { 328 if(lT[k]<cT[k]){lT=cT;break;} 329 if(lT[k]>cT[k]){break;} 330 } 331 } 332 else 333 { 334 if(hc<hl){lT=cT;hl=hilb(lT,1);} 335 } 336 } 337 j=1; 338 attrib(lT,"isSB",1); 339 while(j<=size(T)) 340 { 341 cT=lead(T[j]); 342 attrib(cT,"isSB",1); 343 if((size(reduce(cT,lT))!=0)(size(reduce(lT,cT))!=0)) 344 { 345 T=delete(T,j); 346 if(j==1) { L=L[2..size(L)]; } 347 else 348 { 349 if (j==size(L)) { L=L[1..size(L)1]; } 350 else { L=L[1..j1],L[j+1..size(L)]; } 351 } 352 j; 353 } 354 j++; 355 } 356 return(list(T,L,lT)); 357 } 358 example 359 { "EXAMPLE:"; echo = 2; 360 intvec L=2,3,5,7,11; 361 ring r=0,(y,x),Dp; 362 ideal I1=y2x,y6; 363 ideal I2=yx2,y3x,x5,y6; 364 ideal I3=y2x,x3y,x5,y6; 365 ideal I4=y2x,x3y,x5; 366 ideal I5=y2x,yx3,x5,y6; 367 list T=I1,I2,I3,I4,I5; 368 list TT=deleteUnluckyPrimes(T,L); 369 TT; 370 } 371 /////////////////////////////////////////////////////////////////////////////// 372 proc liftPoly(list T, intvec L) 373 "USAGE: liftPoly(T,L); T list of polys, L intvec of primes 374 RETURN: poly p in Q[x] such that p mod L[i]=T[i] 375 EXAMPLE: example liftPoly; shows an example 376 " 377 { 378 int i; 379 list TT; 380 for(i=size(T);i>0;i) 381 { TT[i]=ideal(T[i]); } 382 T=TT; 383 ideal hh=chinrem(T,L); 384 poly h=hh[1]; 385 poly p=lead(h); 386 poly result; 387 number n; 388 bigint N=L[1]; 389 for(i=size(L);i>1;i) 390 { 391 N=N*L[i]; 392 } 393 while(h!=0) 394 { 395 n=Farey(N,bigint(leadcoef(h))); 396 result=result+n*p; 397 h=hlead(h); 398 p=leadmonom(h); 399 } 400 return(result); 401 } 402 example 403 { "EXAMPLE:"; echo = 2; 404 ring R = 0,(x,y),dp; 405 intvec L=32003,181,241,499; 406 list T=ideal(x2+7000x+13000),ideal(x2+100x+147y+40),ideal(x2+120x+191y+10),ideal(x2+x+67y+100); 407 liftPoly(T,L); 408 } 409 /////////////////////////////////////////////////////////////////////////// 410 proc liftPoly1(list T, intvec L) 411 "USAGE: liftPoly1(T,L); T list of polys, L intvec of primes 412 RETURN: poly p in Q[x] such that p mod L[i]=T[i] 413 EXAMPLE: example liftPoly1; shows an example 414 " 415 { 416 poly result; 417 int i; 418 poly p; 419 list TT; 420 number n; 421 422 bigint N=L[1]; 423 for(i=2;i<=size(L);i++) 424 { 425 N=N*L[i]; 426 } 1179 int i,j; 1180 1181 bigint p = 2134567879; 1182 if(size(#)!=0) { p=#[1]; } 1183 while(!primeTest(I,p)) 1184 { 1185 p = prime(random(2000000000,2134567879)); 1186 } 1187 1188 def R = basering; 1189 module F,PrevG,PrevZ,Z2; 1190 ideal testG,testG1,G1,G2,G3,Gp; 1191 list L = p; 1192 list rl = ringlist(R); 1193 rl[1] = int(p); 1194 1195 def S = ring(rl); 1196 setring S; 1197 option(redSB); 1198 module Z,M,Z2; 1199 ideal I = imap(R,I); 1200 ideal Gp,G1,G2,G3; 1201 Gp,Z = liftstd1(I); 1202 attrib(Gp,"isSB",1); 1203 module ZZ = syz(I); 1204 attrib(ZZ,"isSB",1); 1205 Z = reduce(Z,ZZ); 1206 1207 setring R; 1208 Gp = imap(S,Gp); 1209 PrevZ = imap(S,Z); 1210 PrevG = module(Gp); 1211 F = module(I); 1212 testG = farey(Gp,p); 1213 attrib(testG,"isSB",1); 427 1214 while(1) 428 1215 { 429 p=leadmonom(T[1]); 430 for(i=2;i<=size(T);i++) 431 { 432 if(leadmonom(T[i])>p) 433 { 434 p=leadmonom(T[i]); 435 } 436 } 437 if (p==0) {return(result);} 438 for(i=1;i<=size(T);i++) 439 { 440 if(p==leadmonom(T[i])) 441 { 442 TT[i]=leadcoef(T[i]); 443 T[i]=T[i]lead(T[i]); 444 } 445 else 446 { 447 TT[i]=0; 448 } 449 } 450 n=chineseR(TT,L,N); 451 n=Farey(N,bigint(n)); 452 result=result+n*p; 453 } 454 } 455 example 456 { "EXAMPLE:"; echo = 2; 457 ring R = 0,(x,y),dp; 458 intvec L=32003,181,241,499; 459 list T=x2+7000x+13000,x2+100x+147y+40,x2+120x+191y+10,x2+x+67y+100; 460 liftPoly1(T,L); 461 } 462 /////////////////////////////////////////////////////////////////////////////// 463 proc fareyIdeal(ideal I,intvec L) 464 { 465 poly result,p; 466 int i,j; 467 number n; 468 bigint N=L[1]; 469 for(i=2;i<=size(L);i++) 470 { 471 N=N*L[i]; 472 } 473 474 for(i=1;i<=size(I);i++) 475 { 476 p=I[i]; 477 result=lead(p); 478 while(1) 479 { 480 if (p==0) {break;} 481 p=plead(p); 482 n=Farey(N,bigint(leadcoef(p))); 483 result=result+n*leadmonom(p); 484 } 485 I[i]=result; 486 } 487 return(I); 488 } 489 /////////////////////////////////////////////////////////////////////////////// 490 proc Farey (bigint P, bigint N) 491 "USAGE: Farey (P,N); P, N number; 492 RETURN: a rational number a/b such that a/b=N mod P 493 and a,b<(P/2)^{1/2} 494 " 495 { 496 if (P<0){P=P;} 497 if (N<0){N=N+P;} 498 bigint A,B,C,D,E; 499 E=P; 500 B=1; 501 while (N!=0) 502 { 503 if (2*N^2<P) 504 { 505 return(number(N)/number(B)); 506 } 507 D=E mod N; 508 C=A(E div N)*B; 509 E=N; 510 N=D; 511 A=B; 512 B=C; 513 } 514 return(0); 515 } 516 example 517 { "EXAMPLE:"; echo = 2; 518 ring R = 0,x,dp; 519 Farey(32003,12345); 520 } 521 /////////////////////////////////////////////////////////////////////////////// 522 proc chineseR(list T,intvec L,number N) 523 "USAGE: chineseR(T,L,N); 524 RETURN: x such that x = T[i] mod L[i], N=product(L[i]) 525 NOTE: chinese remainder theorem 526 EXAMPLE:example chineseR; shows an example 527 " 528 { 529 number x; 530 if(size(L)==1) 531 { 532 x=T[1] mod L[1]; 533 return(x); 534 } 535 int i; 536 int n=size(L); 537 list M; 538 for(i=1;i<=n;i++) 539 { 540 M[i]=N/L[i]; 541 } 542 list S=eexgcdN(M); 543 for(i=1;i<=n;i++) 544 { 545 x=x+S[i]*M[i]*T[i]; 546 } 547 x=x mod N; 548 return(x); 549 } 550 example 551 { "EXAMPLE:"; echo = 2; 552 ring R = 0,x,dp; 553 chineseR(list(24,15,7),intvec(2,3,5),30); 554 } 555 556 /////////////////////////////////////////////////////////////////////////////// 557 proc pStd(int p,ideal i) 558 "USAGE: pStd(p,i);p integer, i ideal; 559 RETURN: an ideal G which is the groebner base for i 560 EXAMPLE: example pStd; shows an example 561 " 562 { 563 def r=basering; 564 list rl=ringlist(r); 565 rl[1]=p; 566 def r1=ring(rl); 567 setring r1; 568 option(redSB); 569 ideal j=fetch(r,i); 570 ideal GP=groebner(j); 571 setring r; 572 ideal G=fetch(r1,GP); 573 attrib(G,"isSB",1); 574 matrix Z=transmat(p,i,G); 575 matrix G1=gstrich1(p,Z,i,G); 576 ideal g1=G1; 577 ideal g22=reduce(g1,G); 578 matrix G22=transpose(matrix(g22)); 579 matrix M=redmat(G,G1,G22); 580 matrix Z2=M*Z; 581 kill r1; 582 number c=p; 583 bigint cb=p; 584 matrix G0=transpose(matrix(G)); 585 G0= MmodN(G0+ (c)* G22,cb^2); 586 matrix GF=fareyMatrix(G0,cb^2); 587 Z=MmodN(Z+(c)*Z2,cb^2); 588 matrix C=transpose(G); 589 int n=3; 590 while(GF<>C) 591 { 592 C=GF; 593 G1= gstrich2(c,Z,i,G0,n); 594 g1=G1; 595 g22=reduce(g1,G); 596 G22=transpose(matrix(g22)); 597 M=redmat(G,G1,G22); 598 Z2=M*Z; 599 Z=MmodN(Z+(c^(n1))*Z2,cb^n); 600 G0= MmodN(G0+ (c^(n1))* G22,cb^n); 601 GF=fareyMatrix(G0,cb^n); 602 n++; 603 } 604 return(ideal(GF)); 605 } 606 example 607 { "EXAMPLE:"; echo = 2; 608 ring r=0,(x,y,z),dp; 609 ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4; 610 ideal J=pStd(32003,I); 611 J; 612 } 613 /////////////////////////////////////////////////////////////////////////// 614 proc transmat(int p,ideal i,ideal G) 615 "USAGE: transmat(p,I,G); p integer, I,G ideal; 616 RETURN: the transformationmatrix Z for the ideal i mod p and the groebner base for i mod p 617 EXAMPLE: example transmat; shows an example 618 " 619 { 620 def r=basering; 621 int n=nvars(r); 622 list rl=ringlist(r); 623 rl[1]=p; 624 def r1=ring(rl); 625 setring r1; 626 ideal i=fetch(r,i); 627 ideal G=fetch(r,G); 628 attrib(G,"isSB",1); 629 ring rhelp=p,x(1..n),dp; 630 list lhelp=ringlist(rhelp); 631 list l=lhelp[3]; 632 setring r; 633 rl[3]=l; 634 def r2=ring(rl); 635 setring r2; 636 ideal i=fetch(r,i); 637 option(redSB); 638 ideal j=std(i); 639 matrix T=lift(i,j); 640 setring r1; 641 matrix T=fetch(r2,T); 642 ideal j=fetch(r2,j); 643 matrix M=lift(j,G); 644 matrix Z=transpose(T*M); 645 setring r; 646 matrix Z=fetch(r1,Z); 647 return(Z); 648 } 649 example 650 { "EXAMPLE:"; echo = 2; 651 ring r=0,(x,y,z),dp; 652 ideal i=3x3+x2+1,11y5+y3+2,5z4+z2+4; 653 ideal g=x360x260, z436z2+37, y5+33y3+66; 654 int p=181; 655 matrix Z=transmat(p,i,g); 656 Z; 657 } 658 659 /////////////////////////////////////////////////////////////////////////// 660 proc gstrich1(int p, matrix Z, ideal i, ideal gp) 661 "USAGE: gstrich1 (p,Z,i,gp); p integer, Z matrix, i,gp ideals; 662 RETURN: a matrix G such that (Z*FGP)/p, where F and GP are the matrices of the ideals i and gp 663 " 664 { 665 matrix F=transpose(matrix(i)); 666 matrix GP=transpose(matrix(gp)); 667 matrix G=(Z*FGP)/p; 668 return(G); 669 } 670 /////////////////////////////////////////////////////////////////////////// 671 proc gstrich2(number p, matrix Z, ideal i, ideal gp, int n) 672 "USAGE: gstrich2 (p,Z,i,gp,n); p,n integer, Z matrix, i,gp ideals; 673 RETURN: a matrix G such that (Z*FGP)/(p^(n1)), where F and GP are the matrices of the ideals i and gp 674 " 675 { 676 matrix F=transpose(matrix(i)); 677 matrix GP=transpose(matrix(gp)); 678 matrix G=(Z*FGP)/(p^(n1)); 679 return(G); 680 } 681 /////////////////////////////////////////////////////////////////////////// 682 proc redmat(ideal i, matrix h, matrix g) 683 "USAGE: redmat(i,h,g); i ideal , h,g matrices; 684 RETURN: a matrix M such that i=M*h+g 685 " 686 { 687 matrix c=hg; 688 ideal f=transpose(c); 689 matrix N=lift(i,f); 690 matrix M=transpose(N); 691 return(M); 692 } 693 /////////////////////////////////////////////////////////////////////////// 694 proc fareyMatrix(matrix m,bigint N) 695 "USAGE: fareyMatrix(m,y); m matrix, y integer; 696 RETURN: a matrix k of the matrix m with Farey rational numbers a/b as coefficients 697 EXAMPLE: example fareyMatrix; shows an example 698 " 699 { 700 ideal I=m; 701 poly result,p; 702 int i,j; 703 number n; 704 for(i=1;i<=size(I);i++) 705 { 706 p=I[i]; 707 result=lead(p); 708 while(1) 709 { 710 if (p==0) {break;} 711 p=plead(p); 712 n=Farey(N,bigint(leadcoef(p))); 713 result=result+n*leadmonom(p); 714 } 715 I[i]=result; 716 } 717 matrix k=transpose(I); 718 return(k); 719 } 720 example 721 {"EXAMPLE:"; echo = 2; 722 ring r=0,(x,y,z),dp; 723 matrix m[3][1]=x3+682794673x2+682794673,z4+204838402z2+819353608, y5+186216729y3+372433458; 724 int p=32003; 725 matrix b=fareyMatrix(m,p^2); 726 b; 727 } 728 /////////////////////////////////////////////////////////////////////////// 729 proc MmodN(matrix Z,bigint N) 730 "USAGE: MmodN(Z,N);Z matrix, N bigint; 731 RETURN: the matrix Z mod N 732 EXAMPLE: example MmodN; 733 " 734 { 735 int i,j,k; 736 poly m,p; 737 number c; 738 for(i=1;i<=nrows(Z);i++) 739 { 740 for(j=1;j<=ncols(Z);j++) 741 { 742 for(k=1;k<=size(Z[i,j]);k++) 743 { 744 m=leadmonom(Z[i,j][k]); 745 c=bigint(leadcoef(Z[i,j][k])) mod N; 746 p=p+c*m; 747 } 748 Z[i,j]=p; 749 p=0; 750 } 751 } 752 return(Z); 1216 i++; 1217 G1 = ideal(1/(p^i) * sum(F*PrevZ,(1)*PrevG)); 1218 setring S; 1219 G1 = imap(R,G1); 1220 G2 = reduce(G1,Gp); 1221 G3 = sum(G1,(1)*G2); 1222 M = lift(Gp,G3); 1223 Z2 = (1)*Z*M; 1224 1225 setring R; 1226 G2 = imap(S,G2); 1227 Z2 = imap(S,Z2); 1228 PrevG = sum(PrevG, module(p^i*G2)); 1229 PrevZ = sum(PrevZ, multiply(poly(p^i),Z2)); 1230 testG1 = farey(ideal(PrevG),p^(i+1)); 1231 attrib(testG1,"isSB",1); 1232 if(size(reduce(testG1,testG)) == 0) 1233 { 1234 if(size(reduce(I,testG1)) == 0) // I is in testG1 1235 { 1236 if(pTestSB(I,testG1,L,2)) 1237 { 1238 G3 = std(testG1); // testG1 is SB 1239 if(size(reduce(G3,testG1)) == 0) 1240 { 1241 return(G3); 1242 } 1243 } 1244 } 1245 } 1246 testG = testG1; 1247 attrib(testG,"isSB",1); 1248 } 753 1249 } 754 1250 example 755 1251 { "EXAMPLE:"; echo = 2; 756 1252 ring r = 0,(x,y,z),dp; 757 matrix m[3][1]= x3+10668x2+10668, z412801z2+12802, y58728y3+14547; 758 bigint p=32003; 759 matrix b=MmodN(m,p^2); 760 b; 761 } 762 /////////////////////////////////////////////////////////////////////////////// 1253 ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4; 1254 ideal J = modHenselStd(I); 1255 J; 1256 } 1257 1258 //////////////////////////////////////////////////////////////////////////////// 1259 1260 static proc sum(list #) 1261 { 1262 if(typeof(#[1])=="ideal") 1263 { 1264 ideal M; 1265 } 1266 else 1267 { 1268 module M; 1269 } 1270 1271 int i; 1272 for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; } 1273 return(M); 1274 } 1275 1276 //////////////////////////////////////////////////////////////////////////////// 1277 1278 static proc multiply(poly p, list #) 1279 { 1280 if(typeof(#[1])=="ideal") 1281 { 1282 ideal M; 1283 } 1284 else 1285 { 1286 module M; 1287 } 1288 1289 int i; 1290 for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; } 1291 return(M); 1292 } 1293 1294 1295 ////////////////////////////// further examples //////////////////////////////// 1296 763 1297 /* 764 ring r =0,(x,y,z),lp;1298 ring r = 0, (x,y,z), lp; 765 1299 poly s1 = 5x3y2z+3y3x2z+7xy2z2; 766 1300 poly s2 = 3xy2z2+x5+11y2z2; … … 769 1303 ideal i = s1, s2, s3, s4; 770 1304 771 ring r =0,(x,y,z),lp;1305 ring r = 0, (x,y,z), lp; 772 1306 poly s1 = 2xy4z2+x3y2zx2y3z+2xyz2+7y3+7; 773 1307 poly s2 = 2x2y4z+x2yz2xy2z2+2x2yz12x+12y; … … 776 1310 ideal i = s1, s2, s3, s4; 777 1311 778 ring r =0,(x,y,z),lp;1312 ring r = 0, (x,y,z), lp; 779 1313 poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz; 780 1314 poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4; 781 1315 poly s3 = 8x3 + 12y3 + xz2 + 3; 782 1316 poly s4 = 7x2y4 + 18xy3z2 + y3z3; 783 ideal i = 1317 ideal i = s1, s2, s3, s4; 784 1318 785 1319 int n = 6; 786 1320 ring r = 0,(x(1..n)),lp; 787 1321 ideal i = cyclic(n); 788 ring s =0,(x(1..n),t),lp;789 ideal i =imap(r,i);790 i =homog(i,t);791 792 ring r =0,(x(1..4),s),(dp(4),dp);793 poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);1322 ring s = 0, (x(1..n),t), lp; 1323 ideal i = imap(r,i); 1324 i = homog(i,t); 1325 1326 ring r = 0, (x(1..4),s), (dp(4),dp); 1327 poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4); 794 1328 poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4); 795 1329 poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4); … … 798 1332 ideal i = s1, s2, s3, s4, s5; 799 1333 800 ring r=0,(x,y,z),ds; 801 int a =16; 802 int b =15; 803 int c =4; 804 int t =1; 805 poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c1)+x^(c1)*y^(c1)*z3+x^(c2)*y^c*(y2+t*x)^2; 806 ideal i= jacob(f); 807 808 ring r=0,(x,y,z),ds; 809 int a =25; 810 int b =25; 811 int c =5; 812 int t =1; 813 poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c1)+x^(c1)*y^(c1)*z3+x^(c2)*y^c*(y2+t*x)^2; 814 ideal i= jacob(f),f; 815 816 ring r=0,(x,y,z),ds; 817 int a=10; 818 poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a; 819 ideal i= jacob(f); 820 821 ring r=0,(x,y,z),ds; 822 int a =6; 823 int b =8; 824 int c =10; 825 int alpha =5; 826 int beta= 5; 827 int t= 1; 828 poly f =x^a+y^b+z^c+x^alpha*y^(beta5)+x^(alpha2)*y^(beta3)+x^(alpha3)*y^(beta4)*z^2+x^(alpha4)*y^(beta4)*(y^2+t*x)^2; 829 ideal i= jacob(f); 830 1334 ring r = 0, (x,y,z), ds; 1335 int a = 16; 1336 int b = 15; 1337 int c = 4; 1338 int t = 1; 1339 poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c1)+x^(c1)*y^(c1)*z3+x^(c2)*y^c*(y2+t*x)^2; 1340 ideal i = jacob(f); 1341 1342 ring r = 0, (x,y,z), ds; 1343 int a = 25; 1344 int b = 25; 1345 int c = 5; 1346 int t = 1; 1347 poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c1)+x^(c1)*y^(c1)*z3+x^(c2)*y^c*(y2+t*x)^2; 1348 ideal i = jacob(f),f; 1349 1350 ring r = 0, (x,y,z), ds; 1351 int a = 10; 1352 poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a; 1353 ideal i = jacob(f); 1354 1355 ring r = 0, (x,y,z), ds; 1356 int a = 6; 1357 int b = 8; 1358 int c = 10; 1359 int alpha = 5; 1360 int beta = 5; 1361 int t = 1; 1362 poly f = x^a+y^b+z^c+x^alpha*y^(beta5)+x^(alpha2)*y^(beta3)+x^(alpha3)*y^(beta4)*z^2+x^(alpha4)*y^(beta4)*(y^2+t*x)^2; 1363 ideal i = jacob(f); 831 1364 */ 832 1365 833 /*834 ring r=0,(x,y,z),lp;835 poly s1 = 5x3y2z+3y3x2z+7xy2z2;836 poly s2 = 3xy2z2+x5+11y2z2;837 poly s3 = 4xyz+7x3+12y3+1;838 poly s4 = 3x34y3+yz2;839 ideal i = s1, s2, s3, s4;840 841 ring r=0,(x,y,z),lp;842 poly s1 = 2xy4z2+x3y2zx2y3z+2xyz2+7y3+7;843 poly s2 = 2x2y4z+x2yz2xy2z2+2x2yz12x+12y;844 poly s3 = 2y5z+x2y2zxy3zxy3+y4+2y2z;845 poly s4 = 3xy4z3+x2y2zxy3z+4y3z2+3xyz3+4z2x+y;846 ideal i = s1, s2, s3, s4;847 848 ring r=0,(x,y,z),lp;849 poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;850 poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;851 poly s3 = 8x3 + 12y3 + xz2 + 3;852 poly s4 = 7x2y4 + 18xy3z2 + y3z3;853 ideal i = s1, s2, s3, s4;854 855 int n = 6;856 ring r = 0,(x(1..n)),lp;857 ideal i = cyclic(n);858 ring s=0,(x(1..n),t),lp;859 ideal i=imap(r,i);860 i=homog(i,t);861 862 ring r=0,(x(1..4),s),(dp(4),dp);863 poly s1 =1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);864 poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);865 poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);866 poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4);867 poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);868 ideal i = s1, s2, s3, s4, s5;869 870 ring r=0,(x,y,z),ds;871 int a =16;872 int b =15;873 int c =4;874 int t =1;875 poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c1)+x^(c1)*y^(c1)*z3+x^(c2)*y^c*(y2+t*x)^2;876 ideal i= jacob(f);877 878 ring r=0,(x,y,z),ds;879 int a =25;880 int b =25;881 int c =5;882 int t =1;883 poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c1)+x^(c1)*y^(c1)*z3+x^(c2)*y^c*(y2+t*x)^2;884 ideal i= jacob(f),f;885 886 ring r=0,(x,y,z),ds;887 int a=10;888 poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;889 ideal i= jacob(f);890 891 ring r=0,(x,y,z),ds;892 int a =6;893 int b =8;894 int c =10;895 int alpha =5;896 int beta= 5;897 int t= 1;898 poly f =x^a+y^b+z^c+x^alpha*y^(beta5)+x^(alpha2)*y^(beta3)+x^(alpha3)*y^(beta4)*z^2+x^(alpha4)*y^(beta4)*(y^2+t*x)^2;899 ideal i= jacob(f);900 901 */902
Note: See TracChangeset
for help on using the changeset viewer.