1 | //////////////////////////////////////////////////////////////////// |
---|
2 | version="version ehv.lib 4.1.2.0 Feb_2019 "; // $Id$ |
---|
3 | category="Commutative Algebra"; |
---|
4 | |
---|
5 | info=" |
---|
6 | LIBRARY: EHV.lib PROCEDURES FOR PRIMARY DECOMPOSITION OF IDEALS |
---|
7 | AUTHORS: Kai Dehmann, dehmann@mathematik.uni-kl.de; |
---|
8 | |
---|
9 | OVERVIEW: |
---|
10 | Algorithms for primary decomposition and radical-computation |
---|
11 | based on the ideas of Eisenbud, Huneke, and Vasconcelos. |
---|
12 | |
---|
13 | PROCEDURES: |
---|
14 | equiMaxEHV(I); equidimensional part of I |
---|
15 | removeComponent(I,e); intersection of the primary components |
---|
16 | of I of dimension >= e |
---|
17 | AssOfDim(I,e); an ideal such that the associated primes |
---|
18 | are exactly the associated primes of I |
---|
19 | having dimension e |
---|
20 | equiRadEHV(I [,Strategy]); equidimensional radical of I |
---|
21 | radEHV(I [,Strategy]); radical of I |
---|
22 | IntAssOfDim1(I,e); intersection of the associated primes of I |
---|
23 | having dimension e |
---|
24 | IntAssOfDim2(I,e); another way of computing the intersection |
---|
25 | of the associated primes of I |
---|
26 | having dimension e |
---|
27 | decompEHV(I); decomposition of a zero-dimensional |
---|
28 | radical ideal I |
---|
29 | AssEHV(I [,Strategy]); associated primes of I |
---|
30 | minAssEHV(I [,Strategy]); minimal associated primes of I |
---|
31 | localize(I,P,l); the contraction of the ideal generated by I |
---|
32 | in the localization w.r.t P |
---|
33 | componentEHV(I,P,L [,Strategy]); a P-primary component for I |
---|
34 | primdecEHV(I [,Strategy]); a minimal primary decomposition of I |
---|
35 | primDecsAreEquivalent(L, K); procedure for comparing the output of |
---|
36 | primary decomposition algorithms (checks |
---|
37 | if the computed associated primes coincide) |
---|
38 | "; |
---|
39 | |
---|
40 | LIB "ring.lib"; |
---|
41 | LIB "general.lib"; |
---|
42 | LIB "elim.lib"; |
---|
43 | LIB "poly.lib"; |
---|
44 | LIB "random.lib"; |
---|
45 | LIB "inout.lib"; |
---|
46 | LIB "matrix.lib"; |
---|
47 | LIB "algebra.lib"; |
---|
48 | LIB "normal.lib"; |
---|
49 | |
---|
50 | |
---|
51 | ///////////////////////////////////////////////////////////////////// |
---|
52 | // // |
---|
53 | // G E N E R A L A L G O R I T H M S // |
---|
54 | // // |
---|
55 | ///////////////////////////////////////////////////////////////////// |
---|
56 | |
---|
57 | ///////////////////////////////////////////////////////////////////// |
---|
58 | static proc AnnExtEHV(int n,list re) |
---|
59 | "USAGE: AnnExtEHV(n,re); n integer, re resolution |
---|
60 | RETURN: ideal, the annihilator of Ext^n(R/I,R) with given |
---|
61 | resolution re of I" |
---|
62 | { |
---|
63 | if(printlevel > 2){"Entering AnnExtEHV.";} |
---|
64 | |
---|
65 | if(n < 0) |
---|
66 | { |
---|
67 | ideal ann = ideal(1); |
---|
68 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
69 | return(ann); |
---|
70 | } |
---|
71 | int l = size(re); |
---|
72 | |
---|
73 | if(n < l) |
---|
74 | { |
---|
75 | matrix f = transpose(re[n+1]); |
---|
76 | if(n == 0) |
---|
77 | { |
---|
78 | matrix g = 0*gen(ncols(f)); |
---|
79 | } |
---|
80 | else |
---|
81 | { |
---|
82 | matrix g = transpose(re[n]); |
---|
83 | } |
---|
84 | module k = syz(f); |
---|
85 | ideal ann = quotient1(g,k); |
---|
86 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
87 | return(ann); |
---|
88 | } |
---|
89 | |
---|
90 | if(n == l) |
---|
91 | { |
---|
92 | ideal ann = Ann(transpose(re[n])); |
---|
93 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
94 | return(ann); |
---|
95 | } |
---|
96 | |
---|
97 | ideal ann = ideal(1); |
---|
98 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
99 | return(ann); |
---|
100 | } |
---|
101 | |
---|
102 | |
---|
103 | ///////////////////////////////////////////////////////////////////// |
---|
104 | //static |
---|
105 | proc isSubset(ideal I,ideal J) |
---|
106 | "USAGE: isSubset(I,J); I, J ideals |
---|
107 | RETURN: integer, 1 if I is a subset of J and 0 otherwise |
---|
108 | " |
---|
109 | { |
---|
110 | if ( attrib(J,"isSB") == 0) { J = groebner(J); } |
---|
111 | int s = ncols(I); |
---|
112 | for(int i=1; i<=s; i++) |
---|
113 | { |
---|
114 | if(reduce(I[i],J)!=0) |
---|
115 | { |
---|
116 | return(0); |
---|
117 | } |
---|
118 | } |
---|
119 | return(1); |
---|
120 | } |
---|
121 | |
---|
122 | |
---|
123 | ///////////////////////////////////////////////////////////////////// |
---|
124 | // // |
---|
125 | // T H E E Q U I D I M E N S I O N A L P A R T // |
---|
126 | // // |
---|
127 | ///////////////////////////////////////////////////////////////////// |
---|
128 | |
---|
129 | ///////////////////////////////////////////////////////////////////// |
---|
130 | proc equiMaxEHV(ideal I) |
---|
131 | "USAGE: equiMaxEHV(I); I ideal |
---|
132 | RETURN: ideal, the equidimensional part of I. |
---|
133 | NOTE: Uses algorithm of Eisenbud, Huneke, and Vasconcelos. |
---|
134 | EXAMPLE: example equiMaxEHV; shows an example |
---|
135 | " |
---|
136 | { |
---|
137 | if(printlevel > 2){"Entering equiMaxEHV.";} |
---|
138 | if(attrib(basering,"global")==0) |
---|
139 | { |
---|
140 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
141 | } |
---|
142 | ideal J = groebner(I); |
---|
143 | int cod = nvars(basering)-dim(J); |
---|
144 | |
---|
145 | //If I is the entire ring... |
---|
146 | if(cod > nvars(basering)) |
---|
147 | { |
---|
148 | //...then return the ideal generated by 1. |
---|
149 | return(ideal(1)); |
---|
150 | } |
---|
151 | |
---|
152 | //Compute a resolution of I. |
---|
153 | if(printlevel > 2){"Computing resolution.";} |
---|
154 | if(homog(I)==1) |
---|
155 | { |
---|
156 | list re = sres(J,cod+1); |
---|
157 | re = minres(re); |
---|
158 | } |
---|
159 | else |
---|
160 | { |
---|
161 | list re = mres(I,cod+1); |
---|
162 | } |
---|
163 | if(printlevel > 2){"Finished computing resolution.";} |
---|
164 | |
---|
165 | //Compute the annihilator of the cod-th EXT-module. |
---|
166 | ideal ann = AnnExtEHV(cod,re); |
---|
167 | attrib(ann,"isEquidimensional",1); |
---|
168 | if(printlevel > 2){"Leaving equiMaxEHV.";} |
---|
169 | return(ann); |
---|
170 | } |
---|
171 | |
---|
172 | example |
---|
173 | { |
---|
174 | "EXAMPLE:"; |
---|
175 | echo = 2; |
---|
176 | ring r = 0,(x,y,z),dp; |
---|
177 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
178 | equiMaxEHV(I); |
---|
179 | } |
---|
180 | |
---|
181 | |
---|
182 | ///////////////////////////////////////////////////////////////////// |
---|
183 | proc removeComponent(ideal I, int e) |
---|
184 | "USAGE: removeComponent(I,e); I ideal, e integer |
---|
185 | RETURN: ideal, the intersection of the primary components |
---|
186 | of I of dimension >= e |
---|
187 | EXAMPLE: example removeComponent; shows an example" |
---|
188 | { |
---|
189 | if(attrib(basering,"global")==0) |
---|
190 | { |
---|
191 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
192 | } |
---|
193 | |
---|
194 | ideal J = groebner(I); |
---|
195 | |
---|
196 | //Compute a resolution of I |
---|
197 | if(homog(I)==1) |
---|
198 | { |
---|
199 | list re = sres(J,0); |
---|
200 | re = minres(re); |
---|
201 | } |
---|
202 | else |
---|
203 | { |
---|
204 | list re = mres(I,0); |
---|
205 | } |
---|
206 | |
---|
207 | int f = nvars(basering); |
---|
208 | int cod; |
---|
209 | ideal ann; |
---|
210 | int g = nvars(basering) - e; |
---|
211 | while(f > g) |
---|
212 | { |
---|
213 | ann = AnnExtEHV(f,re); |
---|
214 | cod = nvars(basering) - dim(groebner(ann)); |
---|
215 | if( cod == f ) |
---|
216 | { |
---|
217 | I = quotient(I,ann); |
---|
218 | } |
---|
219 | f = f-1; |
---|
220 | } |
---|
221 | return(I); |
---|
222 | } |
---|
223 | |
---|
224 | example |
---|
225 | { |
---|
226 | "EXAMPLE:"; |
---|
227 | echo = 2; |
---|
228 | ring r = 0,(x,y,z),dp; |
---|
229 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
230 | removeComponent(I,1); |
---|
231 | } |
---|
232 | |
---|
233 | |
---|
234 | ///////////////////////////////////////////////////////////////////// |
---|
235 | proc AssOfDim(ideal I, int e) |
---|
236 | "USAGE: AssOfDim(I,e); I ideal, e integer |
---|
237 | RETURN: ideal, such that the associated primes are exactly |
---|
238 | the associated primes of I having dimension e |
---|
239 | EXAMPLE: example AssOfDim; shows an example" |
---|
240 | { |
---|
241 | if(attrib(basering,"global")==0) |
---|
242 | { |
---|
243 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
244 | } |
---|
245 | int g = nvars(basering) - e; |
---|
246 | |
---|
247 | //Compute a resolution of I. |
---|
248 | ideal J = std(I); |
---|
249 | if(homog(I)==1) |
---|
250 | { |
---|
251 | list re = sres(J,g+1); |
---|
252 | re = minres(re); |
---|
253 | } |
---|
254 | else |
---|
255 | { |
---|
256 | list re = mres(I,g+1); |
---|
257 | } |
---|
258 | |
---|
259 | ideal ann = AnnExtEHV(g,re); |
---|
260 | int cod = nvars(basering) - dim(std(ann)); |
---|
261 | |
---|
262 | //If the codimension of I_g:=Ann(Ext^g(R/I,R)) equals g... |
---|
263 | if(cod == g) |
---|
264 | { |
---|
265 | //...then return the equidimensional part of I_g... |
---|
266 | ann = equiMaxEHV(ann); |
---|
267 | attrib(ann,"isEquidimensional",1); |
---|
268 | return(ann); |
---|
269 | } |
---|
270 | //...otherwise... |
---|
271 | else |
---|
272 | { |
---|
273 | //...I has no associated primes of dimension e. |
---|
274 | return(ideal(1)); |
---|
275 | } |
---|
276 | } |
---|
277 | |
---|
278 | example |
---|
279 | { |
---|
280 | "EXAMPLE:"; |
---|
281 | echo = 2; |
---|
282 | ring r = 0,(x,y,z),dp; |
---|
283 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
284 | AssOfDim(I,1); |
---|
285 | } |
---|
286 | |
---|
287 | ///////////////////////////////////////////////////////////////////// |
---|
288 | // // |
---|
289 | // T H E R A D I C A L // |
---|
290 | // // |
---|
291 | ///////////////////////////////////////////////////////////////////// |
---|
292 | |
---|
293 | ///////////////////////////////////////////////////////////////////// |
---|
294 | static proc aJacob(ideal I, int a) |
---|
295 | "USAGE: aJacob(I,a); I ideal, a integer |
---|
296 | RETURN: ideal, the ath-Jacobian ideal of I" |
---|
297 | { |
---|
298 | matrix M = jacob(I); |
---|
299 | int n = nvars(basering); |
---|
300 | if(n-a <= 0) |
---|
301 | { |
---|
302 | return(ideal(1)); |
---|
303 | } |
---|
304 | if(n-a > nrows(M) or n-a > ncols(M)) |
---|
305 | { |
---|
306 | return(ideal(0)); |
---|
307 | } |
---|
308 | ideal J = minor(M,n-a); |
---|
309 | return(J); |
---|
310 | } |
---|
311 | |
---|
312 | |
---|
313 | ///////////////////////////////////////////////////////////////////// |
---|
314 | proc equiRadEHV(ideal I, list #) |
---|
315 | "USAGE: equiRadEHV(I [,Strategy]); I ideal, Strategy list |
---|
316 | RETURN: ideal, the equidimensional radical of I, |
---|
317 | i.e. the intersection of the minimal associated primes of I |
---|
318 | having the same dimension as I |
---|
319 | NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos, |
---|
320 | Works only in characteristic 0 or p large. |
---|
321 | The (optional) second argument determines the strategy used: |
---|
322 | Strategy[1] > strategy for the equidimensional part |
---|
323 | = 0 : uses equiMaxEHV |
---|
324 | = 1 : uses equidimMax |
---|
325 | Strategy[2] > strategy for the radical |
---|
326 | = 0 : combination of strategy 1 and 2 |
---|
327 | = 1 : computation of the radical just with the |
---|
328 | help of regular sequences |
---|
329 | = 2 : does not try to find a regular sequence |
---|
330 | If no second argument is given then Strategy=(0,0,0) is used. |
---|
331 | EXAMPLE: example equiRadEHV; shows an example" |
---|
332 | { |
---|
333 | if(printlevel > 2){"Entering equiRadEHV.";} |
---|
334 | if(attrib(basering,"global")==0) |
---|
335 | { |
---|
336 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
337 | } |
---|
338 | if((char(basering)<100)&&(char(basering)!=0)) |
---|
339 | { |
---|
340 | "WARNING: The characteristic is too small, the result may be wrong"; |
---|
341 | } |
---|
342 | |
---|
343 | //Define the Strategy to be used. |
---|
344 | if(size(#) > 0) |
---|
345 | { |
---|
346 | if(#[1]!=1) |
---|
347 | { |
---|
348 | int equStr = 0; |
---|
349 | } |
---|
350 | else |
---|
351 | { |
---|
352 | int equStr = 1; |
---|
353 | } |
---|
354 | if(size(#) > 1) |
---|
355 | { |
---|
356 | if(#[2]!=1 and #[2]!=2) |
---|
357 | { |
---|
358 | int strategy = 0; |
---|
359 | } |
---|
360 | else |
---|
361 | { |
---|
362 | int strategy = #[2]; |
---|
363 | } |
---|
364 | } |
---|
365 | else |
---|
366 | { |
---|
367 | int strategy = 0; |
---|
368 | } |
---|
369 | } |
---|
370 | else |
---|
371 | { |
---|
372 | int equStr = 0; |
---|
373 | int strategy = 0; |
---|
374 | } |
---|
375 | |
---|
376 | ideal J,I0,radI0,L,radI1,I2,radI2; |
---|
377 | int l,n; |
---|
378 | intvec op = option(get); |
---|
379 | matrix M; |
---|
380 | |
---|
381 | option(redSB); |
---|
382 | list m = mstd(I); |
---|
383 | option(set,op); |
---|
384 | |
---|
385 | int d = dim(m[1]); |
---|
386 | if(d==-1) |
---|
387 | { |
---|
388 | return(ideal(1)); |
---|
389 | } |
---|
390 | |
---|
391 | if(strategy != 2) |
---|
392 | { |
---|
393 | ///////////////////////////////////////////// |
---|
394 | // Computing the equidimensional radical // |
---|
395 | // via regular sequenves // |
---|
396 | ///////////////////////////////////////////// |
---|
397 | |
---|
398 | if(printlevel > 2){"Trying to find a regular sequence.";} |
---|
399 | int cod = nvars(basering)-d; |
---|
400 | |
---|
401 | //Complete intersection case: |
---|
402 | if(cod==size(m[2])) |
---|
403 | { |
---|
404 | J = aJacob(m[2],d); |
---|
405 | if(printlevel > 2){"Leaving equiRadEHV.";} |
---|
406 | return(quotient(m[2],J)); |
---|
407 | } |
---|
408 | |
---|
409 | //First codim elements of I are a complete intersection: |
---|
410 | for(l=1; l<=cod; l++) |
---|
411 | { |
---|
412 | I0[l] = m[2][l]; |
---|
413 | } |
---|
414 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
415 | |
---|
416 | //Last codim elements of I are a complete intersection: |
---|
417 | if(n!=0) |
---|
418 | { |
---|
419 | for(l=1; l<=cod; l++) |
---|
420 | { |
---|
421 | I0[l] = m[2][size(m[2])-l+1]; |
---|
422 | } |
---|
423 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
424 | } |
---|
425 | |
---|
426 | //Taking a generic linear combination of the input: |
---|
427 | if(n!=0) |
---|
428 | { |
---|
429 | M = transpose(sparsetriag(size(m[2]),cod,95,1)); |
---|
430 | I0 = ideal(M*transpose(m[2])); |
---|
431 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
432 | } |
---|
433 | |
---|
434 | //Taking a more generic linear combination of the input: |
---|
435 | if(n!=0) |
---|
436 | { |
---|
437 | while(strategy == 1 and n!=0) |
---|
438 | { |
---|
439 | M = transpose(sparsetriag(size(m[2]),cod,0,100)); |
---|
440 | I0 = ideal(M*transpose(m[2])); |
---|
441 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
442 | } |
---|
443 | } |
---|
444 | |
---|
445 | if(n==0) |
---|
446 | { |
---|
447 | J = aJacob(I0,d); |
---|
448 | if(printlevel > 2){"1st quotient.";} |
---|
449 | radI0 = quotient(I0,J); |
---|
450 | if(printlevel > 2){"2nd quotient.";} |
---|
451 | L = quotient(radI0,m[2]); |
---|
452 | if(printlevel > 2){"3rd quotient.";} |
---|
453 | radI1 = quotient(radI0,L); |
---|
454 | attrib(radI1,"isEquidimensional",1); |
---|
455 | attrib(radI1,"isRadical",1); |
---|
456 | if(printlevel > 2){"Leaving equiRadEHV.";} |
---|
457 | return(radI1); |
---|
458 | } |
---|
459 | } |
---|
460 | |
---|
461 | //////////////////////////////////////////////////// |
---|
462 | // Computing the equidimensional radical directly // |
---|
463 | //////////////////////////////////////////////////// |
---|
464 | |
---|
465 | if(printlevel > 2){"Computing the equidimensional radical directly";} |
---|
466 | |
---|
467 | //Compute the equidimensional part depending on the chosen strategy |
---|
468 | if(equStr == 0) |
---|
469 | { |
---|
470 | I = equiMaxEHV(I); |
---|
471 | } |
---|
472 | if(equStr == 1) |
---|
473 | { |
---|
474 | I = equidimMax(I); |
---|
475 | } |
---|
476 | int a = nvars(basering)-1; |
---|
477 | |
---|
478 | while(a > d) |
---|
479 | { |
---|
480 | if(printlevel > 2){"While-Loop: "+string(a);} |
---|
481 | J = aJacob(I,a); |
---|
482 | while(dim(groebner(J+I))==d) |
---|
483 | { |
---|
484 | if(printlevel > 2){"Quotient-Computation.";} |
---|
485 | I = quotient(I,J); |
---|
486 | if(printlevel > 2){"Computing the a-th Jacobian";} |
---|
487 | J = aJacob(I,a); |
---|
488 | } |
---|
489 | a = a-1; |
---|
490 | } |
---|
491 | if(printlevel > 2){"We left While-Loop.";} |
---|
492 | if(printlevel > 2){"Computing the a-th Jacobian";} |
---|
493 | J = aJacob(I,d); |
---|
494 | if(printlevel > 2){"Quotient-Computation.";} |
---|
495 | I = quotient(I,J); |
---|
496 | attrib(I,"isEquidimensional",1); |
---|
497 | attrib(I,"isRadical",1); |
---|
498 | if(printlevel > 2){"Leaving equiRadEHV.";} |
---|
499 | return(I); |
---|
500 | } |
---|
501 | |
---|
502 | example |
---|
503 | { |
---|
504 | "EXAMPLE:"; |
---|
505 | echo = 2; |
---|
506 | ring r = 0,(x,y,z),dp; |
---|
507 | poly p = z2+1; |
---|
508 | poly q = z3+2; |
---|
509 | ideal i = p*q^2,y-z2; |
---|
510 | ideal pr= equiRadEHV(i); |
---|
511 | pr; |
---|
512 | } |
---|
513 | |
---|
514 | ///////////////////////////////////////////////////////////////////// |
---|
515 | proc radEHV(ideal I, list #) |
---|
516 | "USAGE: radEHV(I [,Strategy]); ideal I, Strategy list |
---|
517 | RETURN: ideal, the radical of I |
---|
518 | NOTE: uses the algorithm of Eisenbud/Huneke/Vasconcelos |
---|
519 | Works only in characteristic 0 or p large. |
---|
520 | The (optional) second argument determines the strategy used: |
---|
521 | Strategy[1] > strategy for the equidimensional part |
---|
522 | = 0 : uses equiMaxEHV |
---|
523 | = 1 : uses equidimMax |
---|
524 | Strategy[2] > strategy for the radical |
---|
525 | = 0 : combination of strategy 1 and 2 |
---|
526 | = 1 : computation of the radical just with the |
---|
527 | help of regular sequences |
---|
528 | = 2 : does not try to find a regular sequence |
---|
529 | Strategy[3] > strategy for the computation of ideal quotients |
---|
530 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
531 | If no second argument is given then Strategy=(0,0,0) is used. |
---|
532 | EXAMPLE: example radEHV; shows an example" |
---|
533 | { |
---|
534 | if(printlevel > 2){"Entering radEHV.";} |
---|
535 | if(attrib(basering,"global")==0) |
---|
536 | { |
---|
537 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
538 | } |
---|
539 | if (size(I) == 0) { return( ideal(0) ); } |
---|
540 | //Compute the equidimensional radical J of I. |
---|
541 | ideal J = equiRadEHV(I,#); |
---|
542 | |
---|
543 | //If I is the entire ring... |
---|
544 | if(deg(J[1]) <= 0) |
---|
545 | { |
---|
546 | //...then return the ideal generated by 1... |
---|
547 | return(ideal(1)); |
---|
548 | } |
---|
549 | |
---|
550 | //...else remove the maximal dimensional components and |
---|
551 | //compute the radical K of the lower dimensional components ... |
---|
552 | ideal K = radEHV(sat(I,J)[1],#); |
---|
553 | |
---|
554 | //..and intersect it with J. |
---|
555 | K = intersect(J,K); |
---|
556 | attrib(K,"isRadical",1); |
---|
557 | return(K); |
---|
558 | } |
---|
559 | |
---|
560 | example |
---|
561 | { |
---|
562 | "EXAMPLE:"; |
---|
563 | echo = 2; |
---|
564 | ring r = 0,(x,y,z),dp; |
---|
565 | poly p = z2+1; |
---|
566 | poly q = z3+2; |
---|
567 | ideal i = p*q^2,y-z2; |
---|
568 | ideal pr= radical(i); |
---|
569 | pr; |
---|
570 | } |
---|
571 | |
---|
572 | ///////////////////////////////////////////////////////////////////// |
---|
573 | proc IntAssOfDim1(ideal I, int e) |
---|
574 | "USAGE: IntAssOfDim1(I,e); I idea, e integer |
---|
575 | RETURN: ideal, the intersection of the associated primes of I having dimension e |
---|
576 | EXAMPLE: example IntAssOfDim1; shows an example" |
---|
577 | { |
---|
578 | if(attrib(basering,"global")==0) |
---|
579 | { |
---|
580 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
581 | } |
---|
582 | int g = nvars(basering) - e; |
---|
583 | |
---|
584 | //Compute a resolution of I. |
---|
585 | ideal J = groebner(I); |
---|
586 | if(homog(I)==1) |
---|
587 | { |
---|
588 | list re = sres(J,g+1); |
---|
589 | re = minres(re); |
---|
590 | } |
---|
591 | else |
---|
592 | { |
---|
593 | list re = mres(I,g+1); |
---|
594 | } |
---|
595 | |
---|
596 | ideal ann = AnnExtEHV(g,re); |
---|
597 | int cod = nvars(basering) - dim(groebner(ann)); |
---|
598 | //If the codimension of I_g:=Ann(Ext^g(R/I,I)) equals g... |
---|
599 | if(cod == g) |
---|
600 | { |
---|
601 | //...then return the equidimensional radical of I_g... |
---|
602 | ann = equiRadEHV(ann); |
---|
603 | attrib(ann,"isEquidimensional",1); |
---|
604 | attrib(ann,"isRadical",1); |
---|
605 | return(ann); |
---|
606 | } |
---|
607 | else |
---|
608 | { |
---|
609 | //...else I has no associated primes of dimension e. |
---|
610 | return(ideal(1)); |
---|
611 | } |
---|
612 | } |
---|
613 | |
---|
614 | example |
---|
615 | { |
---|
616 | "EXAMPLE:"; |
---|
617 | echo = 2; |
---|
618 | ring r = 0,(x,y,z),dp; |
---|
619 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
620 | IntAssOfDim1(I,1); |
---|
621 | } |
---|
622 | |
---|
623 | ///////////////////////////////////////////////////////////////////// |
---|
624 | proc IntAssOfDim2(ideal I, int e) |
---|
625 | "USAGE: IntAssOfDim2(I,e); I ideal, e integer |
---|
626 | RETURN: ideal, the intersection of the associated primes of I having dimension e |
---|
627 | EXAMPLE: example IntAssOfDim2; shows an example" |
---|
628 | { |
---|
629 | if(attrib(basering,"global")==0) |
---|
630 | { |
---|
631 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
632 | } |
---|
633 | ideal I1 = removeComponent(I,e); |
---|
634 | ideal I2 = removeComponent(I,e+1); |
---|
635 | ideal b = quotient(I1,I2); |
---|
636 | b = equiRadEHV(b); |
---|
637 | return(b); |
---|
638 | } |
---|
639 | |
---|
640 | example |
---|
641 | { |
---|
642 | "EXAMPLE:"; |
---|
643 | echo = 2; |
---|
644 | ring r = 0,(x,y,z),dp; |
---|
645 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
646 | IntAssOfDim2(I,1); |
---|
647 | } |
---|
648 | |
---|
649 | |
---|
650 | ///////////////////////////////////////////////////////////////////// |
---|
651 | // // |
---|
652 | // Z E R O - D I M E N S I O N A L D E C O M P O S I T I O N // |
---|
653 | // // |
---|
654 | ///////////////////////////////////////////////////////////////////// |
---|
655 | |
---|
656 | ///////////////////////////////////////////////////////////////////// |
---|
657 | proc decompEHV(ideal I) |
---|
658 | "USAGE: decompEHV(I); I zero-dimensional radical ideal |
---|
659 | RETURN: list, the associated primes of I |
---|
660 | EXAMPLE: example decompEHV; shows an example" |
---|
661 | { |
---|
662 | if(printlevel > 2){"Entering decompEHV.";} |
---|
663 | if(attrib(basering,"global")==0) |
---|
664 | { |
---|
665 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
666 | } |
---|
667 | |
---|
668 | int e,m,vd,nfact; |
---|
669 | list l,k,L; |
---|
670 | poly f,h; |
---|
671 | ideal J; |
---|
672 | def base = basering; |
---|
673 | |
---|
674 | if( attrib(I,"isSB")!=1 ) |
---|
675 | { |
---|
676 | I = groebner(I); |
---|
677 | } |
---|
678 | |
---|
679 | while(1) |
---|
680 | { |
---|
681 | //Choose a random polynomial f from R. |
---|
682 | e = random(0,100); |
---|
683 | f = sparsepoly(e); |
---|
684 | |
---|
685 | //Check if f lies not in I. |
---|
686 | if(reduce(f,I)!=0) |
---|
687 | { |
---|
688 | J = quotient1(I,f); |
---|
689 | |
---|
690 | //If f is a zerodivisor modulo I... |
---|
691 | if( isSubset(J,I) == 0 ) |
---|
692 | { |
---|
693 | //...then use recursion... |
---|
694 | if(printlevel > 2){"We found a zerodivisor -- recursion";} |
---|
695 | l = decompEHV(J) + decompEHV(I+f); |
---|
696 | return(l); |
---|
697 | } |
---|
698 | if(printlevel > 2){"We found a non-zero-divisor.";} |
---|
699 | |
---|
700 | //...else compute the vectorspace dimension vd of I and... |
---|
701 | vd = vdim(I); |
---|
702 | |
---|
703 | //...compute m minimal such that 1,f,f^2,...,f^m are linearly dependent. |
---|
704 | qring Q = I; |
---|
705 | poly g = fetch(base,f); |
---|
706 | k = algDependent(g); |
---|
707 | def R = k[2]; |
---|
708 | setring R; |
---|
709 | if(size(ker)!=1) |
---|
710 | { |
---|
711 | //Calculate a generator for ker. |
---|
712 | ker = mstd(ker)[2]; |
---|
713 | } |
---|
714 | poly g = ker[1]; |
---|
715 | m = deg(g); |
---|
716 | |
---|
717 | //If m and vd coincide... |
---|
718 | if(m==vd) |
---|
719 | { |
---|
720 | if(printlevel > 2){"We have a good candidate.";} |
---|
721 | //...then factorize g. |
---|
722 | if(printlevel > 2){"Factorizing.";} |
---|
723 | L = factorize(g,2); //returns non-constant factors and multiplicities |
---|
724 | nfact = size(L[1]); |
---|
725 | |
---|
726 | //If g is irreducible... |
---|
727 | if(nfact==1 and L[2][1]==1) |
---|
728 | { |
---|
729 | if(printlevel > 2){"The element is irreducible.";} |
---|
730 | setring base; |
---|
731 | //..then I is a maximal ideal... |
---|
732 | l[1] = I; |
---|
733 | kill R,Q; |
---|
734 | return(l); |
---|
735 | } |
---|
736 | //...else... |
---|
737 | else |
---|
738 | { |
---|
739 | |
---|
740 | if(printlevel > 2){"The element is not irreducible -- recursion.";} |
---|
741 | //...take a non-trivial factor g1 of g ... |
---|
742 | poly g1 = L[1][1]; |
---|
743 | //..and insert f in g1... |
---|
744 | ring newR = create_ring(ringlist(R)[1], "(y(1))", "(" + ordstr(R) + ")", "no_minpoly"); |
---|
745 | poly g2 = imap(R,g1); |
---|
746 | setring base; |
---|
747 | h = fetch(newR,g2); |
---|
748 | h = subst(h,var(1),f); |
---|
749 | //...and use recursion... |
---|
750 | kill R,Q,newR; |
---|
751 | l = l + decompEHV(quotient1(I,h)); |
---|
752 | l = l + decompEHV(I+h); |
---|
753 | return(l); |
---|
754 | } |
---|
755 | } |
---|
756 | setring base; |
---|
757 | kill R,Q; |
---|
758 | } |
---|
759 | } |
---|
760 | } |
---|
761 | |
---|
762 | example |
---|
763 | { |
---|
764 | "EXAMPLE:"; |
---|
765 | echo = 2; |
---|
766 | ring r = 32003,(x,y),dp; |
---|
767 | ideal i = x2+y2-10,x2+xy+2y2-16; |
---|
768 | decompEHV(i); |
---|
769 | } |
---|
770 | |
---|
771 | |
---|
772 | ///////////////////////////////////////////////////////////////////// |
---|
773 | // // |
---|
774 | // A S S O C I A T E D P R I M E S // |
---|
775 | // // |
---|
776 | ///////////////////////////////////////////////////////////////////// |
---|
777 | |
---|
778 | |
---|
779 | ///////////////////////////////////////////////////////////////////// |
---|
780 | static proc idempotent(ideal I) |
---|
781 | "USAGE: idempotent(I); ideal I (weighted) homogeneous radical ideal, |
---|
782 | I intersected K[x(1),...,x(k)] zero-dimensional |
---|
783 | where deg(x(i))=0 for all i <= k and deg(x(i))>0 for all i>k. |
---|
784 | RETURN: a list of ideals I(1),...,I(t) such that |
---|
785 | K[x]/I = K[x]/I(1) x ... x K[x]/I(t)" |
---|
786 | { |
---|
787 | if(printlevel > 2){"Entering idempotent.";} |
---|
788 | if(attrib(basering,"global")==0) |
---|
789 | { |
---|
790 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
791 | } |
---|
792 | int n = nvars(basering); |
---|
793 | poly f,g; |
---|
794 | string j; |
---|
795 | int i,k,splits; |
---|
796 | list l; |
---|
797 | |
---|
798 | //Collect all variables of degree 0. |
---|
799 | for(i=1; i<=n; i++) |
---|
800 | { |
---|
801 | if(deg(var(i)) > 0) |
---|
802 | { |
---|
803 | f = f*var(i); |
---|
804 | } |
---|
805 | else |
---|
806 | { |
---|
807 | splits = 1; |
---|
808 | j = j + string(var(i)) + "," ; |
---|
809 | } |
---|
810 | } |
---|
811 | //If there are no variables of degree 0 |
---|
812 | //then there are no idempotents and we are done... |
---|
813 | if(splits == 0) |
---|
814 | { |
---|
815 | l[1]=I; |
---|
816 | return(l); |
---|
817 | } |
---|
818 | |
---|
819 | //...else compute J = I intersected K[x(1),...,x(k)]... |
---|
820 | ideal J = eliminate(I,f); |
---|
821 | def base = basering; |
---|
822 | j = j[1,size(j)-1]; |
---|
823 | j = "ring @r = (" + charstr(basering) + "),(" + j + "),(" + ordstr(basering) + ");"; |
---|
824 | execute(j); |
---|
825 | ideal J = imap(base,J); |
---|
826 | |
---|
827 | //...and compute the associated primes of the zero-dimensional ideal J. |
---|
828 | list L = decompEHV(J); |
---|
829 | int s = size(L); |
---|
830 | ideal K,Z; |
---|
831 | poly g; |
---|
832 | //For each associated prime ideal P_i of J... |
---|
833 | for(i=1; i<=s; i++) |
---|
834 | { |
---|
835 | K = ideal(1); |
---|
836 | //...comnpute the intersection K of the other associated prime ideals... |
---|
837 | for(k=1; k<=s; k++) |
---|
838 | { |
---|
839 | if(i!=k) |
---|
840 | { |
---|
841 | K = intersect(K,L[k]); |
---|
842 | } |
---|
843 | } |
---|
844 | |
---|
845 | //...and find an element that lies in K but not in P_i... |
---|
846 | g = randomid(K,1)[1]; |
---|
847 | Z = L[i]; |
---|
848 | Z = std(Z); |
---|
849 | while(reduce(g,Z)==0) |
---|
850 | { |
---|
851 | g = randomid(K,1)[1]; |
---|
852 | } |
---|
853 | setring base; |
---|
854 | g = imap(@r,g); |
---|
855 | //...and compute the corresponding ideal I(i) |
---|
856 | l[i] = quotient(I,g); |
---|
857 | setring @r; |
---|
858 | |
---|
859 | } |
---|
860 | setring base; |
---|
861 | return(l); |
---|
862 | } |
---|
863 | |
---|
864 | |
---|
865 | ///////////////////////////////////////////////////////////////////// |
---|
866 | static proc equiAssEHV(ideal I) |
---|
867 | "USAGE: equiAssEHV(I); I equidimensional, radical, and homogeneous ideal |
---|
868 | RETURN: a list, the associated prime ideals of I" |
---|
869 | { |
---|
870 | if(printlevel > 2){"Entering equiAssEHV.";} |
---|
871 | if(attrib(basering,"global")==0) |
---|
872 | { |
---|
873 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
874 | } |
---|
875 | list L; |
---|
876 | def base = basering; |
---|
877 | int n = nvars(basering); |
---|
878 | int i,j; |
---|
879 | |
---|
880 | //Compute the normalization of I. |
---|
881 | if(printlevel > 2){"Entering Normalization.";} |
---|
882 | list norOut = normal(I, "noDeco"); |
---|
883 | list K = norOut[1]; |
---|
884 | if(printlevel > 2){"Leaving Normalisation.";} |
---|
885 | |
---|
886 | //The normalization algorithm returns k factors. |
---|
887 | int k = size(K); |
---|
888 | if(printlevel > 1){"Normalization algorithm splits ideal in " + string(k) + " factors.";} |
---|
889 | |
---|
890 | //Examine each factor of the normalization. |
---|
891 | def P; |
---|
892 | for(i=1; i<=k; i++) |
---|
893 | { |
---|
894 | P = K[i]; |
---|
895 | setring P; |
---|
896 | |
---|
897 | //Use procedure idempotent to split the i-th factor of |
---|
898 | //the normalization in a product of integral domains. |
---|
899 | if(printlevel > 1){"Examining " + string(i) +". factor.";} |
---|
900 | list l = idempotent(norid); |
---|
901 | int leng = size(l); |
---|
902 | if(printlevel > 1){"Idempotent algorithm splits factor " + string(i) + " in " + string(leng) + " factors.";} |
---|
903 | |
---|
904 | //Intersect the minimal primes corresponding to the |
---|
905 | //integral domains obtained from idempotent with the groundring, |
---|
906 | //i.e. compute the preimages w.r.t. the corresponding normalization map |
---|
907 | ideal J; |
---|
908 | for(j=1; j<=leng; j++) |
---|
909 | { |
---|
910 | J = l[j]; |
---|
911 | setring base; |
---|
912 | L[size(L)+1] = preimage(P,normap,J); |
---|
913 | setring P; |
---|
914 | } |
---|
915 | kill l, leng, J; |
---|
916 | setring base; |
---|
917 | } |
---|
918 | |
---|
919 | return(L); |
---|
920 | } |
---|
921 | |
---|
922 | |
---|
923 | ///////////////////////////////////////////////////////////////////// |
---|
924 | proc AssEHV(ideal I, list #) |
---|
925 | "USAGE: AssEHV(I [,Strategy]); I Ideal, Strategy list |
---|
926 | RETURN: a list, the associated prime ideals of I |
---|
927 | NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos. |
---|
928 | The (optional) second argument determines the strategy used: |
---|
929 | Strategy[1] > strategy for the equidimensional part |
---|
930 | = 0 : uses equiMaxEHV |
---|
931 | = 1 : uses equidimMax |
---|
932 | Strategy[2] > strategy for the equidimensional radical |
---|
933 | = 0 : uses equiRadEHV |
---|
934 | = 1 : uses equiRadical |
---|
935 | Strategy[3] > strategy for equiRadEHV |
---|
936 | = 0 : combination of strategy 1 and 2 |
---|
937 | = 1 : computation of the radical just with the |
---|
938 | help of regular sequences |
---|
939 | = 2 : does not try to find a regular sequence |
---|
940 | Strategy[4] > strategy for the computation of ideal quotients |
---|
941 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
942 | If no second argument is given, Strategy=(0,0,0,0) is used. |
---|
943 | EXAMPLE: example AssEHV; shows an example" |
---|
944 | { |
---|
945 | if(attrib(basering,"global")==0) |
---|
946 | { |
---|
947 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
948 | } |
---|
949 | if(printlevel > 2){"Entering AssEHV";} |
---|
950 | |
---|
951 | //Specify the strategy to be used. |
---|
952 | if(size(#)==0) |
---|
953 | { |
---|
954 | # = 0,0,0,0; |
---|
955 | } |
---|
956 | if(size(#)==1) |
---|
957 | { |
---|
958 | # = #[1],0,0,0; |
---|
959 | } |
---|
960 | if(size(#)==2) |
---|
961 | { |
---|
962 | # = #[1],#[2],0,0; |
---|
963 | } |
---|
964 | if(size(#)==3) |
---|
965 | { |
---|
966 | # = #[1],#[2],#[3],0; |
---|
967 | } |
---|
968 | |
---|
969 | list L; |
---|
970 | ideal K; |
---|
971 | def base = basering; |
---|
972 | int m,j; |
---|
973 | ideal J = groebner(I); |
---|
974 | int n = nvars(basering); |
---|
975 | int d = dim(J); |
---|
976 | |
---|
977 | //Compute a resolution of I. |
---|
978 | if(printlevel > 2){"Computing resolution.";} |
---|
979 | if(homog(I)==1) |
---|
980 | { |
---|
981 | list re = sres(J,0); |
---|
982 | re = minres(re); |
---|
983 | } |
---|
984 | else |
---|
985 | { |
---|
986 | list re = mres(I,0); |
---|
987 | } |
---|
988 | ideal ann; |
---|
989 | int cod; |
---|
990 | |
---|
991 | //For 0<=i<= dim(I) compute the intersection of the i-dimensional associated primes of I |
---|
992 | for(int i=0; i<=d; i++) |
---|
993 | { |
---|
994 | if(printlevel > 1){"Are there components of dimension " + string(i) + "?";} |
---|
995 | ann = AnnExtEHV(n-i,re); |
---|
996 | cod = n - dim(groebner(ann)); |
---|
997 | |
---|
998 | //If there are associated primes of dimension i... |
---|
999 | if(cod == n-i) |
---|
1000 | { |
---|
1001 | if(printlevel > 1){"Yes. There are components of dimension " + string(i) + ".";} |
---|
1002 | //...then compute the intersection K of all associated primes of I of dimension i |
---|
1003 | if(#[2]==0) |
---|
1004 | { |
---|
1005 | if(size(#) > 3) |
---|
1006 | { |
---|
1007 | K = equiRadEHV(ann,#[1],#[3],#[4]); |
---|
1008 | } |
---|
1009 | if(size(#) > 2) |
---|
1010 | { |
---|
1011 | K = equiRadEHV(ann,#[1],#[3]); |
---|
1012 | } |
---|
1013 | else |
---|
1014 | { |
---|
1015 | K = equiRadEHV(ann,#[1]); |
---|
1016 | } |
---|
1017 | } |
---|
1018 | if(#[2]==1) |
---|
1019 | { |
---|
1020 | K = equiRadical(ann); |
---|
1021 | } |
---|
1022 | attrib(K,"isEquidimensional",1); |
---|
1023 | attrib(K,"isRadical",1); |
---|
1024 | |
---|
1025 | //If K is already homogeneous then use equiAssEHV to recover the associated primes of K... |
---|
1026 | if(homog(K)==1) |
---|
1027 | { |
---|
1028 | if(printlevel > 2){"Input already homogeneous.";} |
---|
1029 | L = L + equiAssEHV(K); |
---|
1030 | } |
---|
1031 | |
---|
1032 | //...else... |
---|
1033 | else |
---|
1034 | { |
---|
1035 | //...homogenize K w.r.t. t,... |
---|
1036 | if(printlevel > 2){"Input not homogeneous; must homogenize.";} |
---|
1037 | def homoR=changeord(list(list("dp",1:nvars(basering)))); |
---|
1038 | setring homoR; |
---|
1039 | ideal homoJ = fetch(base,K); |
---|
1040 | homoJ = groebner(homoJ); |
---|
1041 | list l2; |
---|
1042 | for (int ii = 1; ii <= n; ii++) |
---|
1043 | { |
---|
1044 | l2[i] = "x("+string(ii)+")"; |
---|
1045 | } |
---|
1046 | l2[n+1] = "t"; |
---|
1047 | ring newR = create_ring(ringlist(base)[1], l2, "dp", "no_minpoly"); |
---|
1048 | ideal homoK = fetch(homoR,homoJ); |
---|
1049 | homoK = homog(homoK,t); |
---|
1050 | attrib(homoK,"isEquidimensional",1); |
---|
1051 | attrib(homoK,"isRadical",1); |
---|
1052 | |
---|
1053 | //...compute the associated primes of the homogenization using equiAssEHV,... |
---|
1054 | list l = equiAssEHV(homoK); |
---|
1055 | |
---|
1056 | //...and set t=1 in the generators of the associated primes just computed. |
---|
1057 | ideal Z; |
---|
1058 | for(j=1; j<=size(l); j++) |
---|
1059 | { |
---|
1060 | Z = subst(l[j],t,1); |
---|
1061 | setring base; |
---|
1062 | L[size(L)+1] = fetch(newR,Z); |
---|
1063 | setring newR; |
---|
1064 | } |
---|
1065 | setring base; |
---|
1066 | kill homoR; |
---|
1067 | kill newR; |
---|
1068 | } |
---|
1069 | } |
---|
1070 | else |
---|
1071 | { |
---|
1072 | if(printlevel > 1){"No. There are no components of dimension " + string(i) + ".";} |
---|
1073 | } |
---|
1074 | } |
---|
1075 | return(L); |
---|
1076 | } |
---|
1077 | |
---|
1078 | example |
---|
1079 | { |
---|
1080 | "EXAMPLE:"; |
---|
1081 | echo = 2; |
---|
1082 | ring r = 0,(x,y,z),dp; |
---|
1083 | poly p = z2+1; |
---|
1084 | poly q = z3+2; |
---|
1085 | ideal i = p*q^2,y-z2; |
---|
1086 | list pr = AssEHV(i); |
---|
1087 | pr; |
---|
1088 | } |
---|
1089 | |
---|
1090 | ///////////////////////////////////////////////////////////////////// |
---|
1091 | proc minAssEHV(ideal I, list #) |
---|
1092 | "USAGE: minAssEHV(I [,Strategy]); I ideal, Strategy list |
---|
1093 | RETURN: a list, the minimal associated prime ideals of I |
---|
1094 | NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos. |
---|
1095 | The (optional) second argument determines the strategy used: |
---|
1096 | Strategy[1] > strategy for the equidimensional part |
---|
1097 | = 0 : uses equiMaxEHV |
---|
1098 | = 1 : uses equidimMax |
---|
1099 | Strategy[2] > strategy for the equidimensional radical |
---|
1100 | = 0 : uses equiRadEHV, resp. radicalEHV |
---|
1101 | = 1 : uses equiRadical, resp. radical |
---|
1102 | Strategy[3] > strategy for equiRadEHV |
---|
1103 | = 0 : combination of strategy 1 and 2 |
---|
1104 | = 1 : computation of the radical just with the |
---|
1105 | help of regular sequences |
---|
1106 | = 2 : does not try to find a regular sequence |
---|
1107 | Strategy[4] > strategy for the computation of ideal quotients |
---|
1108 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
1109 | If no second argument is given, Strategy=(0,0,0,0) is used. |
---|
1110 | EXAMPLE: example minAssEHV; shows an example" |
---|
1111 | { |
---|
1112 | if(attrib(basering,"global")==0) |
---|
1113 | {ERROR("// Not implemented for this ordering, please change to global ordering.");} |
---|
1114 | |
---|
1115 | //Specify the strategy to be used. |
---|
1116 | if(size(#)==0) |
---|
1117 | { |
---|
1118 | # = 0,0,0,0; |
---|
1119 | } |
---|
1120 | if(size(#)==1) |
---|
1121 | { |
---|
1122 | # = #[1],0,0,0; |
---|
1123 | } |
---|
1124 | if(size(#)==2) |
---|
1125 | { |
---|
1126 | # = #[1],#[2],0,0; |
---|
1127 | } |
---|
1128 | if(size(#)==3) |
---|
1129 | { |
---|
1130 | # = #[1],#[2],#[3],0; |
---|
1131 | } |
---|
1132 | |
---|
1133 | //Compute the radical of I. |
---|
1134 | if(#[2]==0) |
---|
1135 | { |
---|
1136 | I = radEHV(I,#); |
---|
1137 | } |
---|
1138 | if(#[2]==1) |
---|
1139 | { |
---|
1140 | I = radical(I); |
---|
1141 | } |
---|
1142 | |
---|
1143 | //Compute the associated primes of the radical. |
---|
1144 | return(AssEHV(I,#)); |
---|
1145 | } |
---|
1146 | |
---|
1147 | example |
---|
1148 | { |
---|
1149 | "EXAMPLE:"; |
---|
1150 | echo = 2; |
---|
1151 | ring r = 0,(x,y,z),dp; |
---|
1152 | poly p = z2+1; |
---|
1153 | poly q = z3+2; |
---|
1154 | ideal i = p*q^2,y-z2; |
---|
1155 | list pr = minAssEHV(i); |
---|
1156 | pr; |
---|
1157 | } |
---|
1158 | |
---|
1159 | |
---|
1160 | ///////////////////////////////////////////////////////////////////// |
---|
1161 | // // |
---|
1162 | // P R I M A R Y D E C O M P O S I T I O N // |
---|
1163 | // // |
---|
1164 | ///////////////////////////////////////////////////////////////////// |
---|
1165 | |
---|
1166 | |
---|
1167 | ///////////////////////////////////////////////////////////////////// |
---|
1168 | proc localize(ideal I, ideal P, list l) |
---|
1169 | "USAGE: localize(I,P,l); I ideal, P an associated prime ideal of I, |
---|
1170 | l list of all associated primes of I |
---|
1171 | RETURN: ideal, the contraction of the ideal generated by I |
---|
1172 | in the localization w.r.t P |
---|
1173 | EXAMPLE: example localize; shows an example" |
---|
1174 | { |
---|
1175 | if(attrib(basering,"global")==0) |
---|
1176 | { |
---|
1177 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
1178 | } |
---|
1179 | |
---|
1180 | ideal Intersection = ideal(1); |
---|
1181 | int s = size(l); |
---|
1182 | if(attrib(P,"isSB")!=1) |
---|
1183 | { |
---|
1184 | P = groebner(P); |
---|
1185 | } |
---|
1186 | |
---|
1187 | //Compute the intersection of all associated primes of I that are not contained in P. |
---|
1188 | for(int i=1; i<=s; i++) |
---|
1189 | { |
---|
1190 | if(isSubset(l[i],P)!=1) |
---|
1191 | { |
---|
1192 | Intersection = intersect(Intersection,l[i]); |
---|
1193 | } |
---|
1194 | } |
---|
1195 | Intersection = groebner(Intersection); |
---|
1196 | |
---|
1197 | //If the intersection is the entire ring... |
---|
1198 | if(reduce(1,Intersection)==0) |
---|
1199 | { |
---|
1200 | //...then return I... |
---|
1201 | return(I); |
---|
1202 | } |
---|
1203 | //...else try to find an element f that lies in the intersection but outside P... |
---|
1204 | poly f = 0; |
---|
1205 | while(reduce(f,P) == 0) |
---|
1206 | { |
---|
1207 | f = randomid(Intersection,1)[1]; |
---|
1208 | } |
---|
1209 | |
---|
1210 | //...and saturate I w.r.t. f. |
---|
1211 | I = sat(I,f)[1]; |
---|
1212 | return(I); |
---|
1213 | } |
---|
1214 | |
---|
1215 | ///////////////////////////////////////////////////////////////////// |
---|
1216 | proc componentEHV(ideal I, ideal P, list L, list #) |
---|
1217 | "USAGE: componentEHV(I,P,L [,Strategy]); |
---|
1218 | I ideal, P associated prime of I, |
---|
1219 | L list of all associated primes of I, Strategy list |
---|
1220 | RETURN: ideal, a P-primary component Q for I |
---|
1221 | NOTE: The (optional) second argument determines the strategy used: |
---|
1222 | Strategy[1] > strategy for equidimensional part |
---|
1223 | = 0 : uses equiMaxEHV |
---|
1224 | = 1 : uses equidimMax |
---|
1225 | If no second argument is given then Strategy=0 is used. |
---|
1226 | EXAMPLE: example componentEHV; shows an example" |
---|
1227 | { |
---|
1228 | if(printlevel > 2){"Entering componentEHV.";} |
---|
1229 | if(attrib(basering,"global")==0) |
---|
1230 | { |
---|
1231 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
1232 | } |
---|
1233 | |
---|
1234 | //If no strategy is specified use standard strategy. |
---|
1235 | if(size(#)==0) |
---|
1236 | { |
---|
1237 | # = 0; |
---|
1238 | } |
---|
1239 | |
---|
1240 | ideal T = P; |
---|
1241 | ideal Q; |
---|
1242 | |
---|
1243 | //Compute the localization of I at P... |
---|
1244 | ideal IP = groebner(localize(I,P,L)); |
---|
1245 | |
---|
1246 | //...and compute the saturation of the localization w.r.t. P. |
---|
1247 | ideal IP2 = sat(IP,P)[1]; |
---|
1248 | |
---|
1249 | //As long as we have not found a primary component... |
---|
1250 | int isPrimaryComponent = 0; |
---|
1251 | while(isPrimaryComponent!=1) |
---|
1252 | { |
---|
1253 | //...compute the equidimensional part Q of I+P^n... |
---|
1254 | if(#[1]==0) |
---|
1255 | { |
---|
1256 | Q = equiMaxEHV(I+T); |
---|
1257 | } |
---|
1258 | if(#[1]==1) |
---|
1259 | { |
---|
1260 | Q = equidimMax(I+T); |
---|
1261 | } |
---|
1262 | //...and check if it is a primary component for P. |
---|
1263 | if(isSubset(intersect(IP2,Q),IP)==1) |
---|
1264 | { |
---|
1265 | isPrimaryComponent = 1; |
---|
1266 | } |
---|
1267 | else |
---|
1268 | { |
---|
1269 | T = T*P; |
---|
1270 | } |
---|
1271 | } |
---|
1272 | if(printlevel > 2){"Leaving componentEHV.";} |
---|
1273 | return(Q); |
---|
1274 | } |
---|
1275 | |
---|
1276 | example |
---|
1277 | { |
---|
1278 | "EXAMPLE:"; |
---|
1279 | echo = 2; |
---|
1280 | ring r = 0,(x,y,z),dp; |
---|
1281 | poly p = z2+1; |
---|
1282 | poly q = z3+2; |
---|
1283 | ideal i = p*q^2,y-z2; |
---|
1284 | list pr = AssEHV(i); |
---|
1285 | componentEHV(i,pr[1],pr); |
---|
1286 | } |
---|
1287 | |
---|
1288 | |
---|
1289 | ///////////////////////////////////////////////////////////////////// |
---|
1290 | proc primdecEHV(ideal I, list #) |
---|
1291 | "USAGE: primdecEHV(I [,Strategy]); I ideal, Strategy list |
---|
1292 | RETURN: a list pr of primary ideals and their associated primes: |
---|
1293 | pr[i][1] the i-th primary component, |
---|
1294 | pr[i][2] the i-th prime component. |
---|
1295 | NOTE: Algorithm of Eisenbud/Huneke/Vasconcelos. |
---|
1296 | The (optional) second argument determines the strategy used: |
---|
1297 | Strategy[1] > strategy for equidimensional part |
---|
1298 | = 0 : uses equiMaxEHV |
---|
1299 | = 1 : uses equidimMax |
---|
1300 | Strategy[2] > strategy for equidimensional radical |
---|
1301 | = 0 : uses equiRadEHV, resp. radicalEHV |
---|
1302 | = 1 : uses equiRadical, resp. radical |
---|
1303 | Strategy[3] > strategy for equiRadEHV |
---|
1304 | = 0 : combination of strategy 1 and 2 |
---|
1305 | = 1 : computation of the radical just with the |
---|
1306 | help of regular sequences |
---|
1307 | = 2 : does not try to find a regular sequence |
---|
1308 | Strategy[4] > strategy for the computation of ideal quotients |
---|
1309 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
1310 | If no second argument is given then Strategy=(0,0,0,0) is used. |
---|
1311 | EXAMPLE: example primdecEHV; shows an example" |
---|
1312 | { |
---|
1313 | if(printlevel > 2){"Entering primdecEHV.";} |
---|
1314 | if(attrib(basering,"global")==0) |
---|
1315 | { |
---|
1316 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
1317 | } |
---|
1318 | list L,K; |
---|
1319 | |
---|
1320 | //Specify the strategy to be used. |
---|
1321 | if(size(#)==0) |
---|
1322 | { |
---|
1323 | # = 0,0,0,0; |
---|
1324 | } |
---|
1325 | if(size(#)==1) |
---|
1326 | { |
---|
1327 | # = #[1],0,0,0; |
---|
1328 | } |
---|
1329 | if(size(#)==2) |
---|
1330 | { |
---|
1331 | # = #[1],#[2],0,0; |
---|
1332 | } |
---|
1333 | if(size(#)==3) |
---|
1334 | { |
---|
1335 | # = #[1],#[2],#[3],0; |
---|
1336 | } |
---|
1337 | |
---|
1338 | //Compute the associated primes of I... |
---|
1339 | L = AssEHV(I,#); |
---|
1340 | if(printlevel > 0){"We have " + string(size(L)) + " prime components.";} |
---|
1341 | |
---|
1342 | //...and compute for each associated prime of I a corresponding primary component. |
---|
1343 | int l = size(L); |
---|
1344 | for(int i=1; i<=l; i++) |
---|
1345 | { |
---|
1346 | K[i] = list(); |
---|
1347 | K[i][2] = L[i]; |
---|
1348 | K[i][1] = componentEHV(I,L[i],L,#[1]); |
---|
1349 | } |
---|
1350 | if(printlevel > 2){"Leaving primdecEHV.";} |
---|
1351 | return(K); |
---|
1352 | |
---|
1353 | } |
---|
1354 | |
---|
1355 | example |
---|
1356 | { |
---|
1357 | "EXAMPLE:"; |
---|
1358 | echo = 2; |
---|
1359 | ring r = 0,(x,y,z),dp; |
---|
1360 | poly p = z2+1; |
---|
1361 | poly q = z3+2; |
---|
1362 | ideal i = p*q^2,y-z2; |
---|
1363 | list pr = primdecEHV(i); |
---|
1364 | pr; |
---|
1365 | } |
---|
1366 | |
---|
1367 | |
---|
1368 | ///////////////////////////////////////////////////////////////////// |
---|
1369 | // // |
---|
1370 | // A L G O R I T H M S F O R T E S T I N G // |
---|
1371 | // // |
---|
1372 | ///////////////////////////////////////////////////////////////////// |
---|
1373 | proc removeRedundantComponents(list primdecResult) |
---|
1374 | "USAGE: removeRedundantComponents(L); L a primary decomposition |
---|
1375 | RETURN: an irredundant primary decomposition. |
---|
1376 | EXAMPLE: example removeRedundantComponents; shows an example" |
---|
1377 | proc removeRedundantComponents(list primdecResult) |
---|
1378 | { |
---|
1379 | int i,j; |
---|
1380 | |
---|
1381 | i = 1; |
---|
1382 | j = 1; |
---|
1383 | while ( i<=size(primdecResult) ) |
---|
1384 | { |
---|
1385 | ASSUME(1, i>0 && i<=size(primdecResult) ); |
---|
1386 | ASSUME(1, j>0 && j<=size(primdecResult) ); |
---|
1387 | |
---|
1388 | if (i!=j) |
---|
1389 | { |
---|
1390 | if ( isSubset( std(primdecResult[j][1]) , std(primdecResult[i][1]) ) ) |
---|
1391 | { |
---|
1392 | primdecResult[i] = primdecResult[j]; |
---|
1393 | primdecResult = delete(primdecResult,j); |
---|
1394 | j = j - 1; /*correct j due to deletion */ |
---|
1395 | if ( j<i ) { i = i -1; } /*correct i */ |
---|
1396 | } |
---|
1397 | } |
---|
1398 | j = j + 1; |
---|
1399 | if ( j > size(primdecResult) ) |
---|
1400 | { |
---|
1401 | i = i + 1; |
---|
1402 | j = 1; |
---|
1403 | } |
---|
1404 | } |
---|
1405 | |
---|
1406 | return (primdecResult); |
---|
1407 | } |
---|
1408 | example |
---|
1409 | { |
---|
1410 | "EXAMPLE:"; |
---|
1411 | ring rng = 0,(x,y),dp; |
---|
1412 | def L1 = list(list(ideal(x-1),ideal(x-1)),list(ideal(x-1),ideal(x-1)),list(ideal(y-2),ideal(y-2))); |
---|
1413 | def L2 = list(list(ideal(x-1),ideal(x-1)),list(ideal(y-2),ideal(y-2))); |
---|
1414 | L1; |
---|
1415 | L2; |
---|
1416 | ASSUME(0, primDecsAreEquivalent( L1, L2 ) ); |
---|
1417 | } |
---|
1418 | |
---|
1419 | ///////////////////////////////////////////////////////////////////// |
---|
1420 | proc primDecsAreEquivalent(list L, list K) |
---|
1421 | "USAGE: primDecsAreEquivalent(L,K); L,K list of ideals |
---|
1422 | RETURN: integer, 1 if the lists are the same up to ordering and redundant components; 0 otherwise |
---|
1423 | EXAMPLE: example primDecsAreEquivalent; shows an example" |
---|
1424 | { |
---|
1425 | L = removeRedundantComponents(L); |
---|
1426 | K = removeRedundantComponents(K); |
---|
1427 | int s1 = size(L); |
---|
1428 | int s2 = size(K); |
---|
1429 | if(s1!=s2) |
---|
1430 | { |
---|
1431 | return(0); |
---|
1432 | } |
---|
1433 | list L1, K1; |
---|
1434 | int i,j,t; |
---|
1435 | list N; |
---|
1436 | for(i=1; i<=s1; i++) |
---|
1437 | { |
---|
1438 | L1[i]=std(L[i][2]); |
---|
1439 | K1[i]=std(K[i][2]); |
---|
1440 | } |
---|
1441 | for(i=1; i<=s1; i++) |
---|
1442 | { |
---|
1443 | for(j=1; j<=s1; j++) |
---|
1444 | { |
---|
1445 | if(isSubset(L1[i],K1[j])) |
---|
1446 | { |
---|
1447 | if(isSubset(K1[j],L1[i])) |
---|
1448 | { |
---|
1449 | for(t=1; t<=size(N); t++) |
---|
1450 | { |
---|
1451 | if(N[t]==j) |
---|
1452 | { |
---|
1453 | return(0); |
---|
1454 | } |
---|
1455 | } |
---|
1456 | N[size(N)+1]=j; |
---|
1457 | } |
---|
1458 | |
---|
1459 | } |
---|
1460 | } |
---|
1461 | } |
---|
1462 | if ( size(N) != s1 ) |
---|
1463 | { |
---|
1464 | return(0); |
---|
1465 | } |
---|
1466 | return(1); |
---|
1467 | } |
---|
1468 | |
---|
1469 | example |
---|
1470 | { |
---|
1471 | "EXAMPLE:"; |
---|
1472 | echo = 2; |
---|
1473 | ring r = 0,(x,y),dp; |
---|
1474 | ideal i = x2,xy; |
---|
1475 | list L1 = primdecGTZ(i); |
---|
1476 | list L2 = primdecEHV(i); |
---|
1477 | primDecsAreEquivalent(L1,L2); |
---|
1478 | } |
---|