Changeset 5ebd453 in git
- Timestamp:
- May 25, 2009, 6:29:04 PM (15 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- 8512090f711dc406701bd340e3fd7e96f79a294c
- Parents:
- 5b2685e36d3336e4237ec537cdbf20e8dd4f57e3
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/normal.lib
r5b2685e r5ebd453 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: normal.lib,v 1.5 3 2009-04-15 11:15:56 seelischExp $";2 version="$Id: normal.lib,v 1.54 2009-05-25 16:29:04 Singular Exp $"; 3 3 category="Commutative Algebra"; 4 4 info=" … … 166 166 if ( find(method,"withDelta") or find(method,"wd") or find(method,"withdelta")) 167 167 { 168 if((decomp == 0) or (decomp == 3)){ 168 if((decomp == 0) or (decomp == 3)) 169 { 169 170 withDelta = 1; 170 } else { 171 } 172 else 173 { 171 174 "WARNING: the delta invariants cannot be computed with an equidimensional"; 172 175 " decomposition."; … … 191 194 // R is the ring where computations will be done 192 195 193 if((useRing == 1) and (isGlobal == 1)){ 196 if((useRing == 1) and (isGlobal == 1)) 197 { 194 198 def globR = basering; 195 } else { 199 } 200 else 201 { 196 202 // We change to dp ordering. 197 203 list rl = ringlist(origR); … … 206 212 //------------------------ trivial checkings ------------------------ 207 213 id = groebner(id); 208 if((size(id) == 0) or (id[1] == 1)){ 214 if((size(id) == 0) or (id[1] == 1)) 215 { 209 216 // The original ring R/I was normal. Nothing to do. 210 217 // We define anyway a new ring, equal to R, to be able to return it. … … 218 225 export normap; 219 226 setring origR; 220 if(withDelta){ 227 if(withDelta) 228 { 221 229 result = list(list(ideal(1), poly(1), ROut), 0); 222 230 } else { … … 228 236 //------------------------ preliminary decomposition----------------------- 229 237 list prim; 230 if(decomp == 2){ 238 if(decomp == 2) 239 { 231 240 dbprint(dbg, "// Computing the equidimensional decomposition..."); 232 241 prim = equidim(id); 233 242 } 234 if((decomp == 0) or (decomp == 1)){ 243 if((decomp == 0) or (decomp == 1)) 244 { 235 245 prim = id; 236 246 } 237 if(decomp == 3){ 247 if(decomp == 3) 248 { 238 249 dbprint(dbg, "// Computing the minimal associated primes..."); 239 250 if( noFac ) … … 253 264 //----------------- back to the original ring if required ------------------ 254 265 // if ring was not global and useRing is on, we go back to the original ring 255 if((useRing == 1) and (isGlobal != 1)){ 266 if((useRing == 1) and (isGlobal != 1)) 267 { 256 268 setring origR; 257 269 def R = basering; 258 270 list prim = fetch(globR, prim); 259 } else { 271 } 272 else 273 { 260 274 def R = basering; 261 if(useRing == 0){ 275 if(useRing == 0) 276 { 262 277 ideal U; 263 278 poly c; … … 283 298 { 284 299 if(dbg>=2){pause();} 285 if(dbg>=1){ 300 if(dbg>=1) 301 { 286 302 "// start computation of component",i; 287 303 " --------------------------------"; 288 304 } 289 if(groebner(prim[i])[1] != 1){ 290 if(dbg>=2){ 305 if(groebner(prim[i])[1] != 1) 306 { 307 if(dbg>=2) 308 { 291 309 "We compute the normalization in the ring"; basering; 292 310 } … … 294 312 norComp = normalM(prim[i], decomp, withDelta); 295 313 printlevel = printlevel - 1; 296 for(j = 1; j <= size(norComp); j++){ 314 for(j = 1; j <= size(norComp); j++) 315 { 297 316 newR = norComp[j][3]; 298 317 newRList = ringlist(newR); 299 318 U = norComp[j][1]; 300 319 c = norComp[j][2]; 301 if(withDelta){ 320 if(withDelta) 321 { 302 322 delt = norComp[j][4]; 303 if((delt >= 0) and (deltI >= 0)){ 323 if((delt >= 0) and (deltI >= 0)) 324 { 304 325 deltI = deltI + delt; 305 } else { 326 } 327 else 328 { 306 329 deltI = -1; 307 330 } 308 331 } 309 332 // -- incorporate result for this component to the list of results --- 310 if(useRing == 0){ 333 if(useRing == 0) 334 { 311 335 // We go back to the original ring. 312 336 setring origR; … … 315 339 newRListO = imap(R, newRList); 316 340 // We change the ordering in the new ring. 317 if(nvars(newR) > nvars(origR)){ 341 if(nvars(newR) > nvars(origR)) 342 { 318 343 newRListO[3]=insert(origOrd, newRListO[3][1]); 319 } else { 344 } 345 else 346 { 320 347 newRListO[3] = origOrd; 321 348 } … … 329 356 totalComps++; 330 357 result[totalComps] = list(U, c, newROrigOrd); 331 if(withDelta){ 358 if(withDelta) 359 { 332 360 result[totalComps] = insert(result[totalComps], delt, 3); 333 361 } 334 362 setring R; 335 } else { 363 } 364 else 365 { 336 366 setring R; 337 367 totalComps++; … … 343 373 344 374 // -------------------------- delta computation ---------------------------- 345 if(withDelta == 1){ 375 if(withDelta == 1) 376 { 346 377 // Intersection multiplicities of list prim, sp=size(prim). 347 378 if ( dbg >= 1 ) … … 371 402 list MG; // Module generators 372 403 intvec DV; // Vector of delta's of each component 373 for(i = 1; i <= size(result); i++){ 404 for(i = 1; i <= size(result); i++) 405 { 374 406 RL[i] = result[i][3]; 375 407 MG[i] = lineUpLast(result[i][1], result[i][2]); 376 if(withDelta){ 408 if(withDelta) 409 { 377 410 DV[i] = result[i][4]; 378 411 } … … 380 413 list resultNew; 381 414 382 if(withDelta){ 415 if(withDelta) 416 { 383 417 resultNew = list(RL, MG, list(DV, deltI)); 384 418 } else { … … 391 425 { 392 426 ""; 393 if(!withDelta){ 427 if(!withDelta) 428 { 394 429 "// 'normal' created a list, say nor, of two elements."; 395 430 } else { … … 410 445 "// 1/ci * Ui, with Ui the i-th ideal of nor[2] and ci the last element"; 411 446 "// Ui[size(Ui)] of Ui."; 412 if(withDelta){ 447 if(withDelta) 448 { 413 449 "// * nor[3] is a list of an intvec of size r, the delta invariants "; 414 450 "// of the r components, and an integer, the total delta invariant "; … … 484 520 poly p = Li[4]; 485 521 int noRed = 0; 486 if(size(Li) > 4){ 487 if(Li[5] == 1){ 522 if(size(Li) > 4) 523 { 524 if(Li[5] == 1) 525 { 488 526 noRed = 1; 489 527 } … … 732 770 733 771 734 if(noRed == 0){ 772 if(noRed == 0) 773 { 735 774 L = substpart(endid,endphi,homo,rw); 736 775 def lastRing=L[1]; … … 3248 3287 proc normalP(ideal id,list #) 3249 3288 "USAGE: normalP(id [,choose]); id a radical ideal, choose a comma separated 3250 list of optional strings: \"withRing\" or \"isPrim\" or \"noFac\".@* 3251 If choose = \"noFac\", factorization is avoided during the computation 3252 of the minimal associated primes; choose = \"isPrim\" disables the 3253 computation of the minimal associated primes (should only be used 3254 if the ideal is known to be prime).@* 3289 list of optional strings \"withRing\", \"isPrim\", \"noFac\", 3290 \"noRed\", where@* 3291 - \"noFac\", factorization is avoided during the computation 3292 of the minimal associated primes.@* 3293 - \"isPrim\", disables the computation of the minimal associated 3294 primes (should only be used if the ideal is known to be prime).@* 3295 - \"withRing\", the ring structure of the normalization is 3296 computed. The number of variables in the new ring is reduced as much 3297 as possible.@* 3298 - \"noRed\", when computing the ring structure no reduction on the 3299 number of variables is done, it creates one new variable for each 3300 of the new generators of the integral closure in the quotient field.@* 3255 3301 ASSUME: The characteristic of the ground field must be positive. The ideal 3256 3302 id is assumed to be radical. However, if choose does not contain … … 3301 3347 " 3302 3348 { 3303 int i,j,y, wring, isPrim, noFac, sr, del, co;;3349 int i,j,y, sr, del, co; 3304 3350 intvec deli; 3305 3351 list resu, Resu, prim, Gens, mstdid; 3306 3352 ideal gens; 3353 3354 // Default options 3355 int wring = 0; // The ring structure is not computed. 3356 int noRed = 0; // No reduction is done in the ring structure 3357 int isPrim = 0; // Ideal is not assumed to be prime 3358 int noFac = 0; // Use facstd when computing the decomposition 3359 3360 3307 3361 y = printlevel-voice+2; 3308 3362 … … 3337 3391 wring=1; 3338 3392 } 3393 if ( find(method,"noRed") or find(method,"nored") ) 3394 { 3395 noRed=1; 3396 } 3339 3397 if ( find(method,"isPrim") or find(method,"isprim") ) 3340 3398 { 3341 3399 isPrim=1; 3342 3400 } 3343 if ( find(method,"noFac") )3401 if ( find(method,"noFac") or find(method,"nofac")) 3344 3402 { 3345 3403 noFac=1; … … 3404 3462 { 3405 3463 //------ 2-nd main subprocedure: compute ring structure ----------- 3406 resu = list(computeRing(II,prim[i])) + resu; 3464 if(noRed == 0){ 3465 resu = list(computeRing(II,prim[i])) + resu; 3466 } else { 3467 resu = list(computeRing(II,prim[i], "noRed")) + resu; 3468 } 3407 3469 } 3408 3470 printlevel = printlevel-1; … … 3642 3704 list g3 = gnirlist[3]; //contains blocks of orderings 3643 3705 int n3 = size(g3); 3706 3644 3707 //----------------- first identify module ordering ------------------ 3645 3708 if ( g3[n3][1]== "c" or g3[n3][1] == "C" ) … … 3656 3719 } 3657 3720 //---- delete variables which were substituted and weights -------- 3658 //gnirlist;"";3659 3721 intvec V; 3660 3722 int n(0); … … 3663 3725 for ( ii=1; ii<=n3-1; ii++ ) 3664 3726 { 3665 V = V,g3[ii][2]; //copy weights for ordering in each block 3666 if ( ii==1 ) //into one intvector 3667 { 3668 V = V[2..size(V)]; 3669 } 3727 // If order is a matrix ordering, it is replaced by dp ordering. 3728 // TODO: replace it only when some of the original 3729 // variables are eliminated. 3730 if(g3[ii][1] == "M"){ 3731 g3[ii][1] = "dp"; 3732 g3[ii][2] = (1..sqroot(size(g3[ii][2])))*0+1; 3733 } 3734 V = V,g3[ii][2]; //copy weights for ordering in each block 3735 if ( ii==1 ) //into one intvector 3736 { 3737 V = V[2..size(V)]; 3738 } 3670 3739 // int n(ii) = size(g3[ii][2]); 3671 3740 int n(ii) = size(V); … … 3674 3743 for ( jj = n(ii-1)+1; jj<=n(ii); jj++) 3675 3744 { 3676 //"ii, jj", ii,jj;3677 3745 if( Le[4][jj] !=0 ) //jj=index of var which was not substituted 3678 3746 { … … 3680 3748 newg2[kk] = g2[jj]; //not substituted var from varlist 3681 3749 V(ii)=V(ii),V[jj]; //weight of not substituted variable 3682 //"V(",ii,")", V(ii);3683 3750 } 3684 3751 } … … 4009 4076 def R = basering; 4010 4077 4011 // ---------------- computation of the singular locus --------------------- 4012 // We compute the radical of the ideal of minors modulo the original ideal. 4013 // This is done only in the first step. 4014 4078 //------------------------ Groebner bases and dimension of I----------------- 4015 4079 if(isGlobal == 1){ 4016 4080 list IM = mstd(I); … … 4025 4089 int d = dim(I); 4026 4090 4091 // ---------------- computation of the singular locus --------------------- 4092 // We compute the radical of the ideal of minors modulo the original ideal. 4093 // This is done only in the first step. 4027 4094 qring Q = I; // We work in the quotient by the groebner base of the ideal I 4028 4095 option("redSB"); … … 4033 4100 4034 4101 dbprint(dbg, "Computing the jacobian ideal..."); 4102 4035 4103 ideal J = minor(jacob(IMin), nvars(basering) - d, I); 4036 4104 J = groebner(J); 4105 4106 //------------------ We check if the singular locus is empty ------------- 4107 if(J[1] == 1){ 4108 // The original ring R/I was normal. Nothing to do. 4109 // We define anyway a new ring, equal to R, to be able to return it. 4110 setring R; 4111 list lR = ringlist(R); 4112 def ROut = ring(lR); 4113 setring ROut; 4114 ideal norid = fetch(R, I); 4115 ideal normap = maxideal(1); 4116 export norid; 4117 export normap; 4118 setring R; 4119 if(withDelta){ 4120 list output = ideal(1), poly(1), ROut, 0; 4121 } else { 4122 list output = ideal(1), poly(1), ROut; 4123 } 4124 return(list(output)); 4125 } 4126 4037 4127 4038 4128 // -------------------- election of the conductor ------------------------- … … 4148 4238 // I = Id1 \cap Id2 4149 4239 printlevel = printlevel + 1; 4240 4150 4241 list nor1 = normalM(Id1, decomp, withDelta)[1]; 4151 4242 list nor2 = normalM(Id2, decomp, withDelta)[1]; … … 4210 4301 ideal oldU = 1; 4211 4302 4212 if(dbg >= 2){ 4303 if(dbg >= 2) 4304 { 4213 4305 "The quotient is"; U; 4214 4306 } … … 4217 4309 // We check if the equality in Grauert - Remmert criterion is satisfied. 4218 4310 isNormal = checkInclusions(D*oldU, U); 4219 if(isNormal == 0){ 4220 if(dbg >= 1){ 4311 if(isNormal == 0) 4312 { 4313 if(dbg >= 1) 4314 { 4221 4315 "In this step, we have the ring 1/c * U, with c =", c; 4222 4316 "and U = "; U; 4223 4317 } 4224 } else { 4318 } 4319 else 4320 { 4225 4321 // The original ring R/I was normal. Nothing to do. 4226 4322 // We define anyway a new ring, equal to R, to be able to return it. … … 4229 4325 def ROut = ring(lR); 4230 4326 setring ROut; 4231 ideal norid = 0;4327 ideal norid = fetch(R, I); 4232 4328 ideal normap = maxideal(1); 4233 4329 export norid; 4234 4330 export normap; 4235 4331 setring R; 4236 if(withDelta){ 4332 if(withDelta) 4333 { 4237 4334 list output = ideal(1), poly(1), ROut, 0; 4238 } else { 4335 } 4336 else 4337 { 4239 4338 list output = ideal(1), poly(1), ROut; 4240 4339 } … … 4243 4342 4244 4343 // ----- computation of the chain of ideals A1 c A2 c ... c An ------------ 4245 while(isNormal == 0){ 4344 while(isNormal == 0) 4345 { 4246 4346 step++; 4247 if(dbg >= 1){ 4347 if(dbg >= 1) 4348 { 4248 4349 ""; 4249 4350 "Step ", step, " begins."; … … 4270 4371 cJMod = getGenerators(cJ, U, c); 4271 4372 4272 if(dbg >= 2){ 4373 if(dbg >= 2) 4374 { 4273 4375 "The test ideal in this step is "; 4274 4376 cJMod; … … 4287 4389 isNormal = checkInclusions(D*oldU, U); 4288 4390 4289 if(isNormal == 1){ 4391 if(isNormal == 1) 4392 { 4290 4393 // We go one step back. In the last step we didnt get antyhing new, 4291 4394 // we just verified that the ring was already normal. … … 4293 4396 dbprint(dbg, ""); 4294 4397 U = oldU; 4295 } else { 4398 } 4399 else 4400 { 4296 4401 // ------------- preparation for next iteration ---------------------- 4297 4402 // We have to go on. 4298 4403 // The new denominator is chosen. 4299 4404 c = D * c; 4300 if(deg(c) > deg(condu)){ 4405 if(deg(c) > deg(condu)) 4406 { 4301 4407 U = changeDenominatorQ(U, c, condu); 4302 4408 c = condu; 4303 4409 } 4304 if(dbg >= 1){ 4410 if(dbg >= 1) 4411 { 4305 4412 "In this step, we have the ring 1/c * U, with c =", c; 4306 4413 "and U = "; … … 4312 4419 4313 4420 // ------------------------- delta computation ---------------------------- 4314 if(withDelta){ 4421 if(withDelta) 4422 { 4315 4423 ideal UD = groebner(U); 4316 4424 delt = vdim(std(modulo(UD, c))); … … 4340 4448 int i; 4341 4449 ideal newU; 4342 for (i = 1; i <= ncols(U); i++){ 4343 if(U[i] != c){ 4344 if(size(newU) == 0){ 4450 for (i = 1; i <= ncols(U); i++) 4451 { 4452 if(U[i] != c) 4453 { 4454 if(size(newU) == 0) 4455 { 4345 4456 newU = U[i]; 4346 } else { 4457 } 4458 else 4459 { 4347 4460 newU = newU, U[i]; 4348 4461 } 4349 4462 } 4350 4463 } 4351 if(size(newU) == 0){ 4464 if(size(newU) == 0) 4465 { 4352 4466 newU = c; 4353 } else { 4467 } 4468 else 4469 { 4354 4470 newU = newU, c; 4355 4471 } … … 4364 4480 int i; 4365 4481 ideal newU = c; 4366 for (i = 1; i <= ncols(U); i++){ 4367 if(U[i] != c){ 4482 for (i = 1; i <= ncols(U); i++) 4483 { 4484 if(U[i] != c) 4485 { 4368 4486 newU = newU, U[i]; 4369 4487 } … … 4374 4492 /////////////////////////////////////////////////////////////////////////////// 4375 4493 4376 static proc getSmallest(ideal J){ 4494 static proc getSmallest(ideal J) 4495 { 4377 4496 // Computes the polynomial of smallest degree of J. 4378 4497 // If there are more than one, it chooses the one with smallest number … … 4382 4501 int d = deg(p); 4383 4502 int di; 4384 for(i = 2; i <= ncols(J); i++){ 4385 if(J[i] != 0){ 4503 for(i = 2; i <= ncols(J); i++) 4504 { 4505 if(J[i] != 0) 4506 { 4386 4507 di = deg(J[i]); 4387 if(di < d){ 4508 if(di < d) 4509 { 4388 4510 p = J[i]; 4389 4511 d = di; 4390 } else { 4391 if(di == d){ 4392 if(size(J[i]) < size(p)){ 4512 } 4513 else 4514 { 4515 if(di == d) 4516 { 4517 if(size(J[i]) < size(p)) 4518 { 4393 4519 p = J[i]; 4394 4520 } … … 4402 4528 /////////////////////////////////////////////////////////////////////////////// 4403 4529 4404 static proc getGenerators(ideal J, ideal U, poly c){ 4530 static proc getGenerators(ideal J, ideal U, poly c) 4531 { 4405 4532 // Computes the generators of J as an ideal in the original ring, 4406 4533 // where J is given by generators in the new ring. … … 4414 4541 4415 4542 if(dbg>1){"Checking for new generators...";} 4416 for(i = 1; i <= ncols(J); i++){ 4417 for(j = 1; j <= ncols(U); j++){ 4543 for(i = 1; i <= ncols(J); i++) 4544 { 4545 for(j = 1; j <= ncols(U); j++) 4546 { 4418 4547 p = lift(c, J[i]*U[j])[1,1]; 4419 4548 p = reduce(p, JGr); 4420 if(p != 0){ 4421 if(dbg>1){ 4549 if(p != 0) 4550 { 4551 if(dbg>1) 4552 { 4422 4553 "New polynoial added:", p; 4423 4554 if(dbg>4){pause();} … … 4434 4565 /////////////////////////////////////////////////////////////////////////////// 4435 4566 4436 static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D){ 4567 static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D) 4568 { 4437 4569 // Internal procedure, used in normalM. 4438 4570 // Computes the test ideal in the new ring. … … 4461 4593 int isGlobal = ord_test(origEre); // Checks if the original ring has 4462 4594 // global ordering. 4463 if(isGlobal != 1){ 4595 if(isGlobal != 1) 4596 { 4464 4597 list rl = ringlist(origEre); 4465 4598 list origOrd = rl[3]; … … 4470 4603 setring ere; 4471 4604 ideal norid = imap(origEre, norid); 4472 } else { 4605 } 4606 else 4607 { 4473 4608 def ere = origEre; 4474 4609 } … … 4488 4623 J = groebner(J); 4489 4624 //J = interred(J); 4490 if(dbg > 1){ 4625 if(dbg > 1) 4626 { 4491 4627 "The radical in the generated ring is"; 4492 4628 J; … … 4517 4653 // ---------------------- step 1 of the mapping --------------------------- 4518 4654 intvec degs; 4519 for(i = 1; i<=ncols(J); i++){ 4655 for(i = 1; i<=ncols(J); i++) 4656 { 4520 4657 degs[i] = degSubring(J[i], @v); 4521 4658 } 4522 if(dbg > 1){ 4659 if(dbg > 1) 4660 { 4523 4661 "The degrees with respect to the new variables are"; 4524 4662 degs; … … 4532 4670 // ---------------------- step 3 of the mapping --------------------------- 4533 4671 ideal z; // The variables of the original ring in order. 4534 for(i = 1; i<=nvars(R); i++){ 4672 for(i = 1; i<=nvars(R); i++) 4673 { 4535 4674 z[i] = var(i); 4536 4675 } 4537 4676 4538 4677 map f = ere, U[2..ncols(U)], z[1..ncols(z)]; // The map to the original ring. 4539 if(dbg > 1){ 4678 if(dbg > 1) 4679 { 4540 4680 "The map is "; 4541 4681 f; … … 4543 4683 } 4544 4684 4545 if(dbg > 1){ 4685 if(dbg > 1) 4686 { 4546 4687 "Computing the map..."; 4547 4688 } 4548 4689 4549 4690 J = f(mapJ); 4550 if(dbg > 1){ 4691 if(dbg > 1) 4692 { 4551 4693 "The ideal J mapped back (before lifting) is"; 4552 4694 J; … … 4558 4700 ideal J = imap(R, J); 4559 4701 poly c = imap(R, c); 4560 for(i = 1; i<=ncols(J); i++){ 4561 if(degs[i]>1){ 4702 for(i = 1; i<=ncols(J); i++) 4703 { 4704 if(degs[i]>1) 4705 { 4562 4706 J[i] = lift(c^(degs[i]-1), J[i])[1,1]; 4563 } else { 4564 if(degs[i]==0){ 4707 } 4708 else 4709 { 4710 if(degs[i]==0) 4711 { 4565 4712 J[i] = c*J[i]; 4566 4713 } … … 4568 4715 } 4569 4716 4570 if(dbg > 1){ 4717 if(dbg > 1) 4718 { 4571 4719 "The ideal J lifted is"; 4572 4720 J; … … 4585 4733 /////////////////////////////////////////////////////////////////////////////// 4586 4734 4587 static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I){ 4735 static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I) 4736 { 4588 4737 // Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the 4589 4738 // ring is 1/c2 * U2. … … 4603 4752 /////////////////////////////////////////////////////////////////////////////// 4604 4753 4605 static proc changeDenominatorQ(ideal U1, poly c1, poly c2){ 4754 static proc changeDenominatorQ(ideal U1, poly c1, poly c2) 4755 { 4606 4756 // Given a ring in the form 1/c1 * U, it computes a new U2 st the ring 4607 4757 // is 1/c2 * U2. … … 4619 4769 /////////////////////////////////////////////////////////////////////////////// 4620 4770 4621 static proc checkInclusions(ideal U1, ideal U2){ 4771 static proc checkInclusions(ideal U1, ideal U2) 4772 { 4622 4773 // Checks if the identity A = Hom(J, J) of Grauert-Remmert criterion is 4623 4774 // satisfied. … … 4636 4787 // The following check should always be satisfied. 4637 4788 // This is only used for debugging. 4638 if(dbg > 1){ 4789 if(dbg > 1) 4790 { 4639 4791 "and the inclusion A c Hom(J, J): (this should always be satisfied)"; 4640 4792 // This interred is used only because a bug in groebner! … … 4645 4797 "Something went wrong... (this inclusion should always be satisfied)"; 4646 4798 ~; 4647 } else { 4799 } 4800 else 4801 { 4648 4802 if(dbg>4){pause();} 4649 4803 } 4650 4804 } 4651 4805 4652 if(size(reduction1[1]) == 0){ 4806 if(size(reduction1[1]) == 0) 4807 { 4653 4808 // We are done! The ring computed in the last step was normal. 4654 4809 return(1); 4655 } else { 4810 } 4811 else 4812 { 4656 4813 return(0); 4657 4814 } … … 4660 4817 /////////////////////////////////////////////////////////////////////////////// 4661 4818 4662 static proc degSubring(poly p, intvec @v){ 4819 static proc degSubring(poly p, intvec @v) 4820 { 4663 4821 // Computes the degree of a polynomial taking only some variables as variables 4664 4822 // and the others as parameters. … … 4668 4826 int d = 0; // The degree 4669 4827 int e; // Degree (auxiliar variable) 4670 for(i = 1; i <= size(p); i++){ 4828 for(i = 1; i <= size(p); i++) 4829 { 4671 4830 e = sum(leadexp(p[i]), @v); 4672 4831 if(e > d){d = e;} … … 4677 4836 /////////////////////////////////////////////////////////////////////////////// 4678 4837 4679 static proc mapBackIdeal(ideal I, poly c, intvec @v){ 4838 static proc mapBackIdeal(ideal I, poly c, intvec @v) 4839 { 4680 4840 // Modifies all polynomials in I so that a map x(i) -> y(i)/c can be 4681 4841 // carried out. … … 4684 4844 4685 4845 int i; // counter 4686 for(i = 1; i <= ncols(I); i++){ 4846 for(i = 1; i <= ncols(I); i++) 4847 { 4687 4848 I[i] = mapBackPoly(I[i], c, @v); 4688 4849 } … … 4692 4853 /////////////////////////////////////////////////////////////////////////////// 4693 4854 4694 static proc mapBackPoly(poly p, poly c, intvec @v){ 4855 static proc mapBackPoly(poly p, poly c, intvec @v) 4856 { 4695 4857 // Multiplies each monomial of p by a power of c so that a map x(i) -> y(i)/c 4696 4858 // can be carried out. … … 4701 4863 int d = degSubring(p, @v); 4702 4864 poly g = 0; 4703 for(i = 1; i <= size(p); i++){ 4865 for(i = 1; i <= size(p); i++) 4866 { 4704 4867 e = sum(leadexp(p[i]), @v); 4705 4868 g = g + p[i] * c^(d-e); … … 5567 5730 return(result); 5568 5731 } 5569 5570 5732 example 5571 5733 { "EXAMPLE:"; … … 5617 5779 int n=size(L); 5618 5780 for (int i=1;i<=n;i++) 5781 { 5782 if (defined(R(i))) 5619 5783 { 5620 if (defined(R(i))) { 5621 string s="Fixed name R("+string(i)+") leads to conflict with existing " 5622 +"object having this name"; 5623 ERROR(s); 5624 } 5625 def R(i)=L[i]; 5626 export R(i); 5627 } 5628 5784 string s="Fixed name R("+string(i)+") leads to conflict with existing " 5785 +"object having this name"; 5786 ERROR(s); 5787 } 5788 def R(i)=L[i]; 5789 export R(i); 5790 } 5629 5791 return(); 5630 5792 } … … 5662 5824 // ("primCl" option must also be set). 5663 5825 5664 proc timeNormal(ideal I, list #){ 5826 proc timeNormal(ideal I, list #) 5827 { 5828 def r = basering; 5829 5665 5830 //--------------------------- define the method --------------------------- 5666 5831 int isPrim, useRing; … … 5700 5865 5701 5866 int tt; 5702 if(norM){ 5867 if(norM) 5868 { 5703 5869 tt = timer; 5704 if(decomp == 0){ 5870 if(decomp == 0) 5871 { 5705 5872 "Running normal(useRing, isPrim)..."; 5706 5873 list a1 = normal(I, "useRing", "isPrim"); 5707 5874 "Time normal(useRing, isPrim): ", timer - tt; 5708 } else { 5875 } 5876 else 5877 { 5709 5878 "Running normal(useRing)..."; 5710 5879 list a1 = normal(I, "useRing"); … … 5713 5882 ""; 5714 5883 } 5715 if(norP){ 5884 if(norP) 5885 { 5716 5886 tt = timer; 5717 if(decomp == 0){ 5887 if(decomp == 0) 5888 { 5718 5889 "Running normalP(isPrim)..."; 5719 5890 list a2 = normalP(I, "isPrim"); 5720 5891 "Time normalP(isPrim): ", timer - tt; 5721 } else { 5892 } 5893 else 5894 { 5722 5895 "Running normalP()..."; 5723 5896 list a2 = normalP(I); … … 5727 5900 } 5728 5901 5729 if(norC){ 5902 if(norC) 5903 { 5730 5904 tt = timer; 5731 if(decomp == 0){ 5905 if(decomp == 0) 5906 { 5732 5907 "Running normalC(isPrim)..."; 5733 5908 list a3 = normalC(I, "isPrim"); 5734 5909 "Time normalC(isPrim): ", timer - tt; 5735 } else { 5910 } 5911 else 5912 { 5736 5913 "Running normalC()..."; 5737 5914 list a3 = normalC(I); … … 5741 5918 } 5742 5919 5743 if(primCl){ 5920 if(primCl) 5921 { 5744 5922 tt = timer; 5745 if(decomp == 0){ 5923 if(decomp == 0) 5924 { 5746 5925 "Running normalC(withGens, isPrim)..."; 5747 5926 list a4 = normalC(I, "isPrim", "withGens"); 5748 5927 "Time normalC(withGens, isPrim): ", timer - tt; 5749 } else { 5928 } 5929 else 5930 { 5750 5931 "Running normalC(withGens)..."; 5751 5932 list a4 = normalC(I, "withGens"); … … 5755 5936 } 5756 5937 5757 if(check111 and norM){ 5938 if(check111 and norM) 5939 { 5758 5940 "Checking output with norTest..."; 5759 5941 "WARNING: this checking only works if the original ideal was prime."; … … 5762 5944 } 5763 5945 5764 if(checkP and norP and norM){ 5946 if(checkP and norP and norM) 5947 { 5765 5948 "Comparing with normalP output..."; 5766 if(size(a2) > 0){ 5949 if(size(a2) > 0) 5950 { 5767 5951 "WARNING: this checking only works if the original ideal was prime."; 5768 5952 U1 = a1[2][1]; … … 5775 5959 ideal W = fetch(r, W); 5776 5960 ch = 0; 5777 if(size(reduce(U2, groebner(W))) == 0){ 5961 if(size(reduce(U2, groebner(W))) == 0) 5962 { 5778 5963 "U2 c U1"; 5779 5964 ch = 1; 5780 5965 } 5781 if(size(reduce(W, groebner(U2))) == 0){ 5966 if(size(reduce(W, groebner(U2))) == 0) 5967 { 5782 5968 "U1 c U2"; 5783 5969 ch = ch + 1; 5784 5970 } 5785 if(ch == 2){ 5971 if(ch == 2) 5972 { 5786 5973 "Output of normalP is equal."; 5787 } else { 5974 } 5975 else 5976 { 5788 5977 "ERROR: Output of normalP is different."; 5789 5978 } 5790 5979 setring r; 5791 5980 kill q; 5792 } else { 5981 } 5982 else 5983 { 5793 5984 "normalP returned no output. Comparison is not possible."; 5794 5985 } … … 5796 5987 } 5797 5988 5798 if(checkPC and norM and primCl){ 5989 if(checkPC and norM and primCl) 5990 { 5799 5991 "Comparing with primeClosure output..."; 5800 if(size(a4) > 0){ 5992 if(size(a4) > 0) 5993 { 5801 5994 "WARNING: this checking only works if the original ideal was prime."; 5802 5995 // primeClosure check … … 5810 6003 ideal W = fetch(r, W); 5811 6004 ch = 0; 5812 if(size(reduce(U2, groebner(W))) == 0){ 6005 if(size(reduce(U2, groebner(W))) == 0) 6006 { 5813 6007 "U2 c U1"; 5814 6008 ch = 1; 5815 6009 } 5816 if(size(reduce(W, groebner(U2))) == 0){ 6010 if(size(reduce(W, groebner(U2))) == 0) 6011 { 5817 6012 "U1 c U2"; 5818 6013 ch = ch + 1; 5819 6014 } 5820 if(ch == 2){ 6015 if(ch == 2) 6016 { 5821 6017 "Output of normalC(withGens) is equal."; 5822 } else { 6018 } 6019 else 6020 { 5823 6021 "ERROR: Output of normalC(withGens) is different."; 5824 6022 } 5825 6023 setring r; 5826 6024 kill q; 5827 } else { 6025 } 6026 else 6027 { 5828 6028 "normalC(withGens) returned no output. Comparison is not possible."; 5829 6029 } 5830 6030 ""; 5831 6031 } 6032 } 6033 6034 /////////////////////////////////////////////////////////////////////////// 6035 static proc sqroot(int n); 6036 { 6037 int s = 1; 6038 while(s*s < n) 6039 { 6040 s++; 6041 } 6042 return(s); 5832 6043 } 5833 6044
Note: See TracChangeset
for help on using the changeset viewer.