///////////////////////////////////////////////////////////////////// version="$Id$"; category="Commutative Algebra"; info=" LIBRARY: EHV.lib PROCEDURES FOR PRIMARY DECOMPOSITION OF IDEALS AUTHORS: Kai Dehmann, dehmann@mathematik.uni-kl.de; OVERVIEW: Algorithms for primary decomposition and radical-computation based on the ideas of Eisenbud, Huneke, and Vasconcelos. PROCEDURES: equiMaxEHV(I); equidimensional part of I removeComponent(I,e); intersection of the primary components of I of dimension >= e AssOfDim(I,e); an ideal such that the associated primes are exactly the associated primes of I having dimension e equiRadEHV(I [,Strategy]); equidimensional radical of I radEHV(I [,Strategy]); radical of I IntAssOfDim1(I,e); intersection of the associated primes of I having dimension e IntAssOfDim2(I,e); another way of computing the intersection of the associated primes of I having dimension e decompEHV(I); decomposition of a zero-dimensional radical ideal I AssEHV(I [,Strategy]); associated primes of I minAssEHV(I [,Strategy]); minimal associated primes of I localize(I,P,l); the contraction of the ideal generated by I in the localization w.r.t P componentEHV(I,P,L [,Strategy]); a P-primary component for I primdecEHV(I [,Strategy]); a minimal primary decomposition of I compareLists(L, K); procedure for comparing the output of primary decomposition algorithms (checks if the computed associated primes coincide) "; LIB "ring.lib"; LIB "general.lib"; LIB "elim.lib"; LIB "poly.lib"; LIB "random.lib"; LIB "inout.lib"; LIB "matrix.lib"; LIB "algebra.lib"; LIB "normal.lib"; ///////////////////////////////////////////////////////////////////// // // // G E N E R A L A L G O R I T H M S // // // ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// static proc AnnExtEHV(int n,list re) "USAGE: AnnExtEHV(n,re); n integer, re resolution RETURN: ideal, the annihilator of Ext^n(R/I,R) with given resolution re of I" { if(printlevel > 2){"Entering AnnExtEHV.";} if(n < 0) { ideal ann = ideal(1); if(printlevel > 2){"Leaving AnnExtEHV.";} return(ann); } int l = size(re); if(n < l) { matrix f = transpose(re[n+1]); if(n == 0) { matrix g = 0*gen(ncols(f)); } else { matrix g = transpose(re[n]); } module k = syz(f); ideal ann = quotient1(g,k); if(printlevel > 2){"Leaving AnnExtEHV.";} return(ann); } if(n == l) { ideal ann = Ann(transpose(re[n])); if(printlevel > 2){"Leaving AnnExtEHV.";} return(ann); } ideal ann = ideal(1); if(printlevel > 2){"Leaving AnnExtEHV.";} return(ann); } ///////////////////////////////////////////////////////////////////// //static proc isSubset(ideal I,ideal J) "USAGE: isSubset(I,J); I, J ideals RETURN: integer, 1 if I is a subset of J and 0 otherwise NOTE: if J is not a standard basis the result may be wrong " { int s = size(I); for(int i=1; i<=s; i++) { if(reduce(I[i],J)!=0) { return(0); } } return(1); } ///////////////////////////////////////////////////////////////////// // // // T H E E Q U I D I M E N S I O N A L P A R T // // // ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// proc equiMaxEHV(ideal I) "USAGE: equiMaxEHV(I); I ideal RETURN: ideal, the equidimensional part of I. NOTE: Uses algorithm of Eisenbud, Huneke, and Vasconcelos. EXAMPLE: example equiMaxEHV; shows an example " { if(printlevel > 2){"Entering equiMaxEHV.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } ideal J = groebner(I); int cod = nvars(basering)-dim(J); //If I is the entire ring... if(cod > nvars(basering)) { //...then return the ideal generated by 1. return(ideal(1)); } //Compute a resolution of I. if(printlevel > 2){"Computing resolution.";} if(homog(I)==1) { list re = sres(J,cod+1); re = minres(re); } else { list re = mres(I,cod+1); } if(printlevel > 2){"Finished computing resolution.";} //Compute the annihilator of the cod-th EXT-module. ideal ann = AnnExtEHV(cod,re); attrib(ann,"isEquidimensional",1); if(printlevel > 2){"Leaving equiMaxEHV.";} return(ann); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); equiMaxEHV(I); } ///////////////////////////////////////////////////////////////////// proc removeComponent(ideal I, int e) "USAGE: removeComponent(I,e); I ideal, e integer RETURN: ideal, the intersection of the primary components of I of dimension >= e EXAMPLE: example removeComponent; shows an example" { if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } ideal J = groebner(I); //Compute a resolution of I if(homog(I)==1) { list re = sres(J,0); re = minres(re); } else { list re = mres(I,0); } int f = nvars(basering); int cod; ideal ann; int g = nvars(basering) - e; while(f > g) { ann = AnnExtEHV(f,re); cod = nvars(basering) - dim(groebner(ann)); if( cod == f ) { I = quotient(I,ann); } f = f-1; } return(I); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); removeComponent(I,1); } ///////////////////////////////////////////////////////////////////// proc AssOfDim(ideal I, int e) "USAGE: AssOfDim(I,e); I ideal, e integer RETURN: ideal, such that the associated primes are exactly the associated primes of I having dimension e EXAMPLE: example AssOfDim; shows an example" { if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } int g = nvars(basering) - e; //Compute a resolution of I. ideal J = std(I); if(homog(I)==1) { list re = sres(J,g+1); re = minres(re); } else { list re = mres(I,g+1); } ideal ann = AnnExtEHV(g,re); int cod = nvars(basering) - dim(std(ann)); //If the codimension of I_g:=Ann(Ext^g(R/I,R)) equals g... if(cod == g) { //...then return the equidimensional part of I_g... ann = equiMaxEHV(ann); attrib(ann,"isEquidimensional",1); return(ann); } //...otherwise... else { //...I has no associated primes of dimension e. return(ideal(1)); } } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); AssOfDim(I,1); } ///////////////////////////////////////////////////////////////////// // // // T H E R A D I C A L // // // ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// static proc aJacob(ideal I, int a) "USAGE: aJacob(I,a); I ideal, a integer RETURN: ideal, the ath-Jacobian ideal of I" { matrix M = jacob(I); int n = nvars(basering); if(n-a <= 0) { return(ideal(1)); } if(n-a > nrows(M) or n-a > ncols(M)) { return(ideal(0)); } ideal J = minor(M,n-a); return(J); } ///////////////////////////////////////////////////////////////////// proc equiRadEHV(ideal I, list #) "USAGE: equiRadEHV(I [,Strategy]); I ideal, Strategy list RETURN: ideal, the equidimensional radical of I, i.e. the intersection of the minimal associated primes of I having the same dimension as I NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos, Works only in characteristic 0 or p large. The (optional) second argument determines the strategy used: Strategy[1] > strategy for the equidimensional part = 0 : uses equiMaxEHV = 1 : uses equidimMax Strategy[2] > strategy for the radical = 0 : combination of strategy 1 and 2 = 1 : computation of the radical just with the help of regular sequences = 2 : does not try to find a regular sequence Strategy[3] > strategy for the computation of ideal quotients = n : uses quot(.,.,n) for the ideal quotient computations If no second argument is given then Strategy=(0,0,0) is used. EXAMPLE: example equiRadEHV; shows an example" { if(printlevel > 2){"Entering equiRadEHV.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } if((char(basering)<100)&&(char(basering)!=0)) { "WARNING: The characteristic is too small, the result may be wrong"; } //Define the Strategy to be used. if(size(#) > 0) { if(#[1]!=1) { int equStr = 0; } else { int equStr = 1; } if(size(#) > 1) { if(#[2]!=1 and #[2]!=2) { int strategy = 0; } else { int strategy = #[2]; } if(size(#) > 2) { int quoStr = #[3]; } else { int quoStr = 0; } } else { int strategy = 0; int quoStr = 0; } } else { int equStr = 0; int strategy = 0; int quoStr = 0; } ideal J,I0,radI0,L,radI1,I2,radI2; int l,n; intvec op = option(get); matrix M; option(redSB); list m = mstd(I); option(set,op); int d = dim(m[1]); if(d==-1) { return(ideal(1)); } if(strategy != 2) { ///////////////////////////////////////////// // Computing the equidimensional radical // // via regular sequenves // ///////////////////////////////////////////// if(printlevel > 2){"Trying to find a regular sequence.";} int cod = nvars(basering)-d; //Complete intersection case: if(cod==size(m[2])) { J = aJacob(m[2],d); if(printlevel > 2){"Leaving equiRadEHV.";} return(quot1(m[2],J,quoStr)); } //First codim elements of I are a complete intersection: for(l=1; l<=cod; l++) { I0[l] = m[2][l]; } n = dim(groebner(I0))+cod-nvars(basering); //Last codim elements of I are a complete intersection: if(n!=0) { for(l=1; l<=cod; l++) { I0[l] = m[2][size(m[2])-l+1]; } n = dim(groebner(I0))+cod-nvars(basering); } //Taking a generic linear combination of the input: if(n!=0) { M = transpose(sparsetriag(size(m[2]),cod,95,1)); I0 = ideal(M*transpose(m[2])); n = dim(groebner(I0))+cod-nvars(basering); } //Taking a more generic linear combination of the input: if(n!=0) { while(strategy == 1 and n!=0) { M = transpose(sparsetriag(size(m[2]),cod,0,100)); I0 = ideal(M*transpose(m[2])); n = dim(groebner(I0))+cod-nvars(basering); } } if(n==0) { J = aJacob(I0,d); if(printlevel > 2){"1st quotient.";} radI0 = quot1(I0,J,quoStr); if(printlevel > 2){"2nd quotient.";} L = quot1(radI0,m[2],quoStr); if(printlevel > 2){"3rd quotient.";} radI1 = quot1(radI0,L,quoStr); attrib(radI1,"isEquidimensional",1); attrib(radI1,"isRadical",1); if(printlevel > 2){"Leaving equiRadEHV.";} return(radI1); } } //////////////////////////////////////////////////// // Computing the equidimensional radical directly // //////////////////////////////////////////////////// if(printlevel > 2){"Computing the equidimensional radical directly";} //Compute the equidimensional part depending on the chosen strategy if(equStr == 0) { I = equiMaxEHV(I); } if(equStr == 1) { I = equidimMax(I); } int a = nvars(basering)-1; while(a > d) { if(printlevel > 2){"While-Loop: "+string(a);} J = aJacob(I,a); while(dim(groebner(J+I))==d) { if(printlevel > 2){"Quotient-Computation.";} I = quot1(I,J,quoStr); if(printlevel > 2){"Computing the a-th Jacobian";} J = aJacob(I,a); } a = a-1; } if(printlevel > 2){"We left While-Loop.";} if(printlevel > 2){"Computing the a-th Jacobian";} J = aJacob(I,d); if(printlevel > 2){"Quotient-Computation.";} I = quot1(I,J,quoStr); attrib(I,"isEquidimensional",1); attrib(I,"isRadical",1); if(printlevel > 2){"Leaving equiRadEHV.";} return(I); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; ideal pr= equiRadEHV(i); pr; } ///////////////////////////////////////////////////////////////////// proc radEHV(ideal I, list #) "USAGE: radEHV(I [,Strategy]); ideal I, Strategy list RETURN: ideal, the radical of I NOTE: uses the algorithm of Eisenbud/Huneke/Vasconcelos Works only in characteristic 0 or p large. The (optional) second argument determines the strategy used: Strategy[1] > strategy for the equidimensional part = 0 : uses equiMaxEHV = 1 : uses equidimMax Strategy[2] > strategy for the radical = 0 : combination of strategy 1 and 2 = 1 : computation of the radical just with the help of regular sequences = 2 : does not try to find a regular sequence Strategy[3] > strategy for the computation of ideal quotients = n : uses quot(.,.,n) for the ideal quotient computations If no second argument is given then Strategy=(0,0,0) is used. EXAMPLE: example radEHV; shows an example" { if(printlevel > 2){"Entering radEHV.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } //Compute the equidimensional radical J of I. ideal J = equiRadEHV(I,#); //If I is the entire ring... if(deg(J[1]) <= 0) { //...then return the ideal generated by 1... return(ideal(1)); } //...else remove the maximal dimensional components and //compute the radical K of the lower dimensional components ... ideal K = radEHV(sat(I,J)[1],#); //..and intersect it with J. K = intersect(J,K); attrib(K,"isRadical",1); return(K); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; ideal pr= radical(i); pr; } ///////////////////////////////////////////////////////////////////// proc IntAssOfDim1(ideal I, int e) "USAGE: IntAssOfDim1(I,e); I idea, e integer RETURN: ideal, the intersection of the associated primes of I having dimension e EXAMPLE: example IntAssOfDim1; shows an example" { if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } int g = nvars(basering) - e; //Compute a resolution of I. ideal J = groebner(I); if(homog(I)==1) { list re = sres(J,g+1); re = minres(re); } else { list re = mres(I,g+1); } ideal ann = AnnExtEHV(g,re); int cod = nvars(basering) - dim(groebner(ann)); //If the codimension of I_g:=Ann(Ext^g(R/I,I)) equals g... if(cod == g) { //...then return the equidimensional radical of I_g... ann = equiRadEHV(ann); attrib(ann,"isEquidimensional",1); attrib(ann,"isRadical",1); return(ann); } else { //...else I has no associated primes of dimension e. return(ideal(1)); } } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); IntAssOfDim1(I,1); } ///////////////////////////////////////////////////////////////////// proc IntAssOfDim2(ideal I, int e) "USAGE: IntAssOfDim2(I,e); I ideal, e integer RETURN: ideal, the intersection of the associated primes of I having dimension e EXAMPLE: example IntAssOfDim2; shows an example" { if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } ideal I1 = removeComponent(I,e); ideal I2 = removeComponent(I,e+1); ideal b = quotient(I1,I2); b = equiRadEHV(b); return(b); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); IntAssOfDim2(I,1); } ///////////////////////////////////////////////////////////////////// // // // 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 // // // ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// proc decompEHV(ideal I) "USAGE: decompEHV(I); I zero-dimensional radical ideal RETURN: list, the associated primes of I EXAMPLE: example decompEHV; shows an example" { if(printlevel > 2){"Entering decompEHV.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } int e,m,vd,nfact; list l,k,L; poly f,h; ideal J; def base = basering; if( attrib(I,"isSB")!=1 ) { I = groebner(I); } while(1) { //Choose a random polynomial f from R. e = random(0,100); f = sparsepoly(e); //Check if f lies not in I. if(reduce(f,I)!=0) { J = quotient1(I,f); //If f is a zerodivisor modulo I... if( isSubset(J,I) == 0 ) { //...then use recursion... if(printlevel > 2){"We found a zerodivisor -- recursion";} l = decompEHV(J) + decompEHV(I+f); return(l); } if(printlevel > 2){"We found a non-zero-divisor.";} //...else compute the vectorspace dimension vd of I and... vd = vdim(I); //...compute m minimal such that 1,f,f^2,...,f^m are linearly dependent. qring Q = I; poly g = fetch(base,f); k = algDependent(g); def R = k[2]; setring R; if(size(ker)!=1) { //Calculate a generator for ker. ker = mstd(ker)[2]; } poly g = ker[1]; m = deg(g); //If m and vd coincide... if(m==vd) { if(printlevel > 2){"We have a good candidate.";} //...then factorize g. if(printlevel > 2){"Factorizing.";} L = factorize(g,2); //returns non-constant factors and multiplicities nfact = size(L[1]); //If g is irreducible... if(nfact==1 and L[2][1]==1) { if(printlevel > 2){"The element is irreducible.";} setring base; //..then I is a maximal ideal... l[1] = I; kill R,Q; return(l); } //...else... else { if(printlevel > 2){"The element is not irreducible -- recursion.";} //...take a non-trivial factor g1 of g ... poly g1 = L[1][1]; //..and insert f in g1... execute("ring newR = (" + charstr(R) + "),(y(1)),(" + ordstr(R) + ");"); poly g2 = imap(R,g1); setring base; h = fetch(newR,g2); h = subst(h,var(1),f); //...and use recursion... kill R,Q,newR; l = l + decompEHV(quotient1(I,h)); l = l + decompEHV(I+h); return(l); } } setring base; kill R,Q; } } } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y),dp; ideal i = x2+y2-10,x2+xy+2y2-16; decompEHV(i); } ///////////////////////////////////////////////////////////////////// // // // A S S O C I A T E D P R I M E S // // // ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// static proc idempotent(ideal I) "USAGE: idempotent(I); ideal I (weighted) homogeneous radical ideal, I intersected K[x(1),...,x(k)] zero-dimensional where deg(x(i))=0 for all i <= k and deg(x(i))>0 for all i>k. RETURN: a list of ideals I(1),...,I(t) such that K[x]/I = K[x]/I(1) x ... x K[x]/I(t)" { if(printlevel > 2){"Entering idempotent.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } int n = nvars(basering); poly f,g; string j; int i,k,splits; list l; //Collect all variables of degree 0. for(i=1; i<=n; i++) { if(deg(var(i)) > 0) { f = f*var(i); } else { splits = 1; j = j + string(var(i)) + "," ; } } //If there are no variables of degree 0 //then there are no idempotents and we are done... if(splits == 0) { l[1]=I; return(l); } //...else compute J = I intersected K[x(1),...,x(k)]... ideal J = eliminate(I,f); def base = basering; j = j[1,size(j)-1]; j = "ring @r = (" + charstr(basering) + "),(" + j + "),(" + ordstr(basering) + ");"; execute(j); ideal J = imap(base,J); //...and compute the associated primes of the zero-dimensional ideal J. list L = decompEHV(J); int s = size(L); ideal K,Z; poly g; //For each associated prime ideal P_i of J... for(i=1; i<=s; i++) { K = ideal(1); //...comnpute the intersection K of the other associated prime ideals... for(k=1; k<=s; k++) { if(i!=k) { K = intersect(K,L[k]); } } //...and find an element that lies in K but not in P_i... g = randomid(K,1)[1]; Z = L[i]; Z = std(Z); while(reduce(g,Z)==0) { g = randomid(K,1)[1]; } setring base; g = imap(@r,g); //...and compute the corresponding ideal I(i) l[i] = quotient(I,g); setring @r; } setring base; return(l); } ///////////////////////////////////////////////////////////////////// static proc equiAssEHV(ideal I) "USAGE: equiAssEHV(I); I equidimensional, radical, and homogeneous ideal RETURN: a list, the associated prime ideals of I" { if(printlevel > 2){"Entering equiAssEHV.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } list L; def base = basering; int n = nvars(basering); int i,j; //Compute the normalization of I. if(printlevel > 2){"Entering Normalization.";} list norOut = normal(I, "noDeco"); list K = norOut[1]; if(printlevel > 2){"Leaving Normalisation.";} //The normalization algorithm returns k factors. int k = size(K); if(printlevel > 1){"Normalization algorithm splits ideal in " + string(k) + " factors.";} //Examine each factor of the normalization. def P; for(i=1; i<=k; i++) { P = K[i]; setring P; //Use procedure idempotent to split the i-th factor of //the normalization in a product of integral domains. if(printlevel > 1){"Examining " + string(i) +". factor.";} list l = idempotent(norid); int leng = size(l); if(printlevel > 1){"Idempotent algorithm splits factor " + string(i) + " in " + string(leng) + " factors.";} //Intersect the minimal primes corresponding to the //integral domains obtained from idempotent with the groundring, //i.e. compute the preimages w.r.t. the corresponding normalization map ideal J; for(j=1; j<=leng; j++) { J = l[j]; setring base; L[size(L)+1] = preimage(P,normap,J); setring P; } kill l, leng, J; setring base; } return(L); } ///////////////////////////////////////////////////////////////////// proc AssEHV(ideal I, list #) "USAGE: AssEHV(I [,Strategy]); I Ideal, Strategy list RETURN: a list, the associated prime ideals of I NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos. The (optional) second argument determines the strategy used: Strategy[1] > strategy for the equidimensional part = 0 : uses equiMaxEHV = 1 : uses equidimMax Strategy[2] > strategy for the equidimensional radical = 0 : uses equiRadEHV = 1 : uses equiRadical Strategy[3] > strategy for equiRadEHV = 0 : combination of strategy 1 and 2 = 1 : computation of the radical just with the help of regular sequences = 2 : does not try to find a regular sequence Strategy[4] > strategy for the computation of ideal quotients = n : uses quot(.,.,n) for the ideal quotient computations If no second argument is given, Strategy=(0,0,0,0) is used. EXAMPLE: example AssEHV; shows an example" { if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } if(printlevel > 2){"Entering AssEHV";} //Specify the strategy to be used. if(size(#)==0) { # = 0,0,0,0; } if(size(#)==1) { # = #[1],0,0,0; } if(size(#)==2) { # = #[1],#[2],0,0; } if(size(#)==3) { # = #[1],#[2],#[3],0; } list L; ideal K; def base = basering; int m,j; ideal J = groebner(I); int n = nvars(basering); int d = dim(J); //Compute a resolution of I. if(printlevel > 2){"Computing resolution.";} if(homog(I)==1) { list re = sres(J,0); re = minres(re); } else { list re = mres(I,0); } ideal ann; int cod; //For 0<=i<= dim(I) compute the intersection of the i-dimensional associated primes of I for(int i=0; i<=d; i++) { if(printlevel > 1){"Are there components of dimension " + string(i) + "?";} ann = AnnExtEHV(n-i,re); cod = n - dim(groebner(ann)); //If there are associated primes of dimension i... if(cod == n-i) { if(printlevel > 1){"Yes. There are components of dimension " + string(i) + ".";} //...then compute the intersection K of all associated primes of I of dimension i if(#[2]==0) { if(size(#) > 3) { K = equiRadEHV(ann,#[1],#[3],#[4]); } if(size(#) > 2) { K = equiRadEHV(ann,#[1],#[3]); } else { K = equiRadEHV(ann,#[1]); } } if(#[2]==1) { K = equiRadical(ann); } attrib(K,"isEquidimensional",1); attrib(K,"isRadical",1); //If K is already homogeneous then use equiAssEHV to recover the associated primes of K... if(homog(K)==1) { if(printlevel > 2){"Input already homogeneous.";} L = L + equiAssEHV(K); } //...else... else { //...homogenize K w.r.t. t,... if(printlevel > 2){"Input not homogeneous; must homogenize.";} def homoR=changeord(list(list("dp",1:nvars(basering)))); setring homoR; ideal homoJ = fetch(base,K); homoJ = groebner(homoJ); execute("ring newR = (" + charstr(base) + "),(x(1..n),t),dp;"); ideal homoK = fetch(homoR,homoJ); homoK = homog(homoK,t); attrib(homoK,"isEquidimensional",1); attrib(homoK,"isRadical",1); //...compute the associated primes of the homogenization using equiAssEHV,... list l = equiAssEHV(homoK); //...and set t=1 in the generators of the associated primes just computed. ideal Z; for(j=1; j<=size(l); j++) { Z = subst(l[j],t,1); setring base; L[size(L)+1] = fetch(newR,Z); setring newR; } setring base; kill homoR; kill newR; } } else { if(printlevel > 1){"No. There are no components of dimension " + string(i) + ".";} } } return(L); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; list pr = AssEHV(i); pr; } ///////////////////////////////////////////////////////////////////// proc minAssEHV(ideal I, list #) "USAGE: minAssEHV(I [,Strategy]); I ideal, Strategy list RETURN: a list, the minimal associated prime ideals of I NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos. The (optional) second argument determines the strategy used: Strategy[1] > strategy for the equidimensional part = 0 : uses equiMaxEHV = 1 : uses equidimMax Strategy[2] > strategy for the equidimensional radical = 0 : uses equiRadEHV, resp. radicalEHV = 1 : uses equiRadical, resp. radical Strategy[3] > strategy for equiRadEHV = 0 : combination of strategy 1 and 2 = 1 : computation of the radical just with the help of regular sequences = 2 : does not try to find a regular sequence Strategy[4] > strategy for the computation of ideal quotients = n : uses quot(.,.,n) for the ideal quotient computations If no second argument is given, Strategy=(0,0,0,0) is used. EXAMPLE: example minAssEHV; shows an example" { if(attrib(basering,"global")==0) {ERROR("// Not implemented for this ordering, please change to global ordering.");} //Specify the strategy to be used. if(size(#)==0) { # = 0,0,0,0; } if(size(#)==1) { # = #[1],0,0,0; } if(size(#)==2) { # = #[1],#[2],0,0; } if(size(#)==3) { # = #[1],#[2],#[3],0; } //Compute the radical of I. if(#[2]==0) { I = radEHV(I,#); } if(#[2]==1) { I = radical(I); } //Compute the associated primes of the radical. return(AssEHV(I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; list pr = minAssEHV(i); pr; } ///////////////////////////////////////////////////////////////////// // // // P R I M A R Y D E C O M P O S I T I O N // // // ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// proc localize(ideal I, ideal P, list l) "USAGE: localize(I,P,l); I ideal, P an associated prime ideal of I, l list of all associated primes of I RETURN: ideal, the contraction of the ideal generated by I in the localization w.r.t P EXAMPLE: example localize; shows an example" { if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } ideal Intersection = ideal(1); int s = size(l); if(attrib(P,"isSB")!=1) { P = groebner(P); } //Compute the intersection of all associated primes of I that are not contained in P. for(int i=1; i<=s; i++) { if(isSubset(l[i],P)!=1) { Intersection = intersect(Intersection,l[i]); } } Intersection = groebner(Intersection); //If the intersection is the entire ring... if(reduce(1,Intersection)==0) { //...then return I... return(I); } //...else try to find an element f that lies in the intersection but outside P... poly f = 0; while(reduce(f,P) == 0) { f = randomid(Intersection,1)[1]; } //...and saturate I w.r.t. f. I = sat(I,f)[1]; return(I); } ///////////////////////////////////////////////////////////////////// proc componentEHV(ideal I, ideal P, list L, list #) "USAGE: componentEHV(I,P,L [,Strategy]); I ideal, P associated prime of I, L list of all associated primes of I, Strategy list RETURN: ideal, a P-primary component Q for I NOTE: The (optional) second argument determines the strategy used: Strategy[1] > strategy for equidimensional part = 0 : uses equiMaxEHV = 1 : uses equidimMax If no second argument is given then Strategy=0 is used. EXAMPLE: example componentEHV; shows an example" { if(printlevel > 2){"Entering componentEHV.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } //If no strategy is specified use standard strategy. if(size(#)==0) { # = 0; } ideal T = P; ideal Q; //Compute the localization of I at P... ideal IP = groebner(localize(I,P,L)); //...and compute the saturation of the localization w.r.t. P. ideal IP2 = sat(IP,P)[1]; //As long as we have not found a primary component... int isPrimaryComponent = 0; while(isPrimaryComponent!=1) { //...compute the equidimensional part Q of I+P^n... if(#[1]==0) { Q = equiMaxEHV(I+T); } if(#[1]==1) { Q = equidimMax(I+T); } //...and check if it is a primary component for P. if(isSubset(intersect(IP2,Q),IP)==1) { isPrimaryComponent = 1; } else { T = T*P; } } if(printlevel > 2){"Leaving componentEHV.";} return(Q); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; list pr = AssEHV(i); componentEHV(i,pr[1],pr); } ///////////////////////////////////////////////////////////////////// proc primdecEHV(ideal I, list #) "USAGE: primdecEHV(I [,Strategy]); I ideal, Strategy list RETURN: a list pr of primary ideals and their associated primes: pr[i][1] the i-th primary component, pr[i][2] the i-th prime component. NOTE: Algorithm of Eisenbud/Huneke/Vasconcelos. The (optional) second argument determines the strategy used: Strategy[1] > strategy for equidimensional part = 0 : uses equiMaxEHV = 1 : uses equidimMax Strategy[2] > strategy for equidimensional radical = 0 : uses equiRadEHV, resp. radicalEHV = 1 : uses equiRadical, resp. radical Strategy[3] > strategy for equiRadEHV = 0 : combination of strategy 1 and 2 = 1 : computation of the radical just with the help of regular sequences = 2 : does not try to find a regular sequence Strategy[4] > strategy for the computation of ideal quotients = n : uses quot(.,.,n) for the ideal quotient computations If no second argument is given then Strategy=(0,0,0,0) is used. EXAMPLE: example primdecEHV; shows an example" { if(printlevel > 2){"Entering primdecEHV.";} if(attrib(basering,"global")==0) { ERROR("// Not implemented for this ordering, please change to global ordering."); } list L,K; //Specify the strategy to be used. if(size(#)==0) { # = 0,0,0,0; } if(size(#)==1) { # = #[1],0,0,0; } if(size(#)==2) { # = #[1],#[2],0,0; } if(size(#)==3) { # = #[1],#[2],#[3],0; } //Compute the associated primes of I... L = AssEHV(I,#); if(printlevel > 0){"We have " + string(size(L)) + " prime components.";} //...and compute for each associated prime of I a corresponding primary component. int l = size(L); for(int i=1; i<=l; i++) { K[i] = list(); K[i][2] = L[i]; K[i][1] = componentEHV(I,L[i],L,#[1]); } if(printlevel > 2){"Leaving primdecEHV.";} return(K); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; list pr = primdecEHV(i); pr; } ///////////////////////////////////////////////////////////////////// // // // A L G O R I T H M S F O R T E S T I N G // // // ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// proc compareLists(list L, list K) "USAGE: checkLists(L,K); L,K list of ideals RETURN: integer, 1 if the lists are the same up to ordering and 0 otherwise EXAMPLE: example checkLists; shows an example" { int s1 = size(L); int s2 = size(K); if(s1!=s2) { return(0); } list L1, K1; int i,j,t; list N; for(i=1; i<=s1; i++) { L1[i]=std(L[i][2]); K1[i]=std(K[i][2]); } for(i=1; i<=s1; i++) { for(j=1; j<=s1; j++) { if(isSubset(L1[i],K1[j])) { if(isSubset(K1[j],L1[i])) { for(t=1; t<=size(N); t++) { if(N[t]==j) { return(0); } } N[size(N)+1]=j; } } } } return(1); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; ideal i = x2,xy; list L1 = primdecGTZ(i); list L2 = primdecEHV(i); compareLists(L1,L2); }