//////////////////////////////////////////////////////////////////////////// version="version primdec.lib 4.1.1.0 Dec_2017 "; // $Id$ category="Commutative Algebra"; info=" LIBRARY: primdec.lib Primary Decomposition and Radical of Ideals AUTHORS: Gerhard Pfister, pfister@mathematik.uni-kl.de (GTZ)@* Wolfram Decker, decker@math.uni-sb.de (SY)@* Hans Schoenemann, hannes@mathematik.uni-kl.de (SY)@* Santiago Laplagne, slaplagn@dm.uba.ar (GTZ) OVERVIEW: Algorithms for primary decomposition based on the ideas of Gianni, Trager and Zacharias (implementation by Gerhard Pfister), respectively based on the ideas of Shimoyama and Yokoyama (implementation by Wolfram Decker and Hans Schoenemann).@* The procedures are implemented to be used in characteristic 0.@* They also work in positive characteristic >> 0.@* In small characteristic and for algebraic extensions, primdecGTZ may not terminate.@* Algorithms for the computation of the radical based on the ideas of Krick, Logar, Laplagne and Kemper (implementation by Gerhard Pfister and Santiago Laplagne). They work in any characteristic.@* Baserings must have a global ordering and no quotient ideal. Exceptions: primdecGTZ, absPrimdecGTZ, minAssGTZ, primdecSY, minAssChar, radical accept non-global ordering. PROCEDURES: Ann(M); annihilator of R^n/M, R=basering, M in R^n primdecGTZ(I); complete primary decomposition via Gianni,Trager,Zacharias primdecGTZE(I); complete primary decomposition via Gianni,Trager,Zacharias. Returns empty list for the unit ideal primdecSY(I...); complete primary decomposition via Shimoyama-Yokoyama primdecSYE(I,..); complete primary decomposition via Shimoyama-Yokoyama. Returns empty list for the unit ideal minAssGTZ(I); the minimal associated primes via Gianni,Trager,Zacharias (with modifications by Laplagne) minAssGTZE(I); the minimal associated primes via Gianni,Trager,Zacharias. Returns empty list for unit ideal minAssChar(I...); the minimal associated primes using characteristic sets minAssCharE(I..); the minimal associated primes using characteristic sets. Returns empty list for unit ideal testPrimary(L,k); tests the result of the primary decomposition testPrimaryE(L,k); tests the result of the primary decomposition. Handles also empty list L. radical(I); computes the radical of I via Krick/Logar (with modifications by Laplagne) and Kemper radicalEHV(I); computes the radical of I via Eisenbud,Huneke,Vasconcelos equiRadical(I); the radical of the equidimensional part of the ideal I prepareAss(I); list of radicals of the equidimensional components of I equidim(I); weak equidimensional decomposition of I equidimMax(I); equidimensional locus of I equidimMaxEHV(I); equidimensional locus of I via Eisenbud,Huneke,Vasconcelos zerodec(I); zerodimensional decomposition via Monico absPrimdecGTZ(I); the absolute prime components of I absPrimdecGTZE(I); the absolute prime components of I. Assumes I is not unit ideal. sep(f,k); the separabel part of f as polynomial in Fp(t1,...,tm) SEE ALSO: primdecint_lib KEYWORDS: primary decomposition "; LIB "general.lib"; LIB "elim.lib"; LIB "poly.lib"; LIB "random.lib"; LIB "inout.lib"; LIB "matrix.lib"; LIB "triang.lib"; LIB "absfact.lib"; LIB "ring.lib"; /////////////////////////////////////////////////////////////////////////////// // // Gianni/Trager/Zacharias // /////////////////////////////////////////////////////////////////////////////// static proc sat1 (ideal id, poly p) "USAGE: sat1(id,j); id ideal, j polynomial RETURN: saturation of id with respect to j (= union_(k=1...) of id:j^k) NOTE: result is a std basis in the basering " { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int @k; ideal inew=std(id); ideal iold; intvec op=option(get); option(returnSB); while(specialIdealsEqual(iold,inew)==0 ) { iold=inew; inew=quotient(iold,p); @k++; } @k--; option(set,op); list L =inew,p^@k; return (L); } /////////////////////////////////////////////////////////////////////////////// static proc sat2 (ideal id, ideal h) "USAGE: sat2(id,j); id ideal, j polynomial RETURN: saturation of id with respect to j (= union_(k=1...) of id:j^k) NOTE: result is a std basis in the basering " { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int @k,@i; def @P= basering; if(ordstr(basering)[1,2]!="dp") { def @Phelp=changeordTo(basering,"dp"); setring @Phelp; ideal inew=std(imap(@P,id)); ideal @h=imap(@P,h); } else { ideal @h=h; ideal inew=std(id); } ideal fac; for(@i=1;@i<=ncols(@h);@i++) { if(deg(@h[@i])>0) { fac=fac+factorize(@h[@i],1); } } fac=simplify(fac,6); poly @f=1; if(deg(fac[1])>0) { ideal iold; for(@i=1;@i<=size(fac);@i++) { @f=@f*fac[@i]; } intvec op = option(get); option(returnSB); while(specialIdealsEqual(iold,inew)==0 ) { iold=inew; if(deg(iold[size(iold)])!=1) { inew=quotient(iold,@f); } else { inew=iold; } @k++; } option(set,op); @k--; } if(ordstr(@P)[1,2]!="dp") { setring @P; ideal inew=std(imap(@Phelp,inew)); poly @f=imap(@Phelp,@f); } list L =inew,@f^@k; return (L); } /////////////////////////////////////////////////////////////////////////////// proc minSat(ideal inew, ideal h) { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; int i,k; poly f=1; ideal iold,fac; list quotM,l; for(i=1;i<=ncols(h);i++) { if(deg(h[i])>0) { fac=fac+factorize(h[i],1); } } fac=simplify(fac,6); if(size(fac)==0) { l=inew,1; return(l); } fac=sort(fac)[1]; for(i=1;i<=size(fac);i++) { f=f*fac[i]; } quotM[1]=inew; quotM[2]=fac; quotM[3]=f; f=1; intvec op = option(get); option(returnSB); while(specialIdealsEqual(iold,quotM[1])==0) { if(k>0) { f=f*quotM[3]; } iold=quotM[1]; quotM=quotMin(quotM); k++; } option(set,op); l=quotM[1],f; return(l); } static proc quotMin(list tsil) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int i,j,k,action; ideal verg; list l; poly g; ideal laedi=tsil[1]; ideal fac=tsil[2]; poly f=tsil[3]; ideal star=quotient(laedi,f); if(specialIdealsEqual(star,laedi)) { l=star,fac,f; return(l); } action=1; while(action==1) { if(size(fac)==1) { action=0; break; } for(i=1;i<=size(fac);i++) { g=1; verg=laedi; for(j=1;j<=size(fac);j++) { if(i!=j) { g=g*fac[j]; } } verg=quotient(laedi,g); if(specialIdealsEqual(verg,star)==1) { f=g; fac[i]=0; fac=simplify(fac,2); break; } if(i==size(fac)) { action=0; } } } l=star,fac,f; return(l); } /////////////////////////////////////////////////////////////////////////////// static proc testFactor(list act,poly p) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; poly keep=p; int i; poly q=act[1][1]^act[2][1]; for(i=2;i<=size(act[1]);i++) { q=q*act[1][i]^act[2][i]; } q=1/leadcoef(q)*q; p=1/leadcoef(p)*p; if(p-q!=0) { "ERROR IN FACTOR, please inform the authors"; } } /////////////////////////////////////////////////////////////////////////////// static proc factor(poly p) "USAGE: factor(p) p poly RETURN: list=; NOTE: EXAMPLE: example factor; shows an example " { ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; ideal @i; list @l; intvec @v,@w; int @j,@k,@n; @l=factorize(p); for(@j=1;@j<=size(@l[1]);@j++) { if(leadcoef(@l[1][@j])==@l[1][@j]) { @n++; } } if(@n>0) { if(@n==size(@l[1])) { @l[1]=ideal(1); @v=1; @l[2]=@v; } else { @k=0; int pleh; for(@j=1;@j<=size(@l[1]);@j++) { if(leadcoef(@l[1][@j])!=@l[1][@j]) { @k++; @i=@i+ideal(@l[1][@j]); if(size(@i)==pleh) { "//factorization error"; @l; @k--; @v[@k]=@v[@k]+@l[2][@j]; } else { pleh++; @v[@k]=@l[2][@j]; } } } @l[1]=@i; @l[2]=@v; } } // } return(@l); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = (x+y)^2*(y-z)^3; list l = factor(p); l; ring r1 =(0,b,d,f,g),(a,c,e),lp; poly p =(1*d)*e^2+(1*d*f^2*g); list l = factor(p); l; ring r2 =(0,b,f,g),(a,c,e,d),lp; poly p =(1*d)*e^2+(1*d*f^2*g); list l = factor(p); l; } /////////////////////////////////////////////////////////////////////////////// proc idealsEqual( ideal k, ideal j) { return(stdIdealsEqual(std(k),std(j))); } static proc specialIdealsEqual( ideal k1, ideal k2) { int j; if(size(k1)==size(k2)) { for(j=1;j<=size(k1);j++) { if(leadexp(k1[j])!=leadexp(k2[j])) { return(0); } } return(1); } return(0); } static proc stdIdealsEqual( ideal k1, ideal k2) { int j; if(size(k1)==size(k2)) { for(j=1;j<=size(k1);j++) { if(leadexp(k1[j])!=leadexp(k2[j])) { return(0); } } attrib(k2,"isSB",1); if(size(reduce(k1,k2,5))==0) { return(1); } } return(0); } /////////////////////////////////////////////////////////////////////////////// proc primaryTest (ideal i, poly p) { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; if(i[1]==1){return(ideal(1));} int m=1; int n=nvars(basering); int e,f; poly t; ideal h; list act; ideal prm=p; attrib(prm,"isSB",1); while (n>1) { n--; m++; //search for i[m] which has a power of var(n) as leading term if (n==1) { m=size(i); } else { while (lead(i[m])/var(n-1)==0) { m++; } m--; } //check whether i[m] =(c*var(n)+h)^e modulo prm for some //h in K[var(n+1),...,var(nvars(basering))], c in K //if not (0) is returned, else var(n)+h is added to prm e=deg(lead(i[m])); if(char(basering)!=0) { f=1; if(e mod char(basering)==0) { if ( voice >=16 ) { "// WARNING: The characteristic is perhaps too small to use"; "// the algorithm of Gianni/Trager/Zacharias."; "// This may result in an infinte loop"; "// loop in primaryTest, voice:",voice;""; } while(e mod char(basering)==0) { f=f*char(basering); e=e div char(basering); } } t=leadcoef(i[m])*e*var(n)^f+(i[m]-lead(i[m]))/var(n)^((e-1)*f); i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m]; if (reduce(i[m]-t^e,prm,5) !=0) { return(ideal(0)); } if(f>1) { act=factorize(t); if(size(act[1])>2) { return(ideal(0)); } if(deg(act[1][2])>1) { return(ideal(0)); } t=act[1][2]; } } else { t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1); i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m]; if (reduce(i[m]-t^e,prm,5) !=0) { return(ideal(0)); } } h=interred(t); t=h[1]; prm = prm,t; attrib(prm,"isSB",1); } return(prm); } /////////////////////////////////////////////////////////////////////////////// proc gcdTest(ideal act) { ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; int i,j; if(size(act)<=1) { return(0); } for (i=1;i0) { return(0); } } } return(1); } /////////////////////////////////////////////////////////////////////////////// static proc splitPrimary(list l,ideal ser,int @wr,list sact) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int i,j,k,s,r,w; list keepresult,act,keepprime; poly @f; int sl=size(l); for(i=sl div 2;i>=1;i--) { if(sact[2][i]>1) { keepprime[i]=l[2*i-1]+ideal(sact[1][i]); } else { keepprime[i]=l[2*i-1]; } } i=0; while(i0)&&(size(reduce(ser,l[2*i-1],5))==0)) { l[2*i-1]=ideal(1); l[2*i]=ideal(1); continue; } if(size(l[2*i])==0) { if(homog(l[2*i-1])==1) { l[2*i]=maxideal(1); continue; } j=0; /* if(i<=sl div 2) { j=1; } */ while(j1)) { keepprime[i]=interred(keepprime[i]+ideal(act[1][1])); if(homog(keepprime[i])==1) { l[2*i]=maxideal(1); break; } } if(gcdTest(act[1])==1) { for(k=2;k<=r;k++) { keepprime[size(l) div 2+k-1]=interred(keepprime[i]+ideal(act[1][k])); } keepprime[i]=interred(keepprime[i]+ideal(act[1][1])); for(k=1;k<=r;k++) { if(@wr==0) { keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]); } else { keepresult[k]=std(l[2*i-1],act[1][k]); } } l[2*i-1]=keepresult[1]; if(vdim(keepresult[1])==deg(act[1][1])) { l[2*i]=keepresult[1]; } if((homog(keepresult[1])==1)||(homog(keepprime[i])==1)) { l[2*i]=maxideal(1); } s=size(l)-2; for(k=2;k<=r;k++) { l[s+2*k-1]=keepresult[k]; keepprime[s div 2+k]=interred(keepresult[k]+ideal(act[1][k])); if(vdim(keepresult[k])==deg(act[1][k])) { l[s+2*k]=keepresult[k]; } else { l[s+2*k]=ideal(0); } if((homog(keepresult[k])==1)||(homog(keepprime[s div 2+k])==1)) { l[s+2*k]=maxideal(1); } } i--; break; } if(r>=2) { s=size(l); @f=act[1][1]; act=sat1(l[2*i-1],act[1][1]); if(deg(act[1][1])>0) { l[s+2]=ideal(0); l[s+1]=std(l[2*i-1],act[2]); if(homog(l[s+1])==1) { l[s+2]=maxideal(1); } //else { l[s+2]=ideal(0); } keepprime[s div 2+1]=interred(keepprime[i]+ideal(@f)); if(homog(keepprime[s div 2+1])==1) { l[s+2]=maxideal(1); } keepprime[i]=act[1]; l[2*i-1]=act[1]; attrib(l[2*i-1],"isSB",1); if(homog(l[2*i-1])==1) { l[2*i]=maxideal(1); } i--; break; } } } } } if(sl==size(l)) { return(l); } for(i=1;i<=size(l) div 2;i++) { attrib(l[2*i-1],"isSB",1); if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],5))==0)&&(deg(l[2*i-1][1])>0)) { "Achtung in split"; l[2*i-1]=ideal(1); l[2*i]=ideal(1); } if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1)) { keepprime[i]=std(keepprime[i]); if(homog(keepprime[i])==1) { l[2*i]=maxideal(1); } else { act=zero_decomp(keepprime[i],ideal(0),@wr,1); if(size(act)==2) { l[2*i]=act[2]; } } } } return(l); } example { "EXAMPLE:"; echo=2; ring r = 32003,(x,y,z),lp; ideal i1=x*(x+1),yz,(z+1)*(z-1); ideal i2=xy,yz,(x-2)*(x+3); list l=i1,ideal(0),i2,ideal(0),i2,ideal(1); list l1=splitPrimary(l,ideal(0),0); l1; } /////////////////////////////////////////////////////////////////////////////// static proc splitCharp(list l) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; if((char(basering)==0)||(npars(basering)>0)) { return(l); } def op = option(get); def P=basering; int i,j,k,m,q,d,o; int n = nvars(basering); ideal s,t,u,sact; poly ni; string minp,gnir,va; list sa,keep,rp,keep1; for(i=1;i<=size(l) div 2;i++) { if(size(l[2*i])==0) { if(deg(l[2*i-1][1])==vdim(l[2*i-1])) { l[2*i]=l[2*i-1]; } } } for(i=1;i<=size(l) div 2;i++) { if(size(l[2*i])==0) { s=factorize(l[2*i-1][1],1); //vermeiden!!! t=l[2*i-1]; m=size(t); ni=s[1]; if(deg(ni)>1) { va=varstr(P); j=size(va); while(va[j]!=","){j--;} va=va[1..j-1]; gnir="ring RL=("+string(char(P))+","+string(var(n))+"),("+va+"),lp;"; execute(gnir); minpoly=leadcoef(imap(P,ni)); ideal act; ideal t=imap(P,t); for(k=2;k<=m;k++) { act=factorize(t[k],1); if(size(act)>1){break;} } setring P; sact=imap(RL,act); if(size(sact)>1) { sa=sat1(l[2*i-1],sact[1]); keep[size(keep)+1]=std(l[2*i-1],sa[2]); if(sa[1][1]==l[2*i-1][1]) { l[2*i-1]=std(sa[1]); l[2*i]=primaryTest(sa[1],s[1]); } else { l[2*i-1]=std(sa[1]); l[2*i]=primaryTest(sa[1],factorize(sa[1][1],1)[1]); } } if((size(sact)==1)&&(m==2)) { l[2*i]=std(l[2*i-1],sact[1]); } if((size(sact)==1)&&(m>2)) { setring RL; option(redSB); t=std(t); list sp=zero_decomp(t,0,0); setring P; rp=imap(RL,sp); for(o=1;o<=size(rp);o++) { rp[o]=interred(simplify(rp[o],1)+ideal(ni)); } l[2*i-1]=rp[1]; l[2*i]=rp[2]; rp=delete(rp,1); rp=delete(rp,1); keep1=keep1+rp; option(set,op); } kill RL; } } } if(size(keep)>0) { for(i=1;i<=size(keep);i++) { if(deg(keep[i][1])>0) { l[size(l)+1]=keep[i]; l[size(l)+1]=primaryTest(keep[i],factorize(keep[i][1],1)[1]); } } } l=l+keep1; option(set,op); return(l); } /////////////////////////////////////////////////////////////////////////////// proc zero_decomp (ideal j,ideal ser,int @wr,list #) "USAGE: zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1 (@wr=0 for primary decomposition, @wr=1 for computation of associated primes) RETURN: list = list of primary ideals and their radicals (at even positions in the list) if the input is zero-dimensional and a standardbases with respect to lex-ordering If ser!=(0) and ser is contained in j or if j is not zero-dimen- sional then ideal(1),ideal(1) is returned NOTE: Algorithm of Gianni/Trager/Zacharias EXAMPLE: example zero_decomp; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; def @P = basering; int uytrewq; int nva = nvars(basering); int @k,@s,@n,@k1,@zz; list primary,lres0,lres1,act,@lh,@wh; map phi,psi,phi1,psi1; ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1; intvec @vh,@hilb; string @ri; poly @f; if (dim(j)>0) { primary[1]=ideal(1); primary[2]=ideal(1); return(primary); } intvec save=option(get); option(redSB); j=interred(j); attrib(j,"isSB",1); if(vdim(j)==deg(j[1])) { act=factor(j[1]); for(@k=1;@k<=size(act[1]);@k++) { @qh=j; if(@wr==0) { @qh[1]=act[1][@k]^act[2][@k]; } else { @qh[1]=act[1][@k]; } primary[2*@k-1]=interred(@qh); @qh=j; @qh[1]=act[1][@k]; primary[2*@k]=interred(@qh); attrib( primary[2*@k-1],"isSB",1); if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],5))==0)) { primary[2*@k-1]=ideal(1); primary[2*@k]=ideal(1); } } option(set,save); return(primary); } option(set,save); if(homog(j)==1) { primary[1]=j; if((size(ser)>0)&&(size(reduce(ser,j,5))==0)) { primary[1]=ideal(1); primary[2]=ideal(1); return(primary); } if(dim(j)==-1) { primary[1]=ideal(1); primary[2]=ideal(1); } else { primary[2]=maxideal(1); } return(primary); } //the first element in the standardbase is factorized if(deg(j[1])>0) { act=factor(j[1]); testFactor(act,j[1]); } else { primary[1]=ideal(1); primary[2]=ideal(1); return(primary); } //with the factors new ideals (hopefully the primary decomposition) //are created if(size(act[1])>1) { if(size(#)>1) { primary[1]=ideal(1); primary[2]=ideal(1); primary[3]=ideal(1); primary[4]=ideal(1); return(primary); } for(@k=1;@k<=size(act[1]);@k++) { if(@wr==0) { primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]); } else { primary[2*@k-1]=std(j,act[1][@k]); } if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k]))) { primary[2*@k] = primary[2*@k-1]; } else { primary[2*@k] = primaryTest(primary[2*@k-1],act[1][@k]); } } } else { primary[1]=j; if((size(#)>0)&&(act[2][1]>1)) { act[2]=1; primary[1]=std(primary[1],act[1][1]); } if(@wr!=0) { primary[1]=std(j,act[1][1]); } if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1]))) { primary[2]=primary[1]; } else { primary[2]=primaryTest(primary[1],act[1][1]); } } if(size(#)==0) { primary=splitPrimary(primary,ser,@wr,act); } if((voice>=7)&&(char(basering)<=181)&&(@wr==1)) { primary=splitCharp(primary); } /* if((@wr==2)&&(npars(basering)>0)&&(voice>=7)&&(char(basering)>0)) { //the prime decomposition of Yokoyama in characteristic p list ke,ek; @k=0; while(@k=9){primary=extF(primary);}; //test whether all ideals in the decomposition are primary and //in general position //if not after a random coordinate transformation of the last //variable the corresponding ideal is decomposed again. if((npars(basering)>0)&&(voice>=9)) { poly randp; for(@zz=1;@zz0)&&(voice>=11)) { jmap[size(jmap)]=randp; } for(@n=2;@n<=size(primary[2*@k-1]);@n++) { if(deg(lead(primary[2*@k-1][@n]))==1) { for(@zz=1;@zz<=nva;@zz++) { if(lead(primary[2*@k-1][@n])/var(@zz)!=0) { jmap1[@zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n] +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]); jmap2[@zz]=primary[2*@k-1][@n]; @qht[@n]=var(@zz); } } jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0); } } if(size(subst(jmap[nva],var(1),0)-var(nva))!=0) { // jmap[nva]=subst(jmap[nva],var(1),0); //hier geaendert +untersuchen!!!!!!!!!!!!!! } phi1=@P,jmap1; phi=@P,jmap; for(@n=1;@n<=nva;@n++) { jmap[@n]=-(jmap[@n]-2*var(@n)); } psi=@P,jmap; psi1=@P,jmap2; @qh=phi(@qht); //=================== the new part ============================ if (npars(basering)>1) { @qh=groebner(@qh,"par2var"); } else { @qh=groebner(@qh); } //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // ideal @qh=homog(imap(@P,@qht),@t); // // ideal @qh1=std(@qh); // @hilb=hilb(@qh1,1); // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // ideal @qh=homog(imap(@P,@qh),@t); // kill @Phelp; // @qh=std(@qh,@hilb); // @qh=subst(@qh,@t,1); // setring @P; // @qh=imap(@Phelp1,@qh); // kill @Phelp1; // @qh=clearSB(@qh); // attrib(@qh,"isSB",1); //============================================================= ser1=phi1(ser); @lh=zero_decomp (@qh,phi(ser1),@wr); kill lres0; list lres0; if((size(@lh)==2)&&(@lh[1]!=1)) { helpprim=@lh[2]; lres0[1]=primary[2*@k-1]; attrib(lres0[1],"isSB",1); ser1=psi(helpprim); lres0[2]=psi1(ser1); if(size(reduce(lres0[2],lres0[1],5))==0) { primary[2*@k]=primary[2*@k-1]; continue; } } else { lres1=psi(@lh); lres0=psi1(lres1); } //=================== the new part ============================ primary=delete(primary,2*@k-1); primary=delete(primary,2*@k-1); @k--; if(size(lres0)==2) { lres0[2]=groebner(lres0[2]); } else { for(@n=1;@n<=size(lres0) div 2;@n++) { lres0[2*@n-1]=groebner(lres0[2*@n-1]); lres0[2*@n]=groebner(lres0[2*@n]); } } primary=primary+lres0; //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // list @lvec; // list @lr=imap(@P,lres0); // ideal @lr1; // // if(size(@lr)==2) // { // @lr[2]=homog(@lr[2],@t); // @lr1=std(@lr[2]); // @lvec[2]=hilb(@lr1,1); // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lvec[2*@n]=@lvec[2*@n-1]; // } // else // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lr[2*@n]=homog(@lr[2*@n],@t); // @lr1=std(@lr[2*@n]); // @lvec[2*@n]=hilb(@lr1,1); // // } // } // } // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // list @lr=imap(@Phelp,@lr); // // kill @Phelp; // if(size(@lr)==2) // { // @lr[2]=std(@lr[2],@lvec[2]); // @lr[2]=subst(@lr[2],@t,1); // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=@lr[2*@n-1]; // attrib(@lr[2*@n],"isSB",1); // } // else // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]); // @lr[2*@n]=subst(@lr[2*@n],@t,1); // } // } // } // kill @lvec; // setring @P; // lres0=imap(@Phelp1,@lr); // kill @Phelp1; // for(@n=1;@n<=size(lres0);@n++) // { // lres0[@n]=clearSB(lres0[@n]); // attrib(lres0[@n],"isSB",1); // } // // primary[2*@k-1]=lres0[1]; // primary[2*@k]=lres0[2]; // @s=size(primary) div 2; // for(@n=1;@n<=size(lres0) div 2-1;@n++) // { // primary[2*@s+2*@n-1]=lres0[2*@n+1]; // primary[2*@s+2*@n]=lres0[2*@n+2]; // } // @k--; //============================================================= } } return(primary); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; i=std(i); list pr= zero_decomp(i,ideal(0),0); pr; } /////////////////////////////////////////////////////////////////////////////// proc extF(list l,list #) { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; //zero_dimensional primary decomposition after finite field extension def R=basering; int p=char(R); if((p==0)||(p>13)||(npars(R)>0)){return(l);} int ex=3; if(size(#)>0){ex=#[1];} list peek,peek1; while(size(l)>0) { if(size(l[2])==0) { peek[size(peek)+1]=l[1]; } else { peek1[size(peek1)+1]=l[1]; peek1[size(peek1)+1]=l[2]; } l=delete(l,1); l=delete(l,1); } if(size(peek)==0){return(peek1);} string gnir="ring RH=("+string(p)+"^"+string(ex)+",@a),("+varstr(R)+"),lp;"; execute(gnir); string mp="minpoly="+string(minpoly)+";"; gnir="ring RL=("+string(p)+",@a),("+varstr(R)+"),lp;"; execute(gnir); execute(mp); list L=imap(R,peek); list pr, keep; int i; for(i=1;i<=size(L);i++) { attrib(L[i],"isSB",1); pr=zero_decomp(L[i],0,0); keep=keep+pr; } for(i=1;i<=size(keep);i++) { keep[i]=simplify(keep[i],1); } mp="poly pp="+string(minpoly)+";"; string gnir1="ring RS="+string(p)+",("+varstr(R)+",@a),lp;"; execute(gnir1); execute(mp); list L=imap(RL,keep); for(i=1;i<=size(L);i++) { L[i]=eliminate(L[i]+ideal(pp),@a); } i=0; int j; while(i0) { re[size(re)+1]=L; ke[size(ke)+1]=K; } } } qe=re; re=em; keep=ke; ke=em; } setring R; list qe=imap(P,keep); list pe=imap(P,qe); for(l=1;l<=size(qe);l++) { qe[l]=simplify(qe[l],2); } list rr=pe,qe; return(rr); } /////////////////////////////////////////////////////////////////////////////// proc zeroSepClos(ideal I,ideal F) { //computes the separable closure of the special ideal I //F is the set of special elements of I //returns the separable closure sc(I) of I and an intvec v //such that sc(I)=preimage(frobenius definde by v) //i.e. var(i)----->var(i)^(p^v[i]) ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; if(homog(I)==1){return(maxideal(1));} //assume F[i] irreducible in I and depending only on var(i) def R=basering; int n=nvars(R); int p=char(R); intvec v; v[n]=0; int i,k; list l; for(i=1;i<=n;i++) { l[i]=sep(F[i],i); F[i]=l[i][1]; if(l[i][2]>k){k=l[i][2];} } if(k==0){return(list(I,v));} //the separable case ideal m; for(i=1;i<=n;i++) { m[i]=var(i)^(p^l[i][2]); v[i]=l[i][2]; } map phi=R,m; ideal J=preimage(R,phi,I); return(list(J,v)); } /////////////////////////////////////////////////////////////////////////////// proc insepDecomp_i(int patchPrimaryDecomposition, ideal i) { //decomposes i into special ideals //computes the prime decomposition of the special ideals //and transforms it back to a decomposition of i // if patchPrimaryDecomposition=1, drop unit ideal in the decomposition, // since the unit ideal it is not prime! ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; def R=basering; list pr=zeroSp(i); int l,k; list re,wo,qr; ideal m=maxideal(1); ideal K; map phi=R,m; int p=char(R); intvec op=option(get); for(l=1;l<=size(pr[1]);l++) { wo=zeroSepClos(pr[1][l],pr[2][l]); for(k=1;k<=nvars(basering);k++) { m[k]=var(k)^(p^wo[2][k]); } phi=R,m; qr = decomp_i(patchPrimaryDecomposition,wo[1],2); option(redSB); for(k=1;k<=size(qr) div 2;k++) { K=qr[2*k]; K=phi(K); K=groebner(K); re[size(re)+1]=zeroRad(K); } option(set,op); } option(set,op); return(re); } /////////////////////////////////////////////////////////////////////////////// static proc clearSB (ideal i,list #) "USAGE: clearSB(i); i ideal which is SB ordered by monomial ordering RETURN: ideal = minimal SB NOTE: EXAMPLE: example clearSB; shows an example " { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int k,j; poly m; int c=size(i); if(size(#)==0) { for(j=1;j0) { m=lead(i[j]); for(k=j+1;k<=c;k++) { if(size(lead(i[k])/m)>0) { i[k]=0; } } } } } else { j=0; while(j0) { m=lead(i[j]); for(k=j+1;k<=c;k++) { if(size(lead(i[k])/m)>0) { if((leadexp(m)!=leadexp(i[k]))||(#[j]<=#[k])) { i[k]=0; } else { i[j]=0; break; } } } } } } return(simplify(i,2)); } example { "EXAMPLE:"; echo = 2; ring r = (0,a,b),(x,y,z),dp; ideal i=ax2+y,a2x+y,bx; list l=1,2,1; ideal j=clearSB(i,l); j; } /////////////////////////////////////////////////////////////////////////////// static proc clearSBNeu (ideal i,list #) "USAGE: clearSB(i); i ideal which is SB ordered by monomial ordering RETURN: ideal = minimal SB NOTE: EXAMPLE: example clearSB; shows an example " { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int k,j; intvec m,n,v,w; int c=size(i); w=leadexp(0); v[size(i)]=0; j=0; while(j=0) { m=leadexp(i[j]); for(k=j+1;k<=c;k++) { n=leadexp(i[k]); if(n!=w) { if(((m==n)&&(#[j]>#[k]))||((teilt(n,m))&&(n!=m))) { i[j]=0; v[j]=1; break; } if(((m==n)&&(#[j]<=#[k]))||((teilt(m,n))&&(n!=m))) { i[k]=0; v[k]=1; } } } } } return(v); } static proc teilt(intvec a, intvec b) { int i; for(i=1;i<=size(a);i++) { if(a[i]>b[i]){return(0);} } return(1); } /////////////////////////////////////////////////////////////////////////////// static proc independSet (ideal j) "USAGE: independentSet(i); i ideal RETURN: list = new varstring with the independent set at the end, ordstring with the corresponding block ordering, the integer where the independent set starts in the varstring NOTE: EXAMPLE: example independentSet; shows an example " { int n,k,di; list resu,hilf; string var1,var2; list v=indepSet(j,1); for(n=1;n<=size(v);n++) { di=0; var1=""; var2=""; for(k=1;k<=size(v[n]);k++) { if(v[n][k]!=0) { di++; var2=var2+"var("+string(k)+"),"; } else { var1=var1+"var("+string(k)+"),"; } } if(di>0) { var1=var1+var2; var1=var1[1..size(var1)-1]; hilf[1]=var1; hilf[2]="lp"; //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")"; hilf[3]=di; resu[n]=hilf; } else { resu[n]=varstr(basering),ordstr(basering),0; } } return(resu); } example { "EXAMPLE:"; echo = 2; ring s1=(0,x,y),(a,b,c,d,e,f,g),lp; ideal i=ea-fbg,fa+be,ec-fdg,fc+de; i=std(i); list l=independSet(i); l; i=i,g; l=independSet(i); l; ring s=0,(x,y,z),lp; ideal i=z,yx; list l=independSet(i); l; } /////////////////////////////////////////////////////////////////////////////// static proc maxIndependSet (ideal j) "USAGE: maxIndependentSet(i); i ideal RETURN: list = new varstring with the maximal independent set at the end, ordstring with the corresponding block ordering, the integer where the independent set starts in the varstring NOTE: EXAMPLE: example maxIndependentSet; shows an example " { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int n,k,di; list resu,hilf; string var1,var2; list v=indepSet(j,0); for(n=1;n<=size(v);n++) { di=0; var1=""; var2=""; for(k=1;k<=size(v[n]);k++) { if(v[n][k]!=0) { di++; var2=var2+"var("+string(k)+"),"; } else { var1=var1+"var("+string(k)+"),"; } } if(di>0) { var1=var1+var2; var1=var1[1..size(var1)-1]; hilf[1]=var1; hilf[2]="lp"; hilf[3]=di; resu[n]=hilf; } else { resu[n]=varstr(basering),ordstr(basering),0; } } return(resu); } example { "EXAMPLE:"; echo = 2; ring s1=(0,x,y),(a,b,c,d,e,f,g),lp; ideal i=ea-fbg,fa+be,ec-fdg,fc+de; i=std(i); list l=maxIndependSet(i); l; i=i,g; l=maxIndependSet(i); l; ring s=0,(x,y,z),lp; ideal i=z,yx; list l=maxIndependSet(i); l; } /////////////////////////////////////////////////////////////////////////////// static proc prepareQuotientring (int nnp,string order) "USAGE: prepareQuotientring(nnp, order); nnp int, order string RETURN: Kvar(nnp+1),...,var(nvars)[..rest ] EXAMPLE: example prepareQuotientring; shows an example " { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; list rl=ringlist(basering); if (typeof(rl[1])=="int") { int p=rl[1]; list rl2=rl[2]; rl[1]=list(p, list(rl2[nnp+1..nvars(basering)]), list(list("lp",1:(nvars(basering)-nnp))), ideal(0)); rl[2]=list(rl2[1..nnp]); rl[3]=list(list(order,1:nnp),list("C",0)); } else { if (typeof(rl[1])=="list") { list rl1=rl[1]; list rl2=rl[2]; rl1=list(rl1[1][1], rl[1][2]+list(rl2[nnp+1..nvars(basering)]), list(list("lp",1:(size(rl[1][2])+nvars(basering)-nnp))), ideal(0)); rl[1]=rl1; rl[2]=list(rl2[1..nnp]); rl[3]=list(list(order,1:nnp),list("C",0)); } else { ERROR("Unexpected case in prepareQuotientring. Please inform the authors"); } } def quotring=ring(rl); return(quotring); } example { "EXAMPLE:"; echo = 2; ring s1=(0,x),(a,b,c,d,e,f,g),lp; def Q= prepareQuotientring(3,"lp"); Q; } /////////////////////////////////////////////////////////////////////////////// static proc cleanPrimary(list l) { int i,j; list lh; for(i=1;i<=size(l) div 2;i++) { if(deg(l[2*i-1][1])>0) { j++; lh[j]=l[2*i-1]; j++; lh[j]=l[2*i]; } } return(lh); } /////////////////////////////////////////////////////////////////////////////// proc minAssPrimesoldE(ideal I, list #) "USAGE: minAssPrimesoldE(I); I ideal minAssPrimesold(I,1); I ideal (to use also the factorizing Groebner) RETURN: list = the minimal associated prime ideals of I EXAMPLE: example minAssPrimesoldE; shows an example " { return(minAssPrimesold_i(1,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr= minAssPrimesoldE(i); pr; minAssPrimesoldE(i,1); } proc minAssPrimesold(ideal I, list #) "USAGE: minAssPrimesold(I); I ideal minAssPrimesold(i,1); I ideal (to use also the factorizing Groebner) RETURN: list = the minimal associated prime ideals of I. In case I is unit ideal, returns list(ideal(1)); EXAMPLE: example minAssPrimesold; shows an example " { return(minAssPrimesold_i(0,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr= minAssPrimesold(i); pr; minAssPrimesold(i,1); } static proc minAssPrimesold_i(int patchPrimaryDecomposition, ideal i, list #) { // // parameter patchPrimaryDecomposition : if = 1, patch the decomposition( drop unit ideal in the decomposition), // : if = 0, taken no special action in case the unit ideal is in the decomposition // for other parameters see minAssPrimesold, minAssPrimesoldE ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; def @P=basering; if(size(i)==0) { return(list(ideal(0))); } list qr=simplifyIdeal(i); map phi=@P,qr[2]; i=qr[1]; def gnir=ring(ringlist(@P)); setring gnir; ideal i=fetch(@P,i); if(size(#)==0) { int @wr; list tluser,@res; list primary=decomp_i(patchPrimaryDecomposition,i,2); @res[1]=primary; tluser=union(@res); setring @P; if (size(tluser)>0) { list @res=imap(gnir,tluser); return(phi(@res)); } else { return(tluser); } } list @res,empty; ideal ser; def op = option( get ); option( redSB ); list @pr=facstd(i); //if(size(@pr)==1) // { // attrib(@pr[1],"isSB",1); // if((dim(@pr[1])==0)&&(homog(@pr[1])==1)) // { // setring @P; // list @res=maxideal(1); // return(phi(@res)); // } // if(dim(@pr[1])>1) // { // setring @P; // // kill gnir; // execute ("ring gnir1 = ("+charstr(basering)+"), // ("+varstr(basering)+"),(C,lp);"); // ideal i=fetch(@P,i); // list @pr=facstd(i); // // ideal ser; // setring gnir; // @pr=fetch(gnir1,@pr); // kill gnir1; // } // } // option( noredSB ); option( set, op ); int j,k,odim,ndim,count; attrib(@pr[1],"isSB",1); if(#[1]==77) { odim=dim(@pr[1]); count=1; intvec pos; pos[size(@pr)]=0; for(j=2;j<=size(@pr);j++) { attrib(@pr[j],"isSB",1); ndim=dim(@pr[j]); if(ndim>odim) { for(k=count;k0) { if(deg(i[j])>1){break;} j--; } if(j==0){return(1);} if(deg(i[j])==vdim(i)){return(1);} return(0); } static proc minAssPrimesE(ideal I, list #) "USAGE: minAssPrimesE(I); I ideal Optional parameters in list #: (can be entered in any order) 0, "facstd" -> uses facstd to first decompose the ideal 1, "noFacstd" -> does not use facstd (default) "SL" -> the new algorithm is used (default) "GTZ" -> the old algorithm is used RETURN: list = the minimal associated prime ideals of I EXAMPLE: example minAssPrimesE; shows an example " { return(minAssPrimes_i(1,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal I = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr= minAssPrimesE(I); pr; minAssPrimesE(I,1); } static proc minAssPrimes(ideal I, list #) "USAGE: minAssPrimes(I); I ideal Optional parameters in list #: (can be entered in any order) 0, "facstd" -> uses facstd to first decompose the ideal 1, "noFacstd" -> does not use facstd (default) "SL" -> the new algorithm is used (default) "GTZ" -> the old algorithm is used RETURN: list = the minimal associated prime ideals of I. If I is the unit ideal returns list(ideal(1)) ; EXAMPLE: example minAssPrimes; shows an example " { return(minAssPrimes_i(0,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr= minAssPrimes(i); pr; minAssPrimes(i,1); } static proc minAssPrimes_i(int patchPrimaryDecomposition, ideal i, list #) { // parameter patchPrimaryDecomposition: 1 to patch( remove unit ideal from the decomposition) , // 0 for no special action on unit ideal. // for other parameters see 'minAssPrimes', 'minAssPrimesE' ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; if (size(i) == 0) { return(list(ideal(0))); } intvec origOp = option(get); string algorithm; // Algorithm to be used string facstdOption; // To uses proc facstd int j; // Counter def P0 = basering; list Pl=ringlist(P0); intvec dp_w; for(j=nvars(P0);j>0;j--) {dp_w[j]=1;} Pl[3]=list(list("dp",dp_w),list("C",0)); def P=ring(Pl); setring P; ideal i=imap(P0,i); // Set input parameters algorithm = "SL"; // Default: SL algorithm facstdOption = "Facstd"; // Default: facstd is not used int subsystem=0; if(size(#) > 0) { int valid; for(j = 1; j <= size(#); j++) { valid = 0; if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) { if (#[j] == 0) {facstdOption = "noFacstd"; valid = 1;} // If #[j] == 0, facstd is not used. if (#[j] == 1) {facstdOption = "facstd"; valid = 1;} // If #[j] == 1, facstd is used. } if(typeof(#[j]) == "string") { if(#[j] == "GTZ" || #[j] == "SL") { algorithm = #[j]; valid = 1; } if(#[j] == "noFacstd" || #[j] == "facstd") { facstdOption = #[j]; valid = 1; } if(#[j] == "noSubsystem" || #[j] == "subsystem") { subsystem = (#[j]=="subsystem"); valid = 1; } } if(valid == 0) { dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); } } } i=simplify(i,2); intvec op = option(get); list q = simplifyIdeal(i); list re = maxideal(1); int a; map phi = P,q[2]; i=q[1]; if (subsystem==0) { int k; list result; if(npars(P) == 0) { option(redSB); } if (attrib(i,"isSB")!=1) { i=groebner(q[1]); } else { for(j=1;j<=nvars(basering);j++) { if (q[2][j]!=var(j)) {k=1;break;} } if(k) { i=groebner(q[1]); } } if(size(i)==1) { if ( deg(lead(i[1]))==0 ) // we have the unit ideal. { setring P0; option( set,origOp ); if (patchPrimaryDecomposition==1) { return( list() ); } else { return( list(ideal(1)) ); } } } if( dim(i) == -1 ) { option( set,op ); setring P0; option( set,origOp ); return( ideal(1) ); } if( (dim(i) == 0 ) && ( npars(P) == 0) ) { int di = vdim(i); def gnir=changeordTo(P,"lp"); setring gnir; ideal J = std(imap(P,i)); attrib(J, "isSB", 1); if(vdim(J) != di) { J = fglm(P, i); } // list pr = triangMH(J,2); HIER KOENNEN verschiedene Mengen zu gleichen // asoziierten Primidealen fuehren // Aenderung list pr = triangMH(J,2); list qr, re; for(k = 1; k <= size(pr); k++) { if(primT(pr[k])&&(0)) { re[size(re) + 1] = pr[k]; } else { attrib(pr[k], "isSB", 1); // Lines changed if (algorithm == "GTZ") { qr = decomp_i(patchPrimaryDecomposition,pr[k], 2); } else { qr = minAssSL(pr[k]); } for(j = 1; j <= size(qr) div 2; j++) { re[size(re) + 1] = std(qr[2 * j]); } } } setring P; re = imap(gnir, re); re=phi(re); option(set, op); setring(P0); list re=imap(P,re); option( set,origOp ); return(re); } if ((facstdOption == "noFacstd") || (dim(i) == 0)) { if (algorithm == "GTZ") { re[1] = decomp_i(patchPrimaryDecomposition,i, 2); } else { re[1] = minAssSL(i); } re = union(re); option(set, op); re=phi(re); setring(P0); option( set,origOp ); list re=imap(P,re); return(re); } } list qq; if (subsystem) { qq = prep_decomp(i); } else { qq=list(i); } q=list(); for(int ll=size(qq);ll>0;ll--) { q=q+facstd(qq[ll]); } kill qq; kill ll; option(set,op); // Debug dbprint(printlevel - voice, "Components returned by facstd", size(q), q); for(j = 1; j <= size(q); j++) { if(a == 0){attrib(q[j], "isSB", 1);} // Debug dbprint(printlevel - voice, "We compute the decomp of component", j); // Lines changed if (algorithm == "GTZ") { re[j] = decomp_i(patchPrimaryDecomposition,q[j], 2); } else { re[j] = minAssSL(q[j]); } // Debug dbprint(printlevel - voice, "Number of components obtained for this component:", size(re[j]) div 2); dbprint(printlevel - voice, "re[j]:", re[j]); } re = union(re); re=phi(re); setring(P0); list re=imap(P,re); option( set,origOp ); return(re); } static proc union(list li) { int i,j,k; def P=basering; int liSize=size(li); int li1Size=0; if (size(li)>0) { li1Size=size(li[1]); } def ir=changeordTo(basering,"lp"); setring ir; list l; if ( liSize > 0) { if (li1Size > 0) { l = fetch(P,li); } else { ASSUME(1, size(li)==1); l[1] = list(); } } list @erg; for(k=1;k<=size(l);k++) { for(j=1;j<=size(l[k]) div 2;j++) { if(deg(l[k][2*j][1])!=0) { i++; @erg[i]=l[k][2*j]; } } } list @wos; i=0; ideal i1,i2; while(isize(@erg)) { @wos[size(@wos)+1]=@erg[i]; } } } if (size(@erg)>0) { if(deg(@erg[size(@erg)][1])!=0) { @wos[size(@wos)+1]=@erg[size(@erg)]; } } int @wosSize = size(@wos); setring P; list @ser; if (@wosSize>0) { @ser=fetch(ir,@wos); } return(@ser); } proc equidim(ideal i,list #) "USAGE: equidim(I) or equidim(I,1) ; I ideal RETURN: list of equidimensional ideals a[1],...,a[s] with: - a[s] the equidimensional locus of I, i.e. the intersection of the primary ideals of dimension of I, except I is unit ideal. - a[1],...,a[s-1] the lower dimensional equidimensional loci. If I is the unit ideal, a list containing the unit ideal as a[1] is returned. NOTE: An embedded component q (primary ideal) of I can be replaced in the decomposition by a primary ideal q1 with the same radical as q. @* @code{equidim(I,1)} uses the algorithm of Eisenbud/Huneke/Vasconcelos. EXAMPLE:example equidim; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } intvec op ; def P = basering; list eq; intvec w; int n,m; int g=size(i); int a=attrib(i,"isSB"); int homo=homog(i); if(size(#)!=0) { m=1; } if(((homo==1)||(a))&&(find(ordstr(basering),"l")==0) &&(find(ordstr(basering),"s")==0)) { def gnir=ring(ringlist(basering)); setring gnir; ideal i=imap(P,i); ideal j=i; if(a) { attrib(j,"isSB",1); } else { j=groebner(i); } } else { def gnir=changeordTo(basering,"dp"); setring gnir; ideal i=imap(P,i); ideal j=groebner(i); } if(homo==1) { for(n=1;n<=nvars(basering);n++) { w[n]=ord(var(n)); } intvec hil=hilb(j,1,w); } if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1) ||(dim(j)==0)||(dim(j)+g==nvars(basering))) { setring P; eq[1]=i; return(eq); } if(m==0) { ideal k=equidimMax(j); } else { ideal k=equidimMaxEHV(j); } if(size(reduce(k,j,5))==0) { setring P; eq[1]=i; kill gnir; return(eq); } op=option(get); option(returnSB); j=quotient(j,k); option(set,op); list equi=equidim(j); if(deg(equi[size(equi)][1])<=0) { equi[size(equi)]=k; } else { equi[size(equi)+1]=k; } setring P; eq=imap(gnir,equi); kill gnir; return(eq); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),dp; ideal i = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); equidim(i); } /////////////////////////////////////////////////////////////////////////////// proc equidimMax(ideal i) "USAGE: equidimMax(i); i ideal RETURN: ideal of equidimensional locus (of maximal dimension) of i. EXAMPLE: example equidimMax; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(attrib(basering,"global")!=1) { ERROR( "Not implemented for this ordering, please change to a global ordering." ); } def P = basering; ideal eq; intvec w; int n; int g=size(i); int a=attrib(i,"isSB"); int homo=homog(i); if(((homo==1)||(a))&&(find(ordstr(basering),"l")==0) &&(find(ordstr(basering),"s")==0)) { def gnir=ring(ringlist(basering)); setring gnir; ideal i=imap(P,i); ideal j=i; if(a) { attrib(j,"isSB",1); } else { j=groebner(i); } } else { def gnir=changeordTo(basering,"dp"); setring gnir; ideal i=imap(P,i); ideal j=groebner(i); } list indep; ideal equ,equi; if(homo==1) { for(n=1;n<=nvars(basering);n++) { w[n]=ord(var(n)); } intvec hil=hilb(j,1,w); } if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1) ||(dim(j)==0)||(dim(j)+g==nvars(basering))) { setring P; return(i); } indep=maxIndependSet(j); execute("ring gnir1 = ("+charstr(basering)+"),("+indep[1][1]+"),(" +indep[1][2]+");"); if(homo==1) { ideal j=std(imap(gnir,j),hil,w); } else { ideal j=groebner(imap(gnir,j)); } def quotring=prepareQuotientring(nvars(basering)-indep[1][3],"lp"); setring quotring; ideal j=imap(gnir1,j); kill gnir1; j=clearSB(j); ideal h; for(n=1;n<=size(j);n++) { h[n]=leadcoef(j[n]); } setring gnir; ideal h=imap(quotring,h); kill quotring; list l=minSat(j,h); if(deg(l[2])>0) { equ=l[1]; attrib(equ,"isSB",1); j=std(j,l[2]); if(dim(equ)==dim(j)) { equi=equidimMax(j); equ=interred(intersect(equ,equi)); } } else { equ=i; } setring P; eq=imap(gnir,equ); kill gnir; return(eq); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),dp; ideal i = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); equidimMax(i); } /////////////////////////////////////////////////////////////////////////////// static proc islp() { string s=ordstr(basering); int n=find(s,"lp"); if(!n){return(0);} int k=find(s,","); string t=s[k+1..size(s)]; int l=find(t,","); t=s[1..k-1]; int m=find(t,","); if(l+m){return(0);} return(1); } /////////////////////////////////////////////////////////////////////////////// //w=0: GTZ //w=1: SY //w=2: minAss proc algeDecoE(ideal I, int w) {//reduces primery decomposition over algebraic extensions to //the other cases return( algeDeco_i( 1, I, w) ); } //w=0: GTZ //w=1: SY //w=2: minAss // deprecated. use 'algeDecoE()' proc algeDeco(ideal I, int w) {//reduces primery decomposition over algebraic extensions to //the other cases return( algeDeco_i(0, I, w)); } //w=0: GTZ //w=1: SY //w=2: minAss static proc algeDeco_i(int patchPrimaryDecomposition, ideal i, int w) {//reduces primery decomposition over algebraic extensions to //the other cases // if patchPrimaryDecomposition=1, drop unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; // the really needed things: ASSUME(1, typeof(ringlist(basering)[1])=="list"); // in alg. extension //reduces primery decomposition over algebraic extensions to //the other cases def R=basering; int n=nvars(R); intvec op = option(get); //---Anfang Provisorium if((size(i)==2) && (w==2)) { //treats a special case separately which would otherwise take a lot longer in factorization option( redSB ); ideal J = std(i); option( set, op ); if(size(J)==1) { if ( deg(lead(J[1]))==0 ) // we have the unit ideal { if (patchPrimaryDecomposition==1) { return( list() ); } else { return( list( ideal(1) ) ); } } } if ((size(J)==2)&&(deg(J[1])==1)) { // minAssPrimes correspond to factorization of J[2] ideal keep; poly f; int j; for(j=1;j<=nvars(basering);j++) { f=J[2]; while((f/var(j))*var(j)-f==0) { f=f/var(j); keep=keep,var(j); } J[2]=f; } ideal K=factorize(J[2],1); if(deg(K[1])==0){K=0;} K=K+std(keep); ideal L; list resu; for(j=1;j<=size(K);j++) { L=J[1],K[j]; resu[j]=L; } option( set, op ); return(resu); } } //---Ende Provisorium list R_l=ringlist(R); ideal @p=R_l[1][4]; // minpoly R_l[2]=R_l[2]+R_l[1][2]; // vars R_l[1]=R_l[1][1]; // char R_l[3]=list(list("dp",1:size(R_l[2])),list("C",0)); // ord def RH=ring(R_l); kill R_l;setring RH; ideal @pp=imap(R,@p); poly @p=@pp[1]; ideal i=imap(R,i); ideal I=subst(i,var(nvars(basering)),0); int j; for(j=1;j<=ncols(i);j++) { if(i[j]!=I[j]){break;} } if((j>ncols(i))&&(deg(@p)==1)) { setring R; kill RH; // remove extension, set order to dp: list R_l=ringlist(R); R_l[1]=R_l[1][1]; // char R_l[3]=list(list("dp",1:nvars(R)),list("C",0)); // ord def RH=ring(R_l); kill R_l; setring RH; ideal i=imap(R,i); ideal J; } else { i=i,@p; } list pr; if(w==0) { pr=decomp_i(patchPrimaryDecomposition,i); } if(w==1) { pr=prim_dec_i(patchPrimaryDecomposition,i,1); pr=reconvList(pr); } if(w==2) { pr=minAssPrimes_i(patchPrimaryDecomposition,i); } int sizepr = size(pr); if(n0) { list pr=imap(RH,pr); } ideal K; for(j=1;j<=sizepr;j++) { K=groebner(pr[j]); if (size(K)>1) { K = K[2..size(K)]; } pr[j]=K; } setring R; if (sizepr>0) { list pr=imap(RS,pr); } } else { setring R; if (sizepr>0) { list pr=imap(RH,pr); } } list re; if(w==2) { re=pr; } else { re=convList(pr); } option( set, op ); return( re ); } /////////////////////////////////////////////////////////////////////////////// static proc prepare_absprimdec(list primary) { ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; list resu,tempo; string absotto; resu[size(primary) div 2]=list(); for(int ab=1;ab<=size(primary) div 2;ab++) { absotto= absFactorize(primary[2*ab][1],77); tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); resu[ab]=tempo; } return(resu); } /////////////////////////////////////////////////////////////////////////////// static proc decompE(ideal I,list #) "USAGE: decompE(I); I ideal (for primary decomposition) (resp. decompE(I,1); (for the associated primes of dimension of I) ) decompE(I,2); (for the minimal associated primes) ) decompE(I,3); (for the absolute primary decomposition) ) RETURN: list = list of primary ideals and their associated primes (at even positions in the list) (resp. a list of the minimal associated primes) if I is unit ideal, returns emtpy list NOTE: Algorithm of Gianni/Trager/Zacharias EXAMPLE: example decompE; shows an example " { return(decomp_i(1,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal I = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr= decompE(I); pr; testPrimary( pr, I); } static proc decomp(ideal I,list #) "USAGE: decomp(I); I ideal (for primary decomposition) (resp. decomp(I,1); (for the associated primes of dimension of I) ) decomp(I,2); (for the minimal associated primes) ) decomp(I,3); (for the absolute primary decomposition) ) RETURN: list = list of primary ideals and their associated primes (at even positions in the list) (resp. a list of the minimal associated primes) if I is unit ideal, returns list(ideal(1),ideal(1)) ( resp. list(ideal(1))) EXAMPLE: example decomp; shows an example " { return(decomp_i(0,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr= decomp(i); pr; testPrimary( pr, i); } static proc decomp_i(int patchPrimaryDecomposition, ideal i,list #) { // if patchPrimaryDecomposition=1, drop unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. // for other parameters see 'decomp' or 'decompE' ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; intvec op,@vv; list primary,indep,ltras; intvec @vh,isat,@w; int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn; ideal peek=i; ideal ser,tras; if(size(#)>0) { if (typeof(#[1])=="int") { if((#[1]==0)||(#[1]==1)||(#[1]==2)||(#[1]==3)) { if (#[1]!=0) { @wr=#[1];} if(@wr==3){abspri=1;@wr=0;} if(size(#)>1) { seri=1; peek=#[2]; ser=#[3]; } } } } if(abspri) { list absprimary,abskeep,absprimarytmp,abskeeptmp; } //---------------------------------------------------------------- //i is the zero-ideal ? //---------------------------------------------------------------- if(size(i)==0) { primary=ideal(0),ideal(0); if (abspri) { return(prepare_absprimdec(primary));} return(primary); } intvec initialOp = option(get); def @P = basering; int isS=attrib(i,"isSB"); //---------------------------------------------------------------- //i is homgeneous ? //---------------------------------------------------------------- homo=homog(i); if(homo) { if(isS/*attrib(i,"isSB")*/!=1) { //ltras=mstd(i); tras=groebner(i); } else { tras=i; } ltras = tras,tras; attrib( ltras[1], "isSB", 1); if (size(ltras[1])>0) { if ( deg(lead(ltras[1]))==0 ) // we have the unit ideal. { option(set,initialOp); if (patchPrimaryDecomposition==1) { if (abspri) { return(prepare_absprimdec(list())); } return( list() ); } else { primary[1]=ideal(1); primary[2]=ideal(1); if (abspri) { return(prepare_absprimdec(primary));} return( primary ); } } } tras=ltras[1]; attrib(tras,"isSB",1); if((dim(tras)==0) && (!abspri)) { primary[1]=ltras[2]; primary[2]=maxideal(1); option(set,initialOp); if(@wr>0) { list l; l[1]=maxideal(1); l[2]=maxideal(1); return(l); } return(primary); } for(@n=1;@n<=nvars(basering);@n++) { @w[@n]=ord(var(@n)); } intvec @hilb=hilb(tras,1,@w); intvec keephilb=@hilb; } //---------------------------------------------------------------- //pass to the lexicographical ordering and compute a standardbasis //---------------------------------------------------------------- int lp=islp(); def gnir=changeordTo(basering,"lp"); setring gnir; op=option(get); option(redSB); ideal ser=fetch(@P,ser); if(homo==1) { if(!lp) { ideal @j=std(fetch(@P,i),@hilb,@w); } else { ideal @j=fetch(@P,tras); attrib(@j,"isSB",1); } } else { if(lp&&isS) { ideal @j=fetch(@P,i); attrib(@j,"isSB",1); } else { ideal @j=groebner(fetch(@P,i)); } if(size(@j)==1) { if ( deg( lead(@j[1]) )==0 ) // we have the unit ideal. { setring @P; option(set,initialOp); if (patchPrimaryDecomposition==1) { return( list() ); } else { return( list(ideal(1),ideal(1)) ); } } } } option(set,op); if(seri==1) { ideal peek=fetch(@P,peek); attrib(peek,"isSB",1); } else { ideal peek=@j; } if((size(ser)==0)&&(!abspri)) { ideal fried; @n=size(@j); for(@k=1;@k<=@n;@k++) { if(deg(lead(@j[@k]))==1) { fried[size(fried)+1]=@j[@k]; @j[@k]=0; } } if(size(fried)==nvars(basering)) { setring @P; option(set,initialOp); primary[1]=i; primary[2]=i; if (abspri) { return(prepare_absprimdec(primary));} return(primary); } if(size(fried)>0) { string newva; string newma; poly f; for(@k=1;@k<=nvars(basering);@k++) { @n1=0; for(@n=1;@n<=size(fried);@n++) { if(leadmonom(fried[@n])==var(@k)) { @n1=1; break; } } if(@n1==0) { newva=newva+string(var(@k))+","; newma=newma+string(var(@k))+","; } else { newma=newma+string(0)+","; fried[@n]=fried[@n]/leadcoef(fried[@n]); f=fried[@n]-lead(fried[@n]); @j=subst(@j,var(@k),-f); } } newva[size(newva)]=")"; newma[size(newma)]=";"; execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;"); execute("map @kappa=gnir,"+newma); ideal @j= @kappa(@j); @j=std(@j); list pr=decomp_i(patchPrimaryDecomposition,@j,#); if (size(pr)==0) { setring @P; option(set,initialOp); if (abspri) { return(prepare_absprimdec(list()));} return(list()); } setring gnir; list pr=imap(@deirf,pr); for(@k=1;@k<=size(pr);@k++) { @j=pr[@k]+fried; pr[@k]=@j; } setring @P; option(set,initialOp); primary=imap(gnir,pr); if (abspri) { return(prepare_absprimdec(primary));} return(primary); } } //---------------------------------------------------------------- //j is the ring //---------------------------------------------------------------- if (dim(@j)==-1) { setring @P; option(set,initialOp); primary=ideal(1),ideal(1); if (abspri) { return(prepare_absprimdec(primary));} return(primary); } //---------------------------------------------------------------- // the case of one variable //---------------------------------------------------------------- if(nvars(basering)==1) { list fac=factor(@j[1]); list gprimary; for(@k=1;@k<=size(fac[1]);@k++) { if(@wr==0) { gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } else { gprimary[2*@k-1]=ideal(fac[1][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } } setring @P; option(set,initialOp); primary=fetch(gnir,gprimary); if (abspri) { return(prepare_absprimdec(primary));} return(primary); } //------------------------------------------------------------------ //the zero-dimensional case //------------------------------------------------------------------ if (dim(@j)==0) { op=option(get); option(redSB); list gprimary= zero_decomp(@j,ser,@wr); setring @P; primary=fetch(gnir,gprimary); if(size(ser)>0) { primary=cleanPrimary(primary); } if(abspri) { setring gnir; list primary=imap(@P,primary); list resu,tempo; string absotto; map sigma,invsigma; ideal II,jmap; nn=nvars(basering); for(ab=1;ab<=size(primary) div 2;ab++) { II=primary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])==vdim(II)) { absotto= absFactorize(primary[2*ab][1],77); tempo= primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); } else { invsigma=basering,maxideal(1); jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); absotto = absFactorize(II[1],77); II=var(nn); tempo= primary[2*ab-1],primary[2*ab],absotto,string(invsigma(II)); } resu[ab]=tempo; } primary=resu; setring @P; primary=imap(gnir,primary); } option(set,initialOp); return(primary); } poly @gs,@gh,@p; string @va; def quotring; list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep; ideal @h; int jdim=dim(@j); list fett; int lauf,di,newtest; //------------------------------------------------------------------ //search for a maximal independent set indep,i.e. //look for subring such that the intersection with the ideal is zero //j intersected with K[var(indep[3]+1),...,var(nvar] is zero, //indep[1] is the new varstring and indep[2] the string for block-ordering //------------------------------------------------------------------ if(@wr!=1) { allindep=independSet(@j); for(@m=1;@m<=size(allindep);@m++) { if(allindep[@m][3]==jdim) { di++; indep[di]=allindep[@m]; } else { lauf++; restindep[lauf]=allindep[@m]; } } } else { indep=maxIndependSet(@j); } ideal jkeep=@j; if((ordstr(@P)[1]=="w")&&(size(ringlist(@P)[3])==2)) // weighted ordering { list gnir_l=ringlist(gnir); list @P_l=ringlist(@P); gnir_l[3]=@P_l[3]; // ord def @Phelp=ring(gnir_l); kill gnir_l,@P_l; setring @Phelp; } else { def @Phelp=changeordTo(gnir,"dp"); setring @Phelp; } if(homo==1) { if(((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w") ||(ordstr(@P)[3]=="w"))&&(size(ringlist(@P)[3])==2)) { ideal jwork=imap(@P,tras); attrib(jwork,"isSB",1); } else { ideal jwork=std(imap(gnir,@j),@hilb,@w); } } else { ideal jwork=groebner(imap(gnir,@j)); } list hquprimary; poly @p,@q; ideal @h,fac,ser; ideal @Ptest=1; di=dim(jwork); keepdi=di; setring gnir; for(@m=1;@m<=size(indep);@m++) { isat=0; @n2=0; if((indep[@m][1]==varstr(basering))&&(@m==1)) //this is the good case, nothing to do, just to have the same notations //change the ring { def gnir1=ring(ringlist(basering)); setring gnir1; ideal @j=fetch(gnir,@j); attrib(@j,"isSB",1); ideal ser=fetch(gnir,ser); } else { @va=string(maxideal(1)); if(@m==1) { @j=fetch(@P,i); } execute("ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),(" +indep[@m][2]+");"); execute("map phi=gnir,"+@va+";"); op=option(get); option(redSB); ideal @j=groebner(phi(@j)); ideal ser=phi(ser); option(set,op); } if((deg(@j[1])==0)||(dim(@j)K(var(nnpr+1),..,var(nva))[..rest..] //------------------------------------------------------------------------ quotring=prepareQuotientring(nvars(basering)-indep[@m][3],"lp"); //--------------------------------------------------------------------- //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //--------------------------------------------------------------------- ideal @jj=lead(@j); //!! vorn vereinbaren setring quotring; ideal @jj=imap(gnir1,@jj); @vv=clearSBNeu(@jj,fett); //!! vorn vereinbaren setring gnir1; @k=size(@j); for (lauf=1;lauf<=@k;lauf++) { if(@vv[lauf]==1) { @j[lauf]=0; } } @j=simplify(@j,2); setring quotring; // @j considered in the quotientring ideal @j=imap(gnir1,@j); ideal ser=imap(gnir1,ser); kill gnir1; //j is a standardbasis in the quotientring but usually not minimal //here it becomes minimal attrib(@j,"isSB",1); //we need later ggt(h[1],...)=gh for saturation ideal @h; if(deg(@j[1])>0) { for(@n=1;@n<=size(@j);@n++) { @h[@n]=leadcoef(@j[@n]); } //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] op=option(get); option(redSB); list uprimary= zero_decomp(@j,ser,@wr); if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); } else { list uprimary; uprimary[1]=ideal(1); uprimary[2]=ideal(1); } //we need the intersection of the ideals in the list quprimary with the //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal //but fi polynomials, then the intersection of q with the polynomialring //is the saturation of the ideal generated by f1,...,fr with respect to //h which is the lcm of the leading coefficients of the fi considered in //in the quotientring: this is coded in saturn list saturn; ideal hpl; for(@n=1;@n<=size(uprimary);@n++) { uprimary[@n]=interred(uprimary[@n]); // temporary fix hpl=0; for(@n1=1;@n1<=size(uprimary[@n]);@n1++) { hpl=hpl,leadcoef(uprimary[@n][@n1]); } saturn[@n]=hpl; } //-------------------------------------------------------------------- //we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //back to the polynomialring //--------------------------------------------------------------------- setring gnir; collectprimary=imap(quotring,uprimary); lsau=imap(quotring,saturn); @h=imap(quotring,@h); kill quotring; def quotring; @n2=size(quprimary); @n3=@n2; for(@n1=1;@n1<=size(collectprimary) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } if(size(@h)>0) { //--------------------------------------------------------------- //we change to @Phelp to have the ordering dp for saturation //--------------------------------------------------------------- setring @Phelp; @h=imap(gnir,@h); if(@wr!=1) { if(defined(@LL)){kill @LL;} list @LL=minSat(jwork,@h); @Ptest=intersect(@Ptest,@LL[1]); @q=@LL[2]; } else { fac=ideal(0); for(lauf=1;lauf<=ncols(@h);lauf++) { if(deg(@h[lauf])>0) { fac=fac+factorize(@h[lauf],1); } } fac=simplify(fac,6); @q=1; for(lauf=1;lauf<=size(fac);lauf++) { @q=@q*fac[lauf]; } } jwork=std(jwork,@q); keepdi=dim(jwork); if(keepdi0) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); if(@wr==0) { //HIER STATT DURCHSCHNITT SATURIEREN! ideal htest=@Ptest; } else { ideal htest=hquprimary[2]; for (@n1=2;@n1<=size(hquprimary) div 2;@n1++) { htest=intersect(htest,hquprimary[2*@n1]); } } if(size(ser)>0) { ser=intersect(htest,ser); } else { ser=htest; } setring gnir; ser=imap(@Phelp,ser); } if(size(reduce(ser,peek,5))!=0) { for(@m=1;@m<=size(restindep);@m++) { // if(restindep[@m][3]>=keepdi) // { isat=0; @n2=0; if(restindep[@m][1]==varstr(basering)) //the good case, nothing to do, just to have the same notations //change the ring { def gnir1=ring(ringlist(basering)); setring gnir1; ideal @j=fetch(gnir,jkeep); attrib(@j,"isSB",1); } else { @va=string(maxideal(1)); execute("ring gnir1 = ("+charstr(basering)+"),("+ restindep[@m][1]+"),(" +restindep[@m][2]+");"); execute("map phi=gnir,"+@va+";"); op=option(get); option(redSB); if(homo==1) { ideal @j=std(phi(jkeep),keephilb,@w); } else { ideal @j=groebner(phi(jkeep)); } ideal ser=phi(ser); option(set,op); } for (lauf=1;lauf<=size(@j);lauf++) { fett[lauf]=size(@j[lauf]); } //------------------------------------------------------------------ //we have now the following situation: //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may //pass to this quotientring, j is their still a standardbasis, the //leading coefficients of the polynomials there (polynomials in //K[var(nnp+1),..,var(nva)]) are collected in the list h, //we need their ggt, gh, because of the following: //let (j:gh^n)=(j:gh^infinity) then //j*K(var(nnp+1),..,var(nva))[..the rest..] //intersected with K[var(1),...,var(nva)] is (j:gh^n) //on the other hand j=(j,gh^n) intersected with (j:gh^n) //------------------------------------------------------------------ //the arrangement for the quotientring // K(var(nnp+1),..,var(nva))[..the rest..] //and the map phi:K[var(1),...,var(nva)] ----> //--->K(var(nnpr+1),..,var(nva))[..the rest..] //------------------------------------------------------------------ if (defined(quotring)==voice) {kill quotring;} def quotring=prepareQuotientring(nvars(basering)-restindep[@m][3],"lp"); //------------------------------------------------------------------ //we pass to the quotientring K(var(nnp+1),..,var(nva))[..rest..] //------------------------------------------------------------------ setring quotring; // @j considered in the quotientring ideal @j=imap(gnir1,@j); ideal ser=imap(gnir1,ser); kill gnir1; //j is a standardbasis in the quotientring but usually not minimal //here it becomes minimal @j=clearSB(@j,fett); attrib(@j,"isSB",1); //we need later ggt(h[1],...)=gh for saturation ideal @h; for(@n=1;@n<=size(@j);@n++) { @h[@n]=leadcoef(@j[@n]); } //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..] op=option(get); option(redSB); list uprimary= zero_decomp(@j,ser,@wr); if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); //we need the intersection of the ideals in the list quprimary with //the polynomialring, i.e. let q=(f1,...,fr) in the quotientring //such an ideal but fi polynomials, then the intersection of q with //the polynomialring is the saturation of the ideal generated by //f1,...,fr with respect toh which is the lcm of the leading //coefficients of the fi considered in the quotientring: //this is coded in saturn list saturn; ideal hpl; for(@n=1;@n<=size(uprimary);@n++) { hpl=0; for(@n1=1;@n1<=size(uprimary[@n]);@n1++) { hpl=hpl,leadcoef(uprimary[@n][@n1]); } saturn[@n]=hpl; } //------------------------------------------------------------------ //we leave the quotientring K(var(nnp+1),..,var(nva))[..rest..] //back to the polynomialring //------------------------------------------------------------------ setring gnir; collectprimary=imap(quotring,uprimary); lsau=imap(quotring,saturn); @h=imap(quotring,@h); kill quotring; @n2=size(quprimary); @n3=@n2; for(@n1=1;@n1<=size(collectprimary) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } if(@n2>=@n3+2) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(@wr==0) { ser=intersect(ser,hquprimary[2*@n-1]); } else { ser=intersect(ser,hquprimary[2*@n]); } } setring gnir; ser=imap(@Phelp,ser); } // } } if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { if (deg(quprimary[2*ab][1])!=0) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } } quprimary=resu; @wr=3; } if(size(reduce(ser,peek,5))!=0) { htprimary=decomp_i(patchPrimaryDecomposition, @j,@wr,peek,ser); // here we collect now both results primary(sat(j,gh)) // and primary(j,gh^n) @n=size(quprimary); for (@k=1;@k<=size(htprimary);@k++) { quprimary[@n+@k]=htprimary[@k]; } } } } else { if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } quprimary=resu; } } //--------------------------------------------------------------------------- //back to the ring we started with //the final result: primary //--------------------------------------------------------------------------- setring @P; option(set,initialOp); primary=imap(gnir,quprimary); if(!abspri) { primary=cleanPrimary(primary); } if (size(primary)>0) { if (abspri && (typeof(primary[1][1])=="poly")) { return(prepare_absprimdec(primary));} } return(primary); } /////////////////////////////////////////////////////////////////////////////// static proc powerCoeffs(poly f,int e) //computes a polynomial with the same monomials as f but coefficients //the p^e th power of the coefficients of f { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int i; poly g; int ex=char(basering)^e; for(i=1;i<=size(f);i++) { g=g+leadcoef(f[i])^ex*leadmonom(f[i]); } return(g); } /////////////////////////////////////////////////////////////////////////////// proc sep(poly f,int i, list #) "USAGE: input: a polynomial f depending on the i-th variable and optional an integer k considering the polynomial f defined over Fp(t1,...,tm) as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k)) RETURN: the separabel part of f as polynomial in Fp(t1,...,tm) and an integer k to indicate that f should be considerd as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k)) EXAMPLE: example sep; shows an example { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; def R=basering; int k; if(size(#)>0){k=#[1];} poly h=gcd(f,diff(f,var(i))); if((reduce(f,std(h),5)!=0)||(reduce(diff(f,var(i)),std(h),5)!=0)) { ERROR("FEHLER IN GCD"); } poly g1=lift(h,f)[1][1]; // f/h poly h1; while(h!=h1) { h1=h; h=gcd(h,diff(h,var(i))); } if(deg(h1)==0){return(list(g1,k));} //in characteristic 0 we return here k++; ideal ma=maxideal(1); ma[i]=var(i)^char(R); map phi=R,ma; ideal hh=h; //this is technical because preimage works only for ideals poly u=preimage(R,phi,hh)[1]; //h=u(x(i)^p) list g2=sep(u,i,k); //we consider u(t(1)^(p^-1),...,t(m)^(p^-1)) g1=powerCoeffs(g1,g2[2]-k+1); //to have g1 over the same field as g2[1] list g3=sep(g1*g2[1],i,g2[2]); return(g3); } example { "EXAMPLE:"; echo = 2; ring R=(5,t,s),(x,y,z),dp; poly f=(x^25-t*x^5+t)*(x^3+s); sep(f,1); } /////////////////////////////////////////////////////////////////////////////// proc zeroRad(ideal I,list #) "USAGE: zeroRad(I) , I a zero-dimensional ideal RETURN: the radical of I NOTE: Algorithm of Kemper EXAMPLE: example zeroRad; shows an example" { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; if(homog(I)==1){return(maxideal(1));} //I needs to be a reduced standard basis def R=basering; int m=npars(R); int n=nvars(R); int p=char(R); int d=vdim(I); int i,k; list l; if(((p==0)||(p>d))&&(d==deg(I[1]))) { intvec e=leadexp(I[1]); for(i=1;i<=nvars(basering);i++) { if(e[i]!=0) break; } I[1]=sep(I[1],i)[1]; return(interred(I)); } intvec op=option(get); option(redSB); ASSUME(1, dim(I)==0); ideal F=finduni(I);//F[i] generates I intersected with K[var(i)] option(set,op); if(size(#)>0){I=#[1];} for(i=1;i<=n;i++) { l[i]=sep(F[i],i); F[i]=l[i][1]; if(l[i][2]>k){k=l[i][2];} //computation of the maximal k } if((k==0)||(m==0)) //the separable case { intvec save=option(get); option(redSB); I=interred(I+F); option(set,save); return(I); } //I=simplify(I,1); for(i=1;i<=n;i++) //consider all polynomials over { //Fp(t(1)^(p^-k),...,t(m)^(p^-k)) F[i]=powerCoeffs(F[i],k-l[i][2]); } string cR="ring @R="+string(p)+",("+parstr(R)+","+varstr(R)+"),dp;"; execute(cR); ideal F=imap(R,F); string nR="ring @S="+string(p)+",(@y(1..m),"+varstr(R)+","+parstr(R)+"),dp;"; execute(nR); ideal G=fetch(@R,F); //G[i](t(1)^(p^-k),...,t(m)^(p^-k),x(i))=sep(F[i]) ideal I=imap(R,I); ideal J=I+G; poly el=1; k=p^k; for(i=1;i<=m;i++) { J=J,var(i)^k-var(m+n+i); el=el*var(i); } J=eliminate(J,el); setring R; ideal J=imap(@S,J); return(J); } example { "EXAMPLE:"; echo = 2; ring R=(5,t),(x,y),dp; ideal I=x^5-t,y^5-t; zeroRad(I); } /////////////////////////////////////////////////////////////////////////////// static proc algeRad(ideal i) { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; //reduces radical computation over algebraic extensions to the other cases def R=basering; int n=nvars(R); string mp="poly @p="+string(minpoly)+";"; string gnir="ring RH="+string(char(R))+",("+varstr(R)+","+string(par(1)) +"),dp;"; execute(gnir); execute(mp); ideal i=imap(R,i); i=i,@p; ideal ra=radical(i); gnir="ring RS="+string(char(R))+",("+varstr(RH) +"),(dp("+string(n)+"),lp);"; execute(gnir); ideal K=imap(RH,ra); K=std(K); if (ncols(K)>1) { K=K[2..size(K)];} setring R; ideal ra=imap(RS,K); return( ra ); } /////////////////////////////////////////////////////////////////////////////// proc radicalEHV(ideal i) "USAGE: radicalEHV(i); i ideal. RETURN: ideal, the radical of i. NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos, which reduces the computation to the complete intersection case, by taking, in the general case, a generic linear combination of the input. Works only in characteristic 0 or p large. EXAMPLE: example radicalEHV; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } if ( size(i)==0 ) { return(ideal(0)); } if(hasAlgExtensionCoefficient(basering)) { return(algeRad(i)); } if((char(basering)<100)&&(char(basering)!=0)) { "WARNING: The characteristic is too small, the result may be wrong"; } ideal J,I,I0,radI0,L,radI1,I2,radI2; int l,n; intvec op=option(get); matrix M; option(redSB); list m=mstd(i); I=m[2]; option(set,op); if ( dim(m[1])<0 ) { return(ideal(1)); } int cod=nvars(basering)-dim(m[1]); //-------------------complete intersection case:---------------------- if(cod==size(m[2])) { J=minor(jacob(I),cod); return(quotient(I,J)); } //-----first codim elements of I are a complete intersection:--------- for(l=1;l<=cod;l++) { I0[l]=I[l]; } n=dim(std(I0))+cod-nvars(basering); //-----last codim elements of I are a complete intersection:---------- if(n!=0) { for(l=1;l<=cod;l++) { I0[l]=I[size(I)-l+1]; } n=dim(std(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(I)); n=dim(std(I0))+cod-nvars(basering); } //-----taking a more generic linear combination of the input:--------- if(n!=0) { M=transpose(sparsetriag(size(m[2]),cod,0,100)); I0=ideal(M*transpose(I)); n=dim(std(I0))+cod-nvars(basering); } if(n==0) { J=minor(jacob(I0),cod); radI0=quotient(I0,J); L=quotient(radI0,I); radI1=quotient(radI0,L); if(size(reduce(radI1,m[1],5))==0) { return(I); } I2=sat(I,radI1)[1]; if(deg(I2[1])<=0) { return(radI1); } return(intersect(radI1,radicalEHV(I2))); } //---------------------general case------------------------------------- return(radical(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= radicalEHV(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc Ann(module M) "USAGE: Ann(M); M module RETURN: ideal, the annihilator of coker(M) NOTE: The output is the ideal of all elements a of the basering R such that a * R^m is contained in M (m=number of rows of M). EXAMPLE: example Ann; shows an example " { M=prune(M); //to obtain a small embedding ideal ann=quotient1(M,freemodule(nrows(M))); return(ann); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; module M = x2-y2,z3; Ann(M); M = [1,x2],[y,x]; Ann(M); qring Q=std(xy-1); module M=imap(r,M); Ann(M); } /////////////////////////////////////////////////////////////////////////////// //computes the equidimensional part of the ideal i of codimension e static proc int_ass_primary_e(ideal i, int e) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; if(homog(i)!=1) { i=std(i); } list re=sres(i,0); //the resolution re=minres(re); //minimized resolution ideal ann = AnnExt_R(e,re); if ( nvars(basering)-dim(std(ann)) != e ) { return( ideal(1) ); } return(ann); } /////////////////////////////////////////////////////////////////////////////// //computes the annihilator of Ext^n(R/i,R) with given resolution re //n is not necessarily the number of variables // !! borrowed correct code from 'ehv.lib::AnnExtEHV' by Kai Dehmann !! duplicate code!! (jk) static proc AnnExt_R(int n,list re) "USAGE: AnnExt_R(n,re); n integer, re resolution RETURN: ideal, the annihilator of Ext^n(R/I,R) with given resolution re of I " { if(n < 0) { return(ideal(1)); } int l = size(re); if(n < l) { matrix f = transpose(re[n+1]); if(n == 0) { matrix g = matrix(0,1,ncols(f)); } else { matrix g = transpose(re[n]); } module k = syz(f); return(quotient1(g,k)); } if(n == l) { return(Ann(transpose(re[n]))); } return(ideal(1)); } /////////////////////////////////////////////////////////////////////////////// static proc analyze(list pr) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int ii,jj; for(ii=1;ii<=size(pr) div 2;ii++) { dim(std(pr[2*ii])); idealsEqual(pr[2*ii-1],pr[2*ii]); "==========================="; } for(ii=size(pr) div 2;ii>1;ii--) { for(jj=1;jj0; n--) { m=coef(p,var(n)); if(m[1,1]!=1) { p=m[2,1]; break; } } if(deg(p)==0) { p=0; } return(p); } /////////////////////////////////////////////////////// // min_ass_prim_charsets // input: generators of an ideal PS and an integer cho // If cho=0, the given ordering of the variables is used. // Otherwise, the system tries to find an "optimal ordering", // which in some cases may considerably speed up the algorithm // output: the minimal associated primes of PS // algorithm: via characteriostic sets ////////////////////////////////////////////////////// static proc min_ass_prim_charsets_i (int patchPrimaryDecomposition, ideal PS, int cho) { // if patchPrimaryDecomposition=1, drop unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, hasGlobalOrdering(basering) ) ; ASSUME(1, not isQuotientRing(basering) ) ; if((cho<0) or (cho>1)) { ERROR(" must be 0 or 1"); } intvec saveopt=option(get); option(notWarnSB); list L; if(cho==0) { L=min_ass_prim_charsets0_i(patchPrimaryDecomposition,PS); } else { L=min_ass_prim_charsets1_i(patchPrimaryDecomposition,PS); } option(set,saveopt); return(L); } /////////////////////////////////////////////////////// // min_ass_prim_charsets0 // input: generators of an ideal PS // output: the minimal associated primes of PS // algorithm: via characteristic sets // the given ordering of the variables is used ////////////////////////////////////////////////////// static proc min_ass_prim_charsets0_i (int patchPrimaryDecomposition, ideal PS) { // if patchPrimaryDecomposition=1, drop unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; if (size(PS)==0) { return( list(ideal(0))); } intvec op; matrix m=char_series(PS); // We compute an irreducible // characteristic series if ((nrows(m)==1) && (ncols(m)==1) && (m[1,1]==1)) // in case of an empty series: min_ass_prim_charsets1 { return (min_ass_prim_charsets1_i(patchPrimaryDecomposition,PS)); } int i,j,k; list PSI; list PHI; // the ideals given by the characteristic series for(i=nrows(m);i>=1; i--) { PHI[i]=ideal(m[i,1..ncols(m)]); } // We compute the radical of each ideal in PHI ideal I,JS,II; int sizeJS, sizeII; for(i=size(PHI);i>=1; i--) { I=0; for(j=size(PHI[i]);j>0;j--) { I=I+ini_mod(PHI[i][j]); } JS=std(PHI[i]); sizeJS=size(JS); for(j=size(I);j>0;j--) { II=0; sizeII=0; k=0; while(k<=sizeII) // successive saturation { op=option(get); option(returnSB); II=quotient(JS,I[j]); option(set,op); sizeII=size(II); if(sizeII==sizeJS) { for(k=1;k<=sizeII;k++) { if(leadexp(II[k])!=leadexp(JS[k])) break; } } JS=II; sizeJS=sizeII; } } PSI=insert(PSI,JS); } int sizePSI=size(PSI); // We eliminate redundant ideals for(i=1;i=1;i--) { if(size(PSI[i])==0) { PSI=delete(PSI,i); } } if(size(PSI)==1) { if (idealsEqual( PSI[1], ideal(1) )) { if (patchPrimaryDecomposition==1) { return( list() ); } else { return( list(ideal(1)) ); } } } return (PSI); } /////////////////////////////////////////////////////// // min_ass_prim_charsets1 // input: generators of an ideal PS // output: the minimal associated primes of PS // algorithm: via characteristic sets // input: generators of an ideal PS and an integer i // The system tries to find an "optimal ordering" of // the variables ////////////////////////////////////////////////////// static proc min_ass_prim_charsets1_i (int patchPrimaryDecomposition, ideal PS) { // if patchPrimaryDecomposition=1, drop unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; if (size(PS)==0) { return( list(ideal(0))); } intvec op; def oldring=basering; string n=system("neworder",PS); execute("ring r=("+charstr(oldring)+"),("+n+"),dp;"); ideal PS=imap(oldring,PS); matrix m=char_series(PS); // We compute an irreducible // characteristic series // this series may be empty (1x1: 1) int i,j,k,cnt; while ((cnt=1; i--) { PHI[i]=simplify(ideal(m[i,1..ncols(m)]),2); I=0; for(j=ncols(PHI[i]);j>0;j--) { I=I,ini_mod(PHI[i][j]); } if (ncols(I)>1) { I=I[2..ncols(I)]; } ITPHI[i]=I; } setring oldring; matrix m=imap(r,m); list PHI=imap(r,PHI); list ITPHI=imap(r,ITPHI); // We compute the radical of each ideal in PHI ideal I,JS,II; int sizeJS, sizeII; for(i=size(PHI);i>=1; i--) { I=0; for(j=size(PHI[i]);j>0;j--) { I=I+ITPHI[i][j]; } JS=std(PHI[i]); sizeJS=size(JS); for(j=size(I);j>0;j--) { II=0; sizeII=0; k=0; while(k<=sizeII) // successive iteration { op=option(get); option(returnSB); II=quotient(JS,I[j]); option(set,op); //std // II=std(II); sizeII=size(II); if(sizeII==sizeJS) { for(k=1;k<=sizeII;k++) { if(leadexp(II[k])!=leadexp(JS[k])) break; } } JS=II; sizeJS=sizeII; } } PSI=insert(PSI,JS); } int sizePSI=size(PSI); // We eliminate redundant ideals for(i=1;i=1;i--) { if(size(PSI[i])==0) { PSI=delete(PSI,i); } } if(size(PSI)==1) { if (idealsEqual( PSI[1], ideal(1) )) { if (patchPrimaryDecomposition==1) { return( list() ); } else { return( list(ideal(1)) ); } } } return (PSI); } ///////////////////////////////////////////////////// // proc prim_dec // input: generators of an ideal I and an integer choose // If choose=0, min_ass_prim_charsets with the given // ordering of the variables is used. // If choose=1, min_ass_prim_charsets with the "optimized" // ordering of the variables is used. // If choose=2, minAssPrimes from primdec.lib is used // If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used // output: a primary decomposition of I, i.e., a list // of pairs consisting of a standard basis of a primary component // of I and a standard basis of the corresponding associated prime. // To compute the minimal associated primes of a given ideal // min_ass_prim_l is called, i.e., the minimal associated primes // are computed via characteristic sets. // In the homogeneous case, the performance of the procedure // will be improved if I is already given by a minimal set of // generators. Apply minbase if necessary. ////////////////////////////////////////////////////////// static proc prim_dec_i(int patchPrimaryDecomposition, ideal I, int choose) { // if patchPrimaryDecomposition=1, drop unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; if((choose<0) or (choose>3)) { ERROR("ERROR: must be 0 or 1 or 2 or 3"); } ideal H=1; // The intersection of the primary components list U; // the leaves of the decomposition tree, i.e., // pairs consisting of a primary component of I // and the corresponding associated prime list W; // the non-leaf vertices in the decomposition tree. // every entry has 6 components: // 1- the vertex itself , i.e., a standard bais of the // given ideal I (type 1), or a standard basis of a // pseudo-primary component arising from // pseudo-primary decomposition (type 2), or a // standard basis of a remaining component arising from // pseudo-primary decomposition or extraction (type 3) // 2- the type of the vertex as indicated above // 3- the weighted_tree_depth of the vertex // 4- the tester of the vertex // 5- a standard basis of the associated prime // of a vertex of type 2, or 0 otherwise // 6- a list of pairs consisting of a standard // basis of a minimal associated prime ideal // of the father of the vertex and the // irreducible factors of the "minimal // divisor" of the seperator or extractor // corresponding to the prime ideal // as computed by the procedure minsat, // if the vertex is of type 3, or // the empty list otherwise ideal SI=std(I); if(SI[1]==1) // primdecSY(ideal(1)) { ASSUME(1, ncols(SI)==1); if (patchPrimaryDecomposition==1) { return( list() ); } else { return( list(list(ideal(1),ideal(1))) ); } } intvec save=option(get); option(notWarnSB); int ncolsSI=ncols(SI); int ncolsH=1; W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree int weighted_tree_depth; int i,j; int check; list V; // current vertex list VV; // new vertex list QQ; list WI; ideal Qi,SQ,SRest,fac; poly tester; while(1) { i=1; while(1) { while(i<=size(W)) // find vertex V of smallest weighted tree-depth { if (W[i][3]<=weighted_tree_depth) break; i++; } if (i<=size(W)) break; i=1; weighted_tree_depth++; } V=W[i]; W=delete(W,i); // delete V from W // now proceed by type of vertex V if (V[2]==2) // extraction needed { SQ,SRest,fac=extraction(V[1],V[5]); // standard basis of primary component, // standard basis of remaining component, // irreducible factors of // the "minimal divisor" of the extractor // as computed by the procedure minsat, check=0; for(j=1;j<=ncolsH;j++) { if (NF(H[j],SQ,1)!=0) // Q is not redundant { check=1; break; } } if(check==1) // Q is not redundant { QQ=list(); QQ[1]=list(SQ,V[5]); // primary component, associated prime, // i.e., standard bases thereof U=U+QQ; H=intersect(H,SQ); H=std(H); ncolsH=ncols(H); check=0; if(ncolsH==ncolsSI) { for(j=1;j<=ncolsSI;j++) { if(leadexp(H[j])!=leadexp(SI[j])) { check=1; break; } } } else { check=1; } if(check==0) // H==I => U is a primary decomposition { option(set,save); return(U); } } if (SRest[1]!=1) // the remaining component is not // the whole ring { if (rad_con(V[4],SRest)==0) // the new vertex is not the // root of a redundant subtree { VV[1]=SRest; // remaining component VV[2]=3; // pseudoprimdec_special VV[3]=V[3]+1; // weighted depth VV[4]=V[4]; // the tester did not change VV[5]=ideal(0); VV[6]=list(list(V[5],fac)); W=insert(W,VV,size(W)); } } } else { if (V[2]==3) // pseudo_prim_dec_special is needed { QQ,SRest=pseudo_prim_dec_special_charsets_i(patchPrimaryDecomposition,V[1],V[6],choose); // QQ = quadruples: // standard basis of pseudo-primary component, // standard basis of corresponding prime, // seperator, irreducible factors of // the "minimal divisor" of the seperator // as computed by the procedure minsat, // SRest=standard basis of remaining component } else // V is the root, pseudo_prim_dec is needed { QQ,SRest=pseudo_prim_dec_charsets_i(patchPrimaryDecomposition,I,SI,choose); // QQ = quadruples: // standard basis of pseudo-primary component, // standard basis of corresponding prime, // seperator, irreducible factors of // the "minimal divisor" of the seperator // as computed by the procedure minsat, // SRest=standard basis of remaining component } //check for(i=size(QQ);i>=1;i--) //for(i=1;i<=size(QQ);i++) { tester=QQ[i][3]*V[4]; Qi=QQ[i][2]; if(NF(tester,Qi,1)!=0) // the new vertex is not the // root of a redundant subtree { VV[1]=QQ[i][1]; VV[2]=2; VV[3]=V[3]+1; VV[4]=tester; // the new tester as computed above VV[5]=Qi; // QQ[i][2]; VV[6]=list(); W=insert(W,VV,size(W)); } } if (SRest[1]!=1) // the remaining component is not // the whole ring { if (rad_con(V[4],SRest)==0) // the vertex is not the root // of a redundant subtree { VV[1]=SRest; VV[2]=3; VV[3]=V[3]+2; VV[4]=V[4]; // the tester did not change VV[5]=ideal(0); WI=list(); for(i=1;i<=size(QQ);i++) { WI=insert(WI,list(QQ[i][2],QQ[i][4])); } VV[6]=WI; W=insert(W,VV,size(W)); } } } } option(set,save); } ////////////////////////////////////////////////////////////////////////// // proc pseudo_prim_dec_charsets // input: Generators of an arbitrary ideal I, a standard basis SI of I, // and an integer choo // If choo=0, min_ass_prim_charsets with the given // ordering of the variables is used. // If choo=1, min_ass_prim_charsets with the "optimized" // ordering of the variables is used. // If choo=2, minAssPrimes from primdec.lib is used // If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used // output: a pseudo primary decomposition of I, i.e., a list // of pseudo primary components together with a standard basis of the // remaining component. Each pseudo primary component is // represented by a quadrupel: A standard basis of the component, // a standard basis of the corresponding associated prime, the // seperator of the component, and the irreducible factors of the // "minimal divisor" of the seperator as computed by the procedure minsat, // calls proc pseudo_prim_dec_i ////////////////////////////////////////////////////////////////////////// static proc pseudo_prim_dec_charsets_i(int patchPrimaryDecomposition, ideal I, ideal SI, int choo) { // if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; list L; // The list of minimal associated primes, // each one given by a standard basis if((choo==0) or (choo==1)) { L=min_ass_prim_charsets_i(patchPrimaryDecomposition,I,choo); } else { if(choo==2) { L=minAssPrimes_i(patchPrimaryDecomposition,I); } else { L=minAssPrimes_i(patchPrimaryDecomposition,I,1); } for(int i=size(L);i>=1;i--) { L[i]=std(L[i]); } } return (pseudo_prim_dec_i_i(patchPrimaryDecomposition,SI,L)); } //////////////////////////////////////////////////////////////// // proc pseudo_prim_dec_special_charsets // input: a standard basis of an ideal I whose radical is the // intersection of the radicals of ideals generated by one prime ideal // P_i together with one polynomial f_i, the list V6 must be the list of // pairs (standard basis of P_i, irreducible factors of f_i), // and an integer choo // If choo=0, min_ass_prim_charsets with the given // ordering of the variables is used. // If choo=1, min_ass_prim_charsets with the "optimized" // ordering of the variables is used. // If choo=2, minAssPrimes from primdec.lib is used // If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used // output: a pseudo primary decomposition of I, i.e., a list // of pseudo primary components together with a standard basis of the // remaining component. Each pseudo primary component is // represented by a quadrupel: A standard basis of the component, // a standard basis of the corresponding associated prime, the // seperator of the component, and the irreducible factors of the // "minimal divisor" of the seperator as computed by the procedure minsat, // calls proc pseudo_prim_dec_i //////////////////////////////////////////////////////////////// static proc pseudo_prim_dec_special_charsets_i (int patchPrimaryDecomposition, ideal SI,list V6, int choo) { // if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int i,j,l; list m; list L; int sizeL; ideal P,SP; ideal fac; int dimSP; for(l=size(V6);l>=1;l--) // creates a list of associated primes // of I, possibly redundant { P=V6[l][1]; fac=V6[l][2]; for(i=ncols(fac);i>=1;i--) { SP=P+fac[i]; SP=std(SP); if(SP[1]!=1) { if((choo==0) or (choo==1)) { m=min_ass_prim_charsets_i(patchPrimaryDecomposition,SP,choo); // a list of SB } else { if(choo==2) { m=minAssPrimes_i(patchPrimaryDecomposition,SP); } else { m=minAssPrimes_i(patchPrimaryDecomposition,SP,1); } for(j=size(m);j>=1;j--) { m[j]=std(m[j]); } } dimSP=dim(SP); for(j=size(m);j>=1; j--) { if(dim(m[j])==dimSP) { L=insert(L,m[j],size(L)); } } } } } sizeL=size(L); for(i=1;i=1;i--) { if(size(L[i])==0) { L=delete(L,i); } } return (pseudo_prim_dec_i_i(patchPrimaryDecomposition,SI,L)); } //////////////////////////////////////////////////////////////// // proc pseudo_prim_dec_i_i // input: A standard basis of an arbitrary ideal I, and standard bases // of the minimal associated primes of I // output: a pseudo primary decomposition of I, i.e., a list // of pseudo primary components together with a standard basis of the // remaining component. Each pseudo primary component is // represented by a quadrupel: A standard basis of the component Q_i, // a standard basis of the corresponding associated prime P_i, the // seperator of the component, and the irreducible factors of the // "minimal divisor" of the seperator as computed by the procedure minsat, //////////////////////////////////////////////////////////////// static proc pseudo_prim_dec_i_i (int patchPrimaryDecomposition, ideal SI, list L) { // if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; list Q; if (size(L)==1) // one minimal associated prime only // the ideal is already pseudo primary { Q=SI,L[1],1; list QQ; QQ[1]=Q; return (QQ,ideal(1)); } poly f0,f,g; ideal fac; int i,j,k,l; ideal SQi; ideal I'=SI; list QP; int sizeL=size(L); for(i=1;i<=sizeL;i++) { fac=0; for(j=1;j<=sizeL;j++) // compute the seperator sep_i // of the i-th component { if (i!=j) // search g not in L[i], but L[j] { for(k=1;k<=ncols(L[j]);k++) { if(NF(L[j][k],L[i],1)!=0) { break; } } fac=fac+L[j][k]; } } // delete superfluous polynomials fac=simplify(fac,8+2); // saturation SQi,f0,f,fac=minsat_ppd(SI,fac); I'=I',f; QP=SQi,L[i],f0,fac; // the quadrupel: // a standard basis of Q_i, // a standard basis of P_i, // sep_i, // irreducible factors of // the "minimal divisor" of the seperator // as computed by the procedure minsat, Q[i]=QP; } I'=std(I'); return (Q, I'); // I' = remaining component } //////////////////////////////////////////////////////////////// // proc extraction // input: A standard basis of a pseudo primary ideal I, and a standard // basis of the unique minimal associated prime P of I // output: an extraction of I, i.e., a standard basis of the primary // component Q of I with associated prime P, a standard basis of the // remaining component, and the irreducible factors of the // "minimal divisor" of the extractor as computed by the procedure minsat //////////////////////////////////////////////////////////////// static proc extraction (ideal SI, ideal SP) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; list indsets=indepSet(SP,0); poly f; if(size(indsets)!=0) //check, whether dim P != 0 { intvec v; // a maximal independent set of variables // modulo P string U; // the independent variables string A; // the dependent variables int j,k; int a; // the size of A int degf; ideal g; list polys; int sizepolys; list newpoly; def R=basering; //intvec hv=hilb(SI,1); for (k=1;k<=size(indsets);k++) { v=indsets[k]; for (j=1;j<=nvars(R);j++) { if (v[j]==1) { U=U+varstr(j)+","; } else { A=A+varstr(j)+","; a++; } } U[size(U)]=")"; // we compute the extractor of I (w.r.t. U) execute("ring RAU=("+charstr(basering)+"),("+A+U+",(dp("+string(a)+"),dp);"); ideal I=imap(R,SI); //I=std(I,hv); // the standard basis in (R[U])[A] I=std(I); // the standard basis in (R[U])[A] A[size(A)]=")"; execute("ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;"); ideal I=imap(RAU,I); //"std in lokalisierung:"+newline,I; ideal h; for(j=ncols(I);j>=1;j--) { h[j]=leadcoef(I[j]); // consider I in (R(U))[A] } setring R; g=imap(Rloc,h); kill RAU,Rloc; U=""; A=""; a=0; f=lcm(g); newpoly[1]=f; polys=polys+newpoly; newpoly=list(); } f=polys[1]; degf=deg(f); sizepolys=size(polys); for (k=2;k<=sizepolys;k++) { if (deg(polys[k])=1;i--) { f0=f0*fac[i]; } poly f=1; ideal iold; list quotM; quotM[1]=SI; quotM[2]=fac; quotM[3]=f0; // we deal seperately with the first quotient; // factors, which do not contribute to this one, // are omitted iold=quotM[1]; quotM=minquot(quotM); fac=quotM[2]; if(quotM[3]==1) { return(quotM[1],f0,f,fac); } while(special_ideals_equal(iold,quotM[1])==0) { f=f*quotM[3]; iold=quotM[1]; quotM=minquot(quotM); } return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f) } ///////////////////////////////////////////////////// // proc minsat_ppd // input: a standard basis of an ideal I and a polynomial p // output: a standard basis IS of the saturation of I w.r. to p, // the maximal squarefree factor f0 of p, // the "minimal divisor" f of f0 such that the saturation of // I w.r. to f equals the saturation of I w.r. to f0 (which is IS), // the irreducible factors of f ////////////////////////////////////////////////////////// static proc minsat_ppd(ideal SI, ideal fac) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; fac=sort(fac)[1]; int i,k; poly f0=1; for(i=ncols(fac);i>=1;i--) { f0=f0*fac[i]; } poly f=1; ideal iold; list quotM; quotM[1]=SI; quotM[2]=fac; quotM[3]=f0; // we deal seperately with the first quotient; // factors, which do not contribute to this one, // are omitted iold=quotM[1]; quotM=minquot(quotM); fac=quotM[2]; if(quotM[3]==1) { return(quotM[1],f0,f,fac); } while(special_ideals_equal(iold,quotM[1])==0) { f=f*quotM[3]; iold=quotM[1]; quotM=minquot(quotM); k++; } return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f) } ///////////////////////////////////////////////////////////////// // proc minquot // input: a list with 3 components: a standard basis // of an ideal I, a set of irreducible polynomials, and // there product f0 // output: a standard basis of the ideal (I:f0), the irreducible // factors of the "minimal divisor" f of f0 with (I:f0) = (I:f), // the "minimal divisor" f ///////////////////////////////////////////////////////////////// static proc minquot(list tsil) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; intvec op; int i,j,k,action; ideal verg; list l; poly g; ideal laedi=tsil[1]; ideal fac=tsil[2]; poly f=tsil[3]; //std // ideal star=quotient(laedi,f); // star=std(star); op=option(get); option(returnSB); ideal star=quotient(laedi,f); option(set,op); if(special_ideals_equal(laedi,star)==1) { return(laedi,ideal(1),1); } action=1; while(action==1) { if(size(fac)==1) { action=0; break; } for(i=1;i<=size(fac);i++) { g=1; for(j=1;j<=size(fac);j++) { if(i!=j) { g=g*fac[j]; } } //std // verg=quotient(laedi,g); // verg=std(verg); op=option(get); option(returnSB); verg=quotient(laedi,g); option(set,op); if(special_ideals_equal(verg,star)==1) { f=g; fac[i]=0; fac=simplify(fac,2); break; } if(i==size(fac)) { action=0; } } } l=star,fac,f; return(l); } ///////////////////////////////////////////////// // proc special_ideals_equal // input: standard bases of ideal k1 and k2 such that // k1 is contained in k2, or k2 is contained ink1 // output: 1, if k1 equals k2, 0 otherwise ////////////////////////////////////////////////// static proc special_ideals_equal( ideal k1, ideal k2) { int j; if(size(k1)==size(k2)) { for(j=1;j<=size(k1);j++) { if(leadexp(k1[j])!=leadexp(k2[j])) { return(0); } } return(1); } return(0); } /////////////////////////////////////////////////////////////////////////////// static proc convList(list l) { int i; list re,he; for(i=1;i<=size(l) div 2;i++) { he=l[2*i-1],l[2*i]; re[i]=he; } return(re); } /////////////////////////////////////////////////////////////////////////////// static proc reconvList(list l) { int i; list re; for(i=size(l);i>0;i--) { re[2*i-1]=l[i][1]; re[2*i]=l[i][2]; } return(re); } /////////////////////////////////////////////////////////////////////////////// // // The main procedures // /////////////////////////////////////////////////////////////////////////////// proc primdecGTZE(ideal I, list #) "USAGE: primdecGTZE(I); i ideal RETURN: a list pr of primary ideals and their associated primes for a proper ideal, and an empty list for the unit ideal. @format pr[i][1] the i-th primary component, pr[i][2] the i-th prime component. @end format NOTE: - Algorithm of Gianni/Trager/Zacharias. - Designed for characteristic 0, works also in char k > 0, if it terminates (may result in an infinite loop in small characteristic!) - For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring - For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument EXAMPLE: example primdecGTZE; shows an example " { return (primdecGTZ_i(1,I, #)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z3+2; ideal I = p*q^2,y-z2; list pr = primdecGTZE(I); pr; ideal J = 1; list prempty = primdecGTZE(J); prempty; } proc primdecGTZ(ideal I, list #) "USAGE: primdecGTZ(I); I ideal RETURN: a list pr of primary ideals and their associated primes for a proper ideal I, otherwise pr = list( list( ideal(1), ideal(1) ) @format pr[i][1] the i-th primary component, pr[i][2] the i-th prime component. @end format NOTE: - Algorithm of Gianni/Trager/Zacharias. - Designed for characteristic 0, works also in char k > 0, if it terminates (may result in an infinite loop in small characteristic!) - For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring - For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument EXAMPLE: example primdecGTZ; shows an example " { return (primdecGTZ_i(0, I , #)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; list pr = primdecGTZ(i); pr; } static proc primdecGTZ_i(int patchPrimaryDecomposition,ideal i, list #) { // if parameter patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. // For other parameters see 'primdecGTZ' or 'primdecGTZE'. ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(size(#)>0) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=primdecGTZ_i(patchPrimaryDecomposition,i,#); int sizeli = size(li); setring r; if (sizeli==0) { return ( list() ); } list li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k][2]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } if(minpoly!=0) { return(algeDeco_i(patchPrimaryDecomposition,i,0)); ERROR( "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal" ); } return(convList(decomp_i(patchPrimaryDecomposition,i,#))); } /////////////////////////////////////////////////////////////////////////////// proc absPrimdecGTZE(ideal I, list #) "USAGE: absPrimdecGTZE(I); I ideal ASSUME: Ground field has characteristic 0. RETURN: a ring containing two lists: @code{absolute_primes}, the absolute prime components of I, and @code{primary_decomp}, the output of @code{primdecGTZ(I)}. Will fail for unit ideal. The list absolute_primes has to be interpreted as follows: each entry describes a class of conjugated absolute primes, @format absolute_primes[i][1] the absolute prime component, absolute_primes[i][2] the number of conjugates. @end format The first entry of @code{absolute_primes[i][1]} is the minimal polynomial of a minimal finite field extension over which the absolute prime component is defined. For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring. For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument NOTE: Algorithm of Gianni/Trager/Zacharias combined with the @code{absFactorize} command. SEE ALSO: primdecGTZ; absFactorize EXAMPLE: example absPrimdecGTZE; shows an example " { return(absPrimdecGTZ_i(1,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z3+2; ideal I = p*q^2,y-z2; def S = absPrimdecGTZE(I); setring S; absolute_primes; } proc absPrimdecGTZ(ideal I, list #) "USAGE: absPrimdecGTZ(I); I ideal ASSUME: Ground field has characteristic 0. RETURN: a ring containing two lists: @code{absolute_primes}, the absolute prime components of I, and @code{primary_decomp}, the output of @code{primdecGTZ(I)}. The list absolute_primes has to be interpreted as follows: each entry describes a class of conjugated absolute primes, @format absolute_primes[i][1] the absolute prime component, absolute_primes[i][2] the number of conjugates. @end format The first entry of @code{absolute_primes[i][1]} is the minimal polynomial of a minimal finite field extension over which the absolute prime component is defined. For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring. For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument NOTE: Algorithm of Gianni/Trager/Zacharias combined with the @code{absFactorize} command. SEE ALSO: primdecGTZ; absFactorize EXAMPLE: example absPrimdecGTZ; shows an example " { return(absPrimdecGTZ_i(0,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; def S = absPrimdecGTZ(i); setring S; absolute_primes; } static proc absPrimdecGTZ_i(int patchPrimaryDecomposition, ideal I, list #) { // if parameter patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. // For other parameters see 'absPrimdecGTZ' or 'absPrimdecGTZE'. ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if (char(basering) != 0) { ERROR("primdec.lib::absPrimdecGTZ is only implemented for "+ +"characteristic 0"); } if(size(#)>0) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithm automatically passes to the global case // hence prepare to go back to an appropriate new ring def r=basering; ideal max_of_r=maxideal(1); def s=changeord(list(list("dp",1:nvars(basering)))); setring s; def I=imap(r,I); def S=absPrimdecGTZ_i(patchPrimaryDecomposition,I); setring S; ring r1=char(basering),var(nvars(r)+1),dp; def rS=r+r1; // move objects to appropriate ring and clean up setring rS; def max_of_r=imap(r,max_of_r); attrib(max_of_r,"isSB",1); def absolute_primes=imap(S,absolute_primes); def primary_decomp=imap(S,primary_decomp); if(!defined(keep_comp)) { ideal tempid; for(int k=size(absolute_primes);k>=1;k--) { tempid=absolute_primes[k][1]; tempid[1]=0; // ignore minimal polynomial if(size(reduce(lead(tempid),max_of_r,5))!=0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering absolute_primes=delete(absolute_primes,k); } } for(k=size(primary_decomp);k>=1;k--) { if(mindeg(std(lead(primary_decomp[k][2]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering primary_decomp=delete(primary_decomp,k); } } kill tempid; } export(primary_decomp); export(absolute_primes); return(rS); } if(minpoly!=0) { //return(algeDeco_i(patchPrimaryDecomposition,I,0)); ERROR( "Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal" ); } def R=basering; int n=nvars(R); list L=decomp_i(patchPrimaryDecomposition,I,3); if (patchPrimaryDecomposition && size(L)==0 ) { "// will not handle case with unit ideal"; } string newvar=L[1][3]; int k=find(newvar,",",find(newvar,",")+1); newvar=newvar[k+1..size(newvar)]; list lR=ringlist(R); int i,de,ii; intvec vv=1:n; //for(i=1;i<=n;i++){vv[i]=1;} list orst; orst[1]=list("dp",vv); orst[2]=list("dp",intvec(1)); orst[3]=list("C",0); lR[3]=orst; lR[2][n+1] = newvar; def Rz = ring(lR); setring Rz; list L=imap(R,L); list absolute_primes,primary_decomp; ideal I,M,N,K; M=maxideal(1); N=maxideal(1); poly p,q,f,g; map phi,psi; string tvar; for(i=1;i<=size(L);i++) { tvar=L[i][4]; ii=find(tvar,"+"); while(ii) { tvar=tvar[ii+1..size(tvar)]; ii=find(tvar,"+"); } for(ii=1;ii<=nvars(basering);ii++) { if(tvar==string(var(ii))) break; } I=L[i][2]; execute("K="+L[i][3]+";"); p=K[1]; q=K[2]; execute("f="+L[i][4]+";"); g=2*var(ii)-f; M[ii]=f; N[ii]=g; de=deg(p); psi=Rz,M; phi=Rz,N; I=phi(I),p,q; I=std(I); absolute_primes[i]=list(psi(I),de); primary_decomp[i]=list(L[i][1],L[i][2]); } export(primary_decomp); export(absolute_primes); setring R; dbprint( printlevel-voice+4," // 'absPrimdecGTZ' created a ring, in which two lists absolute_primes (the // absolute prime components) and primary_decomp (the primary and prime // components over the current basering) are stored. // To access the list of absolute prime components, type (if the name S was // assigned to the return value): setring S; absolute_primes; "); return(Rz); } /////////////////////////////////////////////////////////////////////////////// proc primdecSYE(ideal I, list #) "USAGE: primdecSYE(I, c); I ideal, c int (optional) RETURN: a list pr of primary ideals and their associated primes: @format pr[i][1] the i-th primary component, pr[i][2] the i-th prime component. @end format If I is the unit ideal returns an empty list. NOTE: Algorithm of Shimoyama/Yokoyama. @format if c=0, the given ordering of the variables is used, if c=1, minAssChar tries to use an optimal ordering (default), if c=2, minAssGTZ is used, if c=3, minAssGTZ and facstd are used. @end format For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring. For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as third argument EXAMPLE: example primdecSY; shows an example " { return (primdecSY_i(1,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z3+2; ideal I = p*q^2,y-z2; list pr = primdecSYE(I); pr; ideal J = x; list prUnit = primdecSYE(J); prUnit; } proc primdecSY( ideal I, list #) "USAGE: primdecSY(I, c); I ideal, c int (optional) RETURN: a list pr of primary ideals and their associated primes for proper ideal I, otherwise pr[1] is list( ideal(1),ideal(1) )' @format pr[i][1] the i-th primary component, pr[i][2] the i-th prime component. @end format NOTE: Algorithm of Shimoyama/Yokoyama. @format if c=0, the given ordering of the variables is used, if c=1, minAssChar tries to use an optimal ordering (default), if c=2, minAssGTZ is used, if c=3, minAssGTZ and facstd are used. @end format For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring. For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as third argument EXAMPLE: example primdecSY; shows an example " { return (primdecSY_i(0,I,#)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; list pr = primdecSY(i); pr; } static proc primdecSY_i(int patchPrimaryDecomposition, ideal i, list #) { // if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. // For other paremetes see 'primdecSY' or 'primdecSYE' ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(size(#)>1) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=primdecSY_i(patchPrimaryDecomposition,i); int sizeli = size(li); setring r; if (sizeli==0) { return ( list() ); } def li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k][2]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } i=simplify(i,2); if ((i[1]==0)||(i[1]==1)) { list L = list(ideal(i[1]), ideal(i[1]) ); return(list(L)); } if(minpoly!=0) { return(algeDeco_i(patchPrimaryDecomposition,i,1)); } if (size(#)!=0) { return(prim_dec_i(patchPrimaryDecomposition,i,#[1])); } else { return(prim_dec_i(patchPrimaryDecomposition,i,1)); } } proc minAssGTZE(ideal I,list #) "USAGE: minAssGTZE(I[, l]); I ideal, l list (optional) @* Optional parameters in list l (can be entered in any order): @* 0, \"facstd\" -> uses facstd to first decompose the ideal (default) @* 1, \"noFacstd\" -> does not use facstd @* \"GTZ\" -> the original algorithm by Gianni, Trager and Zacharias is used @* \"SL\" -> GTZ algorithm with modificiations by Laplagne is used (default) RETURN: a list, the minimal associated prime ideals of I. NOTE: - Designed for characteristic 0, works also in char k > 0 based on an algorithm of Yokoyama - For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring - For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument EXAMPLE: example minAssGTZE; shows an example " { list result = minAssGTZ_i(1,I,#); return(result); } 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 = minAssGTZE(I); pr; ideal J = 1; list prempty = minAssGTZE(J); prempty; } proc minAssGTZ(ideal I, list #) "USAGE: minAssGTZ(I[, l]); I ideal, l list (optional) @* Optional parameters in list l (can be entered in any order): @* 0, \"facstd\" -> uses facstd to first decompose the ideal (default) @* 1, \"noFacstd\" -> does not use facstd @* \"GTZ\" -> the original algorithm by Gianni, Trager and Zacharias is used @* \"SL\" -> GTZ algorithm with modificiations by Laplagne is used (default) RETURN: a list, the minimal associated prime ideals of proper ideal I, otherwise ideal(1) NOTE: - Designed for characteristic 0, works also in char k > 0 based on an algorithm of Yokoyama - For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring - For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument EXAMPLE: example minAssGTZ; shows an example " { list result = minAssGTZ_i(0,I,#); return(result); } 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 = minAssGTZ(i); pr; } /////////////////////////////////////////////////////////////////////////////// static proc minAssGTZ_i(int patchPrimaryDecomposition, ideal i,list #) { // if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. // For other parameters see 'minAssGTZ' or 'minAssGTZE' ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(size(#)>0) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=minAssGTZ_i(patchPrimaryDecomposition,i); int sizeli = size(li); setring r; if (sizeli==0) { return(list()); } def li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } int j; string algorithm; string facstdOption; string subsystem="noSubsystem"; int useFac; // Set input parameters algorithm = "SL"; // Default: SL algorithm facstdOption = "facstd"; if(size(#) > 0) { int valid; for(j = 1; j <= size(#); j++) { valid = 0; if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) { if (#[j] == 1) {facstdOption = "noFacstd"; valid = 1;} // If #[j] == 1, facstd is not used. if (#[j] == 0) {facstdOption = "facstd"; valid = 1;} // If #[j] == 0, facstd is used. } if(typeof(#[j]) == "string") { if((#[j] == "GTZ") || (#[j] == "SL")) { algorithm = #[j]; valid = 1; } if((#[j] == "noFacstd") || (#[j] == "facstd")) { facstdOption = #[j]; valid = 1; } if((#[j] == "noSubsystem") || (#[j] == "subsystem")) { subsystem = #[j]; valid = 1; } } if(valid == 0) { dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); } } } if(minpoly!=0) { return(algeDeco_i(patchPrimaryDecomposition,i,2)); } list result = minAssPrimes_i(patchPrimaryDecomposition,i, facstdOption, algorithm, subsystem); return(result); } /////////////////////////////////////////////////////////////////////////////// proc minAssCharE(ideal I, list #) "USAGE: minAssCharE(I[,c]); i ideal, c int (optional). RETURN: list, the minimal associated prime ideals of I. If I is the unit ideal returns an empty list. NOTE: If c=0, the given ordering of the variables is used. @* Otherwise, the system tries to find an optimal ordering, which in some cases may considerably speed up the algorithm. @* For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as third argument EXAMPLE: example minAssCharE; shows an example " { return(minAssChar_i(1,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 = minAssCharE(I); pr; ideal J = 5; list prempty = minAssCharE(J); prempty; } proc minAssChar(ideal I, list #) "USAGE: minAssChar(I[,c]); i ideal, c int (optional). RETURN: list, the minimal associated prime ideals of I. If I is the unit ideal returns list( ideal(1) ) NOTE: If c=0, the given ordering of the variables is used. @* Otherwise, the system tries to find an optimal ordering, which in some cases may considerably speed up the algorithm. @* For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as third argument EXAMPLE: example minAssChar; shows an example " { return(minAssChar_i(0,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 = minAssChar(i); pr; } proc minAssChar_i(int patchPrimaryDecomposition, ideal i, list #) { // if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. // For other parameters see 'minAssChar' or 'minAssCharE' ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0,size(#)<3); if(size(#)>1) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=minAssChar_i(patchPrimaryDecomposition,i); int sizeli = size(li); setring r; if (sizeli==0) { return(list()); } def li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } if (size(#)>0) { return(min_ass_prim_charsets_i(patchPrimaryDecomposition,i,#[1])); } else { return(min_ass_prim_charsets_i(patchPrimaryDecomposition,i,1)); } } /////////////////////////////////////////////////////////////////////////////// proc equiRadical(ideal i) "USAGE: equiRadical(I); I ideal RETURN: ideal, intersection of associated primes of I of maximal dimension. NOTE: A combination of the algorithms of Krick/Logar (with modifications by Laplagne) and Kemper is used. Works also in positive characteristic (Kempers algorithm). EXAMPLE: example equiRadical; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } return(radical(i, 1)); } 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= equiRadical(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc radical(ideal i, list #) "USAGE: radical(I[, l]); I ideal, l list (optional) @* Optional parameters in list l (can be entered in any order): @* 0, \"fullRad\" -> full radical is computed (default) @* 1, \"equiRad\" -> equiRadical is computed @* \"KL\" -> Krick/Logar algorithm is used @* \"SL\" -> modifications by Laplagne are used (default) @* \"facstd\" -> uses facstd to first decompose the ideal (default for non homogeneous ideals) @* \"noFacstd\" -> does not use facstd (default for homogeneous ideals) RETURN: ideal, the radical of I (or the equiradical if required in the input parameters) NOTE: A combination of the algorithms of Krick/Logar (with modifications by Laplagne) and Kemper is used. Works also in positive characteristic (Kempers algorithm). EXAMPLE: example radical; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; dbprint(printlevel - voice, "Radical, version 2006.05.08"); if(size(i) == 0){return(ideal(0));} if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // compute radical and go back def j=radical(i); setring r; def j=imap(s,j); return(j); } if(hasAlgExtensionCoefficient(basering)) { return(algeRad(i)); } int j; def P0 = basering; list Pl=ringlist(P0); intvec dp_w; for(j=nvars(P0);j>0;j--) {dp_w[j]=1;} Pl[3]=list(list("dp",dp_w),list("C",0)); def @P=ring(Pl); setring @P; ideal i=imap(P0,i); int il; string algorithm; int useFac,useSubsystem; // Set input parameters algorithm = "SL"; // Default: SL algorithm il = 0; // Default: Full radical (not only equiRadical) if (homog(i) == 1) { // Default: facStd is used, except if the ideal is homogeneous. useFac = 0; } else { useFac = 1; } if(size(#) > 0) { int valid; for(j = 1; j <= size(#); j++) { valid = 0; if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) { il = #[j]; // If il == 1, equiRadical is computed valid = 1; } if(typeof(#[j]) == "string") { if(#[j] == "KL") { algorithm = "KL"; valid = 1; } if(#[j] == "SL") { algorithm = "SL"; valid = 1; } if(#[j] == "noFacstd") { useFac = 0; valid = 1; } if(#[j] == "facstd") { useFac = 1; valid = 1; } if(#[j] == "equiRad") { il = 1; valid = 1; } if(#[j] == "fullRad") { il = 0; valid = 1; } if(#[j] == "subsystem") { useSubsystem = 1; valid = 1; } } if(valid == 0) { dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); } } } ideal rad = 1; intvec op = option(get); list qr = simplifyIdeal(i); map phi = @P, qr[2]; option(redSB); i = groebner(qr[1]); option(set, op); int di = dim(i); if(di == 0) { i = zeroRad(i, qr[1]); option(redSB); i=interred(phi(i)); option(set, op); setring(P0); i=imap(@P,i); return(i); } option(redSB); list pr; list qq; if (useSubsystem) { qq=prep_decomp(i); } else { qq=i; } if(useFac == 1) { for(int ll=size(qq);ll>0;ll--) { pr=pr+facstd(qq[ll]); } kill ll; } else { pr = i; } kill qq; option(set, op); int s = size(pr); if(useFac == 1) { dbprint(printlevel - voice, "Number of components returned by facstd: ", s); } for(j = 1; j <= s; j++) { attrib(pr[s + 1 - j], "isSB", 1); if((size(reduce(rad, pr[s + 1 - j], 5)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il)) { // SL Debug messages dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]); dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j])); if(algorithm == "KL") { rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il)); } if(algorithm == "SL") { rad = intersect(rad, radicalSL(pr[s + 1 - j], il)); } } else { // SL Debug dbprint(printlevel-voice, "The radical of this component is not needed."); dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 5))); dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j])); dbprint(printlevel-voice, "il", il); } } rad=interred(phi(rad)); setring(P0); i=imap(@P,rad); 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 = radical(i); pr; } /////////////////////////////////////////////////////////////////////////////// // // Computes the radical of I using KL algorithm. // The only difference with the previous implementation of KL algorithm is // that now it uses block dp instead of lp ordering for the reduction to the // zerodimensional case. // The reduction step has been moved to the new routine radicalReduction, so that it can be // used also by radicalSL procedure. // static proc radicalKL(ideal I, ideal ser, list #) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; // ideal I The ideal for which the radical is computed // ideal ser Used to reduce components already obtained // list # If #[1] = 1, equiradical is computed. // I needs to be a Groebner basis. if (attrib(I, "isSB") != 1) { I = groebner(I); } ideal rad; // The radical int allIndep = 1; // All max independent sets are used list result = radicalReduction(I, ser, allIndep, #); int done = result[3]; rad = result[1]; if (done == 0) { rad = intersect(rad, radicalKL(result[2], ideal(1), #)); } return(rad); } /////////////////////////////////////////////////////////////////////////////// // // Computes the radical of I via Laplagne algorithm, using zerodimensional radical in // the zero dimensional case. // For the reduction to the zerodimensional case, it uses the procedure // radical, with some modifications to avoid the recursion. // static proc radicalSL(ideal I, list #) // Input = I, ideal // #, list. If #[1] = 1, then computes only the equiradical. // Output = (P, primaryDec) where P = rad(I) and primaryDec is the list of the radicals // obtained in intermediate steps. { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; ideal rad = 1; ideal equiRad = 1; list primes; int k; // Counter int il; // If il = 1, only the equiradical is required. int iDim; // The dimension of I int stop = 0; // Checks if the radical has been obtained if (attrib(I, "isSB") != 1) { I = groebner(I); } iDim = dim(I); // Checks if only equiradical is required if (size(#) > 0) { il = #[1]; } while(stop == 0) { dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals."); primes = radicalSLIteration(I, rad); // A list of primes or intersections of primes, not included in P dbprint (printlevel - voice, "// Output of Iteration Step:"); dbprint (printlevel - voice, primes); if (size(primes) > 0) { dbprint (printlevel - voice, "// We intersect P with the ideal just obtained."); for(k = 1; k <= size(primes); k++) { rad = intersect(rad, primes[k]); if (il == 1) { if (attrib(primes[k], "isSB") != 1) { primes[k] = groebner(primes[k]); } if (iDim == dim(primes[k])) { equiRad = intersect(equiRad, primes[k]); } } } } else { stop = 1; } } if (il == 0) { return(rad); } else { return(equiRad); } } ////////////////////////////////////////////////////////////////////////// // Based on radicalKL. // It contains all of old version of proc radicalKL except the recursion call. // // Output: // #1 -> output ideal, the part of the radical that has been computed // #2 -> complementary ideal, the part of the ideal I whose radical remains to be computed // = (I, h) in KL algorithm // This is not used in the new algorithm. It is part of KL algorithm // #3 -> done, 1: output = radical, there is no need to continue // 0: radical = output \cap \sqrt{complementary ideal} // This is not used in the new algorithm. It is part of KL algorithm static proc radicalReduction(ideal I, ideal ser, int allIndep, list #) { // allMaximal 1 -> Indicates that the reduction to the zerodim case // must be done for all indep set of the leading terms ideal // 0 -> Otherwise // ideal ser Only for radicalKL. (Same as in radicalKL) // list # Only for radicalKL (If #[1] = 1, // only equiradical is required. // It is used to set the value of done.) ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; attrib(I, "isSB", 1); // I needs to be a reduced standard basis list indep, fett; intvec @w, @hilb, op; int @wr, @n, @m, lauf, di; ideal fac, @h, collectrad, lsau; poly @q; string @va; def @P = basering; int jdim = dim(I); // Computes the dimension of I int homo = homog(I); // Finds out if I is homogeneous ideal rad = ideal(1); // The unit ideal ideal te = ser; if(size(#) > 0) { @wr = #[1]; } if(homo == 1) { for(@n = 1; @n <= nvars(basering); @n++) { @w[@n] = ord(var(@n)); } @hilb = hilb(I, 1, @w); } // SL 2006.04.11 1 Debug messages dbprint(printlevel-voice, "//Computes the radical of the ideal:", I); // SL 2006.04.11 2 Debug messages //--------------------------------------------------------------------------- //j is the ring //--------------------------------------------------------------------------- if (jdim==-1) { return(ideal(1), ideal(1), 1); } //--------------------------------------------------------------------------- //the zero-dimensional case //--------------------------------------------------------------------------- if (jdim==0) { return(zeroRad(I), ideal(1), 1); } //------------------------------------------------------------------------- //search for a maximal independent set indep,i.e. //look for subring such that the intersection with the ideal is zero //j intersected with K[var(indep[3]+1),...,var(nvar)] is zero, //indep[1] is the new varstring, indep[2] the string for the block-ordering //------------------------------------------------------------------------- // SL 2006-04-24 1 If allIndep = 0, then it only computes one maximal // independent set. // This looks better for the new algorithm but not for KL // algorithm list parameters = allIndep; indep = newMaxIndependSetDp(I, parameters); // SL 2006-04-24 2 for(@m = 1; @m <= size(indep); @m++) { if((indep[@m][1] == varstr(basering)) && (@m == 1)) //this is the good case, nothing to do, just to have the same notations //change the ring { def gnir1=ring(ringlist(basering)); setring gnir1; ideal @j = fetch(@P, I); attrib(@j, "isSB", 1); } else { @va = string(maxideal(1)); execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),(" + indep[@m][2] + ");"); execute("map phi = @P," + @va + ";"); if(homo == 1) { ideal @j = std(phi(I), @hilb, @w); } else { ideal @j = groebner(phi(I)); } } if((deg(@j[1]) == 0) || (dim(@j) < jdim)) { setring @P; break; } for (lauf = 1; lauf <= size(@j); lauf++) { fett[lauf] = size(@j[lauf]); } //------------------------------------------------------------------------ // We have now the following situation: // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass // to this quotientring, j is there still a standardbasis, the // leading coefficients of the polynomials there (polynomials in // K[var(nnp+1),..,var(nva)]) are collected in the list h, // we need their LCM, gh, because of the following: // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..] // intersected with K[var(1),...,var(nva)] is (j:gh^n) // on the other hand j = ((j, gh^n) intersected with (j : gh^n)) //------------------------------------------------------------------------ // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..] // and the map phi:K[var(1),...,var(nva)] -----> // K(var(nnpr+1),..,var(nva))[..the rest..] //------------------------------------------------------------------------ def quotring = prepareQuotientring(nvars(basering) - indep[@m][3],"dp"); //------------------------------------------------------------------------ // We pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //------------------------------------------------------------------------ setring quotring; // @j considered in the quotientring ideal @j = imap(gnir1, @j); kill gnir1; // j is a standardbasis in the quotientring but usually not minimal // here it becomes minimal @j = clearSB(@j, fett); // We need later LCM(h[1],...) = gh for saturation ideal @h; if(deg(@j[1]) > 0) { for(@n = 1; @n <= size(@j); @n++) { @h[@n] = leadcoef(@j[@n]); } op = option(get); option(redSB); @j = std(@j); //to obtain a reduced standardbasis option(set, op); // SL 1 Debug messages dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j))); ideal zero_rad = zeroRad(@j); dbprint(printlevel - voice, "zero_rad passed"); // SL 2 } else { ideal zero_rad = ideal(1); } // We need the intersection of the ideals in the list quprimary with the // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal // but fi polynomials, then the intersection of q with the polynomialring // is the saturation of the ideal generated by f1,...,fr with respect to // h which is the lcm of the leading coefficients of the fi considered in // the quotientring: this is coded in saturn zero_rad = std(zero_rad); ideal hpl; for(@n = 1; @n <= size(zero_rad); @n++) { hpl = hpl, leadcoef(zero_rad[@n]); } //------------------------------------------------------------------------ // We leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] // back to the polynomialring //------------------------------------------------------------------------ setring @P; collectrad = imap(quotring, zero_rad); lsau = simplify(imap(quotring, hpl), 2); @h = imap(quotring, @h); kill quotring; // Here the intersection with the polynomialring // mentioned above is really computed collectrad = sat2(collectrad, lsau)[1]; if(deg(@h[1])>=0) { fac = ideal(0); for(lauf = 1; lauf <= ncols(@h); lauf++) { if(deg(@h[lauf]) > 0) { fac = fac + factorize(@h[lauf], 1); } } fac = simplify(fac, 6); @q = 1; for(lauf = 1; lauf <= size(fac); lauf++) { @q = @q * fac[lauf]; } op = option(get); option(returnSB); option(redSB); I = quotient(I + ideal(@q), rad); attrib(I, "isSB", 1); option(set, op); } if((deg(rad[1]) > 0) && (deg(collectrad[1]) > 0)) { rad = intersect(rad, collectrad); te = intersect(te, collectrad); te = simplify(reduce(te, I, 1), 2); } else { if(deg(collectrad[1]) > 0) { rad = collectrad; te = intersect(te, collectrad); te = simplify(reduce(te, I, 1), 2); } } if((dim(I) < jdim)||(size(te) == 0)) { break; } if(homo==1) { @hilb = hilb(I, 1, @w); } } // SL 2006.04.11 1 Debug messages dbprint (printlevel-voice, "// Part of the Radical already computed:", rad); dbprint (printlevel-voice, "// Dimension:", dim(groebner(rad))); // SL 2006.04.11 2 Debug messages // SL 2006.04.21 1 New variable "done". // It tells if the radical is already computed or // if it still has to be computed the radical of the new ideal I int done; if(((@wr == 1) && (dim(I) 0) { allMaximal = #[1]; } else { allMaximal = 1; } int nMax; if (allMaximal == 1) { nMax = size(v); } else { nMax = 1; } for(n = 1; n <= nMax; n++) // SL 2006.04.21 2 { di = 0; var1 = ""; var2 = ""; for(k = 1; k <= size(v[n]); k++) { if(v[n][k] != 0) { di++; var2 = var2 + "var(" + string(k) + "), "; } else { var1 = var1 + "var(" + string(k) + "), "; } } if(di > 0) { var1 = var1 + var2; var1 = var1[1..size(var1) - 2]; // The "- 2" removes the trailer comma hilf[1] = var1; // SL 2006.21.04 1 The order is now block dp instead of lp hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")"; // SL 2006.21.04 2 hilf[3] = di; resu[n] = hilf; } else { resu[n] = varstr(basering), ordstr(basering), 0; } } return(resu); } example { "EXAMPLE:"; echo = 2; ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp; ideal i = ea - fbg, fa + be, ec - fdg, fc + de; i = std(i); list l = newMaxIndependSetDp(i); l; i = i, g; l = newMaxIndependSetDp(i); l; ring s = 0, (x, y, z), lp; ideal i = z, yx; list l = newMaxIndependSetDp(i); l; } /////////////////////////////////////////////////////////////////////////////// proc prepareAss(ideal i) "USAGE: prepareAss(I); I ideal RETURN: list, the radicals of the maximal dimensional components of I. NOTE: Uses algorithm of Eisenbud/Huneke/Vasconcelos. EXAMPLE: example prepareAss; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } ideal j=std(i); int cod=nvars(basering)-dim(j); int e; list er; ideal ann; if(homog(i)==1) { resolution re=sres(j,0); //the resolution re=minres(re); //minimized resolution } else { list re=mres(i,0); } for(e=cod;e<=nvars(basering);e++) { ann=AnnExt_R(e,re); if(nvars(basering)-dim(std(ann))==e) { er[size(er)+1]=equiRadical(ann); } } return(er); } 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 = prepareAss(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc equidimMaxEHV(ideal i) "USAGE: equidimMaxEHV(I); I ideal RETURN: ideal, the equidimensional component (of maximal dimension) of I. NOTE: Uses algorithm of Eisenbud, Huneke and Vasconcelos. EXAMPLE: example equidimMaxEHV; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } ideal j=groebner(i); int cod=nvars(basering)-dim(j); if(cod > nvars(basering)) { dbprint(printlevel,"//If I is the entire ring..."); dbprint(printlevel,"//...then return the ideal generated by 1."); return(ideal(1)); } int e; ideal ann; if(homog(i)==1) { resolution re=sres(j,0); //the resolution re=minres(re); //minimized resolution } else { resolution re=mres(j,0); } ann = AnnExt_R(cod,re); if( nvars(basering)-dim(std(ann) ) != cod) { return( ideal(1) ); } 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)); equidimMaxEHV(i); } proc testPrimaryE(list pr, ideal k) "USAGE: testPrimaryE(pr,k); pr a list, k an ideal. ASSUME: pr is the result of a primary decomposition and may be empty ( for the unit ideal) RETURN: int, 1 if the intersection of the ideals in pr is k, 0 if not EXAMPLE: example testPrimaryE; shows an example " { return(testPrimary_i(1,pr,k)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),dp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr = primdecGTZ(i); testPrimaryE(pr,i); } proc testPrimary(list pr, ideal k) "USAGE: testPrimary(pr,k); pr a list, k an ideal. ASSUME: pr is the result of primdecGTZ(k) or primdecSY(k). RETURN: int, 1 if the intersection of the ideals in pr is k, 0 if not EXAMPLE: example testPrimary; shows an example " { return(testPrimary_i(0,pr,k)); } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),dp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr = primdecGTZ(i); testPrimary(pr,i); } static proc testPrimary_i(int patchPrimaryDecomposition,list pr, ideal k) { // if patchPrimaryDecomposition=1, handle the case of an empty decomposition list. // For other parameters see 'testPrimary' ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; int i; pr=reconvList(pr); if (patchPrimaryDecomposition==1) { if (idealsEqual( k, ideal(1)) ) { return( size(pr)==0 ); //list expected to be empty. } } ideal j=pr[1]; for (i=2;i<=size(pr) div 2;i++) { j=intersect(j,pr[2*i-1]); } return(idealsEqual(j,k)); } /////////////////////////////////////////////////////////////////////////////// proc zerodec(ideal I) "USAGE: zerodec(I); I ideal ASSUME: I is zero-dimensional, the characteristic of the ground field is 0 RETURN: list of primary ideals, the zero-dimensional decomposition of I NOTE: The algorithm (of Monico), works well only for a small total number of solutions (@code{vdim(std(I))} should be < 100) and without parameters. In practice, it works also in large characteristic p>0 but may fail for small p. @* If printlevel > 0 (default = 0) additional information is displayed. EXAMPLE: example zerodec; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(2, dim(groebner(I))==0 ); if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } def R=basering; poly q; int j,time; matrix m; list re; poly va=var(1); ideal J=groebner(I); ideal ba=kbase(J); int d=vdim(J); dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d)); //------ compute matrix of multiplication on R/I with generic element p ----- int e=nvars(basering); poly p=randomLast(100)[e]+random(-50,50); //the generic element matrix n[d][d]; time = timer; for(j=2;j<=e;j++) { va=va*var(j); } for(j=1;j<=d;j++) { q=reduce(p*ba[j],J); m=coeffs(q,ba,va); n[j,1..d]=m[1..d,1]; } dbprint(printlevel-voice+2, "// time for computing multiplication matrix (with generic element) : "+ string(timer-time)); //---------------- compute characteristic polynomial of matrix -------------- execute("ring P1=("+charstr(R)+"),T,dp;"); matrix n=imap(R,n); time = timer; poly charpol=det(n-var(1)*freemodule(d)); dbprint(printlevel-voice+2,"// time for computing char poly: "+ string(timer-time)); //------------------- factorize characteristic polynomial ------------------- //check first if constant term of charpoly is != 0 (which is true for //sufficiently generic element) if(charpol[size(charpol)]!=0) { time = timer; list fac=factor(charpol); testFactor(fac,charpol); dbprint(printlevel-voice+2,"// time for factorizing char poly: "+ string(timer-time)); int f=size(fac[1]); //--------------------------- the irreducible case -------------------------- if(f==1) { setring R; re=I; return(re); } //---------------------------- the reducible case --------------------------- //if f_i are the irreducible factors of charpoly, mult=ri, then //are the primary components where g_i = f_i(p). However, substituting p in //f_i may result in a huge object although the final result may be small. //Hence it is better to simultaneously reduce with I. For this we need a new //ring. execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);"); list rfac=imap(P1,fac); intvec ov=option(get);; option(redSB); list re1; ideal new = var(1)-imap(R,p),imap(R,J); attrib(new, "isSB",1); //we know that new is a standard basis for(j=1;j<=f;j++) { re1[j]=reduce(rfac[1][j]^rfac[2][j],new); } setring R; re = imap(P,re1); for(j=1;j<=f;j++) { J=I,re[j]; re[j]=interred(J); } option(set,ov); return(re); } else //------------------- choice of generic element failed ------------------- { dbprint(printlevel-voice+2,"// try new generic element!"); setring R; return(zerodec(I)); } } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; ideal i = x2-2,y2-2; list pr = zerodec(i); pr; } static proc newDecompStepE(ideal I, list #) { return(newDecompStep_i(1,I,#)); } static proc newDecompStep(ideal I, list #) { return(newDecompStep_i(0,I,#)); } /////////////////////////////////////////////////////////////////////////////// static proc newDecompStep_i(int patchPrimaryDecomposition, ideal i, list #) "USAGE: newDecompStep_i(patchPrimaryDecomposition, I); I ideal (for primary decomposition) newDecompStep_i(patchPrimaryDecomposition, I,1); (for the associated primes of dimension of i) newDecompStep_i(patchPrimaryDecomposition, I,2); (for the minimal associated primes) newDecompStep_i(patchPrimaryDecomposition, I,3); (for the absolute primary decomposition (not tested!)) "oneIndep"; (for using only one max indep set) "intersect"; (returns alse the intersection of the components founded) RETURN: list = list of primary ideals and their associated primes (at even positions in the list) (resp. a list of the minimal associated primes) NOTE: Algorithm of Gianni/Trager/Zacharias if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, since the unit ideal it is not prime, otherwise take no special action. EXAMPLE: example newDecompStep; shows an example " { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; intvec op@P, op,@vv; def @P = basering; list primary,indep,ltras; intvec @vh,isat,@w; int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn; ideal peek=i; ideal ser,tras; list data; list result; intvec @hilb; int isS=attrib(i,"isSB"); // Debug dbprint(printlevel - voice, "newDecompStep, v2.0"); string indepOption = "allIndep"; string intersectOption = "noIntersect"; if(size(#)>0) { int count = 1; if(typeof(#[count]) == "string") { if ((#[count] == "oneIndep") or (#[count] == "allIndep")) { indepOption = #[count]; count++; } } if(typeof(#[count]) == "string") { if ((#[count] == "intersect") or (#[count] == "noIntersect")) { intersectOption = #[count]; count++; } } if((typeof(#[count]) == "int") or (typeof(#[count]) == "number")) { if ((#[count]==1)||(#[count]==2)||(#[count]==3)) { @wr=#[count]; if(@wr==3){abspri = 1; @wr = 0;} count++; } } if(size(#)>count) { seri=1; peek=#[count + 1]; ser=#[count + 2]; } } if(abspri) { list absprimary,abskeep,absprimarytmp,abskeeptmp; } homo=homog(i); if(homo==1) { if(attrib(i,"isSB")) { ltras=i,i; } else { //ltras=mstd(i); tras=groebner(i); ltras=tras,tras; } attrib(ltras[1],"isSB",1); tras = ltras[1]; attrib(tras,"isSB",1); if(dim(tras)==0) { primary[1]=ltras[2]; primary[2]=maxideal(1); if(@wr>0) { list l; l[2]=maxideal(1); l[1]=maxideal(1); if (intersectOption == "intersect") { return(list(l, maxideal(1))); } else { return(l); } } if (intersectOption == "intersect") { return(list(primary, primary[1])); } else { return(primary); } } for(@n=1;@n<=nvars(basering);@n++) { @w[@n]=ord(var(@n)); } @hilb=hilb(tras,1,@w); intvec keephilb=@hilb; } //---------------------------------------------------------------- //i is the zero-ideal //---------------------------------------------------------------- if(size(i)==0) { primary=i,i; if (intersectOption == "intersect") { return(list(primary, i)); } else { return(primary); } } //---------------------------------------------------------------- //pass to the lexicographical ordering and compute a standardbasis //---------------------------------------------------------------- int lp=islp(); op@P = option(get); def gnir=changeordTo(basering,"lp"); setring gnir; op=option(get); option(redSB); ideal ser=fetch(@P,ser); if(homo==1) { if(!lp) { ideal @j=std(fetch(@P,i),@hilb,@w); } else { ideal @j=fetch(@P,tras); attrib(@j,"isSB",1); } } else { if(lp&&isS) { ideal @j=fetch(@P,i); attrib(@j,"isSB",1); } else { ideal @j=groebner(fetch(@P,i)); } } option(set,op); if(seri==1) { ideal peek=fetch(@P,peek); attrib(peek,"isSB",1); } else { ideal peek=@j; } if((size(ser)==0)&&(!abspri)) { ideal fried; @n=size(@j); for(@k=1;@k<=@n;@k++) { if(deg(lead(@j[@k]))==1) { fried[size(fried)+1]=@j[@k]; @j[@k]=0; } } if(size(fried)==nvars(basering)) { setring @P; option(set,op@P); primary[1]=i; primary[2]=i; if (intersectOption == "intersect") { return(list(primary, i)); } else { return(primary); } } if(size(fried)>0) { string newva; string newma; for(@k=1;@k<=nvars(basering);@k++) { @n1=0; for(@n=1;@n<=size(fried);@n++) { if(leadmonom(fried[@n])==var(@k)) { @n1=1; break; } } if(@n1==0) { newva=newva+string(var(@k))+","; newma=newma+string(var(@k))+","; } else { newma=newma+string(0)+","; } } newva[size(newva)]=")"; newma[size(newma)]=";"; execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;"); execute("map @kappa=gnir,"+newma); ideal @j= @kappa(@j); @j=simplify(@j, 2); attrib(@j,"isSB",1); result = newDecompStep_i(patchPrimaryDecomposition, @j, indepOption, intersectOption, @wr); if (intersectOption == "intersect") { list pr = result[1]; ideal intersection = result[2]; } else { list pr = result; } setring gnir; list pr=imap(@deirf,pr); for(@k=1;@k<=size(pr);@k++) { @j=pr[@k]+fried; pr[@k]=@j; } if (intersectOption == "intersect") { ideal intersection = imap(@deirf, intersection); @j = intersection + fried; intersection = @j; } setring @P; option(set,op@P); if (intersectOption == "intersect") { return(list(imap(gnir,pr), imap(gnir,intersection))); } else { return(imap(gnir,pr)); } } } //---------------------------------------------------------------- //j is the ring //---------------------------------------------------------------- if (dim(@j)==-1) { setring @P; option(set,op@P); primary=ideal(1),ideal(1); if (intersectOption == "intersect") { return(list(primary, ideal(1))); } else { return(primary); } } //---------------------------------------------------------------- // the case of one variable //---------------------------------------------------------------- if(nvars(basering)==1) { list fac=factor(@j[1]); list gprimary; poly generator=1; ideal gIntersection; for(@k=1;@k<=size(fac[1]);@k++) { if(@wr==0) { gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } else { gprimary[2*@k-1]=ideal(fac[1][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } if (intersectOption == "intersect") { generator = generator * fac[1][@k]; } } if (intersectOption == "intersect") { gIntersection = generator; } setring @P; primary=fetch(gnir,gprimary); if (intersectOption == "intersect") { ideal intersection = fetch(gnir,gIntersection); } if(abspri) { list resu,tempo; string absotto; for(ab=1;ab<=size(primary) div 2;ab++) { absotto= absFactorize(primary[2*ab][1],77); tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); resu[ab]=tempo; } primary=resu; intersection = 1; for(ab=1;ab<=size(primary);ab++) { intersection = intersect(intersection, primary[ab][2]); } } if (intersectOption == "intersect") { return(list(primary, intersection)); } else { return(primary); } } //------------------------------------------------------------------ //the zero-dimensional case //------------------------------------------------------------------ if (dim(@j)==0) { op=option(get); option(redSB); list gprimary= newZero_decomp(@j,ser,@wr); setring @P; primary=fetch(gnir,gprimary); if(size(ser)>0) { primary=cleanPrimary(primary); } if(abspri) { list resu,tempo; string absotto; for(ab=1;ab<=size(primary) div 2;ab++) { absotto= absFactorize(primary[2*ab][1],77); tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); resu[ab]=tempo; } primary=resu; } option(set,op@P); if (intersectOption == "intersect") { return(list(primary, fetch(gnir,@j))); } else { return(primary); } } poly @gs,@gh,@p; string @va; list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep; ideal @h; int jdim=dim(@j); list fett; int lauf,di,newtest; //------------------------------------------------------------------ //search for a maximal independent set indep,i.e. //look for subring such that the intersection with the ideal is zero //j intersected with K[var(indep[3]+1),...,var(nvar] is zero, //indep[1] is the new varstring and indep[2] the string for block-ordering //------------------------------------------------------------------ if(@wr!=1) { allindep = newMaxIndependSetLp(@j, indepOption); for(@m=1;@m<=size(allindep);@m++) { if(allindep[@m][3]==jdim) { di++; indep[di]=allindep[@m]; } else { lauf++; restindep[lauf]=allindep[@m]; } } } else { indep = newMaxIndependSetLp(@j, indepOption); } ideal jkeep=@j; if((ordstr(@P)[1]=="w")&&(size(ringlist(@P)[3])==2)) // weighted ordering { def @Phelp=ring(ringlist(gnir)); setring @Phelp; } else { def @Phelp=changeordTo(gnir,"dp"); setring @Phelp; } if(homo==1) { if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w") ||(ordstr(@P)[3]=="w")) { ideal jwork=imap(@P,tras); attrib(jwork,"isSB",1); } else { ideal jwork=std(imap(gnir,@j),@hilb,@w); } } else { ideal jwork=groebner(imap(gnir,@j)); } list hquprimary; poly @p,@q; ideal @h,fac,ser; //Aenderung================ ideal @Ptest=1; //========================= di=dim(jwork); keepdi=di; ser = 1; setring gnir; for(@m=1; @m<=size(indep); @m++) { data[1] = indep[@m]; result = newReduction(@j, ser, @hilb, @w, jdim, abspri, @wr, data); quprimary = quprimary + result[1]; if(abspri) { absprimary = absprimary + result[2]; abskeep = abskeep + result[3]; } @h = result[5]; ser = result[4]; if(size(@h)>0) { //--------------------------------------------------------------- //we change to @Phelp to have the ordering dp for saturation //--------------------------------------------------------------- setring @Phelp; @h=imap(gnir,@h); //Aenderung================================== if(defined(@LL)){kill @LL;} list @LL=minSat(jwork,@h); @Ptest=intersect(@Ptest,@LL[1]); ser = intersect(ser, @LL[1]); //=========================================== if(@wr!=1) { //Aenderung================================== @q=@LL[2]; //=========================================== //@q=minSat(jwork,@h)[2]; } else { fac=ideal(0); for(lauf=1;lauf<=ncols(@h);lauf++) { if(deg(@h[lauf])>0) { fac=fac+factorize(@h[lauf],1); } } fac=simplify(fac,6); @q=1; for(lauf=1;lauf<=size(fac);lauf++) { @q=@q*fac[lauf]; } } jwork = std(jwork,@q); keepdi = dim(jwork); if(keepdi < di) { setring gnir; @j = imap(@Phelp, jwork); ser = imap(@Phelp, ser); break; } if(homo == 1) { @hilb = hilb(jwork, 1, @w); } setring gnir; ser = imap(@Phelp, ser); @j = imap(@Phelp, jwork); } } if((size(quprimary)==0)&&(@wr==1)) { @j=ideal(1); quprimary[1]=ideal(1); quprimary[2]=ideal(1); } if((size(quprimary)==0)) { keepdi = di - 1; quprimary[1]=ideal(1); quprimary[2]=ideal(1); } //--------------------------------------------------------------- //notice that j=sat(j,gh) intersected with (j,gh^n) //we finished with sat(j,gh) and have to start with (j,gh^n) //--------------------------------------------------------------- if((deg(@j[1])!=0)&&(@wr!=1)) { if(size(quprimary)>0) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); if(@wr==0) { //Aenderung==================================================== //HIER STATT DURCHSCHNITT SATURIEREN! ideal htest=@Ptest; /* ideal htest=hquprimary[1]; for (@n1=2;@n1<=size(hquprimary) div 2;@n1++) { htest=intersect(htest,hquprimary[2*@n1-1]); } */ //============================================================= } else { ideal htest=hquprimary[2]; for (@n1=2;@n1<=size(hquprimary) div 2;@n1++) { htest=intersect(htest,hquprimary[2*@n1]); } } if(size(ser)>0) { ser=intersect(htest,ser); } else { ser=htest; } setring gnir; ser=imap(@Phelp,ser); } if(size(reduce(ser,peek,5))!=0) { for(@m=1;@m<=size(restindep);@m++) { // if(restindep[@m][3]>=keepdi) // { isat=0; @n2=0; if(restindep[@m][1]==varstr(basering)) //the good case, nothing to do, just to have the same notations //change the ring { def gnir1=ring(ringlist(basering)); setring gnir1; ideal @j=fetch(gnir,jkeep); attrib(@j,"isSB",1); } else { @va=string(maxideal(1)); execute("ring gnir1 = ("+charstr(basering)+"),("+ restindep[@m][1]+"),(" +restindep[@m][2]+");"); execute("map phi=gnir,"+@va+";"); op=option(get); option(redSB); if(homo==1) { ideal @j=std(phi(jkeep),keephilb,@w); } else { ideal @j=groebner(phi(jkeep)); } ideal ser=phi(ser); option(set,op); } for (lauf=1;lauf<=size(@j);lauf++) { fett[lauf]=size(@j[lauf]); } //------------------------------------------------------------------ //we have now the following situation: //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may //pass to this quotientring, j is their still a standardbasis, the //leading coefficients of the polynomials there (polynomials in //K[var(nnp+1),..,var(nva)]) are collected in the list h, //we need their ggt, gh, because of the following: //let (j:gh^n)=(j:gh^infinity) then //j*K(var(nnp+1),..,var(nva))[..the rest..] //intersected with K[var(1),...,var(nva)] is (j:gh^n) //on the other hand j=(j,gh^n) intersected with (j:gh^n) //------------------------------------------------------------------ //the arrangement for the quotientring // K(var(nnp+1),..,var(nva))[..the rest..] //and the map phi:K[var(1),...,var(nva)] ----> //--->K(var(nnpr+1),..,var(nva))[..the rest..] //------------------------------------------------------------------ def quotring=prepareQuotientring(nvars(basering)-restindep[@m][3],"lp"); //------------------------------------------------------------------ //we pass to the quotientring K(var(nnp+1),..,var(nva))[..rest..] //------------------------------------------------------------------ setring quotring; // @j considered in the quotientring ideal @j=imap(gnir1,@j); ideal ser=imap(gnir1,ser); kill gnir1; //j is a standardbasis in the quotientring but usually not minimal //here it becomes minimal @j=clearSB(@j,fett); attrib(@j,"isSB",1); //we need later ggt(h[1],...)=gh for saturation ideal @h; for(@n=1;@n<=size(@j);@n++) { @h[@n]=leadcoef(@j[@n]); } //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..] op=option(get); option(redSB); list uprimary= newZero_decomp(@j,ser,@wr); if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); //we need the intersection of the ideals in the list quprimary with //the polynomialring, i.e. let q=(f1,...,fr) in the quotientring //such an ideal but fi polynomials, then the intersection of q with //the polynomialring is the saturation of the ideal generated by //f1,...,fr with respect toh which is the lcm of the leading //coefficients of the fi considered in the quotientring: //this is coded in saturn list saturn; ideal hpl; for(@n=1;@n<=size(uprimary);@n++) { hpl=0; for(@n1=1;@n1<=size(uprimary[@n]);@n1++) { hpl=hpl,leadcoef(uprimary[@n][@n1]); } saturn[@n]=hpl; } //------------------------------------------------------------------ //we leave the quotientring K(var(nnp+1),..,var(nva))[..rest..] //back to the polynomialring //------------------------------------------------------------------ setring gnir; collectprimary=imap(quotring,uprimary); lsau=imap(quotring,saturn); @h=imap(quotring,@h); kill quotring; @n2=size(quprimary); //================NEU========================================= if(deg(quprimary[1][1])<=0){ @n2=0; } //============================================================ @n3=@n2; for(@n1=1;@n1<=size(collectprimary) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } if(@n2>=@n3+2) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(@wr==0) { ser=intersect(ser,hquprimary[2*@n-1]); } else { ser=intersect(ser,hquprimary[2*@n]); } } setring gnir; ser=imap(@Phelp,ser); } // } } if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { if (deg(quprimary[2*ab][1])!=0) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } } quprimary=resu; @wr=3; } if(size(reduce(ser,peek,5))!=0) { if(@wr>0) { // The following line was dropped to avoid the recursion step: //htprimary=newDecompStep_i(patchPrimaryDecomposition, @j,@wr,peek,ser); htprimary = list(); } else { // The following line was dropped to avoid the recursion step: //htprimary=newDecompStep_i(patchPrimaryDecomposition,@j,peek,ser); htprimary = list(); } // here we collect now both results primary(sat(j,gh)) // and primary(j,gh^n) @n=size(quprimary); if (deg(quprimary[1][1])<=0) { @n=0; } for (@k=1;@k<=size(htprimary);@k++) { quprimary[@n+@k]=htprimary[@k]; } } } } else { if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } quprimary=resu; } } //--------------------------------------------------------------------------- //back to the ring we started with //the final result: primary //--------------------------------------------------------------------------- setring @P; option(set,op@P); primary=imap(gnir,quprimary); if (intersectOption == "intersect") { return(list(primary, imap(gnir, ser))); } else { return(primary); } } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal I = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; int patchDecomposition = 1; list pr = newDecompStep_i(patchDecomposition, I); pr; testPrimary( pr, I); } // This was part of proc decomp. // In proc newDecompStep, used for the computation of the minimal associated primes, // this part was separated as a soubrutine to make the code more clear. // Also, since the reduction is performed twice in proc newDecompStep, it should use both times this routine. // This is not yet implemented, since the reduction is not exactly the same and some changes should be made. static proc newReduction(ideal @j, ideal ser, intvec @hilb, intvec @w, int jdim, int abspri, int @wr, list data) { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; string @va; def quotring; intvec op; intvec @vv; def gnir = basering; ideal isat=0; int @n; int @n1 = 0; int @n2 = 0; int @n3 = 0; int homo = homog(@j); int lauf; int @k; list fett; int keepdi; list collectprimary; list lsau; list lnew; ideal @h; list indepInfo = data[1]; list quprimary = list(); //if(abspri) //{ int ab; list absprimarytmp,abskeeptmp; list absprimary, abskeep; //} // Debug dbprint(printlevel - voice, "newReduction, v2.0"); if((indepInfo[1]==varstr(basering))) // &&(@m==1) //this is the good case, nothing to do, just to have the same notations //change the ring { def gnir1=ring(ringlist(basering)); setring gnir1; ideal @j = fetch(gnir, @j); attrib(@j,"isSB",1); ideal ser = fetch(gnir, ser); } else { @va=string(maxideal(1)); //Aenderung============== //if(@m==1) //{ // @j=fetch(@P,i); //} //======================= execute("ring gnir1 = ("+charstr(basering)+"),("+indepInfo[1]+"),(" +indepInfo[2]+");"); execute("map phi=gnir,"+@va+";"); op=option(get); option(redSB); if(homo==1) { ideal @j=std(phi(@j),@hilb,@w); } else { ideal @j=groebner(phi(@j)); } ideal ser=phi(ser); option(set,op); } if((deg(@j[1])==0)||(dim(@j)K(var(nnpr+1),..,var(nva))[..rest..] //------------------------------------------------------------------------ quotring=prepareQuotientring(nvars(basering)-indepInfo[3],"lp"); //--------------------------------------------------------------------- //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //--------------------------------------------------------------------- ideal @jj=lead(@j); //!! vorn vereinbaren setring quotring; ideal @jj=imap(gnir1,@jj); @vv=clearSBNeu(@jj,fett); //!! vorn vereinbaren setring gnir1; @k=size(@j); for (lauf=1;lauf<=@k;lauf++) { if(@vv[lauf]==1) { @j[lauf]=0; } } @j=simplify(@j,2); setring quotring; // @j considered in the quotientring ideal @j=imap(gnir1,@j); ideal ser=imap(gnir1,ser); kill gnir1; //j is a standardbasis in the quotientring but usually not minimal //here it becomes minimal attrib(@j,"isSB",1); //we need later ggt(h[1],...)=gh for saturation ideal @h; if(deg(@j[1])>0) { for(@n=1;@n<=size(@j);@n++) { @h[@n]=leadcoef(@j[@n]); } //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] op=option(get); option(redSB); int zeroMinAss = @wr; if (@wr == 2) {zeroMinAss = 1;} list uprimary= newZero_decomp(@j, ser, zeroMinAss); //HIER if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); } else { list uprimary; uprimary[1]=ideal(1); uprimary[2]=ideal(1); } //we need the intersection of the ideals in the list quprimary with the //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal //but fi polynomials, then the intersection of q with the polynomialring //is the saturation of the ideal generated by f1,...,fr with respect to //h which is the lcm of the leading coefficients of the fi considered in //in the quotientring: this is coded in saturn list saturn; ideal hpl; for(@n=1;@n<=size(uprimary);@n++) { uprimary[@n]=interred(uprimary[@n]); // temporary fix hpl=0; for(@n1=1;@n1<=size(uprimary[@n]);@n1++) { hpl=hpl,leadcoef(uprimary[@n][@n1]); } saturn[@n]=hpl; } //-------------------------------------------------------------------- //we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //back to the polynomialring //--------------------------------------------------------------------- setring gnir; collectprimary=imap(quotring,uprimary); lsau=imap(quotring,saturn); @h=imap(quotring,@h); kill quotring; @n2=size(quprimary); @n3=@n2; for(@n1=1;@n1<=size(collectprimary) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } return(quprimary, absprimary, abskeep, ser, @h); } //////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // Based on minAssGTZ proc minAssE(ideal I,list #) "USAGE: minAssE(I[, l]); I ideal, l list (optional) of parameters, same as minAssGTZ RETURN: a list, the minimal associated prime ideals of I. NOTE: Designed for characteristic 0, works also in char k > 0 based on an algorithm of Yokoyama EXAMPLE: example minAssE; shows an example " { return(minAss_i(1,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 = minAssE(i); pr; ideal j = 1; list prempty = minAssE(j); prempty; } proc minAss(ideal I,list #) "USAGE: minAss(I[, l]); I ideal, l list (optional) of parameters, same as minAssGTZ RETURN: a list, the minimal associated prime ideals of I. If I is the unit ideal, returns list(ideal(1)); NOTE: Designed for characteristic 0, works also in char k > 0 based on an algorithm of Yokoyama EXAMPLE: example minAss; shows an example " { return(minAss_i(0,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 = minAss(i); pr; } static proc minAss_i(int patchPrimaryDecomposition,ideal I,list #) { // if patchPrimaryDecomposition=1, drop the unit ideal in the decomposition, // since the unit ideal it is not prime, otherwise take no special action. // For other parameters see 'minAss' or 'minAssE' return(minAssGTZ_i(patchPrimaryDecomposition,I,#)); } /////////////////////////////////////////////////////////////////////////////// // // Computes the minimal associated primes of I via Laplagne algorithm, // using primary decomposition in the zero dimensional case. // For reduction to the zerodimensional case, it uses the procedure // decomp, with some modifications to avoid the recursion. // static proc minAssSL(ideal I) // Input = I, ideal // Output = primaryDec where primaryDec is the list of the minimal // associated primes and the primary components corresponding to these primes. { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; ideal P = 1; list pd = list(); int k; int stop = 0; list primaryDec = list(); while (stop == 0) { // Debug dbprint(printlevel - voice, "// We call minAssSLIteration to find new prime ideals!"); pd = minAssSLIteration(I, P); // Debug dbprint(printlevel - voice, "// Output of minAssSLIteration:"); dbprint(printlevel - voice, pd); if (size(pd[1]) > 0) { primaryDec = primaryDec + pd[1]; // Debug dbprint(printlevel - voice, "// We intersect the prime ideals obtained."); P = intersect(P, pd[2]); // Debug dbprint(printlevel - voice, "// Intersection finished."); } else { stop = 1; } } // Returns only the primary components, not the radical. return(primaryDec); } /////////////////////////////////////////////////////////////////////////////// // Given an ideal I and an ideal P (intersection of some minimal prime ideals // associated to I), it calculates new minimal prime ideals associated to I // which were not used to calculate P. P = 1 represents empty intersetion. // This version uses Primary Decomposition in the zerodimensional case. static proc minAssSLIteration(ideal I, ideal P); { ASSUME(1, hasFieldCoefficient(basering) ); ASSUME(1, not isQuotientRing(basering) ) ; ASSUME(1, hasGlobalOrdering(basering) ) ; int k = 1; int good = 0; list primaryDec = list(); // Debug dbprint (printlevel-voice, "// We search for an element in P - sqrt(I)."); while ((k <= size(P)) and (good == 0)) { good = 1 - rad_con(P[k], I); k++; } k--; if (good == 0) { // Debug dbprint (printlevel - voice, "// No element was found, P = sqrt(I)."); return (list(primaryDec, ideal(0))); } // Debug dbprint (printlevel - voice, "// We found h = ", P[k]); dbprint (printlevel - voice, "// We calculate the saturation of I with respect to the element just founded."); ideal J = sat(I, P[k])[1]; // Uses decomp from primdec, modified to avoid the recursion. // Debug dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via decomp."); primaryDec = newDecompStep_i( 1, J, "oneIndep", "intersect", 2); // Debug dbprint(printlevel - voice, "// Proc decomp has found", size(primaryDec) div 2, "new primary components."); return(primaryDec); } /////////////////////////////////////////////////////////////////////////////////// // Based on maxIndependSet // Added list # as parameter // If the first element of # is 0, the output is only 1 max indep set. // If no list is specified or #[1] = 1, the output is all the max indep set of the // leading terms ideal. This is the original output of maxIndependSet proc newMaxIndependSetLp(ideal j, list #) "USAGE: newMaxIndependentSetLp(i); i ideal (returns all maximal independent sets of the corresponding leading terms ideal) newMaxIndependentSetLp(i, 0); i ideal (returns only one maximal independent set) RETURN: list = #1. new varstring with the maximal independent set at the end, #2. ordstring with the lp ordering, #3. the number of independent variables NOTE: EXAMPLE: example newMaxIndependentSetLp; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; int n, k, di; list resu, hilf; string var1, var2; list v = indepSet(j, 0); // SL 2006.04.21 1 Lines modified to use only one independent Set string indepOption; if (size(#) > 0) { indepOption = #[1]; } else { indepOption = "allIndep"; } int nMax; if (indepOption == "allIndep") { nMax = size(v); } else { nMax = 1; } for(n = 1; n <= nMax; n++) // SL 2006.04.21 2 { di = 0; var1 = ""; var2 = ""; for(k = 1; k <= size(v[n]); k++) { if(v[n][k] != 0) { di++; var2 = var2 + "var(" + string(k) + "), "; } else { var1 = var1 + "var(" + string(k) + "), "; } } if(di > 0) { var1 = var1 + var2; var1 = var1[1..size(var1) - 2]; // The "- 2" removes the trailer comma hilf[1] = var1; // SL 2006.21.04 1 The order is now block dp instead of lp //hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")"; // SL 2006.21.04 2 // For decomp, lp ordering is needed. Nothing is changed. hilf[2] = "lp"; hilf[3] = di; resu[n] = hilf; } else { resu[n] = varstr(basering), ordstr(basering), 0; } } return(resu); } example { "EXAMPLE:"; echo = 2; ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp; ideal i = ea - fbg, fa + be, ec - fdg, fc + de; i = std(i); list l = newMaxIndependSetLp(i); l; i = i, g; l = newMaxIndependSetLp(i); l; ring s = 0, (x, y, z), lp; ideal i = z, yx; list l = newMaxIndependSetLp(i); l; } /////////////////////////////////////////////////////////////////////////////// proc newZero_decomp (ideal j, ideal ser, int @wr, list #) "USAGE: newZero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1 (@wr=0 for primary decomposition, @wr=1 for computation of associated primes) if #[1] = "nest", then #[2] indicates the nest level (number of recursive calls) When the nest level is high it indicates that the computation is difficult, and different methods are applied. RETURN: list = list of primary ideals and their radicals (at even positions in the list) if the input is zero-dimensional and a standardbases with respect to lex-ordering If ser!=(0) and ser is contained in j or if j is not zero-dimen- sional then ideal(1),ideal(1) is returned NOTE: Algorithm of Gianni/Trager/Zacharias EXAMPLE: example newZero_decomp; shows an example " { ASSUME(0, hasFieldCoefficient(basering) ); ASSUME(0, not isQuotientRing(basering) ) ; ASSUME(0, hasGlobalOrdering(basering) ) ; def @P = basering; int uytrewq; int nva = nvars(basering); int @k,@s,@n,@k1,@zz; list primary,lres0,lres1,act,@lh,@wh; map phi,psi,phi1,psi1; ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1; intvec @vh,@hilb; string @ri; poly @f; // Debug dbprint(printlevel - voice, "proc newZero_decomp"); if (dim(j)>0) { ERROR("dim(j)>0 . Please send the example to the authors"); } j=interred(j); attrib(j,"isSB",1); int nestLevel = 0; if (size(#) > 0) { if (typeof(#[1]) == "string") { if (#[1] == "nest") { nestLevel = #[2]; } # = list(); } } if(vdim(j)==deg(j[1])) { act=factor(j[1]); for(@k=1;@k<=size(act[1]);@k++) { @qh=j; if(@wr==0) { @qh[1]=act[1][@k]^act[2][@k]; } else { @qh[1]=act[1][@k]; } primary[2*@k-1]=interred(@qh); @qh=j; @qh[1]=act[1][@k]; primary[2*@k]=interred(@qh); attrib( primary[2*@k-1],"isSB",1); if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],5))==0)) { primary[2*@k-1]=ideal(1); primary[2*@k]=ideal(1); } } return(primary); } if(homog(j)==1) { primary[1]=j; if((size(ser)>0)&&(size(reduce(ser,j,5))==0)) { ERROR("dim(ser/j)==-1 unexpected. Please send the example to the authors"); } if(dim(j)==-1) { ERROR("dim(j)==-1 unexpected. Please send the example to the authors"); } else { primary[2]=maxideal(1); } return(primary); } //the first element in the standardbase is factorized if(deg(j[1])>0) { act=factor(j[1]); testFactor(act,j[1]); } else { ERROR("failure in newZero_decomp. Please send the example to the authors"); } //with the factors new ideals (hopefully the primary decomposition) //are created if(size(act[1])>1) { if(size(#)>1) { ERROR("failure in newZero_decomp. Please send the example to the authors"); } for(@k=1;@k<=size(act[1]);@k++) { if(@wr==0) { primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]); } else { primary[2*@k-1]=std(j,act[1][@k]); } if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k]))) { primary[2*@k] = primary[2*@k-1]; } else { primary[2*@k] = primaryTest(primary[2*@k-1],act[1][@k]); } } } else { primary[1]=j; if((size(#)>0)&&(act[2][1]>1)) { act[2]=1; primary[1]=std(primary[1],act[1][1]); } if(@wr!=0) { primary[1]=std(j,act[1][1]); } if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1]))) { primary[2]=primary[1]; } else { primary[2]=primaryTest(primary[1],act[1][1]); } } if(size(#)==0) { primary=splitPrimary(primary,ser,@wr,act); } if((voice>=7)&&(char(basering)<=181)) { primary=splitCharp(primary); } if((@wr==2)&&(npars(basering)>0)&&(voice>=7)&&(char(basering)>0)) { //the prime decomposition of Yokoyama in characteristic p list ke,ek; @k=0; while(@k 1){primary=extF(primary);} //test whether all ideals in the decomposition are primary and //in general position //if not after a random coordinate transformation of the last //variable the corresponding ideal is decomposed again. if((npars(basering)>0)&&(nestLevel > 1)) { poly randp; for(@zz=1;@zz0)&&(nestLevel > 1)) { jmap[size(jmap)]=randp; } for(@n=2;@n<=size(primary[2*@k-1]);@n++) { if(deg(lead(primary[2*@k-1][@n]))==1) { for(@zz=1;@zz<=nva;@zz++) { if(lead(primary[2*@k-1][@n])/var(@zz)!=0) { jmap1[@zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n] +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]); jmap2[@zz]=primary[2*@k-1][@n]; @qht[@n]=var(@zz); } } jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0); } } if(size(subst(jmap[nva],var(1),0)-var(nva))!=0) { // jmap[nva]=subst(jmap[nva],var(1),0); //hier geaendert +untersuchen!!!!!!!!!!!!!! } phi1=@P,jmap1; phi=@P,jmap; for(@n=1;@n<=nva;@n++) { jmap[@n]=-(jmap[@n]-2*var(@n)); } psi=@P,jmap; psi1=@P,jmap2; @qh=phi(@qht); //=================== the new part ============================ if (npars(basering)>1) { @qh=groebner(@qh,"par2var"); } else { @qh=groebner(@qh); } //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // ideal @qh=homog(imap(@P,@qht),@t); // // ideal @qh1=std(@qh); // @hilb=hilb(@qh1,1); // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // ideal @qh=homog(imap(@P,@qh),@t); // kill @Phelp; // @qh=std(@qh,@hilb); // @qh=subst(@qh,@t,1); // setring @P; // @qh=imap(@Phelp1,@qh); // kill @Phelp1; // @qh=clearSB(@qh); // attrib(@qh,"isSB",1); //============================================================= ser1=phi1(ser); @lh=newZero_decomp (@qh,phi(ser1),@wr, list("nest", nestLevel + 1)); kill lres0; list lres0; if(size(@lh)==2) { helpprim=@lh[2]; lres0[1]=primary[2*@k-1]; attrib(lres0[1],"isSB",1); ser1=psi(helpprim); lres0[2]=psi1(ser1); if(size(reduce(lres0[2],lres0[1],5))==0) { primary[2*@k]=primary[2*@k-1]; continue; } } else { lres1=psi(@lh); lres0=psi1(lres1); } //=================== the new part ============================ primary=delete(primary,2*@k-1); primary=delete(primary,2*@k-1); @k--; if(size(lres0)==2) { if (npars(basering)>1) { lres0[2]=groebner(lres0[2],"par2var"); } else { lres0[2]=groebner(lres0[2]); } } else { for(@n=1;@n<=size(lres0) div 2;@n++) { if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1) { lres0[2*@n-1]=groebner(lres0[2*@n-1]); lres0[2*@n]=lres0[2*@n-1]; attrib(lres0[2*@n],"isSB",1); } else { lres0[2*@n-1]=groebner(lres0[2*@n-1]); lres0[2*@n]=groebner(lres0[2*@n]); } } } primary=primary+lres0; //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // list @lvec; // list @lr=imap(@P,lres0); // ideal @lr1; // // if(size(@lr)==2) // { // @lr[2]=homog(@lr[2],@t); // @lr1=std(@lr[2]); // @lvec[2]=hilb(@lr1,1); // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lvec[2*@n]=@lvec[2*@n-1]; // } // else // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lr[2*@n]=homog(@lr[2*@n],@t); // @lr1=std(@lr[2*@n]); // @lvec[2*@n]=hilb(@lr1,1); // // } // } // } // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // list @lr=imap(@Phelp,@lr); // // kill @Phelp; // if(size(@lr)==2) // { // @lr[2]=std(@lr[2],@lvec[2]); // @lr[2]=subst(@lr[2],@t,1); // // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=@lr[2*@n-1]; // attrib(@lr[2*@n],"isSB",1); // } // else // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]); // @lr[2*@n]=subst(@lr[2*@n],@t,1); // } // } // } // kill @lvec; // setring @P; // lres0=imap(@Phelp1,@lr); // kill @Phelp1; // for(@n=1;@n<=size(lres0);@n++) // { // lres0[@n]=clearSB(lres0[@n]); // attrib(lres0[@n],"isSB",1); // } // // primary[2*@k-1]=lres0[1]; // primary[2*@k]=lres0[2]; // @s=size(primary) div 2; // for(@n=1;@n<=size(lres0) div 2-1;@n++) // { // primary[2*@s+2*@n-1]=lres0[2*@n+1]; // primary[2*@s+2*@n]=lres0[2*@n+2]; // } // @k--; //============================================================= } } return(primary); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; i=std(i); list pr= newZero_decomp(i,ideal(0),0); pr; } /////////////////////////////////////////////////////////////////////////////// static proc check_variables(alias ideal I, int NN) { int i=1; ideal J=0; ideal JJ; poly p; J=0; i=1; while(i1) {"found subsystem:",size(J)," of ",size(I)," vars:",size(variables(J)),size(variables(I)); return(J); } return(J); } static proc prep_decomp(alias ideal I) { int N=size(variables(I)); int NN=2; ideal J; list L; while(NN<=N) { J=check_variables(I,NN); if (size(J)>0) { L=minAssPrimes_i(0,J); if(size(L)>1) { "prep_decomp:",size(L)," NN=",NN," N=",N; for(int i=size(L);i>0;i--) { L[i]=simplify(L[i]+I,4); } return(L); } } NN++; } return(list(I)); } //////////////////////////////////////////////////////////////////////////// /* //Beispiele Wenk-Dipl (in ~/Texfiles/Diplom/Wenk/Examples/) //Zeiten: Singular/Singular/Singular -r123456789 -v :wilde13 (PentiumPro200) //Singular for HPUX-9 version 1-3-8 (2000060214) Jun 2 2000 15:31:26 //(wilde13) //1. vdim=20, 3 Komponenten //zerodec-time:2(1) (matrix:1 charpoly:0 factor:1) //primdecGTZ-time: 1(0) ring rs= 0,(a,b,c),dp; poly f1= a^2*b*c + a*b^2*c + a*b*c^2 + a*b*c + a*b + a*c + b*c; poly f2= a^2*b^2*c + a*b^2*c^2 + a^2*b*c + a*b*c + b*c + a + c; poly f3= a^2*b^2*c^2 + a^2*b^2*c + a*b^2*c + a*b*c + a*c + c + 1; ideal gls=f1,f2,f3; int time=timer; printlevel =1; time=timer; list pr1=zerodec(gls); timer-time;size(pr1); time=timer; list pr =primdecGTZ(gls); timer-time;size(pr); time=timer; ideal ra =radical(gls); timer-time;size(pr); //2.cyclic5 vdim=70, 20 Komponenten //zerodec-time:36(28) (matrix:1(0) charpoly:18(19) factor:17(9) //primdecGTZ-time: 28(5) //radical : 0 ring rs= 0,(a,b,c,d,e),dp; poly f0= a + b + c + d + e + 1; poly f1= a + b + c + d + e; poly f2= a*b + b*c + c*d + a*e + d*e; poly f3= a*b*c + b*c*d + a*b*e + a*d*e + c*d*e; poly f4= a*b*c*d + a*b*c*e + a*b*d*e + a*c*d*e + b*c*d*e; poly f5= a*b*c*d*e - 1; ideal gls= f1,f2,f3,f4,f5; //3. random vdim=40, 1 Komponente //zerodec-time:126(304) (matrix:1 charpoly:115(298) factor:10(5)) //primdecGTZ-time:17 (11) ring rs=0,(x,y,z),dp; poly f1=2*x^2 + 4*x + 3*y^2 + 7*x*z + 9*y*z + 5*z^2; poly f2=7*x^3 + 8*x*y + 12*y^2 + 18*x*z + 3*y^4*z + 10*z^3 + 12; poly f3=3*x^4 + 1*x*y*z + 6*y^3 + 3*x*z^2 + 2*y*z^2 + 4*z^2 + 5; ideal gls=f1,f2,f3; //4. introduction into resultants, sturmfels, vdim=28, 1 Komponente //zerodec-time:4 (matrix:0 charpoly:0 factor:4) //primdecGTZ-time:1 ring rs=0,(x,y),dp; poly f1= x4+y4-1; poly f2= x5y2-4x3y3+x2y5-1; ideal gls=f1,f2; //5. 3 quadratic equations with random coeffs, vdim=8, 1 Komponente //zerodec-time:0(0) (matrix:0 charpoly:0 factor:0) //primdecGTZ-time:1(0) ring rs=0,(x,y,z),dp; poly f1=2*x^2 + 4*x*y + 3*y^2 + 7*x*z + 9*y*z + 5*z^2 + 2; poly f2=7*x^2 + 8*x*y + 12*y^2 + 18*x*z + 3*y*z + 10*z^2 + 12; poly f3=3*x^2 + 1*x*y + 6*y^2 + 3*x*z + 2*y*z + 4*z^2 + 5; ideal gls=f1,f2,f3; //6. 3 polys vdim=24, 1 Komponente // run("ex14",2); //zerodec-time:5(4) (matrix:0 charpoly:3(3) factor:2(1)) //primdecGTZ-time:4 (2) ring rs=0,(x1,x2,x3,x4),dp; poly f1=16*x1^2 + 3*x2^2 + 5*x3^4 - 1 - 4*x4 + x4^3; poly f2=5*x1^3 + 3*x2^2 + 4*x3^2*x4 + 2*x1*x4 - 1 + x4 + 4*x1 + x2 + x3 + x4; poly f3=-4*x1^2 + x2^2 + x3^2 - 3 + x4^2 + 4*x1 + x2 + x3 + x4; poly f4=-4*x1 + x2 + x3 + x4; ideal gls=f1,f2,f3,f4; //7. ex43, PoSSo, caprasse vdim=56, 16 Komponenten //zerodec-time:23(15) (matrix:0 charpoly:16(13) factor:3(2)) //primdecGTZ-time:3 (2) ring rs= 0,(y,z,x,t),dp; ideal gls=y^2*z+2*y*x*t-z-2*x, 4*y^2*z*x-z*x^3+2*y^3*t+4*y*x^2*t-10*y^2+4*z*x+4*x^2-10*y*t+2, 2*y*z*t+x*t^2-2*z-x, -z^3*x+4*y*z^2*t+4*z*x*t^2+2*y*t^3+4*z^2+4*z*x-10*y*t-10*t^2+2; //8. Arnborg-System, n=6 (II), vdim=156, 90 Komponenten //zerodec-time (char32003):127(45)(matrix:2(0) charpoly:106(37) factor:16(7)) //primdecGTZ-time(char32003) :81 (18) //ring rs= 0,(a,b,c,d,x,f),dp; ring rs= 32003,(a,b,c,d,x,f),dp; ideal gls=a+b+c+d+x+f, ab+bc+cd+dx+xf+af, abc+bcd+cdx+d*xf+axf+abf, abcd+bcdx+cd*xf+ad*xf+abxf+abcf, abcdx+bcd*xf+acd*xf+abd*xf+abcxf+abcdf, abcd*xf-1; //9. ex42, PoSSo, Methan6_1, vdim=27, 2 Komponenten //zerodec-time:610 (matrix:10 charpoly:557 factor:26) //primdecGTZ-time: 118 //zerodec-time(char32003):2 //primdecGTZ-time(char32003):4 //ring rs= 0,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp; ring rs= 32003,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp; ideal gls=64*x2*x7-10*x1*x8+10*x7*x9+11*x7*x10-320000*x1, -32*x2*x7-5*x2*x8-5*x2*x10+160000*x1-5000*x2, -x3*x8+x6*x8+x9*x10+210*x6+1300000, -x4*x8+700000, x10^2-2*x5, -x6*x8+x7*x9-210*x6, -64*x2*x7-10*x7*x9-11*x7*x10+320000*x1-16*x7+7000000, -10*x1*x8-10*x2*x8-10*x3*x8-10*x4*x8-10*x6*x8+10*x2*x10+11*x7*x10 +20000*x2+14*x5, x4*x8-x7*x9-x9*x10-410*x9, 10*x2*x8+10*x3*x8+10*x6*x8+10*x7*x9-10*x2*x10-11*x7*x10-10*x9*x10 -10*x10^2+1400*x6-4200*x10; //10. ex33, walk-s7, Diplomarbeit von Tim, vdim=114 //zerfaellt in unterschiedlich viele Komponenten in versch. Charkteristiken: //char32003:30, char0:3(2xdeg1,1xdeg112!), char181:4(2xdeg1,1xdeg28,1xdeg84) //char 0: zerodec-time:10075 (ca 3h) (matrix:3 charpoly:9367, factor:680 // + 24 sec fuer Normalform (anstatt einsetzen), total [29623k]) // primdecGTZ-time: 214 //char 32003:zerodec-time:197(68) (matrix:2(1) charpoly:173(60) factor:15(6)) // primdecGTZ-time:14 (5) //char 181:zerodec-time:(87) (matrix:(1) charpoly:(58) factor:(25)) // primdecGTZ-time:(2) //in char181 stimmen Ergebnisse von zerodec und primdecGTZ ueberein (gecheckt) //ring rs= 0,(a,b,c,d,e,f,g),dp; ring rs= 32003,(a,b,c,d,e,f,g),dp; poly f1= 2gb + 2fc + 2ed + a2 + a; poly f2= 2gc + 2fd + e2 + 2ba + b; poly f3= 2gd + 2fe + 2ca + c + b2; poly f4= 2ge + f2 + 2da + d + 2cb; poly f5= 2fg + 2ea + e + 2db + c2; poly f6= g2 + 2fa + f + 2eb + 2dc; poly f7= 2ga + g + 2fb + 2ec + d2; ideal gls= f1,f2,f3,f4,f5,f6,f7; ~/Singular/Singular/Singular -r123456789 -v LIB"./primdec.lib"; timer=1; int time=timer; printlevel =1; option(prot,mem); time=timer; list pr1=zerodec(gls); timer-time; time=timer; list pr =primdecGTZ(gls); timer-time; time=timer; list pr =primdecSY(gls); timer-time; time=timer; ideal ra =radical(gls); timer-time;size(pr); LIB"all.lib"; ring R=0,(a,b,c,d,e,f),dp; ideal I=cyclic(6); minAssGTZ(I); ring S=(2,a,b),(x,y),lp; ideal I=x8-b,y4+a; minAssGTZ(I); ring S1=2,(x,y,a,b),lp; ideal I=x8-b,y4+a; minAssGTZ(I); ring S2=(2,z),(x,y),dp; minpoly=z2+z+1; ideal I=y3+y+1,x4+x+1; primdecGTZ(I); minAssGTZ(I); ring S3=2,(x,y,z),dp; ideal I=y3+y+1,x4+x+1,z2+z+1; primdecGTZ(I); minAssGTZ(I); ring R1=2,(x,y,z),lp; ideal I=y6+y5+y3+y2+1,x4+x+1,z2+z+1; primdecGTZ(I); minAssGTZ(I); ring R2=(2,z),(x,y),lp; minpoly=z3+z+1; ideal I=y2+y+(z2+z+1),x4+x+1; primdecGTZ(I); minAssGTZ(I); */