Changeset a2e51b in git
- Timestamp:
- May 13, 2014, 4:08:04 PM (10 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- c6b5f62a2d93b0c3cf0e5bd7cb3a8c2eee553e1b
- Parents:
- a93f95f272cf2777251e0e930411dedb45d02ba0
- git-author:
- Andreas Steenpass <steenpass@mathematik.uni-kl.de>2014-05-13 16:08:04+02:00
- git-committer:
- Andreas Steenpass <steenpass@mathematik.uni-kl.de>2014-05-20 08:03:42+02:00
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/modstd.lib
ra93f95f ra2e51b 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="version modstd.lib 4.0.0.0 Dec_2013"; // $Id$3 category ="Commutative Algebra";2 version="version modstd.lib 4.0.0.0 May_2014 "; // $Id$ 3 category="Commutative Algebra"; 4 4 info=" 5 LIBRARY: modstd.lib Groebner bas is of ideals5 LIBRARY: modstd.lib Groebner bases of ideals using modular methods 6 6 7 7 AUTHORS: A. Hashemi Amir.Hashemi@lip6.fr 8 @*G. Pfister pfister@mathematik.uni-kl.de9 @*H. Schoenemann hannes@mathematik.uni-kl.de10 @*A. Steenpass steenpass@mathematik.uni-kl.de11 @*S. Steidel steidel@mathematik.uni-kl.de8 G. Pfister pfister@mathematik.uni-kl.de 9 H. Schoenemann hannes@mathematik.uni-kl.de 10 A. Steenpass steenpass@mathematik.uni-kl.de 11 S. Steidel steidel@mathematik.uni-kl.de 12 12 13 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, 403-419 (2003). 14 A library for computing Groebner bases of ideals in the polynomial ring over 15 the rational numbers using modular methods. 16 17 REFERENCES: 18 E. A. Arnold: Modular algorithms for computing Groebner bases. 19 J. Symb. Comp. 35, 403-419 (2003). 20 21 N. Idrees, G. Pfister, S. Steidel: Parallelization of Modular Algorithms. 22 J. Symb. Comp. 46, 672-684 (2011). 20 23 21 24 PROCEDURES: 22 modStd(I); standard basis of I using modular methods (chinese remainder) 23 modS(I,L); liftings to Q of standard bases of I mod p for p in L 24 modHenselStd(I); standard basis of I using modular methods (hensel lifting) 25 modStd(I); standard basis of I using modular methods 25 26 "; 26 27 27 28 LIB "poly.lib"; 28 LIB "ring.lib"; 29 LIB "parallel.lib"; 30 31 //////////////////////////////////////////////////////////////////////////////// 32 33 static proc mod_init() 34 { 35 newstruct("idealPrimeTest", "ideal Ideal"); 36 } 37 38 //////////////////////////////////////////////////////////////////////////////// 39 40 static proc redFork(ideal I, ideal J, int n) 41 { 42 attrib(J,"isSB",1); 43 return(reduce(I,J,1)); 44 } 45 46 //////////////////////////////////////////////////////////////////////////////// 47 48 proc isIncluded(ideal I, ideal J, list #) 49 "USAGE: isIncluded(I,J); I,J ideals 50 RETURN: 1 if J includes I, 51 @* 0 if there is an element f in I which does not reduce to 0 w.r.t. J. 52 EXAMPLE: example isIncluded; shows an example 53 " 54 { 55 def R = basering; 56 setring R; 57 58 attrib(J,"isSB",1); 59 int i,j,k; 60 61 if(size(#) > 0) 62 { 63 int n = #[1]; 64 if(n >= ncols(I)) { n = ncols(I); } 65 if(n > 1) 66 { 67 for(i = 1; i <= n - 1; i++) 68 { 69 //link l(i) = "MPtcp:fork"; 70 link l(i) = "ssi:fork"; 71 open(l(i)); 72 73 write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1))); 74 } 75 76 int t = timer; 77 if(reduce(I[ncols(I)], J, 1) != 0) 78 { 79 for(i = 1; i <= n - 1; i++) 80 { 81 close(l(i)); 29 LIB "modular.lib"; 30 31 proc modStd(ideal I, list #) 32 "USAGE: modStd(I[, exactness]); I ideal, exactness int 33 RETURN: a standard basis of I 34 NOTE: The procedure computes a standard basis of I (over the rational 35 numbers) by using modular methods. 36 @* An optional parameter 'exactness' can be provided. 37 If exactness = 1, the procedure computes a standard basis of I for 38 sure; if exactness = 0, it computes a standard basis of I 39 with high probability. 40 SEE ALSO: modular 41 EXAMPLE: example modStd; shows an example" 42 { 43 /* read optional parameter */ 44 int exactness = 1; 45 if (size(#) > 0) { 46 /* For compatibility, we only test size(#) > 4. This can be changed to 47 * size(#) > 1 in the future. */ 48 if (size(#) > 4 || typeof(#[1]) != "int") { 49 ERROR("wrong optional parameter"); 50 } 51 exactness = #[1]; 52 } 53 54 /* save options */ 55 intvec opt = option(get); 56 option(redSB); 57 58 /* choose the right command */ 59 string command = "groebner"; 60 if (npars(basering) > 0) { 61 command = "Modstd::groebner_norm"; 62 } 63 64 /* call modular() */ 65 if (exactness) { 66 I = modular(command, list(I), primeTest_std, 67 deleteUnluckyPrimes_std, pTest_std, finalTest_std); 68 } 69 else { 70 I = modular(command, list(I), primeTest_std, 71 deleteUnluckyPrimes_std, pTest_std); 72 } 73 74 /* return the result */ 75 attrib(I, "isSB", 1); 76 option(set, opt); 77 return(I); 78 } 79 example 80 { 81 "EXAMPLE:"; 82 echo = 2; 83 ring R1 = 0, (x,y,z,t), dp; 84 ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; 85 ideal J = modStd(I); 86 J; 87 I = homog(I, t); 88 J = modStd(I); 89 J; 90 91 ring R2 = 0, (x,y,z), ds; 92 ideal I = jacob(x5+y6+z7+xyz); 93 ideal J = modStd(I, 0); 94 J; 95 96 ring R3 = 0, x(1..4), lp; 97 ideal I = cyclic(4); 98 ideal J1 = modStd(I, 1); // default 99 ideal J2 = modStd(I, 0); 100 size(reduce(J1, J2)); 101 size(reduce(J2, J1)); 102 } 103 104 /* compute a normalized GB via groebner() */ 105 static proc groebner_norm(ideal I) 106 { 107 I = simplify(groebner(I), 1); 108 attrib(I, "isSB", 1); 109 return(I); 110 } 111 112 /* test if the prime p is suitable for the input, i.e. it does not divide 113 * the numerator or denominator of any of the coefficients */ 114 static proc primeTest_std(int p, alias list args) 115 { 116 /* erase zero generators */ 117 ideal I = simplify(args[1], 2); 118 119 /* clear denominators and count the terms */ 120 ideal J; 121 ideal K; 122 int n = ncols(I); 123 intvec sizes; 124 number cnt; 125 int i; 126 for(i = n; i > 0; i--) { 127 J[i] = cleardenom(I[i]); 128 cnt = leadcoef(J[i])/leadcoef(I[i]); 129 K[i] = numerator(cnt)*var(1)+denominator(cnt); 130 } 131 sizes = size(J[1..n]); 132 133 /* change to characteristic p */ 134 def br = basering; 135 list lbr = ringlist(br); 136 if (typeof(lbr[1]) == "int") { 137 lbr[1] = p; 138 } 139 else { 140 lbr[1][1] = p; 141 } 142 def rp = ring(lbr); 143 setring(rp); 144 ideal Jp = fetch(br, J); 145 ideal Kp = fetch(br, K); 146 147 /* test if any coefficient is missing */ 148 if (intvec(size(Kp[1..n])) != 2:n) { 149 setring(br); 150 return(0); 151 } 152 if (intvec(size(Jp[1..n])) != sizes) { 153 setring(br); 154 return(0); 155 } 156 setring(br); 157 return(1); 158 } 159 160 /* find entries in modresults which come from unlucky primes. 161 * For this, sort the entries into categories depending on their leading 162 * ideal and return the indices in all but the biggest category. */ 163 static proc deleteUnluckyPrimes_std(alias list modresults) 164 { 165 int size_modresults = size(modresults); 166 167 /* sort results into categories. 168 * each category is represented by three entries: 169 * - the corresponding leading ideal 170 * - the number of elements 171 * - the indices of the elements 172 */ 173 list cat; 174 int size_cat; 175 ideal L; 176 int i; 177 int j; 178 for (i = 1; i <= size_modresults; i++) { 179 L = lead(modresults[i]); 180 attrib(L, "isSB", 1); 181 for (j = 1; j <= size_cat; j++) { 182 if (size(L) == size(cat[j][1]) 183 && size(reduce(L, cat[j][1])) == 0 184 && size(reduce(cat[j][1], L)) == 0) { 185 cat[j][2] = cat[j][2]+1; 186 cat[j][3][cat[j][2]] = i; 187 break; 82 188 } 189 } 190 if (j > size_cat) { 191 size_cat++; 192 cat[size_cat] = list(); 193 cat[size_cat][1] = L; 194 cat[size_cat][2] = 1; 195 cat[size_cat][3] = list(i); 196 } 197 } 198 199 /* find the biggest categories */ 200 int cat_max = 1; 201 int max = cat[1][2]; 202 for (i = 2; i <= size_cat; i++) { 203 if (cat[i][2] > max) { 204 cat_max = i; 205 max = cat[i][2]; 206 } 207 } 208 209 /* return all other indices */ 210 list unluckyIndices; 211 for (i = 1; i <= size_cat; i++) { 212 if (i != cat_max) { 213 unluckyIndices = unluckyIndices + cat[i][3]; 214 } 215 } 216 return(unluckyIndices); 217 } 218 219 /* test if 'command' applied to 'args' in characteristic p is the same as 220 'result' mapped to characteristic p */ 221 static proc pTest_std(string command, list args, ideal result, int p) 222 { 223 /* change to characteristic p */ 224 def br = basering; 225 list lbr = ringlist(br); 226 if (typeof(lbr[1]) == "int") { 227 lbr[1] = p; 228 } 229 else { 230 lbr[1][1] = p; 231 } 232 def rp = ring(lbr); 233 setring(rp); 234 ideal Ip = fetch(br, args)[1]; 235 ideal Gp = fetch(br, result); 236 attrib(Gp, "isSB", 1); 237 238 /* test if Ip is in Gp */ 239 int i; 240 for (i = ncols(Ip); i > 0; i--) { 241 if (reduce(Ip[i], Gp, 1) != 0) { 242 setring(br); 83 243 return(0); 84 } 85 t = timer - t; 86 if(t > 60) { t = 60; } 87 int i_sleep = system("sh", "sleep "+string(t)); 88 89 j = ncols(I) - n; 90 91 while(j >= 0) 92 { 93 for(i = 1; i <= n - 1; i++) 94 { 95 if(status(l(i), "read", "ready")) 96 { 97 if(read(l(i)) != 0) 98 { 99 for(i = 1; i <= n - 1; i++) 100 { 101 close(l(i)); 102 } 103 return(0); 104 } 105 else 106 { 107 if(j >= 1) 108 { 109 write(l(i), quote(redFork(eval(I[j]), eval(J), 1))); 110 j--; 111 } 112 else 113 { 114 k++; 115 close(l(i)); 116 } 117 } 118 } 119 } 120 if(k == n - 1) 121 { 122 j--; 123 } 124 i_sleep = system("sh", "sleep "+string(t)); 125 } 126 return(1); 127 } 128 } 129 130 for(i = ncols(I); i >= 1; i--) 131 { 132 if(reduce(I[i],J,1) != 0){ return(0); } 133 } 134 return(1); 135 } 136 example 137 { "EXAMPLE:"; echo = 2; 138 ring r=0,(x,y,z),dp; 139 ideal I = x+1,x+y+1; 140 ideal J = x+1,y; 141 isIncluded(I,J); 142 isIncluded(J,I); 143 isIncluded(I,J,4); 144 145 ring R = 0, x(1..5), dp; 146 ideal I1 = cyclic(4); 147 ideal I2 = I1,x(5)^2; 148 isIncluded(I1,I2,4); 149 } 150 151 //////////////////////////////////////////////////////////////////////////////// 152 153 proc pTestSB(ideal I, ideal J, list L, int variant, list #) 154 "USAGE: pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int 155 RETURN: 1 (resp. 0) if for a randomly chosen prime p that is not in L 156 J mod p is (resp. is not) a standard basis of I mod p 157 EXAMPLE: example pTestSB; shows an example 158 " 159 { 160 int i,j,k,p; 161 def R = basering; 162 list r = ringlist(R); 163 164 while(!j) 165 { 166 j = 1; 167 p = prime(random(1000000000,2134567879)); 168 for(i = 1; i <= size(L); i++) 169 { 170 if(p == L[i]) { j = 0; break; } 171 } 172 if(j) 173 { 174 for(i = 1; i <= ncols(I); i++) 175 { 176 for(k = 2; k <= size(I[i]); k++) 177 { 178 if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; } 179 } 180 if(!j){ break; } 181 } 182 } 183 if(j) 184 { 185 if(!primeTest(I,p)) { j = 0; } 186 } 187 } 188 r[1] = p; 189 def @R = ring(r); 190 setring @R; 191 ideal I = imap(R,I); 192 ideal J = imap(R,J); 193 attrib(J,"isSB",1); 194 195 int t = timer; 196 j = 1; 197 if(isIncluded(I,J) == 0) { j = 0; } 198 199 if(printlevel >= 11) 200 { 201 "isIncluded(I,J) takes "+string(timer - t)+" seconds"; 202 "j = "+string(j); 203 } 204 205 t = timer; 206 if(j) 207 { 208 if(size(#) > 0) 209 { 210 ideal K = modpStd(I,p,variant,#[1])[1]; 211 } 212 else 213 { 214 ideal K = groebner(I); 215 } 216 t = timer; 217 if(isIncluded(J,K) == 0) { j = 0; } 218 219 if(printlevel >= 11) 220 { 221 "isIncluded(J,K) takes "+string(timer - t)+" seconds"; 222 "j = "+string(j); 223 } 224 } 225 setring R; 226 return(j); 227 } 228 example 229 { "EXAMPLE:"; echo = 2; 230 intvec L = 2,3,5; 231 ring r = 0,(x,y,z),dp; 232 ideal I = x+1,x+y+1; 233 ideal J = x+1,y; 234 pTestSB(I,I,L,2); 235 pTestSB(I,J,L,2); 236 } 237 238 //////////////////////////////////////////////////////////////////////////////// 239 240 proc deleteUnluckyPrimes(list T, list L, int ho, list #) 241 "USAGE: deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer 242 RETURN: lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #), 243 lT ideal 244 NOTE: - if ho = 1, the polynomials in T are homogeneous, else ho = 0, 245 @* - lT is prevalent, i.e. the most appearing leading ideal in T 246 EXAMPLE: example deleteUnluckyPrimes; shows an example 247 " 248 { 249 ho = ((ho)||(ord_test(basering) == -1)); 250 int j,k,c; 251 intvec hl,hc; 252 ideal cT,lT,cK; 253 lT = lead(T[size(T)]); 254 attrib(lT,"isSB",1); 255 if(!ho) 256 { 257 for(j = 1; j < size(T); j++) 258 { 259 cT = lead(T[j]); 260 attrib(cT,"isSB",1); 261 if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0)) 262 { 263 cK = cT; 264 c++; 265 } 266 } 267 if(c > size(T) div 2){ lT = cK; } 268 } 269 else 270 { 271 hl = hilb(lT,1); 272 for(j = 1; j < size(T); j++) 273 { 274 cT = lead(T[j]); 275 attrib(cT,"isSB",1); 276 hc = hilb(cT,1); 277 if(hl == hc) 278 { 279 for(k = 1; k <= size(lT); k++) 280 { 281 if(lT[k] < cT[k]) { lT = cT; c++; break; } 282 if(lT[k] > cT[k]) { c++; break; } 283 } 284 } 285 else 286 { 287 if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; } 288 } 289 } 290 } 291 292 int addList; 293 if(size(#) > 0) { list M = #; addList = 1; } 294 j = 1; 295 attrib(lT,"isSB",1); 296 while((j <= size(T))&&(c > 0)) 297 { 298 cT = lead(T[j]); 299 attrib(cT,"isSB",1); 300 if((size(reduce(cT,lT)) != 0)||(size(reduce(lT,cT)) != 0)) 301 { 302 T = delete(T,j); 303 if(j == 1) 304 { 305 L = L[2..size(L)]; 306 if(addList == 1) { M = M[2..size(M)]; } 307 } 308 else 309 { 310 if(j == size(L)) 311 { 312 L = L[1..size(L)-1]; 313 if(addList == 1) { M = M[1..size(M)-1]; } 314 } 315 else 316 { 317 L = L[1..j-1],L[j+1..size(L)]; 318 if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; } 319 } 320 } 321 j--; 322 } 323 j++; 324 } 325 326 for(j = 1; j <= size(L); j++) 327 { 328 L[j] = bigint(L[j]); 329 } 330 331 if(addList == 0) { return(list(T,L,lT)); } 332 if(addList == 1) { return(list(T,L,M,lT)); } 333 } 334 example 335 { "EXAMPLE:"; echo = 2; 336 list L = 2,3,5,7,11; 337 ring r = 0,(y,x),Dp; 338 ideal I1 = 2y2x,y6; 339 ideal I2 = yx2,y3x,x5,y6; 340 ideal I3 = y2x,x3y,x5,y6; 341 ideal I4 = y2x,11x3y,x5; 342 ideal I5 = y2x,yx3,x5,7y6; 343 list T = I1,I2,I3,I4,I5; 344 deleteUnluckyPrimes(T,L,1); 345 list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5); 346 deleteUnluckyPrimes(T,L,1,P); 347 } 348 349 //////////////////////////////////////////////////////////////////////////////// 350 351 proc primeTest(def II, bigint p) 352 { 353 if(typeof(II) == "string") 354 { 355 ideal I = `II`.Ideal; 356 } 357 else 358 { 359 ideal I = II; 360 } 361 362 I = simplify(I, 2); // erase zero generators 363 364 int i,j; 365 poly f; 366 number cnt; 367 for(i = 1; i <= size(I); i++) 368 { 369 f = cleardenom(I[i]); 370 if(f == 0) { return(0); } 371 cnt = leadcoef(I[i])/leadcoef(f); 372 if((numerator(cnt) mod p) == 0) { return(0); } 373 if((denominator(cnt) mod p) == 0) { return(0); } 374 for(j = size(f); j > 0; j--) 375 { 376 if((leadcoef(f[j]) mod p) == 0) { return(0); } 377 } 378 } 379 return(1); 380 } 381 382 //////////////////////////////////////////////////////////////////////////////// 383 384 proc primeList(ideal I, int n, list #) 385 "USAGE: primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal, 386 n integer 387 RETURN: the intvec of n greatest primes <= 2147483647 (resp. n greatest primes 388 < L[size(L)] union with L) such that none of these primes divides any 389 coefficient occuring in I 390 NOTE: The number of cores to use can be defined by ncores, default is 1. 391 EXAMPLE: example primeList; shows an example 392 " 393 { 394 intvec L; 395 int i,p; 396 int ncores = 1; 397 398 //----------------- Initialize optional parameter ncores --------------------- 399 if(size(#) > 0) 400 { 401 if(size(#) == 1) 402 { 403 if(typeof(#[1]) == "int") 404 { 405 ncores = #[1]; 406 # = list(); 407 } 408 } 409 else 410 { 411 ncores = #[2]; 412 } 413 } 414 415 if(size(#) == 0) 416 { 417 p = 2147483647; 418 while(!primeTest(I,p)) 419 { 420 p = prime(p-1); 421 if(p == 2) { ERROR("no more primes"); } 422 } 423 L[1] = p; 424 } 425 else 426 { 427 L = #[1]; 428 p = prime(L[size(L)]-1); 429 while(!primeTest(I,p)) 430 { 431 p = prime(p-1); 432 if(p == 2) { ERROR("no more primes"); } 433 } 434 L[size(L)+1] = p; 435 } 436 if(p == 2) { ERROR("no more primes"); } 437 if(ncores == 1) 438 { 439 for(i = 2; i <= n; i++) 440 { 441 p = prime(p-1); 442 while(!primeTest(I,p)) 443 { 444 p = prime(p-1); 445 if(p == 2) { ERROR("no more primes"); } 446 } 447 L[size(L)+1] = p; 448 } 449 } 450 else 451 { 452 int neededSize = size(L)+n-1;; 453 list parallelResults; 454 list arguments; 455 int neededPrimes = neededSize-size(L); 456 idealPrimeTest Id; 457 Id.Ideal = I; 458 export(Id); 459 while(neededPrimes > 0) 460 { 461 arguments = list(); 462 for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0)) 463 *ncores; i > 0; i--) 464 { 465 p = prime(p-1); 466 if(p == 2) { ERROR("no more primes"); } 467 arguments[i] = list("Id", p); 468 } 469 parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores); 470 for(i = size(arguments); i > 0; i--) 471 { 472 if(parallelResults[i]) 473 { 474 L[size(L)+1] = arguments[i][2]; 475 } 476 } 477 neededPrimes = neededSize-size(L); 478 } 479 kill Id; 480 if(size(L) > neededSize) 481 { 482 L = L[1..neededSize]; 483 } 484 } 485 return(L); 486 } 487 example 488 { "EXAMPLE:"; echo = 2; 489 ring r = 0,(x,y,z),dp; 490 ideal I = 2147483647x+y, z-181; 491 intvec L = primeList(I,10); 492 size(L); 493 L[1]; 494 L[size(L)]; 495 L = primeList(I,5,L); 496 size(L); 497 L[size(L)]; 498 } 499 500 //////////////////////////////////////////////////////////////////////////////// 501 502 static proc liftstd1(ideal I) 503 { 504 def R = basering; 505 list rl = ringlist(R); 506 list ordl = rl[3]; 507 508 int i; 509 for(i = 1; i <= size(ordl); i++) 510 { 511 if((ordl[i][1] == "C") || (ordl[i][1] == "c")) 512 { 513 ordl = delete(ordl, i); 514 break; 515 } 516 } 517 518 ordl = insert(ordl, list("c", 0)); 519 rl[3] = ordl; 520 def newR = ring(rl); 521 setring newR; 522 ideal I = imap(R,I); 523 524 intvec opt = option(get); 525 option(none); 526 option(prompt); 527 528 module M; 529 for(i = 1; i <= size(I); i++) 530 { 531 M = M + module(I[i]*gen(1) + gen(i+1)); 532 M = M + module(gen(i+1)); 533 } 534 535 module sM = std(M); 536 537 ideal sI; 538 if(attrib(R,"global")) 539 { 540 for(i = size(I)+1; i <= size(sM); i++) 541 { 542 sI[size(sI)+1] = sM[i][1]; 543 } 544 matrix T = submat(sM,2..nrows(sM),size(I)+1..ncols(sM)); 545 } 546 else 547 { 548 //"=========================================================="; 549 //"WARNING: Algorithm is not applicable if ordering is mixed."; 550 //"=========================================================="; 551 for(i = 1; i <= size(sM)-size(I); i++) 552 { 553 sI[size(sI)+1] = sM[i][1]; 554 } 555 matrix T = submat(sM,2..nrows(sM),1..ncols(sM)-size(I)); 556 } 557 558 setring R; 559 option(set, opt); 560 return(imap(newR,sI),imap(newR,T)); 561 } 562 example 563 { "EXAMPLE:"; echo = 2; 564 ring R = 0,(x,y,z),dp; 565 poly f = x3+y7+z2+xyz; 566 ideal i = jacob(f); 567 matrix T; 568 ideal sm = liftstd(i,T); 569 sm; 570 print(T); 571 matrix(sm) - matrix(i)*T; 572 573 574 ring S = 32003, x(1..5), lp; 575 ideal I = cyclic(5); 576 ideal sI; 577 matrix T; 578 sI,T = liftstd1(I); 579 matrix(sI) - matrix(I)*T; 580 } 581 582 //////////////////////////////////////////////////////////////////////////////// 583 584 proc modpStd(ideal I, int p, int variant, list #) 585 "USAGE: modpStd(I,p,variant,#); I ideal, p integer, variant integer 586 ASSUME: If size(#) > 0, then #[1] is an intvec describing the Hilbert series. 587 RETURN: ideal - a standard basis of I mod p, integer - p 588 NOTE: The procedure computes a standard basis of the ideal I modulo p and 589 fetches the result to the basering. If size(#) > 0 the Hilbert driven 590 standard basis computation std(.,#[1]) is used instead of groebner. 591 The standard basis computation modulo p does also vary depending on the 592 integer variant, namely 593 @* - variant = 1: std(.,#[1]) resp. groebner, 594 @* - variant = 2: groebner, 595 @* - variant = 3: homog. - std(.,#[1]) resp. groebner - dehomog., 596 @* - variant = 4: fglm. 597 EXAMPLE: example modpStd; shows an example 598 " 599 { 600 def R0 = basering; 601 list rl = ringlist(R0); 602 rl[1] = p; 603 def @r = ring(rl); 604 setring @r; 605 ideal i = fetch(R0,I); 606 607 option(redSB); 608 609 if(variant == 1) 610 { 611 if(size(#) > 0) 612 { 613 i = std(i, #[1]); 614 } 615 else 616 { 617 i = groebner(i); 618 } 619 } 620 621 if(variant == 2) 622 { 623 i = groebner(i); 624 } 625 626 if(variant == 3) 627 { 628 list rl = ringlist(@r); 629 int nvar@r = nvars(@r); 630 631 int k; 632 intvec w; 633 for(k = 1; k <= nvar@r; k++) 634 { 635 w[k] = deg(var(k)); 636 } 637 w[nvar@r + 1] = 1; 638 639 rl[2][nvar@r + 1] = "homvar"; 640 rl[3][2][2] = w; 641 642 def HomR = ring(rl); 643 setring HomR; 644 ideal i = imap(@r, i); 645 i = homog(i, homvar); 646 647 if(size(#) > 0) 648 { 649 if(w == 1) 650 { 651 i = std(i, #[1]); 652 } 653 else 654 { 655 i = std(i, #[1], w); 656 } 657 } 658 else 659 { 660 i = groebner(i); 661 } 662 663 i = subst(i, homvar, 1); 664 i = simplify(i, 34); 665 666 setring @r; 667 i = imap(HomR, i); 668 i = interred(i); 669 kill HomR; 670 } 671 672 if(variant == 4) 673 { 674 def R1 = changeord(list(list("dp",1:nvars(basering)))); 675 setring R1; 676 ideal i = fetch(@r,i); 677 i = std(i); 678 setring @r; 679 i = fglm(R1,i); 680 } 681 682 setring R0; 683 return(list(fetch(@r,i),p)); 684 } 685 example 686 { "EXAMPLE:"; echo = 2; 687 ring r = 0, x(1..4), dp; 688 ideal I = cyclic(4); 689 int p = 181; 690 list P = modpStd(I,p,2); 691 P; 692 693 ring r2 = 0, x(1..5), lp; 694 ideal I = cyclic(5); 695 int q = 32003; 696 list Q = modpStd(I,q,4); 697 Q; 698 } 699 700 static proc algeModStd(ideal i,list #) 701 { 702 //reduces modStd over algebraic extensions to the one over a polynomial ring 703 list L=#; 704 def R=basering; 705 int n=nvars(R); 706 list rl=ringlist(R); 707 poly p=minpoly; 708 rl[2][n+1]=rl[1][2][1]; 709 rl[1]=rl[1][1]; 710 rl[3][size(rl[3])+1]=rl[3][size(rl[3])]; 711 rl[3][size(rl[3])-1]=list("dp",1); 712 def S=ring(rl); 713 setring S; 714 poly p=imap(R,p); 715 ideal i=imap(R,i); 716 i=i,p; 717 ideal j=modStd(i,L); 718 if(j[1]==p) 719 { 720 j[1]=0; 721 } 722 j=simplify(j,2); 723 setring R; 724 ideal j=imap(S,j); 725 return(j); 726 } 727 ////////////////////////////// main procedures ///////////////////////////////// 728 729 proc modStd(ideal I, list #) 730 "USAGE: modStd(I); I ideal 731 ASSUME: If size(#) > 0, then # contains either 1, 2 or 4 integers such that 732 @* - #[1] is the number of available processors for the computation, 733 @* - #[2] is an optional parameter for the exactness of the computation, 734 if #[2] = 1, the procedure computes a standard basis for sure, 735 @* - #[3] is the number of primes until the first lifting, 736 @* - #[4] is the constant number of primes between two liftings until 737 the computation stops. 738 RETURN: a standard basis of I if no warning appears; 739 NOTE: The procedure computes a standard basis of I (over the rational 740 numbers) by using modular methods. 741 By default the procedure computes a standard basis of I for sure, but 742 if the optional parameter #[2] = 0, it computes a standard basis of I 743 with high probability. 744 The procedure distinguishes between different variants for the standard 745 basis computation in positive characteristic depending on the ordering 746 of the basering, the parameter #[2] and if the ideal I is homogeneous. 747 @* - variant = 1, if I is homogeneous, 748 @* - variant = 2, if I is not homogeneous, 1-block-ordering, 749 @* - variant = 3, if I is not homogeneous, complicated ordering (lp or 750 > 1 block), 751 @* - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0. 752 EXAMPLE: example modStd; shows an example 753 " 754 { 755 int TT = timer; 756 int RT = rtimer; 757 758 def R0 = basering; 759 list rl = ringlist(R0); 760 761 int algebraic; 762 if(size(#)>0) 763 { 764 if(#[1]<=0) 765 { 766 algebraic=1; 767 #[1]=-#[1]; 768 if(#[1]==0){list LK;#=LK;} 769 } 770 } 771 772 if((npars(R0) > 0) || (rl[1][1] > 0)) 773 { 774 if(npars(R0)==1) 775 { 776 if(minpoly!=0) 777 { 778 list LM=#; 779 if(size(LM)==0){LM[1]=0;} 780 LM[1]=-LM[1]; 781 return(algeModStd(I,LM)); 782 } 783 } 784 785 ERROR("Characteristic of basering should be zero, basering should 786 have no parameters."); 787 } 788 789 int index = 1; 790 int i,k,c; 791 int j = 1; 792 int pTest, sizeTest; 793 int en = 2134567879; 794 int an = 1000000000; 795 bigint N; 796 797 //-------------------- Initialize optional parameters ------------------------ 798 if(size(#) > 0) 799 { 800 if(size(#) == 1) 801 { 802 int n1 = #[1]; 803 int exactness = 1; 804 if(n1 >= 10) 805 { 806 int n2 = n1 + 1; 807 int n3 = n1; 808 } 809 else 810 { 811 int n2 = 10; 812 int n3 = 10; 813 } 814 } 815 if(size(#) == 2) 816 { 817 int n1 = #[1]; 818 int exactness = #[2]; 819 if(n1 >= 10) 820 { 821 int n2 = n1 + 1; 822 int n3 = n1; 823 } 824 else 825 { 826 int n2 = 10; 827 int n3 = 10; 828 } 829 } 830 if(size(#) == 4) 831 { 832 int n1 = #[1]; 833 int exactness = #[2]; 834 if(n1 >= #[3]) 835 { 836 int n2 = n1 + 1; 837 } 838 else 839 { 840 int n2 = #[3]; 841 } 842 if(n1 >= #[4]) 843 { 844 int n3 = n1; 845 } 846 else 847 { 848 int n3 = #[4]; 849 } 850 } 851 } 852 else 853 { 854 int n1 = 1; 855 int exactness = 1; 856 int n2 = 10; 857 int n3 = 10; 858 } 859 860 if(printlevel >= 10) 861 { 862 "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3) 863 +", exactness = "+string(exactness); 864 } 865 866 //------------------------- Save current options ----------------------------- 867 intvec opt = option(get); 868 869 option(redSB); 870 871 //-------------------- Initialize the list of primes ------------------------- 872 int tt = timer; 873 int rt = rtimer; 874 intvec L = primeList(I,n2,n1); 875 if(printlevel >= 10) 876 { 877 "CPU-time for primeList: "+string(timer-tt)+" seconds."; 878 "Real-time for primeList: "+string(rtimer-rt)+" seconds."; 879 } 880 L[5] = prime(random(an,en)); 881 882 //--------------------- Decide which variant to take ------------------------- 883 int variant; 884 int h = homog(I); 885 886 tt = timer; 887 rt = rtimer; 888 889 if(!hasMixedOrdering()) 890 { 891 if(h) 892 { 893 variant = 1; 894 if(printlevel >= 10) { "variant = 1"; } 895 896 rl[1] = L[5]; 897 def @r = ring(rl); 898 setring @r; 899 def @s = changeord(list(list("dp",1:nvars(basering)))); 900 setring @s; 901 ideal I = std(fetch(R0,I)); 902 intvec hi = hilb(I,1); 903 setring R0; 904 kill @r,@s; 905 } 906 else 907 { 908 string ordstr_R0 = ordstr(R0); 909 int neg = 1 - attrib(R0,"global"); 910 911 if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg) 912 { 913 variant = 2; 914 if(printlevel >= 10) { "variant = 2"; } 915 } 916 else 917 { 918 string order; 919 if(system("nblocks") <= 2) 920 { 921 if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") 922 + find(ordstr_R0, "rp") <= 0) 923 { 924 order = "simple"; 925 } 926 } 927 928 if((order == "simple") || (size(rl) > 4)) 929 { 930 variant = 2; 931 if(printlevel >= 10) { "variant = 2"; } 932 } 933 else 934 { 935 rl[1] = L[5]; 936 def @r = ring(rl); 937 setring @r; 938 939 def @s = changeord(list(list("dp",1:nvars(basering)))); 940 setring @s; 941 ideal I = std(fetch(R0,I)); 942 if(dim(I) == 0) 943 { 944 variant = 4; 945 if(printlevel >= 10) { "variant = 4"; } 946 } 947 else 948 { 949 variant = 3; 950 if(printlevel >= 10) { "variant = 3"; } 951 952 int nvar@r = nvars(@r); 953 intvec w; 954 for(i = 1; i <= nvar@r; i++) 955 { 956 w[i] = deg(var(i)); 957 } 958 w[nvar@r + 1] = 1; 959 960 list hiRi = hilbRing(fetch(R0,I),w); 961 intvec W = hiRi[2]; 962 @s = hiRi[1]; 963 setring @s; 964 965 Id(1) = std(Id(1)); 966 intvec hi = hilb(Id(1), 1, W); 967 } 968 969 setring R0; 970 kill @r,@s; 971 } 972 } 973 } 974 } 975 else 976 { 977 if(exactness == 1) { return(groebner(I)); } 978 if(h) 979 { 980 variant = 1; 981 if(printlevel >= 10) { "variant = 1"; } 982 rl[1] = L[5]; 983 def @r = ring(rl); 984 setring @r; 985 def @s = changeord(list(list("dp",1:nvars(basering)))); 986 setring @s; 987 ideal I = std(fetch(R0,I)); 988 intvec hi = hilb(I,1); 989 setring R0; 990 kill @r,@s; 991 } 992 else 993 { 994 string ordstr_R0 = ordstr(R0); 995 int neg = 1 - attrib(R0,"global"); 996 997 if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg) 998 { 999 variant = 2; 1000 if(printlevel >= 10) { "variant = 2"; } 1001 } 1002 else 1003 { 1004 string order; 1005 if(system("nblocks") <= 2) 1006 { 1007 if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") 1008 + find(ordstr_R0, "rp") <= 0) 1009 { 1010 order = "simple"; 1011 } 1012 } 1013 1014 if((order == "simple") || (size(rl) > 4)) 1015 { 1016 variant = 2; 1017 if(printlevel >= 10) { "variant = 2"; } 1018 } 1019 else 1020 { 1021 variant = 3; 1022 if(printlevel >= 10) { "variant = 3"; } 1023 1024 rl[1] = L[5]; 1025 def @r = ring(rl); 1026 setring @r; 1027 int nvar@r = nvars(@r); 1028 intvec w; 1029 for(i = 1; i <= nvar@r; i++) 1030 { 1031 w[i] = deg(var(i)); 1032 } 1033 w[nvar@r + 1] = 1; 1034 1035 list hiRi = hilbRing(fetch(R0,I),w); 1036 intvec W = hiRi[2]; 1037 def @s = hiRi[1]; 1038 setring @s; 1039 1040 Id(1) = std(Id(1)); 1041 intvec hi = hilb(Id(1), 1, W); 1042 1043 setring R0; 1044 kill @r,@s; 1045 } 1046 } 1047 } 1048 } 1049 if(algebraic){variant=2;} 1050 1051 list P,T1,T2,T3,LL; 1052 1053 ideal J,K,H; 1054 1055 //----- If there is more than one processor available, we parallelize the ---- 1056 //----- main standard basis computations in positive characteristic ---- 1057 1058 if(n1 > 1) 1059 { 1060 ideal I_for_fork = I; 1061 export(I_for_fork); // I available for each link 1062 1063 //----- Create n1 links l(1),...,l(n1), open all of them and compute --------- 1064 //----- standard basis for the primes L[2],...,L[n1 + 1]. --------- 1065 1066 for(i = 1; i <= n1; i++) 1067 { 1068 //link l(i) = "MPtcp:fork"; 1069 link l(i) = "ssi:fork"; 1070 open(l(i)); 1071 if((variant == 1) || (variant == 3)) 1072 { 1073 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), 1074 eval(variant), eval(hi)))); 1075 } 1076 if((variant == 2) || (variant == 4)) 1077 { 1078 write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]), 1079 eval(variant)))); 1080 } 1081 } 1082 1083 int t = timer; 1084 if((variant == 1) || (variant == 3)) 1085 { 1086 P = modpStd(I_for_fork, L[1], variant, hi); 1087 } 1088 if((variant == 2) || (variant == 4)) 1089 { 1090 P = modpStd(I_for_fork, L[1], variant); 1091 } 1092 t = timer - t; 1093 if(t > 60) { t = 60; } 1094 int i_sleep = system("sh", "sleep "+string(t)); 1095 T1[1] = P[1]; 1096 T2[1] = bigint(P[2]); 1097 index++; 1098 1099 j = j + n1 + 1; 1100 } 1101 1102 //-------------- Main standard basis computations in positive ---------------- 1103 //---------------------- characteristic start here --------------------------- 1104 1105 list arguments_farey, results_farey; 1106 1107 while(1) 1108 { 1109 tt = timer; rt = rtimer; 1110 1111 if(printlevel >= 10) { "size(L) = "+string(size(L)); } 1112 1113 if(n1 > 1) 1114 { 1115 while(j <= size(L) + 1) 1116 { 1117 for(i = 1; i <= n1; i++) 1118 { 1119 //--- ask if link l(i) is ready otherwise sleep for t seconds --- 1120 if(status(l(i), "read", "ready")) 1121 { 1122 //--- read the result from l(i) --- 1123 P = read(l(i)); 1124 T1[index] = P[1]; 1125 T2[index] = bigint(P[2]); 1126 index++; 1127 1128 if(j <= size(L)) 1129 { 1130 if((variant == 1) || (variant == 3)) 1131 { 1132 write(l(i), quote(modpStd(I_for_fork, eval(L[j]), 1133 eval(variant), eval(hi)))); 1134 j++; 1135 } 1136 if((variant == 2) || (variant == 4)) 1137 { 1138 write(l(i), quote(modpStd(I_for_fork, 1139 eval(L[j]), eval(variant)))); 1140 j++; 1141 } 1142 } 1143 else 1144 { 1145 k++; 1146 close(l(i)); 1147 } 1148 } 1149 } 1150 //--- k describes the number of closed links --- 1151 if(k == n1) 1152 { 1153 j++; 1154 } 1155 i_sleep = system("sh", "sleep "+string(t)); 1156 } 1157 } 1158 else 1159 { 1160 while(j <= size(L)) 1161 { 1162 if((variant == 1) || (variant == 3)) 1163 { 1164 P = modpStd(I, L[j], variant, hi); 1165 } 1166 if((variant == 2) || (variant == 4)) 1167 { 1168 P = modpStd(I, L[j], variant); 1169 } 1170 1171 T1[index] = P[1]; 1172 T2[index] = bigint(P[2]); 1173 index++; 1174 j++; 1175 } 1176 } 1177 1178 if(printlevel >= 10) 1179 { 1180 "CPU-time for computing list is "+string(timer - tt)+" seconds."; 1181 "Real-time for computing list is "+string(rtimer - rt)+" seconds."; 1182 } 1183 1184 //------------------------ Delete unlucky primes ----------------------------- 1185 //------------- unlucky if and only if the leading ideal is wrong ------------ 1186 1187 LL = deleteUnluckyPrimes(T1,T2,h); 1188 T1 = LL[1]; 1189 T2 = LL[2]; 1190 1191 //------------------- Now all leading ideals are the same -------------------- 1192 //------------------- Lift results to basering via farey --------------------- 1193 1194 tt = timer; rt = rtimer; 1195 N = T2[1]; 1196 for(i = 2; i <= size(T2); i++) { N = N*T2[i]; } 1197 H = chinrem(T1,T2); 1198 if(n1 == 1) 1199 { 1200 J = farey(H,N); 1201 } 1202 else 1203 { 1204 for(i = ncols(H); i > 0; i--) 1205 { 1206 arguments_farey[i] = list(ideal(H[i]), N); 1207 } 1208 results_farey = parallelWaitAll("farey", arguments_farey, 0, n1); 1209 for(i = ncols(H); i > 0; i--) 1210 { 1211 J[i] = results_farey[i][1]; 1212 } 1213 } 1214 if(printlevel >= 10) 1215 { 1216 "CPU-time for lifting-process is "+string(timer - tt)+" seconds."; 1217 "Real-time for lifting-process is "+string(rtimer - rt)+" seconds."; 1218 } 1219 1220 //---------------- Test if we already have a standard basis of I -------------- 1221 1222 tt = timer; rt = rtimer; 1223 if((variant == 1) || (variant == 3)) 1224 { 1225 pTest = pTestSB(I,J,L,variant,hi); 1226 } 1227 if((variant == 2) || (variant == 4)) 1228 { 1229 pTest = pTestSB(I,J,L,variant); 1230 } 1231 1232 if(printlevel >= 10) 1233 { 1234 "CPU-time for pTest is "+string(timer - tt)+" seconds."; 1235 "Real-time for pTest is "+string(rtimer - rt)+" seconds."; 1236 } 1237 1238 if(pTest) 1239 { 1240 if(printlevel >= 10) 1241 { 1242 "CPU-time for computation without final tests is " 1243 +string(timer - TT)+" seconds."; 1244 "Real-time for computation without final tests is " 1245 +string(rtimer - RT)+" seconds."; 1246 } 1247 1248 attrib(J,"isSB",1); 1249 1250 if(exactness == 0) 1251 { 1252 option(set, opt); 1253 if(n1 > 1) { kill I_for_fork; } 1254 return(J); 1255 } 1256 1257 if(exactness == 1) 1258 { 1259 tt = timer; rt = rtimer; 1260 sizeTest = 1 - isIncluded(I,J,n1); 1261 1262 if(printlevel >= 10) 1263 { 1264 "CPU-time for checking if I subset <G> is " 1265 +string(timer - tt)+" seconds."; 1266 "Real-time for checking if I subset <G> is " 1267 +string(rtimer - rt)+" seconds."; 1268 } 1269 1270 if(sizeTest == 0) 1271 { 1272 tt = timer; rt = rtimer; 1273 K = std(J); 1274 1275 if(printlevel >= 10) 1276 { 1277 "CPU-time for last std-computation is " 1278 +string(timer - tt)+" seconds."; 1279 "Real-time for last std-computation is " 1280 +string(rtimer - rt)+" seconds."; 1281 } 1282 1283 if(size(reduce(K,J)) == 0) 1284 { 1285 option(set, opt); 1286 if(n1 > 1) { kill I_for_fork; } 1287 return(J); 1288 } 1289 } 1290 } 1291 } 1292 1293 //-------------- We do not already have a standard basis of I ---------------- 1294 //----------- Therefore do the main computation for more primes -------------- 1295 1296 T1 = H; 1297 T2 = N; 1298 index = 2; 1299 1300 j = size(L) + 1; 1301 tt = timer; rt = rtimer; 1302 L = primeList(I,n3,L,n1); 1303 if(printlevel >= 10) 1304 { 1305 "CPU-time for primeList: "+string(timer-tt)+" seconds."; 1306 "Real-time for primeList: "+string(rtimer-rt)+" seconds."; 1307 } 1308 1309 if(n1 > 1) 1310 { 1311 for(i = 1; i <= n1; i++) 1312 { 1313 open(l(i)); 1314 if((variant == 1) || (variant == 3)) 1315 { 1316 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]), 1317 eval(variant), eval(hi)))); 1318 } 1319 if((variant == 2) || (variant == 4)) 1320 { 1321 write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]), 1322 eval(variant)))); 1323 } 1324 } 1325 j = j + n1; 1326 k = 0; 1327 } 1328 } 1329 } 1330 example 1331 { "EXAMPLE:"; echo = 2; 1332 ring R1 = 0,(x,y,z,t),dp; 1333 ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; 1334 ideal J = modStd(I); 1335 J; 1336 I = homog(I,t); 1337 J = modStd(I); 1338 J; 1339 1340 ring R2 = 0,(x,y,z),ds; 1341 ideal I = jacob(x5+y6+z7+xyz); 1342 ideal J1 = modStd(I,1,0); 1343 J1; 1344 1345 ring R3 = 0,x(1..4),lp; 1346 ideal I = cyclic(4); 1347 ideal J1 = modStd(I,1); 1348 ideal J2 = modStd(I,1,0); 1349 size(reduce(J1,J2)); 1350 size(reduce(J2,J1)); 1351 } 1352 1353 //////////////////////////////////////////////////////////////////////////////// 1354 1355 proc modS(ideal I, list L, list #) 1356 "USAGE: modS(I,L); I ideal, L intvec of primes 1357 if size(#)>0 std is used instead of groebner 1358 RETURN: an ideal which is with high probability a standard basis 1359 NOTE: This procedure is designed for fast experiments. 1360 It is not tested whether the result is a standard basis. 1361 It is not tested whether the result generates I. 1362 EXAMPLE: example modS; shows an example 1363 " 1364 { 1365 int j; 1366 bigint N = 1; 1367 def R0 = basering; 1368 ideal J; 1369 list T; 1370 list rl = ringlist(R0); 1371 if((npars(R0)>0) || (rl[1]>0)) 1372 { 1373 ERROR("Characteristic of basering should be zero."); 1374 } 1375 for(j = 1; j <= size(L); j++) 1376 { 1377 N = N*L[j]; 1378 rl[1] = L[j]; 1379 def @r = ring(rl); 1380 setring @r; 1381 ideal I = fetch(R0,I); 1382 if(size(#) > 0) 1383 { 1384 I = std(I); 1385 } 1386 else 1387 { 1388 I = groebner(I); 1389 } 1390 setring R0; 1391 T[j] = fetch(@r,I); 1392 kill @r; 1393 } 1394 L = deleteUnluckyPrimes(T,L,homog(I)); 1395 // unlucky if and only if the leading ideal is wrong 1396 J = farey(chinrem(L[1],L[2]),N); 1397 attrib(J,"isSB",1); 1398 return(J); 1399 } 1400 example 1401 { "EXAMPLE:"; echo = 2; 1402 list L = 3,5,11,13,181,32003; 1403 ring r = 0,(x,y,z,t),dp; 1404 ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4; 1405 I = homog(I,t); 1406 ideal J = modS(I,L); 1407 J; 1408 } 1409 1410 //////////////////////////////////////////////////////////////////////////////// 1411 1412 proc modHenselStd(ideal I, list #) 1413 "USAGE: modHenselStd(I); 1414 RETURN: a standard basis of I; 1415 NOTE: The procedure computes a standard basis of I (over the rational 1416 numbers) by using modular computations and Hensellifting. 1417 For further experiments see procedure modS. 1418 EXAMPLE: example modHenselStd; shows an example 1419 " 1420 { 1421 int i,j; 1422 1423 bigint p = 2134567879; 1424 if(size(#)!=0) { p=#[1]; } 1425 while(!primeTest(I,p)) 1426 { 1427 p = prime(random(2000000000,2134567879)); 1428 } 1429 1430 def R = basering; 1431 module F,PrevG,PrevZ,Z2; 1432 ideal testG,testG1,G1,G2,G3,Gp; 1433 list L = p; 1434 list rl = ringlist(R); 1435 rl[1] = int(p); 1436 1437 def S = ring(rl); 1438 setring S; 1439 intvec opt = option(get); 1440 option(redSB); 1441 module Z,M,Z2; 1442 ideal I = imap(R,I); 1443 ideal Gp,G1,G2,G3; 1444 Gp,Z = liftstd1(I); 1445 attrib(Gp,"isSB",1); 1446 module ZZ = syz(I); 1447 attrib(ZZ,"isSB",1); 1448 Z = reduce(Z,ZZ); 1449 1450 setring R; 1451 Gp = imap(S,Gp); 1452 PrevZ = imap(S,Z); 1453 PrevG = module(Gp); 1454 F = module(I); 1455 testG = farey(Gp,p); 1456 attrib(testG,"isSB",1); 1457 while(1) 1458 { 1459 i++; 1460 G1 = ideal(1/(p^i) * sum_id(F*PrevZ,(-1)*PrevG)); 1461 setring S; 1462 G1 = imap(R,G1); 1463 G2 = reduce(G1,Gp); 1464 G3 = sum_id(G1,(-1)*G2); 1465 M = lift(Gp,G3); 1466 Z2 = (-1)*Z*M; 1467 1468 setring R; 1469 G2 = imap(S,G2); 1470 Z2 = imap(S,Z2); 1471 PrevG = sum_id(PrevG, module(p^i*G2)); 1472 PrevZ = sum_id(PrevZ, multiply(poly(p^i),Z2)); 1473 testG1 = farey(ideal(PrevG),p^(i+1)); 1474 attrib(testG1,"isSB",1); 1475 if(size(reduce(testG1,testG)) == 0) 1476 { 1477 if(size(reduce(I,testG1)) == 0) // I is in testG1 1478 { 1479 if(pTestSB(I,testG1,L,2)) 1480 { 1481 G3 = std(testG1); // testG1 is SB 1482 if(size(reduce(G3,testG1)) == 0) 1483 { 1484 option(set, opt); 1485 return(G3); 1486 } 1487 } 1488 } 1489 } 1490 testG = testG1; 1491 attrib(testG,"isSB",1); 1492 } 1493 } 1494 example 1495 { "EXAMPLE:"; echo = 2; 1496 ring r = 0,(x,y,z),dp; 1497 ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4; 1498 ideal J = modHenselStd(I); 1499 J; 1500 } 1501 1502 //////////////////////////////////////////////////////////////////////////////// 1503 1504 static proc sum_id(list #) 1505 { 1506 if(typeof(#[1])=="ideal") 1507 { 1508 ideal M; 1509 } 1510 else 1511 { 1512 module M; 1513 } 1514 1515 int i; 1516 for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; } 1517 return(M); 1518 } 1519 1520 //////////////////////////////////////////////////////////////////////////////// 1521 1522 static proc multiply(poly p, list #) 1523 { 1524 if(typeof(#[1])=="ideal") 1525 { 1526 ideal M; 1527 } 1528 else 1529 { 1530 module M; 1531 } 1532 1533 int i; 1534 for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; } 1535 return(M); 244 } 245 } 246 247 /* compute command(args) */ 248 execute("Ip = "+command+"(Ip);"); 249 250 /* test if Gp is in Ip */ 251 for (i = ncols(Gp); i > 0; i--) { 252 if (reduce(Gp[i], Ip, 1) != 0) { 253 setring(br); 254 return(0); 255 } 256 } 257 setring(br); 258 return(1); 259 } 260 261 /* test if 'result' is a GB of the input ideal */ 262 static proc finalTest_std(string command, alias list args, ideal result) 263 { 264 /* test if args[1] is in result */ 265 attrib(result, "isSB", 1); 266 int i; 267 for (i = ncols(args[1]); i > 0; i--) { 268 if (reduce(args[1][i], result, 1) != 0) { 269 return(0); 270 } 271 } 272 273 /* test if result is a GB */ 274 ideal G = std(result); 275 if (reduce_parallel(G, result)) { 276 return(0); 277 } 278 return(1); 279 } 280 281 /* return 1, if I_reduce is _not_ in G_reduce, 282 * 0, otherwise 283 * (same as size(reduce(I_reduce, G_reduce))). 284 * Uses parallelization. */ 285 static proc reduce_parallel(ideal I_reduce, ideal G_reduce) 286 { 287 exportto(Modstd, I_reduce); 288 exportto(Modstd, G_reduce); 289 int size_I = ncols(I_reduce); 290 int chunks = Modular::par_range(size_I); 291 intvec range; 292 int i; 293 for (i = chunks; i > 0; i--) { 294 range = Modular::par_range(size_I, i); 295 task t(i) = "Modstd::reduce_task", list(range); 296 } 297 startTasks(t(1..chunks)); 298 waitAllTasks(t(1..chunks)); 299 int result = 0; 300 for (i = chunks; i > 0; i--) { 301 if (getResult(t(i))) { 302 result = 1; 303 break; 304 } 305 } 306 kill I_reduce; 307 kill G_reduce; 308 return(result); 309 } 310 311 /* compute a chunk of reductions for reduce_parallel */ 312 static proc reduce_task(intvec range) 313 { 314 int result = 0; 315 int i; 316 for (i = range[1]; i <= range[2]; i++) { 317 if (reduce(I_reduce[i], G_reduce, 1) != 0) { 318 result = 1; 319 break; 320 } 321 } 322 return(result); 1536 323 } 1537 324
Note: See TracChangeset
for help on using the changeset viewer.