Changeset f2f858 in git
 Timestamp:
 Jun 2, 2014, 3:52:56 PM (9 years ago)
 Branches:
 (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
 Children:
 9f1bf0bd7d089909777bd2491baed90ca7fcfeaf
 Parents:
 175eb21f4e380c8139d078ea81f56e1483e54df4
 gitauthor:
 Andreas Steenpass <steenpass@mathematik.unikl.de>20140602 15:52:56+02:00
 gitcommitter:
 Andreas Steenpass <steenpass@mathematik.unikl.de>20140603 10:42:24+02:00
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/assprimeszerodim.lib
r175eb21 rf2f858 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="version assprimeszerodim.lib 4.0.0.0 Jun_201 3"; // $Id$3 category ="Commutative Algebra";2 version="version assprimeszerodim.lib 4.0.0.0 Jun_2014 "; // $Id$ 3 category="Commutative Algebra"; 4 4 info=" 5 5 LIBRARY: assprimeszerodim.lib associated primes of a zerodimensional ideal 6 6 7 7 AUTHORS: N. Idrees nazeranjawwad@gmail.com 8 @* G. Pfister pfister@mathematik.unikl.de 9 @* S. Steidel steidel@mathematik.unikl.de 8 G. Pfister pfister@mathematik.unikl.de 9 A. Steenpass steenpass@mathematik.unikl.de 10 S. Steidel steidel@mathematik.unikl.de 10 11 11 12 OVERVIEW: … … 28 29 29 30 proc zeroRadical(ideal I, list #) 30 "USAGE: zeroRadical(I,[n]); I ideal, optional: n number of processors (for 31 parallel computing) 31 "USAGE: zeroRadical(I[, exactness]); I ideal, exactness int 32 32 ASSUME: I is zerodimensional in Q[variables] 33 33 RETURN: the radical of I 34 NOTE: A final test is applied to the result if exactness != 0 (default), 35 otherwise no final test is done. 34 36 EXAMPLE: example zeroRadical; shows an example 35 37 " 36 38 { 37 return(zeroR(modStd(I,#),#)); 39 /* read optional parameter */ 40 int exactness = 1; 41 if(size(#) > 0) 42 { 43 if(size(#) > 1  typeof(#[1]) != "int") 44 { 45 ERROR("wrong optional parameter"); 46 } 47 exactness = #[1]; 48 } 49 50 /* compute a standard basis if necessary */ 51 if (!attrib(I, "isSB")) 52 { 53 I = modStd(I, exactness); 54 } 55 56 /* call modular() */ 57 // TODO: write deleteUnluckyPrimes_zeroRadical() 58 if(exactness) 59 { 60 ideal F = modular("Assprimeszerodim::zeroRadP", list(I), 61 Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default, 62 pTest_zeroRadical, finalTest_zeroRadical); 63 } 64 else 65 { 66 ideal F = modular("Assprimeszerodim::zeroRadP", list(I), 67 Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default, 68 pTest_zeroRadical); 69 } 70 71 /* compute the squarefree parts */ 72 poly f; 73 int k; 74 int i; 75 for(i = nvars(basering); i > 0; i) 76 { 77 f = gcd(F[i], diff(F[i], var(i))); 78 k = k + deg(f); 79 F[i] = F[i]/f; 80 } 81 82 /* return the result */ 83 if(k == 0) 84 { 85 return(I); 86 } 87 else 88 { 89 return(modStd(I + F, exactness)); 90 } 38 91 } 39 92 example … … 46 99 //////////////////////////////////////////////////////////////////////////////// 47 100 48 static proc zeroR(ideal I, list #) 49 // compute the radical of I provided that I is zerodimensional in Q[variables] 50 // and a standard basis 51 { 52 attrib(I,"isSB",1); 53 int i, k; 54 int j = 1; 55 int index = 1; 56 int crit; 57 58 list CO1, CO2, P; 59 ideal G, F; 60 bigint N; 61 poly f; 62 63 // Initialize optional parameter  64 if(size(#) > 0) 65 { 66 int n1 = #[1]; 67 if(n1 >= 10) 68 { 69 int n2 = n1 + 1; 70 int n3 = n1; 71 } 72 else 73 { 74 int n2 = 10; 75 int n3 = 10; 76 } 101 /* The pTest for zeroRadical(), to be used in modular(). */ 102 static proc pTest_zeroRadical(string command, list args, ideal result, int p) 103 { 104 /* change to characteristic p */ 105 def br = basering; 106 list lbr = ringlist(br); 107 if(typeof(lbr[1]) == "int") 108 { 109 lbr[1] = p; 77 110 } 78 111 else 79 112 { 80 int n1 = 1; 81 int n2 = 10; 82 int n3 = 10; 83 } 84 85 // Initialize the list of primes  86 intvec L = primeList(I,n2); 87 L[5] = prime(random(100000000,536870912)); 88 89 if(n1 > 1) 90 { 91 92 // Create n links l(1),...,l(n1), open all of them and compute  93 // polynomial F for the primes L[2],...,L[n1 + 1].  94 95 for(i = 1; i <= n1; i++) 96 { 97 //link l(i) = "MPtcp:fork"; 98 link l(i) = "ssi:fork"; 99 open(l(i)); 100 write(l(i), quote(zeroRadP(eval(I), eval(L[i + 1])))); 101 } 102 103 int t = timer; 104 P = zeroRadP(I, L[1]); 105 t = timer  t; 106 if(t > 60) { t = 60; } 107 int i_sleep = system("sh", "sleep "+string(t)); 108 CO1[index] = P[1]; 109 CO2[index] = bigint(P[2]); 110 index++; 111 112 j = j + n1 + 1; 113 } 114 115 // Main computations in positive characteristic start here  116 117 while(!crit) 118 { 119 if(n1 > 1) 120 { 121 while(j <= size(L) + 1) 122 { 123 for(i = 1; i <= n1; i++) 124 { 125 if(status(l(i), "read", "ready")) 126 { 127 // read the result from l(i)  128 P = read(l(i)); 129 CO1[index] = P[1]; 130 CO2[index] = bigint(P[2]); 131 index++; 132 133 if(j <= size(L)) 134 { 135 write(l(i), quote(zeroRadP(eval(I), eval(L[j])))); 136 j++; 137 } 138 else 139 { 140 k++; 141 close(l(i)); 142 } 143 } 144 } 145 // k describes the number of closed links  146 if(k == n1) 147 { 148 j++; 149 } 150 // sleep for t seconds  151 i_sleep = system("sh", "sleep "+string(t)); 152 } 153 } 154 else 155 { 156 while(j <= size(L)) 157 { 158 P = zeroRadP(I, L[j]); 159 CO1[index] = P[1]; 160 CO2[index] = bigint(P[2]); 161 index++; 162 j++; 163 } 164 } 165 166 // insert deleteUnluckyPrimes 167 G = chinrem(CO1,CO2); 168 N = CO2[1]; 169 for(i = 2; i <= size(CO2); i++){ N = N*CO2[i]; } 170 F = farey(G,N); 171 172 crit = 1; 173 for(i = 1; i <= nvars(basering); i++) 174 { 175 if(reduce(F[i],I) != 0) { crit = 0; break; } 176 } 177 178 if(!crit) 179 { 180 CO1 = G; 181 CO2 = N; 182 index = 2; 183 184 j = size(L) + 1; 185 L = primeList(I,n3,L); 186 187 if(n1 > 1) 188 { 189 for(i = 1; i <= n1; i++) 190 { 191 open(l(i)); 192 write(l(i), quote(zeroRadP(eval(I), eval(L[j+i1])))); 193 } 194 j = j + n1; 195 k = 0; 196 } 197 } 198 } 199 200 k = 0; 201 for(i = 1; i <= nvars(basering); i++) 202 { 203 f = gcd(F[i],diff(F[i],var(i))); 204 k = k + deg(f); 205 F[i] = F[i]/f; 206 } 207 208 if(k == 0) { return(I); } 209 else { return(modStd(I + F, n1)); } 210 } 211 212 //////////////////////////////////////////////////////////////////////////////// 213 214 proc assPrimes(list #) 215 "USAGE: assPrimes(I,[n],[a]); I ideal or module, 216 optional: int n: number of processors (for parallel computing), int a: 217  a = 1: method of Eisenbud/Hunecke/Vasconcelos 218  a = 2: method of Gianni/Trager/Zacharias 219  a = 3: method of Monico 220 assPrimes(I) chooses n = a = 1 113 lbr[1][1] = p; 114 } 115 def rp = ring(lbr); 116 setring(rp); 117 ideal Ip = fetch(br, args)[1]; 118 ideal Fp_result = fetch(br, result); 119 120 /* run the command and compare with given result */ 121 execute("ideal Fp = "+command+"(Ip);"); 122 int i; 123 for(i = nvars(br); i > 0; i) 124 { 125 if(Fp[i] != Fp_result[i]) 126 { 127 setring(br); 128 return(0); 129 } 130 } 131 setring(br); 132 return(1); 133 } 134 135 //////////////////////////////////////////////////////////////////////////////// 136 137 /* The finalTest for zeroRadical, to be used in modular(). */ 138 static proc finalTest_zeroRadical(string command, list args, ideal F) 139 { 140 int i; 141 for(i = nvars(basering); i > 0; i) 142 { 143 if(reduce(F[i], args[1]) != 0) { return(0); } 144 } 145 return(1); 146 } 147 148 //////////////////////////////////////////////////////////////////////////////// 149 150 proc assPrimes(def I, list #) 151 "USAGE: assPrimes(I[, alg, exactness]); I ideal or module, 152 alg string (optional), exactness int (optional) 153  alg = "GTZ": method of Gianni/Trager/Zacharias (default) 154  alg = "EHV": method of Eisenbud/Hunecke/Vasconcelos 155  alg = "Monico": method of Monico 221 156 ASSUME: I is zerodimensional over Q[variables] 222 RETURN: a list Re of associated primes of I: 157 RETURN: a list of the associated primes of I 158 NOTE: A final test is applied to the result if exactness != 0 (default), 159 otherwise no final test is done. 223 160 EXAMPLE: example assPrimes; shows an example 224 161 " 225 162 { 226 ideal I; 227 if(typeof(#[1]) == "ideal") 228 { 229 I = #[1]; 230 } 231 else 232 { 233 module M = #[1]; 234 I = Ann(M); 235 } 236 237 // Initialize optional parameter  238 if(size(#) > 1) 239 { 240 if(size(#) == 2) 241 { 242 int n1 = #[2]; 243 int alg = 1; 244 if(n1 >= 10) 245 { 246 int n2 = n1 + 1; 247 int n3 = n1; 248 } 249 else 250 { 251 int n2 = 10; 252 int n3 = 10; 253 } 254 } 255 if(size(#) == 3) 256 { 257 int n1 = #[2]; 258 int alg = #[3]; 259 if(n1 >= 10) 260 { 261 int n2 = n1 + 1; 262 int n3 = n1; 263 } 264 else 265 { 266 int n2 = 10; 267 int n3 = 10; 268 } 269 } 270 } 271 else 272 { 273 int n1 = 1; 274 int alg = 1; 275 int n2 = 10; 276 int n3 = 10; 277 } 278 279 if(printlevel >= 10) 280 { 281 "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3); 282 } 283 284 int T = timer; 285 int RT = rtimer; 286 int TT; 287 int t_sleep; 288 289 def SPR = basering; 290 map phi; 291 list H = ideal(0); 292 ideal F; 293 poly F1; 294 163 /* read input */ 164 if(typeof(I) != "ideal") 165 { 166 if(typeof(I) != "module") 167 { 168 ERROR("The first argument must be of type 'ideal' or 'module'."); 169 } 170 module M = I; 171 kill I; 172 ideal I = Ann(M); 173 kill M; 174 } 175 176 /* read optional parameters */ 177 list defaults = list("GTZ", 1); 178 int i; 179 for(i = 1; i <= size(defaults); i++) 180 { 181 if(typeof(#[i]) != typeof(defaults[i])) 182 { 183 # = insert(#, defaults[i], i1); 184 } 185 } 186 if(size(#) != size(defaults)) 187 { 188 ERROR("wrong optional parameters"); 189 } 190 string alg = #[1]; 191 int exactness = #[2]; 192 int a; 193 if(alg == "GTZ") 194 { 195 a = 1; 196 } 197 if(alg == "EHV") 198 { 199 a = 2; 200 } 201 if(alg == "Monico") 202 { 203 a = 3; 204 } 205 if(a == 0) // alg != "GTZ" && alg != "EHV" && alg != "Monico" 206 { 207 ERROR("unknown method"); 208 } 209 210 /* compute a standard basis if necessary */ 295 211 if(printlevel >= 10) { "========== Start modStd =========="; } 296 I = modStd(I); 212 if (!attrib(I, "isSB")) { 213 I = modStd(I, exactness); 214 } 297 215 if(printlevel >= 10) { "=========== End modStd ==========="; } 298 if(printlevel >= 9) { "modStd takes "+string(rtimerRT)+" seconds."; }299 216 int d = vdim(I); 300 217 if(d == 1) { ERROR("Ideal is not zerodimensional."); } 301 218 if(homog(I) == 1) { return(list(maxideal(1))); } 302 poly f = findGen(I); 303 if(printlevel >= 9) { "Coordinate change: "+string(f); } 304 305 if(size(f) == nvars(SPR)) 306 { 307 TT = timer; 308 int spT = pTestRad(d,f,I); 309 if(printlevel >= 9) 310 { 311 "pTestRad(d,f,I) = "+string(spT)+" and takes " 312 +string(timerTT)+" seconds."; 313 } 314 if(!spT) 315 { 316 if(typeof(attrib(#[1],"isRad")) == "int") 317 { 318 if(attrib(#[1],"isRad") == 0) 319 { 320 TT = timer; 321 I = zeroR(I,n1); 322 if(printlevel >= 9) 323 { 324 "zeroR(I,n1) takes "+string(timerTT)+" seconds."; 325 } 326 TT = timer; 327 I = modStd(I); 328 if(printlevel >= 9) 329 { 330 "modStd(I) takes "+string(timerTT)+" seconds."; 331 } 332 d = vdim(I); 333 f = findGen(I); 334 } 335 } 336 else 337 { 338 TT = timer; 339 I = zeroR(I,n1); 340 if(printlevel >= 9) 341 { 342 "zeroR(I,n1) takes "+string(timerTT)+" seconds."; 343 } 344 TT = timer; 345 I = modStd(I); 346 if(printlevel >= 9) 347 { 348 "modStd(I) takes "+string(timerTT)+" seconds."; 349 } 350 d = vdim(I); 351 f = findGen(I); 352 } 353 } 354 } 355 if(printlevel >= 9) 356 { 357 "Realtime for radicalcheck is "+string(rtimer  RT)+" seconds."; 358 "CPUtime for radicalcheck is "+string(timer  T)+" seconds."; 359 } 360 361 export(SPR); 362 poly f_for_fork = f; 363 export(f_for_fork); // f available for each link 364 ideal I_for_fork = I; 365 export(I_for_fork); // I available for each link 366 367 // Initialize the list of primes  368 intvec L = primeList(I,n2); 369 L[5] = prime(random(1000000000,2134567879)); 370 371 list Re; 372 373 ring rHelp = 0,T,dp; 374 list CO1,CO2,P,H; 375 ideal F,G,testF; 376 bigint N; 377 378 list ringL = ringlist(SPR); 379 int i,k,e,int_break,s; 380 int j = 1; 381 int index = 1; 382 383 // If there is more than one processor available, we parallelize the  384 // main standard basis computations in positive characteristic  385 386 if(n1 > 1) 387 { 388 389 // Create n1 links l(1),...,l(n1), open all of them and compute  390 // standard basis for the primes L[2],...,L[n1 + 1].  391 392 for(i = 1; i <= n1; i++) 393 { 394 //link l(i) = "MPtcp:fork"; 395 link l(i) = "ssi:fork"; 396 open(l(i)); 397 write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[i + 1]), 398 eval(alg)))); 399 } 400 401 int t = timer; 402 P = modpSpecialAlgDep(ringL, L[1], alg); 403 t = timer  t; 404 if(t > 60) { t = 60; } 405 int i_sleep = system("sh", "sleep "+string(t)); 406 CO1[index] = P[1]; 407 CO2[index] = bigint(P[2]); 408 index++; 409 410 j = j + n1 + 1; 411 } 412 413 // Main computations in positive characteristic start here  414 415 int tt = timer; 416 int rt = rtimer; 417 418 while(1) 419 { 420 tt = timer; 421 rt = rtimer; 422 423 if(printlevel >= 9) { "size(L) = "+string(size(L)); } 424 425 if(n1 > 1) 426 { 427 while(j <= size(L) + 1) 428 { 429 for(i = 1; i <= n1; i++) 430 { 431 // ask if link l(i) is ready otherwise sleep for t seconds  432 if(status(l(i), "read", "ready")) 433 { 434 // read the result from l(i)  435 P = read(l(i)); 436 CO1[index] = P[1]; 437 CO2[index] = bigint(P[2]); 438 index++; 439 440 if(j <= size(L)) 441 { 442 write(l(i), quote(modpSpecialAlgDep(eval(ringL), 443 eval(L[j]), 444 eval(alg)))); 445 j++; 446 } 447 else 448 { 449 k++; 450 close(l(i)); 451 } 452 } 453 } 454 // k describes the number of closed links  455 if(k == n1) 456 { 457 j++; 458 } 459 i_sleep = system("sh", "sleep "+string(t)); 460 } 461 } 462 else 463 { 464 while(j <= size(L)) 465 { 466 P = modpSpecialAlgDep(ringL, L[j], alg); 467 CO1[index] = P[1]; 468 CO2[index] = bigint(P[2]); 469 index++; 470 j++; 471 } 472 } 473 474 if(printlevel >= 9) 475 { 476 "Realtime for computing list in assPrimes is "+string(rtimer  rt)+ 477 " seconds."; 478 "CPUtime for computing list in assPrimes is "+string(timer  tt)+ 479 " seconds."; 480 } 481 482 // Lift results to basering via farey  483 484 tt = timer; 485 G = chinrem(CO1,CO2); 486 N = CO2[1]; 487 for(j = 2; j <= size(CO2); j++){ N = N*CO2[j]; } 488 F = farey(G,N); 489 if(printlevel >= 10) { "Liftingprocess takes "+string(timer  tt) 490 +" seconds"; } 491 492 if(pTestPoly(F[1], ringL, alg, L)) 493 { 494 F = cleardenom(F[1]); 495 496 e = deg(F[1]); 497 if(e == d) 498 { 499 H = factorize(F[1]); 500 501 s = size(H[1]); 502 for(i = 1; i <= s; i++) 503 { 504 if(H[2][i] != 1) 505 { 506 int_break = 1; 507 } 508 } 509 510 if(int_break == 0) 511 { 512 setring SPR; 513 phi = rHelp,var(nvars(SPR)); 514 H = phi(H); 515 516 if(printlevel >= 9) 517 { 518 "Realtime without test is "+string(rtimer  RT)+" seconds."; 519 "CPUtime without test is "+string(timer  T)+" seconds."; 520 } 521 522 T = timer; 523 RT = rtimer; 524 525 F = phi(F); 526 527 if(n1 > 1) 528 { 529 open(l(1)); 530 write(l(1), quote(quickSubst(eval(F[1]), eval(f), eval(I)))); 531 t_sleep = timer; 532 } 533 else 534 { 535 F1 = quickSubst(F[1],f,I); 536 if(F1 != 0) { int_break = 1; } 537 } 538 539 if(int_break == 0) 540 { 541 for(i = 2; i <= s; i++) 542 { 543 H[1][i] = quickSubst(H[1][i],f,I); 544 Re[i1] = I + ideal(H[1][i]); 545 } 546 547 if(n1 > 1) 548 { 549 t_sleep = timer  t_sleep; 550 if(t_sleep > 5) { t_sleep = 5; } 551 552 while(1) 553 { 554 if(status(l(1), "read", "ready")) 555 { 556 F1 = read(l(1)); 557 close(l(1)); 558 break; 559 } 560 i_sleep = system("sh", "sleep "+string(t_sleep)); 561 } 562 if(F1 != 0) { int_break = 1; } 563 } 564 if(printlevel >= 9) 565 { 566 "Realtime for test is "+string(rtimer  RT)+" seconds."; 567 "CPUtime for test is "+string(timer  T)+" seconds."; 568 } 569 if(int_break == 0) 570 { 571 kill f_for_fork; 572 kill I_for_fork; 573 kill SPR; 574 return(Re); 575 } 576 } 577 } 578 } 579 } 580 581 int_break = 0; 582 setring rHelp; 583 testF = F; 584 CO1 = G; 585 CO2 = N; 586 index = 2; 587 588 j = size(L) + 1; 589 590 setring SPR; 591 L = primeList(I,n3,L); 592 setring rHelp; 593 594 if(n1 > 1) 595 { 596 for(i = 1; i <= n1; i++) 597 { 598 open(l(i)); 599 write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[j+i1]), 600 eval(alg)))); 601 } 602 j = j + n1; 603 k = 0; 604 } 605 } 219 220 /* compute the radical if necessary */ 221 ideal J = I; 222 int isRad; 223 poly f; 224 isRad, f = pTestRad(I, d); 225 while(!isRad) 226 { 227 J = zeroRadical(I, exactness); 228 J = modStd(J, exactness); 229 d = vdim(J); 230 isRad, f = pTestRad(J, d); 231 } 232 I = J; 233 kill J; 234 235 /* call modular() */ 236 if(exactness) 237 { 238 ideal F = modular("Assprimeszerodim::modpSpecialAlgDep", 239 list(I, f, d, a), 240 Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default, 241 pTest_assPrimes, finalTest_assPrimes); 242 } 243 else 244 { 245 ideal F = modular("Assprimeszerodim::modpSpecialAlgDep", 246 list(I, f, d, a), 247 Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default, 248 pTest_assPrimes); 249 } 250 251 /* compute the components */ 252 list result; 253 list H = factorize(F[1]); 254 for(i = size(H[1])1; i > 0; i) 255 { 256 result[i] = I + ideal(quickSubst(H[1][i+1], f, I)); 257 } 258 259 /* return the result */ 260 return(result); 606 261 } 607 262 example … … 620 275 //////////////////////////////////////////////////////////////////////////////// 621 276 622 static proc specialAlgDepEHV(poly p, ideal I) 623 { 624 //=== computes a poly F in Q[T] such that <F>=kernel(Q[T]>basering) 625 //=== mapping T to p 277 /* Computes a poly F in Q[T] such that 278 * <F> = kernel(Q[T] > basering, T > f), 279 * T := last variable in the basering. 280 */ 281 static proc specialAlgDepEHV(ideal I, poly f) 282 { 626 283 def R = basering; 627 execute("ring Rhelp="+charstr(R)+",T,dp;"); 628 setring R; 629 map phi = Rhelp,p; 630 setring Rhelp; 631 ideal F = preimage(R,phi,I); //corresponds to std(I,pT) in dp(n),dp(1) 632 export(F); 633 setring R; 634 list L = Rhelp; 635 return(L); 636 } 637 638 //////////////////////////////////////////////////////////////////////////////// 639 640 static proc specialAlgDepGTZ(poly p, ideal I) 641 { 642 //=== assume I is zerodimensional 643 //=== computes a poly F in Q[T] such that <F>=kernel(Q[T]>basering) 644 //=== mapping T to p 284 execute("ring QT = ("+charstr(R)+"), "+varstr(R, nvars(R))+", dp;"); 285 setring(R); 286 map phi = QT, f; 287 setring QT; 288 ideal F = preimage(R, phi, I); // corresponds to std(I, fT) in dp(n),dp(1) 289 setring(R); 290 ideal F = imap(QT, F); 291 return(F); 292 } 293 294 //////////////////////////////////////////////////////////////////////////////// 295 296 /* Assume I is zerodimensional. 297 * Computes a poly F in Q[T] such that 298 * <F> = kernel(Q[T] > basering, T > f), 299 * T := last variable in the basering. 300 */ 301 static proc specialAlgDepGTZ(ideal I, poly f) 302 { 645 303 def R = basering; 646 execute("ring Rhelp = "+charstr(R)+",T,dp;"); 647 setring R; 648 map phi = Rhelp,p; 649 def Rlp = changeord(list(list("dp",1:(nvars(R)1)),list("dp",1:1))); 650 setring Rlp; 651 poly p = imap(R,p); 304 if(nvars(R) > 1) 305 { 306 def Rlp = changeord(list(list("dp", 1:(nvars(R)1)), list("dp", 1:1))); 307 setring(Rlp); 308 poly f = imap(R, f); 309 ideal I; 310 } 652 311 ideal K = maxideal(1); 653 K[nvars(R)] = 2*var(nvars(R)) p;654 map phi = R, K;655 idealI = phi(I);312 K[nvars(R)] = 2*var(nvars(R))f; 313 map phi = R, K; 314 I = phi(I); 656 315 I = std(I); 657 poly q = subst(I[1],var(nvars(R)),var(1)); 658 setring Rhelp; 659 map psi=Rlp,T; 660 ideal F=psi(q); 661 export(F); 662 setring R; 663 list L=Rhelp; 664 return(L); 665 } 666 667 //////////////////////////////////////////////////////////////////////////////// 668 669 static proc specialAlgDepMonico(poly p, ideal I) 670 { 671 //=== assume I is zerodimensional 672 //=== computes a poly F in Q[T], the characteristic polynomial of the map 673 //=== basering/I > baserng/I defined by the multiplication with p 674 //=== in case I is radical it is the same poly as in specialAlgDepEHV 316 ideal F = I[1]; 317 if(nvars(R) > 1) 318 { 319 setring(R); 320 ideal F = imap(Rlp, F); 321 } 322 return(F); 323 } 324 325 //////////////////////////////////////////////////////////////////////////////// 326 327 /* Assume I is zerodimensional. 328 * Computes a poly F in Q[T], the characteristic polynomial of the map 329 * basering/I > basering/I defined by the multiplication with f, 330 * T := last variable in the basering. 331 * In case I is radical, it is the same polynomial as in specialAlgDepEHV. 332 */ 333 static proc specialAlgDepMonico(ideal I, poly f, int d) 334 { 675 335 def R = basering; 676 execute("ring Rhelp = "+charstr(R)+",T,dp;");677 setring R;678 map phi = Rhelp,p;679 poly q;680 336 int j; 681 matrix m ; 682 poly va = var(1); 337 matrix M[d][d]; 683 338 ideal J = std(I); 684 ideal ba = kbase(J); 685 int d = vdim(J); 686 matrix n[d][d]; 687 for(j = 2; j <= nvars(R); j++) 688 { 689 va = va*var(j); 339 ideal basis = kbase(J); 340 poly vars = var(nvars(R)); 341 for(j = nvars(R)1; j > 0; j) 342 { 343 vars = var(j)*vars; 690 344 } 691 345 for(j = 1; j <= d; j++) 692 346 { 693 q = reduce(p*ba[j],J); 694 m = coeffs(q,ba,va); 695 n[j,1..d] = m[1..d,1]; 696 } 697 setring Rhelp; 698 matrix n = imap(R,n); 699 ideal F = det(nT*freemodule(d)); 700 export(F); 701 setring R; 702 list L = Rhelp; 703 return(L); 347 M[1..d, j] = coeffs(reduce(f*basis[j], J), basis, vars); 348 } 349 execute("ring QT = ("+charstr(R)+"), "+varstr(R, nvars(R))+", dp;"); 350 matrix M = imap(R, M); 351 ideal F = det(Mvar(1)*freemodule(d)); 352 setring(R); 353 ideal F = imap(QT, F); 354 return(F); 704 355 } 705 356 … … 711 362 //=== mapping T to p and test if d=deg(F) 712 363 def R = basering; 713 execute("ring Rhelp ="+charstr(R)+",T,dp;");364 execute("ring Rhelp = ("+charstr(R)+"), T, dp;"); 714 365 setring R; 715 366 map phi = Rhelp,p; … … 718 369 int e=deg(F[1]); 719 370 setring R; 720 return((e==d) );721 } 722 723 //////////////////////////////////////////////////////////////////////////////// 724 725 static proc findGen(ideal J, list #) 726 { 727 //=== try to find a sparse linear form r such that 728 //=== vector space dim(basering/J)=deg(F), 729 //=== F a poly in Q[T] such that <F>=kernel(Q[T]>basering) mapping T to r 730 //=== if not found returns a generic (randomly chosen) r 731 int d = vdim(J); 371 return((e==d), fetch(Rhelp, F)[1]); 372 } 373 374 //////////////////////////////////////////////////////////////////////////////// 375 376 /* Assume d = vector space dim(basering/J). 377 * Tries to find a (sparse) linear form r such that d = deg(F), where 378 * F is a poly in Q[T] such that <F> = kernel(Q[T]>basering) mapping T to r. 379 * If found, returns (1, r, F). If not found, returns (0, 0, 0). 380 */ 381 static proc findGen(ideal J, int d) 382 { 732 383 def R = basering; 733 384 int n = nvars(R); 734 list rl = ringlist(R); 735 if(size(#) > 0) { int p = #[1]; } 736 else { int p = prime(random(1000000000,2134567879)); } 737 rl[1] = p; 738 def @R = ring(rl); 739 setring @R; 740 ideal J = imap(R,J); 385 int okay; 386 poly F; 387 388 /* try trivial transformation */ 741 389 poly r = var(n); 742 int i,k; 743 k = specialTest(d,r,J); 744 if(!k) 390 okay, F = specialTest(d, r, J); 391 if(okay) 392 { 393 return(1, r, F); 394 } 395 396 /* try transformations of the form var(n) + var(i) */ 397 int i; 398 for(i = 1; i < n; i++) 399 { 400 okay, F = specialTest(d, r+var(i), J); 401 if(okay) 402 { 403 return(1, r+var(i), F); 404 } 405 } 406 407 /* try transformations of the form var(n) + \sum var(i) */ 408 if(n > 2) 745 409 { 746 410 for(i = 1; i < n; i++) 747 411 { 748 k = specialTest(d,r+var(i),J);749 if(k){ r = r + var(i); break; }750 }751 }752 if((!k) && (n > 2))753 {754 for(i = 1; i < n; i++)755 {756 412 r = r + var(i); 757 k = specialTest(d,r,J); 758 if(k){ break; } 759 } 760 } 761 setring R; 762 poly r = randomLast(100)[nvars(R)]; 763 if(k){ r = imap(@R,r); } 764 return(r); 765 } 766 767 //////////////////////////////////////////////////////////////////////////////// 768 769 static proc pTestRad(int d, poly p1, ideal I) 770 { 771 //=== computes a poly F in Z/q1[T] such that 772 //=== <F> = kernel(Z/q1[T]>Z/q1[vars(basering)]) 773 //=== mapping T to p1 and test if d=deg(squarefreepart(F)), q1 a prime randomly 774 //=== chosen 775 //=== If not choose randomly another prime q2 and another linear form p2 and 776 //=== computes a poly F in Z/q2[T] such that 777 //=== <F> = kernel(Z/q2[T]>Z/q2[vars(basering)]) 778 //=== mapping T to p2 and test if d=deg(squarefreepart(F)) 779 //=== if the test is positive then I is radical 413 okay, F = specialTest(d, r, J); 414 if(okay) 415 { 416 return(1, r, F); 417 } 418 } 419 } 420 421 /* try random transformations */ 422 int N = 2; // arbitrarily chosen 423 for(i = N; i > 0; i) 424 { 425 r = randomLast(100)[n]; 426 okay, F = specialTest(d, r, J); 427 if(okay) 428 { 429 return(1, r, F); 430 } 431 } 432 433 /* not found */ 434 return(0, 0, 0); 435 } 436 437 //////////////////////////////////////////////////////////////////////////////// 438 439 /* Assume d = vector space dim(basering/I). 440 * Tests if I is radical over F_p, where p is some randomly chosen prime. 441 * If yes, chooses a linear form r such that d = deg(squarefreepart(F)), where 442 * F is a poly in Z/p[T] such that <F> = kernel(Z/p[T]>Z/p[vars(basering)]) 443 * mapping T to r. 444 * Returns (1, r), if I is radical over F_p, and (0, 0) otherwise. 445 */ 446 static proc pTestRad(ideal I, int d) 447 { 448 int N = 2; // Try N random primes. Value of N can be chosen arbitrarily. 780 449 def R = basering; 781 450 list rl = ringlist(R); 782 int q1 = prime(random(100000000,536870912)); 783 rl[1] = q1; 784 ring Shelp1 = q1,T,dp; 785 setring R; 786 def Rhelp1 = ring(rl); 787 setring Rhelp1; 788 poly p1 = imap(R,p1); 789 ideal I = imap(R,I); 790 map phi = Shelp1,p1; 791 setring Shelp1; 792 ideal F = preimage(Rhelp1,phi,I); 793 poly f = gcd(F[1],diff(F[1],var(1))); 794 int e = deg(F[1]/f); 795 setring R; 796 if(e != d) 797 { 798 poly p2 = findGen(I,q1); 799 setring Rhelp1; 800 poly p2 = imap(R,p2); 801 phi = Shelp1,p2; 802 setring Shelp1; 803 F = preimage(Rhelp1,phi,I); 804 f = gcd(F[1],diff(F[1],var(1))); 805 e = deg(F[1]/f); 806 setring R; 807 if(e == d){ return(1); } 808 if(e != d) 809 { 810 int q2 = prime(random(100000000,536870912)); 811 rl[1] = q2; 812 ring Shelp2 = q2,T,dp; 813 setring R; 814 def Rhelp2 = ring(rl); 815 setring Rhelp2; 816 poly p1 = imap(R,p1); 817 ideal I = imap(R,I); 818 map phi = Shelp2,p1; 819 setring Shelp2; 820 ideal F = preimage(Rhelp2,phi,I); 821 poly f = gcd(F[1],diff(F[1],var(1))); 822 e = deg(F[1]/f); 823 setring R; 824 if(e == d){ return(1); } 825 } 826 } 827 return((e==d)); 828 } 829 830 //////////////////////////////////////////////////////////////////////////////// 831 832 static proc zeroRadP(ideal I, int p) 833 { 834 //=== computes F=(F_1,...,F_n) such that <F_i>=IZ/p[x_1,...,x_n] intersected 835 //=== with Z/p[x_i], F_i monic 836 def R0 = basering; 837 list ringL = ringlist(R0); 838 ringL[1] = p; 839 def @r = ring(ringL); 840 setring @r; 841 ideal I = fetch(R0,I); 451 int p; 452 int okay; 453 int i; 454 for(i = N; i > 0; i) 455 { 456 p = prime(random(100000000,536870912)); 457 458 // change to characteristic p 459 if(typeof(rl[1]) == "int") 460 { 461 rl[1] = p; 462 } 463 else 464 { 465 rl[1][1] = p; 466 } 467 def Rp(i) = ring(rl); 468 setring Rp(i); 469 ideal I = imap(R, I); 470 471 // find and test transformation 472 poly r; 473 poly F; 474 okay, r, F = findGen(I, d); 475 if(okay) 476 { 477 poly f = gcd(F, diff(F, var(1))); 478 if(d == deg(F/f)) // F squarefree? 479 { 480 setring(R); 481 return(1, imap(Rp(i), r)); 482 } 483 } 484 setring(R); 485 } 486 return(0, 0); 487 } 488 489 //////////////////////////////////////////////////////////////////////////////// 490 491 /* Computes an ideal F such that ncols(F) = nvars(basering), 492 * < F[i] > = (I intersected with K[var(i)]), and F[i] is monic. 493 */ 494 static proc zeroRadP(ideal I) 495 { 496 intvec opt = option(get); 842 497 option(redSB); 843 498 I = std(I); 844 ideal F = finduni(I); //F[i] generates I intersected with K[var(i)] 845 int i; 846 for(i = 1; i <= size(F); i++){ F[i] = simplify(F[i],1); } 847 setring R0; 848 return(list(fetch(@r,F),p)); 499 ideal F = finduni(I); // F[i] generates I intersected with K[var(i)] 500 F = simplify(F, 1); 501 option(set, opt); 502 return(F); 849 503 } 850 504 … … 883 537 //////////////////////////////////////////////////////////////////////////////// 884 538 885 static proc modpSpecialAlgDep(list ringL, int p, list #) 886 { 887 //=== prepare parallel computing 888 //=== #=1: method of Eisenbud/Hunecke/Vasconcelos 889 //=== #=2: method of Gianni/Trager/Zacharias 890 //=== #=3: method of Monico 891 892 def R0 = basering; 893 894 ringL[1] = p; 895 def @r = ring(ringL); 896 setring @r; 897 poly f = fetch(SPR,f_for_fork); 898 ideal I = fetch(SPR,I_for_fork); 899 if(size(#) > 0) 900 { 901 if(#[1] == 1) { list M = specialAlgDepEHV(f,I); } 902 if(#[1] == 2) { list M = specialAlgDepGTZ(f,I); } 903 if(#[1] == 3) { list M = specialAlgDepMonico(f,I); } 539 /* Simple switch for specialAlgDepEHV, specialAlgDepGTZ, and 540 * specialAlgDepMonico. 541 */ 542 static proc modpSpecialAlgDep(ideal I, poly f, int d, int alg) 543 { 544 ideal F; 545 if(alg == 1) { F = specialAlgDepEHV(I, f); } 546 if(alg == 2) { F = specialAlgDepGTZ(I, f); } 547 if(alg == 3) { F = specialAlgDepMonico(I, f, d); } 548 F = simplify(F, 1); 549 return(F); 550 } 551 552 //////////////////////////////////////////////////////////////////////////////// 553 554 /* The pTest for assPrimes(), to be used in modular(). */ 555 static proc pTest_assPrimes(string command, list args, ideal F, int p) 556 { 557 def br = basering; 558 list lbr = ringlist(br); 559 if(typeof(lbr[1]) == "int") 560 { 561 lbr[1] = p; 904 562 } 905 563 else 906 564 { 907 list M = specialAlgDepEHV(f,I); 908 } 909 def @S = M[1]; 910 911 setring R0; 912 return(list(imap(@S,F),p)); 913 } 914 915 //////////////////////////////////////////////////////////////////////////////// 916 917 static proc pTestPoly(poly testF, list ringL, int alg, list L) 918 { 919 int i,j,p; 920 def R0 = basering; 921 922 while(!i) 923 { 924 i = 1; 925 p = prime(random(1000000000,2134567879)); 926 for(j = 1; j <= size(L); j++) 927 { 928 if(p == L[j]) { i = 0; break; } 929 } 930 } 931 932 ringL[1] = p; 933 def @R = ring(ringL); 934 setring @R; 935 poly f = fetch(SPR,f_for_fork); 936 ideal I = fetch(SPR,I_for_fork); 937 if(alg == 1) { list M = specialAlgDepEHV(f,I); } 938 if(alg == 2) { list M = specialAlgDepGTZ(f,I); } 939 if(alg == 3) { list M = specialAlgDepMonico(f,I); } 940 def @S = M[1]; 941 setring @S; 942 poly testF = fetch(R0,testF); 943 int k = (testF == F); 944 945 setring R0; 565 lbr[1][1] = p; 566 } 567 def rp = ring(lbr); 568 setring(rp); 569 list args_p = fetch(br, args); 570 ideal F = fetch(br, F); 571 execute("ideal Fp = "+command+"(" 572 +Tasks::argsToString("args_p", size(args_p))+");"); 573 int k = (Fp[1] == F[1]); 574 setring br; 946 575 return(k); 576 } 577 578 //////////////////////////////////////////////////////////////////////////////// 579 580 /* The finalTest for assPrimes(), to be used in modular(). */ 581 static proc finalTest_assPrimes(string command, alias list args, ideal F) 582 { 583 F = cleardenom(F[1]); 584 if(deg(F[1]) != args[3]) { return(0); } 585 if(gcd(F[1], diff(F[1], var(nvars(basering)))) != 1) { return(0); }; 586 if(quickSubst(F[1], args[2], args[1]) != 0) { return(0); } 587 return(1); 947 588 } 948 589
Note: See TracChangeset
for help on using the changeset viewer.