Changeset d88470 in git
- Timestamp:
- Nov 13, 2006, 1:57:33 PM (17 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
- Children:
- 1aecaec05d17a34b8ff6a44d71e1684cc74281e4
- Parents:
- 11a4417493a03a79effa983f174d085bb0c77f5e
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/primdec.lib
r11a441 rd88470 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: primdec.lib,v 1.12 7 2006-11-08 13:45:02Singular Exp $";2 version="$Id: primdec.lib,v 1.128 2006-11-13 12:57:33 Singular Exp $"; 3 3 category="Commutative Algebra"; 4 4 info=" … … 2324 2324 " 2325 2325 { 2326 if( ord_test(basering)!=1)2326 if(attrib(basering,"global")!=1) 2327 2327 { 2328 2328 ERROR( … … 2430 2430 " 2431 2431 { 2432 if( ord_test(basering)!=1)2432 if(attrib(basering,"global")!=1) 2433 2433 { 2434 2434 ERROR( … … 3799 3799 " 3800 3800 { 3801 if( ord_test(basering)!=1)3801 if(attrib(basering,"global")!=1) 3802 3802 { 3803 3803 ERROR( … … 5080 5080 " 5081 5081 { 5082 if( ord_test(basering)!=1)5082 if(attrib(basering,"global")!=1) 5083 5083 { 5084 5084 ERROR( … … 5132 5132 } 5133 5133 5134 if( ord_test(basering)!=1)5134 if(attrib(basering,"global")!=1) 5135 5135 { 5136 5136 ERROR( … … 5231 5231 " 5232 5232 { 5233 if( ord_test(basering)!=1)5233 if(attrib(basering,"global")!=1) 5234 5234 { 5235 5235 ERROR( … … 5308 5308 } 5309 5309 5310 if( ord_test(basering)!=1)5310 if(attrib(basering,"global")!=1) 5311 5311 { 5312 5312 ERROR( … … 5345 5345 " 5346 5346 { 5347 if( ord_test(basering)!=1)5347 if(attrib(basering,"global")!=1) 5348 5348 { 5349 5349 ERROR( … … 5374 5374 " 5375 5375 { 5376 if(ord_test(basering)!=1)5377 5378 5379 5380 5381 5382 5376 if(attrib(basering,"global")!=1) 5377 { 5378 ERROR( 5379 "// Not implemented for this ordering, please change to global ordering." 5380 ); 5381 } 5382 return(radical(i, 1)); 5383 5383 } 5384 5384 example … … 5408 5408 " 5409 5409 { 5410 dbprint(printlevel - voice, "Radical, version 2006.05.08"); 5411 if(ord_test(basering)!=1) 5412 { 5413 ERROR( 5414 "// Not implemented for this ordering, please change to global ordering." 5415 ); 5416 } 5417 if(size(i) == 0){return(ideal(0));} 5418 int j; 5419 def P0 = basering; 5420 list Pl=ringlist(P0); 5421 intvec dp_w; 5422 for(j=nvars(P0);j>0;j--) {dp_w[j]=1;} 5423 Pl[3]=list(list("dp",dp_w),list("C",0)); 5424 def @P=ring(Pl); 5425 setring @P; 5426 ideal i=imap(P0,i); 5427 5428 int il; 5429 string algorithm; 5430 int useFac; 5431 5432 // Set input parameters 5433 algorithm = "SL"; // Default: SL algorithm 5434 il = 0; // Default: Full radical (not only equiRadical) 5435 if (homog(i) == 1) { // Default: facStd is used, except if the ideal is homogeneous. 5436 useFac = 0; 5437 } else { 5438 useFac = 1; 5439 } 5440 if(size(#) > 0){ 5441 int valid; 5442 for(j = 1; j <= size(#); j++){ 5443 valid = 0; 5444 if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) 5445 { 5446 il = #[j]; // If il == 1, equiRadical is computed 5447 valid = 1; 5448 } 5449 if(typeof(#[j]) == "string"){ 5450 if(#[j] == "KL") { 5451 algorithm = "KL"; 5452 valid = 1;} 5453 if(#[j] == "SL") { 5454 algorithm = "SL"; 5455 valid = 1;} 5456 if(#[j] == "noFacstd") { 5457 useFac = 0; 5458 valid = 1;} 5459 if(#[j] == "facstd") { 5460 useFac = 1; 5461 valid = 1;} 5462 if(#[j] == "equiRad") { 5463 il = 1; 5464 valid = 1;} 5465 if(#[j] == "fullRad") { 5466 il = 0; 5467 valid = 1;} 5468 } 5469 if(valid == 0) 5470 { 5471 dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); 5472 } 5473 } 5474 } 5475 5476 ideal rad = 1; 5477 intvec op = option(get); 5478 list qr = simplifyIdeal(i); 5479 map phi = @P, qr[2]; 5480 5481 option(redSB); 5482 i = groebner(qr[1]); 5483 option(set, op); 5484 int di = dim(i); 5485 5486 if(di == 0) 5487 { 5488 i = zeroRad(i, qr[1]); 5489 i=interred(phi(i)); 5490 setring(P0); 5491 i=imap(@P,i); 5492 return(i); 5493 } 5494 5495 option(redSB); 5496 list pr; 5497 if(useFac == 1) 5498 { 5499 pr = facstd(i); 5500 } else { 5501 pr = i; 5502 } 5503 option(set, op); 5504 int s = size(pr); 5505 if(useFac == 1) { 5506 dbprint(printlevel - voice, "Number of components returned by facstd: ", s); 5507 } 5508 for(j = 1; j <= s; j++) 5509 { 5510 attrib(pr[s + 1 - j], "isSB", 1); 5511 if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il)) 5512 { 5513 // SL Debug messages 5514 dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]); 5515 dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j])); 5516 5517 if(algorithm == "KL") { 5518 rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il)); 5519 } 5520 if(algorithm == "SL") { 5521 rad = intersect(rad, radicalSL(pr[s + 1 - j], il)); 5522 } 5523 } else 5524 { 5525 // SL Debug 5526 dbprint(printlevel-voice, "The radical of this component is not needed."); 5527 dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", 5528 size(reduce(rad, pr[s + 1 - j], 1))); 5529 dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j])); 5530 dbprint(printlevel-voice, "il", il); 5531 } 5532 } 5533 rad=interred(phi(rad)); 5534 setring(P0); 5535 i=imap(@P,rad); 5536 return(i); 5410 dbprint(printlevel - voice, "Radical, version 2006.05.08"); 5411 if(attrib(basering,"global")!=1) 5412 { 5413 ERROR( 5414 "// Not implemented for this ordering, please change to global ordering." 5415 ); 5416 } 5417 if(size(i) == 0){return(ideal(0));} 5418 int j; 5419 def P0 = basering; 5420 list Pl=ringlist(P0); 5421 intvec dp_w; 5422 for(j=nvars(P0);j>0;j--) {dp_w[j]=1;} 5423 Pl[3]=list(list("dp",dp_w),list("C",0)); 5424 def @P=ring(Pl); 5425 setring @P; 5426 ideal i=imap(P0,i); 5427 5428 int il; 5429 string algorithm; 5430 int useFac; 5431 5432 // Set input parameters 5433 algorithm = "SL"; // Default: SL algorithm 5434 il = 0; // Default: Full radical (not only equiRadical) 5435 if (homog(i) == 1) 5436 { // Default: facStd is used, except if the ideal is homogeneous. 5437 useFac = 0; 5438 } else { 5439 useFac = 1; 5440 } 5441 if(size(#) > 0){ 5442 int valid; 5443 for(j = 1; j <= size(#); j++) 5444 { 5445 valid = 0; 5446 if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) 5447 { 5448 il = #[j]; // If il == 1, equiRadical is computed 5449 valid = 1; 5450 } 5451 if(typeof(#[j]) == "string"){ 5452 if(#[j] == "KL") { 5453 algorithm = "KL"; 5454 valid = 1; 5455 } 5456 if(#[j] == "SL") { 5457 algorithm = "SL"; 5458 valid = 1; 5459 } 5460 if(#[j] == "noFacstd") { 5461 useFac = 0; 5462 valid = 1;} 5463 if(#[j] == "facstd") { 5464 useFac = 1; 5465 valid = 1;} 5466 if(#[j] == "equiRad") { 5467 il = 1; 5468 valid = 1;} 5469 if(#[j] == "fullRad") { 5470 il = 0; 5471 valid = 1;} 5472 } 5473 if(valid == 0) 5474 { 5475 dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); 5476 } 5477 } 5478 } 5479 5480 ideal rad = 1; 5481 intvec op = option(get); 5482 list qr = simplifyIdeal(i); 5483 map phi = @P, qr[2]; 5484 5485 option(redSB); 5486 i = groebner(qr[1]); 5487 option(set, op); 5488 int di = dim(i); 5489 5490 if(di == 0) 5491 { 5492 i = zeroRad(i, qr[1]); 5493 i=interred(phi(i)); 5494 setring(P0); 5495 i=imap(@P,i); 5496 return(i); 5497 } 5498 5499 option(redSB); 5500 list pr; 5501 if(useFac == 1) 5502 { 5503 pr = facstd(i); 5504 } else { 5505 pr = i; 5506 } 5507 option(set, op); 5508 int s = size(pr); 5509 if(useFac == 1) { 5510 dbprint(printlevel - voice, "Number of components returned by facstd: ", s); 5511 } 5512 for(j = 1; j <= s; j++) 5513 { 5514 attrib(pr[s + 1 - j], "isSB", 1); 5515 if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il)) 5516 { 5517 // SL Debug messages 5518 dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]); 5519 dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j])); 5520 5521 if(algorithm == "KL") 5522 { 5523 rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il)); 5524 } 5525 if(algorithm == "SL") { 5526 rad = intersect(rad, radicalSL(pr[s + 1 - j], il)); 5527 } 5528 } 5529 else 5530 { 5531 // SL Debug 5532 dbprint(printlevel-voice, "The radical of this component is not needed."); 5533 dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", 5534 size(reduce(rad, pr[s + 1 - j], 1))); 5535 dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j])); 5536 dbprint(printlevel-voice, "il", il); 5537 } 5538 } 5539 rad=interred(phi(rad)); 5540 setring(P0); 5541 i=imap(@P,rad); 5542 return(i); 5537 5543 } 5538 5544 example … … 6121 6127 " 6122 6128 { 6123 if( ord_test(basering)!=1)6129 if(attrib(basering,"global")!=1) 6124 6130 { 6125 6131 ERROR( … … 6169 6175 " 6170 6176 { 6171 if( ord_test(basering)!=1)6177 if(attrib(basering,"global")!=1) 6172 6178 { 6173 6179 ERROR( … … 6237 6243 " 6238 6244 { 6239 if(ord_test(basering)!=1)6240 6241 6242 6243 6244 6245 6246 6247 6248 6249 6250 6251 6252 6253 6254 6245 if(attrib(basering,"global")!=1) 6246 { 6247 ERROR( 6248 "// Not implemented for this ordering, please change to global ordering." 6249 ); 6250 } 6251 def R=basering; 6252 poly q; 6253 int j,time; 6254 matrix m; 6255 list re; 6256 poly va=var(1); 6257 ideal J=groebner(I); 6258 ideal ba=kbase(J); 6259 int d=vdim(J); 6260 dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d)); 6255 6261 //------ compute matrix of multiplication on R/I with generic element p ----- 6256 6257 6258 6259 6260 6261 6262 6263 6264 6265 6266 6267 6268 6269 6270 6271 6272 6262 int e=nvars(basering); 6263 poly p=randomLast(100)[e]+random(-50,50); //the generic element 6264 matrix n[d][d]; 6265 time = timer; 6266 for(j=2;j<=e;j++) 6267 { 6268 va=va*var(j); 6269 } 6270 for(j=1;j<=d;j++) 6271 { 6272 q=reduce(p*ba[j],J); 6273 m=coeffs(q,ba,va); 6274 n[j,1..d]=m[1..d,1]; 6275 } 6276 dbprint(printlevel-voice+2, 6277 "// time for computing multiplication matrix (with generic element) : "+ 6278 string(timer-time)); 6273 6279 //---------------- compute characteristic polynomial of matrix -------------- 6274 6275 6276 6277 6278 6279 6280 execute("ring P1=("+charstr(R)+"),T,dp;"); 6281 matrix n=imap(R,n); 6282 time = timer; 6283 poly charpol=det(n-T*freemodule(d)); 6284 dbprint(printlevel-voice+2,"// time for computing char poly: "+ 6285 string(timer-time)); 6280 6286 //------------------- factorize characteristic polynomial ------------------- 6281 6287 //check first if constant term of charpoly is != 0 (which is true for 6282 6288 //sufficiently generic element) 6283 6284 6285 6286 6287 6288 6289 6290 6289 if(charpol[size(charpol)]!=0) 6290 { 6291 time = timer; 6292 list fac=factor(charpol); 6293 testFactor(fac,charpol); 6294 dbprint(printlevel-voice+2,"// time for factorizing char poly: "+ 6295 string(timer-time)); 6296 int f=size(fac[1]); 6291 6297 //--------------------------- the irreducible case -------------------------- 6292 6293 6294 6295 6296 6297 6298 if(f==1) 6299 { 6300 setring R; 6301 re=I; 6302 return(re); 6303 } 6298 6304 //---------------------------- the reducible case --------------------------- 6299 6305 //if f_i are the irreducible factors of charpoly, mult=ri, then <I,g_i^ri> … … 6302 6308 //Hence it is better to simultaneously reduce with I. For this we need a new 6303 6309 //ring. 6304 6305 6306 6307 6308 6309 6310 6311 6312 6313 6314 6315 6316 6317 6318 6319 6320 6321 6322 6323 6310 execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);"); 6311 list rfac=imap(P1,fac); 6312 intvec ov=option(get);; 6313 option(redSB); 6314 list re1; 6315 ideal new = T-imap(R,p),imap(R,J); 6316 attrib(new, "isSB",1); //we know that new is a standard basis 6317 for(j=1;j<=f;j++) 6318 { 6319 re1[j]=reduce(rfac[1][j]^rfac[2][j],new); 6320 } 6321 setring R; 6322 re = imap(P,re1); 6323 for(j=1;j<=f;j++) 6324 { 6325 J=I,re[j]; 6326 re[j]=interred(J); 6327 } 6328 option(set,ov); 6329 return(re); 6324 6330 } 6325 6331 else 6326 6332 //------------------- choice of generic element failed ------------------- 6327 6333 { 6328 6329 6330 6334 dbprint(printlevel-voice+2,"// try new generic element!"); 6335 setring R; 6336 return(zerodec(I)); 6331 6337 } 6332 6338 }
Note: See TracChangeset
for help on using the changeset viewer.