// $Id: primdec.lib,v 1.13 1998-05-05 11:55:34 krueger Exp $ /////////////////////////////////////////////////////// // primdec.lib // algorithms for primary decomposition based on // the ideas of Gianni,Trager,Zacharias // written by Gerhard Pfister // // algorithms for primary decomposition based on // the ideas of Shimoyama/Yokoyama // written by Wolfram Decker and Hans Schoenemann ////////////////////////////////////////////////////// version="$Id: primdec.lib,v 1.13 1998-05-05 11:55:34 krueger Exp $"; info=" LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I) minAssPrimes (ideal I) //computes the minimal associated primes of the ideal I primdecGTZ (ideal I) // Computes a complete primary decomposition via Gianni,Trager,Zacharias radical(ideal I) //computes the radical of the ideal I equiRadical(ideal I) //computes the radical of the equidimensional part of the ideal I prepareAss(ideal I) //computes the radicals of the equidimensional parts of I min_ass_prim_charsets (ideal I, int choose) // minimal associated primes via the characteristic set // package written by Michael Messollen. // The integer choose must be either 0 or 1. // If choose=0, the given ordering of the variables is used. // If choose=1, the system tries to find an \"optimal ordering\", // which in some cases may considerably speed up the algorithm // You may also may want to try one of the algorithms for // minimal associated primes in the the library // primdec.lib, written by Gerhard Pfister. // These algorithms are variants of the algorithm // of Gianni-Trager-Zacharias primdecSY (ideal I, int choose) // Computes a complete primary decomposition via // a variant of the pseudoprimary approach of // Shimoyama-Yokoyama. // The integer choose must be either 0, 1, 2 or 3. // 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 "; LIB "general.lib"; LIB "elim.lib"; LIB "poly.lib"; LIB "random.lib"; /////////////////////////////////////////////////////////////////////////////// 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 EXAMPLE: example sat; shows an example " { int @k; ideal inew=std(id); ideal iold; option(returnSB); while(specialIdealsEqual(iold,inew)==0 ) { iold=inew; inew=quotient(iold,p); @k++; } @k--; option(noreturnSB); list L =inew,p^@k; return (L); } /////////////////////////////////////////////////////////////////////////////// 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 EXAMPLE: example sat2; shows an example " { int @k,@i; def @P= basering; if(ordstr(basering)[1,2]!="dp") { execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),(C,dp);"; 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,4); poly @f=1; if(deg(fac[1])>0) { ideal iold; for(@i=1;@i<=size(fac);@i++) { @f=@f*fac[@i]; } option(returnSB); while(specialIdealsEqual(iold,inew)==0 ) { iold=inew; if(deg(iold[size(iold)])!=1) { inew=quotient(iold,@f); } else { inew=iold; } @k++; } option(noreturnSB); @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) { 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,4); 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; option(returnSB); while(specialIdealsEqual(iold,quotM[1])==0) { if(k>0) { f=f*quotM[3]; } iold=quotM[1]; quotM=quotMin(quotM); k++; } option(noreturnSB); l=quotM[1],f; return(l); } proc quotMin(list tsil) { 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); 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); } //////////////////////////////////////////////////////////////////////////////// proc testFactor(list act,poly p) { 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"; basering; act; keep; pause; p; q; pause; } } //////////////////////////////////////////////////////////////////////////////// proc factor(poly p) "USAGE: factor(p) p poly RETURN: list=; NOTE: EXAMPLE: example factor; shows an example " { ideal @i; list @l; intvec @v,@w; int @j,@k,@n; if(deg(p)<=1) { @i=ideal(p); @v=1; @l[1]=@i; @l[2]=@v; return(@l); } if (size(p)==1) { @w=leadexp(p); for(@j=1;@j<=nvars(basering);@j++) { if(@w[@j]!=0) { @k++; @v[@k]=@w[@j]; @i=@i+ideal(var(@j)); } } @l[1]=@i; @l[2]=@v; return(@l); } // @l=factorize(p,2); @l=factorize(p); // if(npars(basering)>0) // { for(@j=1;@j<=size(@l[1]);@j++) { if(deg(@l[1][@j])==0) { @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(deg(@l[1][@j])!=0) { @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))); } 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); } 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,1))==0) { return(1); } } return(0); } //////////////////////////////////////////////////////////////////////////////// proc testPrimary(list pr, ideal k) "USAGE: testPrimary(pr,k) pr list, k ideal; RETURN: int = 1, if the intersection of the ideals in pr is k, 0 if not NOTE: EEXAMPLE: example testPrimary ; shows an example " { int i; ideal j=pr[1]; for (i=2;i<=size(pr)/2;i++) { j=intersect(j,pr[2*i-1]); } return(idealsEqual(j,k)); } example { "EXAMPLE:"; echo = 2; ring s = 0,(x,y,z),lp; ideal i=x3-x2-x+1,xy-y; ideal i1=x-1; ideal i2=x-1; ideal i3=y,x2-2x+1; ideal i4=y,x-1; ideal i5=y,x+1; ideal i6=y,x+1; list pr=i1,i2,i3,i4,i5,i6; testPrimary(pr,i); pr[5]=y+1,x+1; testPrimary(pr,i); } //////////////////////////////////////////////////////////////////////////////// proc printPrimary( list l, list #) "USAGE: printPrimary(l) l list; RETURN: nothing NOTE: EXAMPLE: example printPrimary; shows an example " { if(size(#)>0) { " "; " The primary decomposition of the ideal "; #[1]; " "; " is: "; " "; } int k; for (k=1;k<=size(l)/2;k=k+1) { " "; "primary ideal: "; l[2*k-1]; " "; "associated prime ideal "; l[2*k]; " "; } } example { "EXAMPLE:"; echo = 2; } //////////////////////////////////////////////////////////////////////////////// proc randomLast(int b) "USAGE: randomLast RETURN: ideal = maxideal(1) but the last variable exchanged by a sum of it with a linear random combination of the other variables NOTE: EXAMPLE: example randomLast; shows an example " { ideal i=maxideal(1); int k=size(i); i[k]=0; i=randomid(i,size(i),b); ideal ires=maxideal(1); ires[k]=i[1]+var(k); return(ires); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; ideal i = randomLast(10); i; } //////////////////////////////////////////////////////////////////////////////// proc primaryTest (ideal i, poly p) { int m=1; int n=nvars(basering); int e; poly t; ideal h; ideal prm=p; attrib(prm,"isSB",1); while (n>1) { n=n-1; m=m+1; //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+1; } m=m-1; } //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])); // t=hilfe1(i,prm,m,n); 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,1) !=0) { return(ideal(0)); } h=interred(t); t=h[1]; prm = prm,t; attrib(prm,"isSB",1); } return(prm); } proc hilfe(ideal i,ideal prm,int m) { poly t; int e; if(size(i[m])==1) { t=var(n); } else { e=deg(lead(i[m])); if(deg(poly(e))>=0) { 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]; } else { i[m]=i[m]/leadcoef(i[m]); t=reduce(coef(i[m],var(n))[2,e+1],prm); t=var(n)+factorize(t,1)[1]; } } return(t); } proc hilfe1(ideal i,ideal prm,int m,int n) { poly t; int e; if(size(i[m])==1) { t=var(n); } else { e=deg(lead(i[m])); i[m]=i[m]/leadcoef(i[m]); t=reduce(coeffs(i[m],var(n))[1,1],prm); if(size(t)==0){return(var(n));} t=var(n)+factorize(t,1)[1]; } return(t); } /////////////////////////////////////////////////////////////////////////////// proc splitPrimary(list l,ideal ser,int @wr,list sact) { int i,j,k,s,r,w; list keepresult,act,keepprime; poly @f; int sl=size(l); for(i=1;i<=sl/2;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],1))==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/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)/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/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/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+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/2+1]=interred(keepprime[i]+ideal(@f)); if(homog(keepprime[s/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)/2;i++) { attrib(l[2*i-1],"isSB",1); if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==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; LIB "primdec.lib"; 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; } //////////////////////////////////////////////////////////////////////////////// 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 computaion 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, Traeger, Zacharias EXAMPLE: example zero_decomp; shows an example " { def @P = basering; int nva = nvars(basering); int @k,@s,@n,@k1,zz; list primary,lres,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); } 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],1))==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,1))==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); } //withe 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((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); } //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. @k=0; while(@k<(size(primary)/2)) { @k++; if (size(primary[2*@k])==0) { for(zz=1;zz0) { @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); // @lh=zero_decomp (@qh,psi(ser),@wr); kill lres; list lres; if(size(@lh)==2) { helpprim=@lh[2]; lres[1]=primary[2*@k-1]; ser1=psi(helpprim); lres[2]=psi1(ser1); if(size(reduce(lres[2],lres[1],1))==0) { primary[2*@k]=primary[2*@k-1]; continue; } } else { act=factor(@qh[1]); if(2*size(act[1])==size(@lh)) { for(@n=1;@n<=size(act[1]);@n++) { @f=act[1][@n]^act[2][@n]; ser1=psi(@f); lres[2*@n-1]=interred(primary[2*@k-1]+psi1(ser1)); helpprim=@lh[2*@n]; ser1=psi(helpprim); lres[2*@n]=psi1(ser1); } } else { lres1=psi(@lh); lres=psi1(lres1); } } 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,lres); 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)/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)/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; lres=imap(@Phelp1,@lr); kill @Phelp1; for(@n=1;@n<=size(lres);@n++) { lres[@n]=clearSB(lres[@n]); attrib(lres[@n],"isSB",1); } primary[2*@k-1]=lres[1]; primary[2*@k]=lres[2]; @s=size(primary)/2; for(@n=1;@n<=size(lres)/2-1;@n++) { primary[2*@s+2*@n-1]=lres[2*@n+1]; primary[2*@s+2*@n]=lres[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 ggt (ideal i) "USAGE: ggt(i); i list of polynomials RETURN: poly = ggt(i[1],...,i[size(i)]) NOTE: EXAMPLE: example ggt; shows an example " { int k; poly p=i[1]; if(deg(p)==0) { return(1); } for (k=2;k<=size(i);k++) { if(deg(i[k])==0) { return(1) } p=GCD(p,i[k]); if(deg(p)==0) { return(1); } } return(p); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = (x+y)*(y+z); poly q = (z4+2)*(y+z); ideal l=p,q; poly pr= ggt(l); pr; } /////////////////////////////////////////////////////////////////////////////// proc gcdTest(ideal act) { int i,j; if(size(act)<=1) { return(0); } for (i=1;i<=size(act)-1;i++) { for(j=i+1;j<=size(act);j++) { if(deg(std(ideal(act[i],act[j]))[1])>0) { return(0); } } } return(1); } /////////////////////////////////////////////////////////////////////////////// proc coeffLcm(ideal h) { string @pa=parstr(basering); if(size(@pa)==0) { return(lcmP(h)); } def bsr= basering; string @id=string(h); execute "ring @r=0,("+@pa+","+varstr(bsr)+"),(C,dp);"; execute "ideal @i="+@id+";"; poly @p=lcmP(@i); string @ps=string(@p); setring bsr; execute "poly @p="+@ps+";"; return(@p); } example { "EXAMPLE:"; echo = 2; ring r =( 0,a,b),(x,y,z),lp; poly p = (a+b)*(y-z); poly q = (a+b)*(y+z); ideal l=p,q; poly pr= coeffLcm(l); pr; } /////////////////////////////////////////////////////////////////////////////// proc lcmP(ideal i) "USAGE: lcm(i); i list of polynomials RETURN: poly = lcm(i[1],...,i[size(i)]) NOTE: EXAMPLE: example lcm; shows an example " { int k,j; poly p,q; i=simplify(i,10); for(j=1;j<=size(i);j++) { if(deg(i[j])>0) { p=i[j]; break; } } if(deg(p)==-1) { return(1); } for (k=j+1;k<=size(i);k++) { if(deg(i[k])!=0) { q=GCD(p,i[k]); if(deg(q)==0) { p=p*i[k]; } else { p=p/q; p=p*i[k]; } } } return(p); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = (x+y)*(y+z); poly q = (z4+2)*(y+z); ideal l=p,q; poly pr= lcmP(l); pr; l=1,-1,p,1,-1,q,1; pr=lcmP(l); pr; } /////////////////////////////////////////////////////////////////////////////// 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 " { 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; LIB "primdec.lib"; 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; } /////////////////////////////////////////////////////////////////////////////// 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=system("indsetall",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; } /////////////////////////////////////////////////////////////////////////////// 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 " { int n,k,di; list resu,hilf; string var1,var2; list v=system("indsetall",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; } /////////////////////////////////////////////////////////////////////////////// proc prepareQuotientring (int nnp) "USAGE: prepareQuotientring(nnp); nnp int RETURN: string = to define Kvar(nnp+1),...,var(nvars)[..rest ] NOTE: EXAMPLE: example independentSet; shows an example " { ideal @ih,@jh; int npar=npars(basering); int @n; string quotring= "ring quring = ("+charstr(basering); for(@n=nnp+1;@n<=nvars(basering);@n++) { quotring=quotring+",var("+string(@n)+")"; @ih=@ih+var(@n); } quotring=quotring+"),(var(1)"; @jh=@jh+var(1); for(@n=2;@n<=nnp;@n++) { quotring=quotring+",var("+string(@n)+")"; @jh=@jh+var(@n); } quotring=quotring+"),(C,lp);"; return(quotring); } example { "EXAMPLE:"; echo = 2; ring s1=(0,x),(a,b,c,d,e,f,g),lp; def @Q=basering; list l= prepareQuotientring(3); l; execute l[1]; execute l[2]; basering; phi; setring @Q; } /////////////////////////////////////////////////////////////////////// proc projdim(list l) { int i=size(l)+1; while(i>2) { i--; if((size(l[i])>0)&&(deg(l[i][1])>0)) { return(i); } } } /////////////////////////////////////////////////////////////////////////////// proc cleanPrimary(list l) { int i,j; list lh; for(i=1;i<=size(l)/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 minAssPrimes(ideal i, list #) "USAGE: minAssPrimes(i); i ideal minAssPrimes(i,1); i ideal (to use also the factorizing Groebner) RETURN: list = the minimal associated prime ideals of i NOTE: EXAMPLE: example minAssPrimes; shows an example " { #[1]=1; def @P=basering; list qr=simplifyIdeal(i); map phi=@P,qr[2]; i=qr[1]; execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(" +ordstr(basering)+");"; ideal i=fetch(@P,i); if(size(#)==0) { int @wr; list tluser,@res; list primary=decomp(i,2); @res[1]=primary; tluser=union(@res); setring @P; list @res=imap(gnir,tluser); return(phi(@res)); } list @res,empty; ideal ser; 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); 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;k<=j-1;k++) { pos[k]=1; } count=j; odim=ndim; } if(ndimsize(@erg)) { @wos[size(@wos)+1]=@erg[i]; } } } if(deg(@erg[size(@erg)][1])!=0) { @wos[size(@wos)+1]=@erg[size(@erg)]; } setring P; list @ser=fetch(ir,@wos); return(@ser); } /////////////////////////////////////////////////////////////////////////////// proc radicalOld(ideal i) { list pr=minAssPrimes(i,1); int j; ideal k=pr[1]; for(j=2;j<=size(pr);j++) { k=intersect(k,pr[j]); } return(k); } /////////////////////////////////////////////////////////////////////////////// proc equiRadical(ideal i) { return(radical(i,1)); } /////////////////////////////////////////////////////////////////////////////// proc decomp(ideal i,list #) "USAGE: decomp(i); i ideal (for primary decomposition) (resp. decomp(i,1); (for the minimal associated primes) ) 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, Traeger, Zacharias EXAMPLE: example decomp; shows an example " { def @P = basering; list primary,indep,ltras; intvec @vh,isat; int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi; ideal peek=i; ideal ser,tras; if(size(#)>0) { if((#[1]==1)||(#[1]==2)) { @wr=#[1]; if(size(#)>1) { seri=1; peek=#[2]; ser=#[3]; } } else { seri=1; peek=#[1]; ser=#[2]; } } homo=homog(i); if(homo==1) { if(attrib(i,"isSB")!=1) { ltras=mstd(i); attrib(ltras[1],"isSB",1); } else { ltras=i,i; } tras=ltras[1]; if(dim(tras)==0) { primary[1]=ltras[2]; primary[2]=maxideal(1); if(@wr>0) { list l; l[1]=maxideal(1); l[2]=maxideal(1); return(l); } return(primary); } intvec @hilb=hilb(tras,1); intvec keephilb=@hilb; } //---------------------------------------------------------------- //i is the zero-ideal //---------------------------------------------------------------- if(size(i)==0) { primary=i,i; return(primary); } //---------------------------------------------------------------- //pass to the lexicographical ordering and compute a standardbasis //---------------------------------------------------------------- execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);"; option(redSB); ideal ser=fetch(@P,ser); if(homo==1) { if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)")) { ideal @j=std(fetch(@P,i),@hilb); } else { ideal @j=fetch(@P,tras); attrib(@j,"isSB",1); } } else { ideal @j=std(fetch(@P,i)); } option(noredSB); if(seri==1) { ideal peek=fetch(@P,peek); attrib(peek,"isSB",1); } else { ideal peek=@j; } if(size(ser)==0) { 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)>0) { @j=simplify(@j,2); attrib(@j,"isSB",1); list pr=decomp(@j); for(@k=1;@k<=size(pr);@k++) { @j=pr[@k]+fried; pr[@k]=@j; } setring @P; return(fetch(gnir,pr)); } } //---------------------------------------------------------------- //j is the ring //---------------------------------------------------------------- if (dim(@j)==-1) { setring @P; return(ideal(0)); } //---------------------------------------------------------------- // 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; primary=fetch(gnir,gprimary); return(primary); } //------------------------------------------------------------------ //the zero-dimensional case //------------------------------------------------------------------ if (dim(@j)==0) { option(redSB); list gprimary= zero_decomp(@j,ser,@wr); option(noredSB); setring @P; primary=fetch(gnir,gprimary); if(size(ser)>0) { primary=cleanPrimary(primary); } return(primary); } poly @gs,@gh,@p; string @va,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 the 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") { execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");"; } else { execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);"; } 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); } } else { ideal jwork=std(imap(gnir,@j)); } list hquprimary; poly @p,@q; ideal @h,fac,ser; di=dim(jwork); keepdi=di; setring gnir; for(@m=1;@m<=size(indep);@m++) { isat=0; @n2=0; option(redSB); 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 { execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),(" +ordstr(basering)+");"; ideal @j=fetch(gnir,@j); attrib(@j,"isSB",1); ideal ser=fetch(gnir,ser); } else { @va=string(maxideal(1)); execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),(" +indep[@m][2]+");"; execute "map phi=gnir,"+@va+";"; if(homo==1) { ideal @j=std(phi(@j),@hilb); } else { ideal @j=std(phi(@j)); } ideal ser=phi(ser); } option(noredSB); if((deg(@j[1])==0)||(dim(@j)K(var(nnpr+1),..,var(nva))[..the rest..] //------------------------------------------------------------------------------------- quotring=prepareQuotientring(nvars(basering)-indep[@m][3]); //--------------------------------------------------------------------- //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //--------------------------------------------------------------------- execute 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; 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..] option(redSB); list uprimary= zero_decomp(@j,ser,@wr); option(noredSB); } 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 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))[..the rest..] //back to the polynomialring //--------------------------------------------------------------------- setring gnir; collectprimary=imap(quring,uprimary); lsau=imap(quring,saturn); @h=imap(quring,@h); kill quring; @n2=size(quprimary); @n3=@n2; for(@n1=1;@n1<=size(collectprimary)/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]; } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3/2+1;@n<=@n2/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) { @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,4); @q=1; for(lauf=1;lauf<=size(fac);lauf++) { @q=@q*fac[lauf]; } } jwork=std(jwork,@q); keepdi=dim(jwork); if(keepdi0)) { @j=ideal(1); quprimary[1]=ideal(1); quprimary[2]=ideal(1); } if((size(quprimary)==0)) { keepdi=di-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) { ideal htest=hquprimary[1]; for (@n1=2;@n1<=size(hquprimary)/2;@n1++) { htest=intersect(htest,hquprimary[2*@n1-1]); } } else { ideal htest=hquprimary[2]; for (@n1=2;@n1<=size(hquprimary)/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,1))!=0) { for(@m=1;@m<=size(restindep);@m++) { // if(restindep[@m][3]>=keepdi) // { isat=0; @n2=0; option(redSB); if(restindep[@m][1]==varstr(basering)) //this is the good case, nothing to do, just to have the same notations //change the ring { execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),(" +ordstr(basering)+");"; 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+";"; if(homo==1) { ideal @j=std(phi(jkeep),keephilb); } else { ideal @j=std(phi(jkeep)); } ideal ser=phi(ser); } option(noredSB); 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..] //------------------------------------------------------------------------------------- quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]); //--------------------------------------------------------------------- //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //--------------------------------------------------------------------- execute 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))[..the rest..] option(redSB); list uprimary= zero_decomp(@j,ser,@wr); option(noredSB); //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 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))[..the rest..] //back to the polynomialring //--------------------------------------------------------------------- setring gnir; collectprimary=imap(quring,uprimary); lsau=imap(quring,saturn); @h=imap(quring,@h); kill quring; @n2=size(quprimary); @n3=@n2; for(@n1=1;@n1<=size(collectprimary)/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]; } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3/2+1;@n<=@n2/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/2+1;@n<=@n2/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(size(reduce(ser,peek,1))!=0) { if(@wr>0) { htprimary=decomp(@j,@wr,peek,ser); } else { htprimary=decomp(@j,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]; } } } } //------------------------------------------------------------ //back to the ring we started with //the final result: primary //------------------------------------------------------------ setring @P; primary=imap(gnir,quprimary); 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; LIB "primdec.lib"; list pr= decomp(i); pr; testPrimary( pr, i); } /////////////////////////////////////////////////////////////////////////////// proc primdecGTZ(ideal i) { return(decomp(i)); } /////////////////////////////////////////////////////////////////////////////// proc primdecSY(ideal i,int j) { return(prim_dec(i,j)); } /////////////////////////////////////////////////////////////////////////////// proc radical(ideal i,list #) { def @P=basering; int j,il; if(size(#)>0) { il=#[1]; } ideal re=1; option(redSB); list qr=simplifyIdeal(i); map phi=@P,qr[2]; i=qr[1]; list pr=facstd(i); if(size(pr)==1) { attrib(pr[1],"isSB",1); if((dim(pr[1])==0)&&(homog(pr[1])==1)) { ideal @res=maxideal(1); return(phi(@res)); } if(dim(pr[1])>1) { execute "ring gnir = ("+charstr(basering)+"), ("+varstr(basering)+"),(C,lp);"; ideal i=fetch(@P,i); list @pr=facstd(i); setring @P; pr=fetch(gnir,@pr); } } option(noredSB); int s=size(pr); if(s==1) { i=radicalEHV(i,ideal(1),il); return(phi(i)); } intvec pos; pos[s]=0; if(il==1) { int ndim,k; attrib(pr[1],"isSB",1); int odim=dim(pr[1]); int count=1; for(j=2;j<=s;j++) { attrib(pr[j],"isSB",1); ndim=dim(pr[j]); if(ndim>odim) { for(k=count;k<=j-1;k++) { pos[k]=1; } count=j; odim=ndim; } if(ndim0) { @wr=#[1]; } @j=m[1]; @j1=m[2]; int jdim=dim(@j); if(size(reduce(ser,@j,1))==0) { return(ser); } if(homo==1) { if(jdim==0) { option(noredSB); return(maxideal(1)); } intvec @hilb=hilb(@j,1); } //---------------------------------------------------------------- //j is the ring //---------------------------------------------------------------- if (jdim==-1) { option(noredSB); return(ideal(0)); } //---------------------------------------------------------------- // the case of one variable //---------------------------------------------------------------- if(nvars(basering)==1) { fac=factorize(@j[1],1); @p=1; for(@k=1;@k<=size(fac);@k++) { @p=@p*fac[@k]; } option(noredSB); return(ideal(@p)); } //------------------------------------------------------------------ //the case of a complete intersection //------------------------------------------------------------------ if(jdim+size(@j1)==nvars(basering)) { // ideal jac=minor(jacob(@j1),size(@j1)); // return(quotient(@j1,jac)); } //------------------------------------------------------------------ //the zero-dimensional case //------------------------------------------------------------------ if (jdim==0) { @j1=finduni(@j); for(@k=1;@k<=size(@j1);@k++) { fac=factorize(cleardenom(@j1[@k]),1); @p=fac[1]; for(@n=2;@n<=size(fac);@n++) { @p=@p*fac[@n]; } @j=@j,@p; } @j=std(@j); option(noredSB); return(@j); } //------------------------------------------------------------------ //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 the block-ordering //------------------------------------------------------------------ indep=maxIndependSet(@j); di=dim(@j); 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 { execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),(" +ordstr(basering)+");"; ideal @j=fetch(@P,@j); 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(@j),@hilb); } else { ideal @j=std(phi(@j)); } } if((deg(@j[1])==0)||(dim(@j)K(var(nnpr+1),..,var(nva))[..the rest..] //------------------------------------------------------------------------------------- quotring=prepareQuotientring(nvars(basering)-indep[@m][3]); //--------------------------------------------------------------------- //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //--------------------------------------------------------------------- execute 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); 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..] option(redSB); @j=interred(@j); attrib(@j,"isSB",1); list @mo=@j,@j; ideal zero_rad= radicalKL(@mo,ideal(1)); } 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 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(quring,zero_rad); lsau=simplify(imap(quring,hpl),2); @h=imap(quring,@h); kill quring; //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,4); @q=1; for(lauf=1;lauf<=size(fac);lauf++) { @q=@q*fac[lauf]; } @mu=mstd(quotient(@j+ideal(@q),rad)); @j=@mu[1]; attrib(@j,"isSB",1); } if((deg(rad[1])>0)&&(deg(collectrad[1])>0)) { rad=intersect(rad,collectrad); } else { if(deg(collectrad[1])>0) { rad=collectrad; } } te=simplify(reduce(te*rad,@j),2); if((dim(@j)0) { il=#[1]; } option(redSB); list m=mstd(i); I=m[2]; option(noredSB); if(size(reduce(re,m[1],1))==0) { return(re); } int cod=nvars(basering)-dim(m[1]); if((nvars(basering)<=5)&&(size(m[2])<=5)) { if(cod==size(m[2])) { J=minor(jacob(I),cod); return(quotient(I,J)); } for(l=1;l<=cod;l++) { I0[l]=I[l]; } if(dim(std(I0))+cod==nvars(basering)) { J=minor(jacob(I0),cod); radI0=quotient(I0,J); L=quotient(radI0,I); radI1=quotient(radI0,L); if(size(reduce(radI1,m[1],1))==0) { return(I); } if(il==1) { return(radI1); } I2=sat(I,radI1)[1]; if(deg(I2[1])<=0) { return(radI1); } return(intersect(radI1,radicalEHV(I2,re,il))); } } return(radicalKL(m,re,il)); } proc Ann(module M) { M=prune(M); //to obtain a small embedding ideal ann=quotient1(M,freemodule(nrows(M))); return(ann); } //computes the equidimensional part of the ideal i of codimension e proc int_ass_primary_e(ideal i, int e) { 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 all equidimensional parts of the ideal i proc prepareAss(ideal i) { ideal j=std(i); int cod=nvars(basering)-dim(j); int e; list er; ideal ann; if(homog(i)==1) { list re=sres(i,0); //the resolution re=minres(re); //minimized resolution } else { list re=mres(i,0); //fehler in sres } 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); } //computes the annihilator of Ext^n(R/i,R) with given resolution re //n is not necessarily the number of variables proc AnnExt_R(int n,list re) { if(n1;ii--) { for(jj=1;jj0; n=n-1) { 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 ////////////////////////////////////////////////////// proc min_ass_prim_charsets (ideal PS, int cho) { if((cho<0) and (cho>1)) { "ERROR: must be 0 or 1" return(); } if(system("version")>933) { option(notWarnSB); } if(cho==0) { return(min_ass_prim_charsets0(PS)); } else { return(min_ass_prim_charsets1(PS)); } } /////////////////////////////////////////////////////// // 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 ////////////////////////////////////////////////////// proc min_ass_prim_charsets0 (ideal PS) { matrix m=char_series(PS); // We compute an irreducible // characteristic series 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=j-1) { I=I+ini_mod(PHI[i][j]); } JS=std(PHI[i]); sizeJS=size(JS); for(j=size(I);j>0;j=j-1) { II=0; sizeII=0; k=0; while(k<=sizeII) // successive saturation { option(returnSB); II=quotient(JS,I[j]); option(noreturnSB); //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); } } 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 ////////////////////////////////////////////////////// proc min_ass_prim_charsets1 (ideal PS) { 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 int i,j,k; ideal I; list PSI; list PHI; // the ideals given by the characteristic series list ITPHI; // their initial terms for(i=nrows(m);i>=1; i--) { PHI[i]=ideal(m[i,1..ncols(m)]); I=0; for(j=size(PHI[i]);j>0;j=j-1) { I=I,ini_mod(PHI[i][j]); } 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=j-1) { I=I+ITPHI[i][j]; } JS=std(PHI[i]); sizeJS=size(JS); for(j=size(I);j>0;j=j-1) { II=0; sizeII=0; k=0; while(k<=sizeII) // successive iteration { option(returnSB); II=quotient(JS,I[j]); option(noreturnSB); //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); } } 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. ////////////////////////////////////////////////////////// proc prim_dec(ideal I, int choose) { if((choose<0) or (choose>3)) { "ERROR: must be 0 or 1 or 2 or 3"; return(); } if(system("version")>933) { option(notWarnSB); } 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); 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 { 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(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,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)); } } } } } ////////////////////////////////////////////////////////////////////////// // 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 ////////////////////////////////////////////////////////////////////////// proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo) { 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,choo); } else { if(choo==2) { L=minAssPrimes(I); } else { L=minAssPrimes(I,1); } for(int i=size(L);i>=1;i=i-1) { L[i]=std(L[i]); } } return (pseudo_prim_dec_i(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 //////////////////////////////////////////////////////////////// proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo) { 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(SP,choo); // a list of SB } else { if(choo==2) { m=minAssPrimes(SP); } else { m=minAssPrimes(SP,1); } for(j=size(m);j>=1;j=j-1) { 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(SI,L)); } //////////////////////////////////////////////////////////////// // proc pseudo_prim_dec_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, //////////////////////////////////////////////////////////////// proc pseudo_prim_dec_i (ideal SI, list L) { 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); // 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 //////////////////////////////////////////////////////////////// proc extraction (ideal SI, ideal SP) { list indsets=system("indsetall",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 ////////////////////////////////////////////////////////// proc minsat_ppd(ideal SI, ideal fac) { 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 ///////////////////////////////////////////////////////////////// proc minquot(list tsil) { 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); option(returnSB); ideal star=quotient(laedi,f); option(noreturnSB); 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); option(returnSB); verg=quotient(laedi,g); option(noreturnSB); 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 ////////////////////////////////////////////////// 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); }