/////////////////////////////////////////////////////////////////////////////// version="$Id: primdec.lib,v 1.98.2.12 2002/11/05 13:01:12 Singular Exp $"; 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) 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 and Kemper (implementation by Gerhard Pfister). PROCEDURES: Ann(M); annihilator of R^n/M, R=basering, M in R^n primdecGTZ(I); complete primary decomposition via Gianni,Trager,Zacharias primdecSY(I...); complete primary decomposition via Shimoyama-Yokoyama minAssGTZ(I); the minimal associated primes via Gianni,Trager,Zacharias minAssChar(I...); the minimal associated primes using characteristic sets testPrimary(L,k); tests the result of the primary decomposition radical(I); computes the radical of I via Krick/Logar 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 "; LIB "general.lib"; LIB "elim.lib"; LIB "poly.lib"; LIB "random.lib"; LIB "inout.lib"; LIB "matrix.lib"; LIB "triang.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 EXAMPLE: example sat1; shows an example " { 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 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]; } 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); } /////////////////////////////////////////////////////////////////////////////// static 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; 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) { 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) { 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 " { 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))); } 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,1))==0) { return(1); } } return(0); } /////////////////////////////////////////////////////////////////////////////// proc primaryTest (ideal i, poly p) { 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=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])); if(char(basering)!=0) { f=1; if(e mod char(basering)==0) { if ( voice >=15 ) { "// 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/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,1) !=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,1) !=0) { return(ideal(0)); } } h=interred(t); t=h[1]; prm = prm,t; attrib(prm,"isSB",1); } return(prm); } /////////////////////////////////////////////////////////////////////////////// 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); } /////////////////////////////////////////////////////////////////////////////// static 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; 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) { if((char(basering)==0)||(npars(basering)>0)) { return(l); } 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)/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)/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]); l[2*i-1]=std(sa[1]); l[2*i]=primaryTest(sa[1],sa[1][1]); } if((size(sact)==1)&&(m==2)) { l[2*i]=l[2*i-1]; attrib(l[2*i],"isSB",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(noredSB); } 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],keep[i][1]); } } } l=l+keep1; 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 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/Trager/Zacharias EXAMPLE: example zero_decomp; shows an example " { 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); } 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); } //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((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>=6)&&(char(basering)<=181)) { primary=splitCharp(primary); } if((@wr==2)&&(npars(basering)>0)&&(voice>=6)&&(char(basering)>0)) { //the prime decomposition of Yokoyama in characteristic p list ke,ek; @k=0; while(@k=8){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>=8)) { poly randp; for(zz=1;zz0)&&(voice>=10)) { 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 ============================ @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) { helpprim=@lh[2]; lres0[1]=primary[2*@k-1]; ser1=psi(helpprim); lres0[2]=psi1(ser1); if(size(reduce(lres0[2],lres0[1],1))==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)/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)/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; // 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)/2; // for(@n=1;@n<=size(lres0)/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 #) { //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]) 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(ideal i) { //decomposes i into special ideals //computes the prime decomposition of the special ideals //and transforms it back to a decomposition of i 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(wo[1],2); option(redSB); for(k=1;k<=size(qr)/2;k++) { K=qr[2*k]; K=phi(K); K=groebner(K); re[size(re)+1]=zeroRad(K); } option(noredSB); } 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 " { 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 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 " { 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) "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; } /////////////////////////////////////////////////////////////////////////////// static 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 minAssPrimesold(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 EXAMPLE: example minAssPrimes; shows an example " { 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(ndim0) { if(deg(i[j])>1){break;} j--; } if(j==0){return(1);} if(deg(i[j])==vdim(i)){return(1);} return(0); } static 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 EXAMPLE: example minAssPrimes; shows an example " { def P=basering; list q=simplifyIdeal(i); list re=maxideal(1); int j,a,k; intvec op=option(get); map phi=P,q[2]; if(npars(P)==0){option(redSB);} i=std(q[1]); if((dim(i)==0)&&(npars(P)==0)) { int di=vdim(i); execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;"); ideal J=imap(P,i); attrib(J,"isSB",1); if(vdim(J)!=di) { J=fglm(P,i); } list pr=triangMH(J,2); list qr,re; for(k=1;k<=size(pr);k++) { if(primT(pr[k])) { re[size(re)+1]=pr[k]; } else { attrib(pr[k],"isSB",1); qr=decomp(pr[k],2); for(j=1;j<=size(qr)/2;j++) { re[size(re)+1]=qr[2*j]; } } } setring P; re=imap(gnir,re); option(set,op); return(phi(re)); } if((size(#)==0)||(dim(i)==0)) { re[1]=decomp(i,2); re=union(re); option(set,op); return(phi(re)); } q=facstd(i); if((size(q)==1)&&(dim(i)>1)) { execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;"); list p=facstd(fetch(P,i)); if(size(p)>1) { a=1; setring P; q=fetch(gnir,p); } else { setring P; } kill gnir; } option(set,op); for(j=1;j<=size(q);j++) { if(a==0){attrib(q[j],"isSB",1);} re[j]=decomp(q[j],2); } re=union(re); return(phi(re)); } 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 union(list li) { int i,j,k; def P=basering; execute("ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);"); list l=fetch(P,li); list @erg; for(k=1;k<=size(l);k++) { for(j=1;j<=size(l[k])/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(deg(@erg[size(@erg)][1])!=0) { @wos[size(@wos)+1]=@erg[size(@erg)]; } setring P; list @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 - a[1],...,a[s-1] the lower dimensional equidimensional loci. 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 " { if(ord_test(basering)!=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==1))&&(find(ordstr(basering),"l")==0) &&(find(ordstr(basering),"s")==0)) { execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(" +ordstr(basering)+");"); ideal i=imap(P,i); ideal j=i; if(a==1) { attrib(j,"isSB",1); } else { j=groebner(i); } } else { execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;"); 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,1))==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 " { if(ord_test(basering)!=1) { ERROR( "// Not implemented for this ordering, please change to 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==1))&&(find(ordstr(basering),"l")==0) &&(find(ordstr(basering),"s")==0)) { execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(" +ordstr(basering)+");"); ideal i=imap(P,i); ideal j=i; if(a==1) { attrib(j,"isSB",1); } else { j=groebner(i); } } else { execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;"); 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)); } string quotring=prepareQuotientring(nvars(basering)-indep[1][3]); execute(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(quring,h); kill quring; 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); } /////////////////////////////////////////////////////////////////////////////// proc algeDeco(ideal i, int w) { //reduces primery decomposition 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; if(w==2) { list pr=minAssPrimes(i); } else { list pr=decomp(i); } int j; gnir="ring RS="+string(char(R))+",("+varstr(RH) +"),(dp("+string(n)+"),lp);"; execute(gnir); list pr=imap(RH,pr); ideal K; for(j=1;j<=size(pr);j++) { K=groebner(pr[j]); K=K[2..size(K)]; pr[j]=K; } setring R; list pr=imap(RS,pr); list re; if(w==2) { re=pr; } else { re=convList(pr); } return(re); } /////////////////////////////////////////////////////////////////////////////// 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) ) 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 EXAMPLE: example decomp; shows an example " { intvec op; def @P = basering; list primary,indep,ltras; intvec @vh,isat,@w; int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi; ideal peek=i; ideal ser,tras; int isS=(attrib(i,"isSB")==1); 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; 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[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; } //---------------------------------------------------------------- //i is the zero-ideal //---------------------------------------------------------------- if(size(i)==0) { primary=i,i; return(primary); } //---------------------------------------------------------------- //pass to the lexicographical ordering and compute a standardbasis //---------------------------------------------------------------- int lp=islp(); execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);"); 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) { 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; primary[1]=i; primary[2]=i; 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); list pr=decomp(@j); setring gnir; list pr=imap(@deirf,pr); for(@k=1;@k<=size(pr);@k++) { @j=pr[@k]+fried; pr[@k]=@j; } setring @P; return(imap(gnir,pr)); } } //---------------------------------------------------------------- //j is the ring //---------------------------------------------------------------- if (dim(@j)==-1) { setring @P; primary=ideal(1),ideal(1); 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; primary=fetch(gnir,gprimary); return(primary); } //------------------------------------------------------------------ //the zero-dimensional case //------------------------------------------------------------------ if (dim(@j)==0) { op=option(get); option(redSB); list gprimary= zero_decomp(@j,ser,@wr); option(set,op); 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 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,@w); } } else { ideal jwork=groebner(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; 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+";"); 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)-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..] op=option(get); option(redSB); list uprimary= zero_decomp(@j,ser,@wr); 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(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) { 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; if(restindep[@m][1]==varstr(basering)) //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+";"); 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..] //------------------------------------------------------------------ quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]); //------------------------------------------------------------------ //we pass to the quotientring K(var(nnp+1),..,var(nva))[..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))[..rest..] op=option(get); option(redSB); list uprimary= zero_decomp(@j,ser,@wr); 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(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; list pr= decomp(i); pr; testPrimary( pr, i); } /////////////////////////////////////////////////////////////////////////////// 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 { 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); } /////////////////////////////////////////////////////////////////////////////// static 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 { def R=basering; int k; if(size(#)>0){k=#[1];} poly h=gcd(f,diff(f,var(i))); 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 { 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 i,k; list l; intvec op=option(get); option(redSB); 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)){return(interred(I+F));} //the separable case 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*y(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 radicalKL (ideal i,ideal ser,list #) { 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,quotring; def @P = basering; int jdim=dim(i); int homo=homog(i); ideal rad=ideal(1); 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); } //--------------------------------------------------------------------------- //j is the ring //--------------------------------------------------------------------------- if (jdim==-1) { return(ideal(1)); } //--------------------------------------------------------------------------- //the zero-dimensional case //--------------------------------------------------------------------------- if (jdim==0) { return(zeroRad(i)); } //------------------------------------------------------------------------- //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 //------------------------------------------------------------------------- indep=maxIndependSet(i); 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,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) //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); //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]); } op=option(get); option(redSB); @j=interred(@j); //to obtain a reduced standardbasis attrib(@j,"isSB",1); option(set,op); ideal zero_rad= zeroRad(@j); } 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(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]; } 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)1;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 ////////////////////////////////////////////////////// static 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 ////////////////////////////////////////////////////// static proc min_ass_prim_charsets0 (ideal PS) { intvec op; 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 { 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); } } 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 (ideal PS) { 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 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 { 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); } } 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(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); if(SI[1]==1) // primdecSY(ideal(1)) { return(list()); } 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 ////////////////////////////////////////////////////////////////////////// static 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 //////////////////////////////////////////////////////////////// static 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, //////////////////////////////////////////////////////////////// static 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 //////////////////////////////////////////////////////////////// static proc extraction (ideal SI, ideal SP) { 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) { 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) { 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)/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=1;i<=size(l);i++) { re[2*i-1]=l[i][1]; re[2*i]=l[i][2]; } return(re); } /////////////////////////////////////////////////////////////////////////////// // // The main procedures // /////////////////////////////////////////////////////////////////////////////// proc primdecGTZ(ideal i) "USAGE: primdecGTZ(i); i ideal 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 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!) EXAMPLE: example primdecGTZ; shows an example " { if(ord_test(basering)!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } if(minpoly!=0) { return(algeDeco(i,0)); ERROR( "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal" ); } return(convList(decomp(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; } /////////////////////////////////////////////////////////////////////////////// proc primdecSY(ideal i, list #) "USAGE: primdecSY(i); i ideal, c int 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 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, if c=2, minAssGTZ is used, if c=3, minAssGTZ and facstd are used. @end format EXAMPLE: example primdecSY; shows an example " { if(ord_test(basering)!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } i=simplify(i,2); if (i[1]==0) { list L=list(ideal(0),ideal(0)); return(list(L)); } if (size(#)==1) { return(prim_dec(i,#[1])); } else { return(prim_dec(i,1)); } } 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; } /////////////////////////////////////////////////////////////////////////////// proc minAssGTZ(ideal i,list #) "USAGE: minAssGTZ(i); i ideal minAssGTZ(i,1); i ideal does not use the factorizing Groebner 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 minAssGTZ; shows an example " { if(ord_test(basering)!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } if(minpoly!=0) { return(algeDeco(i,2)); ERROR( "// Not implemented for algebraic extensions. Simulate the ring extension by adding the minpoly to the ideal" ); } if(size(#)==0) { return(minAssPrimes(i,1)); } return(minAssPrimes(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 = minAssGTZ(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc minAssChar(ideal i, list #) "USAGE: minAssChar(i[,c]); i ideal, c int. RETURN: list, the minimal associated prime ideals of i. 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. @* Due to a bug in the factorization, the result may be not completely decomposed in small characteristic. EXAMPLE: example minAssChar; shows an example " { if(ord_test(basering)!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } if (size(#)==1) { return(min_ass_prim_charsets(i,#[1])); } else { return(min_ass_prim_charsets(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; list pr = minAssChar(i); pr; } /////////////////////////////////////////////////////////////////////////////// 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 and Kemper is used. Works also in positive characteristic (Kempers algorithm). EXAMPLE: example equiRadical; shows an example " { if(ord_test(basering)!=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); i ideal. RETURN: ideal, the radical of i. NOTE: A combination of the algorithms of Krick/Logar and Kemper is used. Works also in positive characteristic (Kempers algorithm). EXAMPLE: example radical; shows an example " { if(ord_test(basering)!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } def @P=basering; int j,il; if(size(#)>0){il=#[1];} if(size(i)==0){return(ideal(0));} ideal re=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]); return(interred(phi(i))); } option(redSB); list pr=i; if (!homog(i)) { pr=facstd(i); } option(set,op); int s=size(pr); for(j=1;j<=s;j++) { attrib(pr[s+1-j],"isSB",1); if((size(reduce(re,pr[s+1-j],1))!=0)&&((dim(pr[s+1-j])==di)||!il)) { re=intersect(re,radicalKL(pr[s+1-j],re,il)); } } return(interred(phi(re))); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; ideal pr= radical(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc 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 " { if(ord_test(basering)!=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) { list 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 " { if(ord_test(basering)!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } ideal j=groebner(i); int cod=nvars(basering)-dim(j); int e; ideal ann; if(homog(i)==1) { list re=sres(j,0); //the resolution re=minres(re); //minimized resolution } else { list re=mres(i,0); } ann=AnnExt_R(cod,re); 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 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 " { int i; pr=reconvList(pr); 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 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); } /////////////////////////////////////////////////////////////////////////////// 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 " { if(ord_test(basering)!=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-T*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 = T-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; } //////////////////////////////////////////////////////////////////////////// /* //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); */