/////////////////////////////////////////////////////////////////////////////// version="$Id$"; category="Commutative Algebra"; info=" LIBRARY: primdec.lib Primary Decomposition and Radical of Ideals AUTHORS: Gerhard Pfister, pfister@mathematik.uni-kl.de (GTZ)@* Wolfram Decker, decker@math.uni-sb.de (SY)@* Hans Schoenemann, hannes@mathematik.uni-kl.de (SY)@* Santiago Laplagne, slaplagn@dm.uba.ar (GTZ) OVERVIEW: Algorithms for primary decomposition based on the ideas of Gianni, Trager and Zacharias (implementation by Gerhard Pfister), respectively based on the ideas of Shimoyama and Yokoyama (implementation by Wolfram Decker and Hans Schoenemann).@* The procedures are implemented to be used in characteristic 0.@* They also work in positive characteristic >> 0.@* In small characteristic and for algebraic extensions, primdecGTZ may not terminate.@* Algorithms for the computation of the radical based on the ideas of Krick, Logar, Laplagne and Kemper (implementation by Gerhard Pfister and Santiago Laplagne). They work in any characteristic.@* Baserings must have a global ordering and no quotient ideal. 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 (with modifications by Laplagne) 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 (with modifications by Laplagne) and Kemper radicalEHV(I); computes the radical of I via Eisenbud,Huneke,Vasconcelos equiRadical(I); the radical of the equidimensional part of the ideal I prepareAss(I); list of radicals of the equidimensional components of I equidim(I); weak equidimensional decomposition of I equidimMax(I); equidimensional locus of I equidimMaxEHV(I); equidimensional locus of I via Eisenbud,Huneke,Vasconcelos zerodec(I); zerodimensional decomposition via Monico absPrimdecGTZ(I); the absolute prime components of I "; LIB "general.lib"; LIB "elim.lib"; LIB "poly.lib"; LIB "random.lib"; LIB "inout.lib"; LIB "matrix.lib"; LIB "triang.lib"; LIB "absfact.lib"; LIB "ring.lib"; /////////////////////////////////////////////////////////////////////////////// // // Gianni/Trager/Zacharias // /////////////////////////////////////////////////////////////////////////////// static proc sat1 (ideal id, poly p) "USAGE: sat1(id,j); id ideal, j polynomial RETURN: saturation of id with respect to j (= union_(k=1...) of id:j^k) NOTE: result is a std basis in the basering " { 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 " { 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,6); poly @f=1; if(deg(fac[1])>0) { ideal iold; for(@i=1;@i<=size(fac);@i++) { @f=@f*fac[@i]; } intvec op = option(get); option(returnSB); while(specialIdealsEqual(iold,inew)==0 ) { iold=inew; if(deg(iold[size(iold)])!=1) { inew=quotient(iold,@f); } else { inew=iold; } @k++; } option(set,op); @k--; } if(ordstr(@P)[1,2]!="dp") { setring @P; ideal inew=std(imap(@Phelp,inew)); poly @f=imap(@Phelp,@f); } list L =inew,@f^@k; return (L); } /////////////////////////////////////////////////////////////////////////////// proc minSat(ideal inew, ideal h) { int i,k; poly f=1; ideal iold,fac; list quotM,l; for(i=1;i<=ncols(h);i++) { if(deg(h[i])>0) { fac=fac+factorize(h[i],1); } } fac=simplify(fac,6); if(size(fac)==0) { l=inew,1; return(l); } fac=sort(fac)[1]; for(i=1;i<=size(fac);i++) { f=f*fac[i]; } quotM[1]=inew; quotM[2]=fac; quotM[3]=f; f=1; intvec op = option(get); option(returnSB); while(specialIdealsEqual(iold,quotM[1])==0) { if(k>0) { f=f*quotM[3]; } iold=quotM[1]; quotM=quotMin(quotM); k++; } option(set,op); l=quotM[1],f; return(l); } static proc quotMin(list tsil) { 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; @l=factorize(p); for(@j=1;@j<=size(@l[1]);@j++) { if(leadcoef(@l[1][@j])==@l[1][@j]) { @n++; } } if(@n>0) { if(@n==size(@l[1])) { @l[1]=ideal(1); @v=1; @l[2]=@v; } else { @k=0; int pleh; for(@j=1;@j<=size(@l[1]);@j++) { if(leadcoef(@l[1][@j])!=@l[1][@j]) { @k++; @i=@i+ideal(@l[1][@j]); if(size(@i)==pleh) { "//factorization error"; @l; @k--; @v[@k]=@v[@k]+@l[2][@j]; } else { pleh++; @v[@k]=@l[2][@j]; } } } @l[1]=@i; @l[2]=@v; } } // } return(@l); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = (x+y)^2*(y-z)^3; list l = factor(p); l; ring r1 =(0,b,d,f,g),(a,c,e),lp; poly p =(1*d)*e^2+(1*d*f^2*g); list l = factor(p); l; ring r2 =(0,b,f,g),(a,c,e,d),lp; poly p =(1*d)*e^2+(1*d*f^2*g); list l = factor(p); l; } /////////////////////////////////////////////////////////////////////////////// proc idealsEqual( ideal k, ideal j) { return(stdIdealsEqual(std(k),std(j))); } static proc specialIdealsEqual( ideal k1, ideal k2) { int j; if(size(k1)==size(k2)) { for(j=1;j<=size(k1);j++) { if(leadexp(k1[j])!=leadexp(k2[j])) { return(0); } } return(1); } return(0); } static proc stdIdealsEqual( ideal k1, ideal k2) { int j; if(size(k1)==size(k2)) { for(j=1;j<=size(k1);j++) { if(leadexp(k1[j])!=leadexp(k2[j])) { return(0); } } attrib(k2,"isSB",1); if(size(reduce(k1,k2,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--; m++; //search for i[m] which has a power of var(n) as leading term if (n==1) { m=size(i); } else { while (lead(i[m])/var(n-1)==0) { m++; } m--; } //check whether i[m] =(c*var(n)+h)^e modulo prm for some //h in K[var(n+1),...,var(nvars(basering))], c in K //if not (0) is returned, else var(n)+h is added to prm e=deg(lead(i[m])); if(char(basering)!=0) { f=1; if(e mod char(basering)==0) { if ( voice >=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 div char(basering); } } t=leadcoef(i[m])*e*var(n)^f+(i[m]-lead(i[m]))/var(n)^((e-1)*f); i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m]; if (reduce(i[m]-t^e,prm,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 div 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 div 2) { j=1; } */ while(j1)) { keepprime[i]=interred(keepprime[i]+ideal(act[1][1])); if(homog(keepprime[i])==1) { l[2*i]=maxideal(1); break; } } if(gcdTest(act[1])==1) { for(k=2;k<=r;k++) { keepprime[size(l) div 2+k-1]=interred(keepprime[i]+ideal(act[1][k])); } keepprime[i]=interred(keepprime[i]+ideal(act[1][1])); for(k=1;k<=r;k++) { if(@wr==0) { keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]); } else { keepresult[k]=std(l[2*i-1],act[1][k]); } } l[2*i-1]=keepresult[1]; if(vdim(keepresult[1])==deg(act[1][1])) { l[2*i]=keepresult[1]; } if((homog(keepresult[1])==1)||(homog(keepprime[i])==1)) { l[2*i]=maxideal(1); } s=size(l)-2; for(k=2;k<=r;k++) { l[s+2*k-1]=keepresult[k]; keepprime[s div 2+k]=interred(keepresult[k]+ideal(act[1][k])); if(vdim(keepresult[k])==deg(act[1][k])) { l[s+2*k]=keepresult[k]; } else { l[s+2*k]=ideal(0); } if((homog(keepresult[k])==1)||(homog(keepprime[s div 2+k])==1)) { l[s+2*k]=maxideal(1); } } i--; break; } if(r>=2) { s=size(l); @f=act[1][1]; act=sat1(l[2*i-1],act[1][1]); if(deg(act[1][1])>0) { l[s+1]=std(l[2*i-1],act[2]); if(homog(l[s+1])==1) { l[s+2]=maxideal(1); } else { l[s+2]=ideal(0); } keepprime[s div 2+1]=interred(keepprime[i]+ideal(@f)); if(homog(keepprime[s div 2+1])==1) { l[s+2]=maxideal(1); } keepprime[i]=act[1]; l[2*i-1]=act[1]; attrib(l[2*i-1],"isSB",1); if(homog(l[2*i-1])==1) { l[2*i]=maxideal(1); } i--; break; } } } } } if(sl==size(l)) { return(l); } for(i=1;i<=size(l) div 2;i++) { attrib(l[2*i-1],"isSB",1); if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],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) div 2;i++) { if(size(l[2*i])==0) { if(deg(l[2*i-1][1])==vdim(l[2*i-1])) { l[2*i]=l[2*i-1]; } } } for(i=1;i<=size(l) div 2;i++) { if(size(l[2*i])==0) { s=factorize(l[2*i-1][1],1); //vermeiden!!! t=l[2*i-1]; m=size(t); ni=s[1]; if(deg(ni)>1) { va=varstr(P); j=size(va); while(va[j]!=","){j--;} va=va[1..j-1]; gnir="ring RL=("+string(char(P))+","+string(var(n))+"),("+va+"),lp;"; execute(gnir); minpoly=leadcoef(imap(P,ni)); ideal act; ideal t=imap(P,t); for(k=2;k<=m;k++) { act=factorize(t[k],1); if(size(act)>1){break;} } setring P; sact=imap(RL,act); if(size(sact)>1) { sa=sat1(l[2*i-1],sact[1]); keep[size(keep)+1]=std(l[2*i-1],sa[2]); 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 computation of associated primes) RETURN: list = list of primary ideals and their radicals (at even positions in the list) if the input is zero-dimensional and a standardbases with respect to lex-ordering If ser!=(0) and ser is contained in j or if j is not zero-dimen- sional then ideal(1),ideal(1) is returned NOTE: Algorithm of Gianni/Trager/Zacharias EXAMPLE: example zero_decomp; shows an example " { def @P = basering; int uytrewq; int nva = nvars(basering); int @k,@s,@n,@k1,zz; list primary,lres0,lres1,act,@lh,@wh; map phi,psi,phi1,psi1; ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1; intvec @vh,@hilb; string @ri; poly @f; if (dim(j)>0) { primary[1]=ideal(1); primary[2]=ideal(1); return(primary); } intvec save=option(get); option(redSB); j=interred(j); attrib(j,"isSB",1); if(vdim(j)==deg(j[1])) { act=factor(j[1]); for(@k=1;@k<=size(act[1]);@k++) { @qh=j; if(@wr==0) { @qh[1]=act[1][@k]^act[2][@k]; } else { @qh[1]=act[1][@k]; } primary[2*@k-1]=interred(@qh); @qh=j; @qh[1]=act[1][@k]; primary[2*@k]=interred(@qh); attrib( primary[2*@k-1],"isSB",1); if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0)) { primary[2*@k-1]=ideal(1); primary[2*@k]=ideal(1); } } option(set,save); return(primary); } option(set,save); if(homog(j)==1) { primary[1]=j; if((size(ser)>0)&&(size(reduce(ser,j,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(@wr!=0) { primary[1]=std(j,act[1][1]); } if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1]))) { primary[2]=primary[1]; } else { primary[2]=primaryTest(primary[1],act[1][1]); } } if(size(#)==0) { primary=splitPrimary(primary,ser,@wr,act); } if((voice>=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 ============================ if (npars(basering)>1) { @qh=groebner(@qh,"par2var"); } else { @qh=groebner(@qh); } //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // ideal @qh=homog(imap(@P,@qht),@t); // // ideal @qh1=std(@qh); // @hilb=hilb(@qh1,1); // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // ideal @qh=homog(imap(@P,@qh),@t); // kill @Phelp; // @qh=std(@qh,@hilb); // @qh=subst(@qh,@t,1); // setring @P; // @qh=imap(@Phelp1,@qh); // kill @Phelp1; // @qh=clearSB(@qh); // attrib(@qh,"isSB",1); //============================================================= ser1=phi1(ser); @lh=zero_decomp (@qh,phi(ser1),@wr); kill lres0; list lres0; if(size(@lh)==2) { 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) div 2;@n++) { if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1) { lres0[2*@n-1]=groebner(lres0[2*@n-1]); lres0[2*@n]=lres0[2*@n-1]; attrib(lres0[2*@n],"isSB",1); } else { lres0[2*@n-1]=groebner(lres0[2*@n-1]); lres0[2*@n]=groebner(lres0[2*@n]); } } } primary=primary+lres0; //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // list @lvec; // list @lr=imap(@P,lres0); // ideal @lr1; // // if(size(@lr)==2) // { // @lr[2]=homog(@lr[2],@t); // @lr1=std(@lr[2]); // @lvec[2]=hilb(@lr1,1); // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lvec[2*@n]=@lvec[2*@n-1]; // } // else // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lr[2*@n]=homog(@lr[2*@n],@t); // @lr1=std(@lr[2*@n]); // @lvec[2*@n]=hilb(@lr1,1); // // } // } // } // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // list @lr=imap(@Phelp,@lr); // // kill @Phelp; // if(size(@lr)==2) // { // @lr[2]=std(@lr[2],@lvec[2]); // @lr[2]=subst(@lr[2],@t,1); // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=@lr[2*@n-1]; // attrib(@lr[2*@n],"isSB",1); // } // else // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]); // @lr[2*@n]=subst(@lr[2*@n],@t,1); // } // } // } // kill @lvec; // setring @P; // lres0=imap(@Phelp1,@lr); // kill @Phelp1; // for(@n=1;@n<=size(lres0);@n++) // { // lres0[@n]=clearSB(lres0[@n]); // attrib(lres0[@n],"isSB",1); // } // // primary[2*@k-1]=lres0[1]; // primary[2*@k]=lres0[2]; // @s=size(primary) div 2; // for(@n=1;@n<=size(lres0) div 2-1;@n++) // { // primary[2*@s+2*@n-1]=lres0[2*@n+1]; // primary[2*@s+2*@n]=lres0[2*@n+2]; // } // @k--; //============================================================= } } return(primary); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; i=std(i); list pr= 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) div 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 clearSBNeu (ideal i,list #) "USAGE: clearSB(i); i ideal which is SB ordered by monomial ordering RETURN: ideal = minimal SB NOTE: EXAMPLE: example clearSB; shows an example " { int k,j; intvec m,n,v,w; int c=size(i); w=leadexp(0); v[size(i)]=0; j=0; while(j=0) { m=leadexp(i[j]); for(k=j+1;k<=c;k++) { n=leadexp(i[k]); if(n!=w) { if(((m==n)&&(#[j]>#[k]))||((teilt(n,m))&&(n!=m))) { i[j]=0; v[j]=1; break; } if(((m==n)&&(#[j]<=#[k]))||((teilt(m,n))&&(n!=m))) { i[k]=0; v[k]=1; } } } } } return(v); } static proc teilt(intvec a, intvec b) { int i; for(i=1;i<=size(a);i++) { if(a[i]>b[i]){return(0);} } return(1); } /////////////////////////////////////////////////////////////////////////////// static proc independSet (ideal j) "USAGE: independentSet(i); i ideal RETURN: list = new varstring with the independent set at the end, ordstring with the corresponding block ordering, the integer where the independent set starts in the varstring NOTE: EXAMPLE: example independentSet; shows an example " { int n,k,di; list resu,hilf; string var1,var2; list v=indepSet(j,1); for(n=1;n<=size(v);n++) { di=0; var1=""; var2=""; for(k=1;k<=size(v[n]);k++) { if(v[n][k]!=0) { di++; var2=var2+"var("+string(k)+"),"; } else { var1=var1+"var("+string(k)+"),"; } } if(di>0) { var1=var1+var2; var1=var1[1..size(var1)-1]; hilf[1]=var1; hilf[2]="lp"; //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")"; hilf[3]=di; resu[n]=hilf; } else { resu[n]=varstr(basering),ordstr(basering),0; } } return(resu); } example { "EXAMPLE:"; echo = 2; ring s1=(0,x,y),(a,b,c,d,e,f,g),lp; ideal i=ea-fbg,fa+be,ec-fdg,fc+de; i=std(i); list l=independSet(i); l; i=i,g; l=independSet(i); l; ring s=0,(x,y,z),lp; ideal i=z,yx; list l=independSet(i); l; } /////////////////////////////////////////////////////////////////////////////// static proc maxIndependSet (ideal j) "USAGE: maxIndependentSet(i); i ideal RETURN: list = new varstring with the maximal independent set at the end, ordstring with the corresponding block ordering, the integer where the independent set starts in the varstring NOTE: EXAMPLE: example maxIndependentSet; shows an example " { 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) div 2;i++) { if(deg(l[2*i-1][1])>0) { j++; lh[j]=l[2*i-1]; j++; lh[j]=l[2*i]; } } return(lh); } /////////////////////////////////////////////////////////////////////////////// proc 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; if(size(i)==0){return(list(ideal(0)));} 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 Optional parameters in list #: (can be entered in any order) 0, "facstd" -> uses facstd to first decompose the ideal 1, "noFacstd" -> does not use facstd (default) "SL" -> the new algorithm is used (default) "GTZ" -> the old algorithm is used RETURN: list = the minimal associated prime ideals of i EXAMPLE: example minAssPrimes; shows an example " { if(size(i) == 0){return(list(ideal(0)));} string algorithm; // Algorithm to be used string facstdOption; // To uses proc facstd int j; // Counter def P0 = basering; list Pl=ringlist(P0); intvec dp_w; for(j=nvars(P0);j>0;j--) {dp_w[j]=1;} Pl[3]=list(list("dp",dp_w),list("C",0)); def P=ring(Pl); setring P; ideal i=imap(P0,i); // Set input parameters algorithm = "SL"; // Default: SL algorithm facstdOption = "Facstd"; // Default: facstd is not used if(size(#) > 0) { int valid; for(j = 1; j <= size(#); j++) { valid = 0; if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) { if (#[j] == 0) {facstdOption = "noFacstd"; valid = 1;} // If #[j] == 0, facstd is not used. if (#[j] == 1) {facstdOption = "facstd"; valid = 1;} // If #[j] == 1, facstd is used. } if(typeof(#[j]) == "string") { if(#[j] == "GTZ" || #[j] == "SL") { algorithm = #[j]; valid = 1; } if(#[j] == "noFacstd" || #[j] == "facstd") { facstdOption = #[j]; valid = 1; } } if(valid == 0) { dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); } } } list q = simplifyIdeal(i); list re = maxideal(1); int a, k; intvec op = option(get); map phi = P,q[2]; list result; if(npars(P) == 0){option(redSB);} if(attrib(i,"isSB")!=1) { i=groebner(q[1]); } else { for(j=1;j<=nvars(basering);j++) { if(q[2][j]!=var(j)){k=1;break;} } if(k) { i=groebner(q[1]); } } if(dim(i) == -1){setring P0;return(ideal(1));} if((dim(i) == 0) && (npars(P) == 0)) { int di = vdim(i); execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;"); ideal J = std(imap(P,i)); attrib(J, "isSB", 1); if(vdim(J) != di) { J = fglm(P, i); } // list pr = triangMH(J,2); HIER KOENNEN verschiedene Mengen zu gleichen // asoziierten Primidealen fuehren // Aenderung list pr = triangMH(J,2); list qr, re; for(k = 1; k <= size(pr); k++) { if(primT(pr[k])&&(0)) { re[size(re) + 1] = pr[k]; } else { attrib(pr[k], "isSB", 1); // Lines changed if (algorithm == "GTZ") { qr = decomp(pr[k], 2); } else { qr = minAssSL(pr[k]); } for(j = 1; j <= size(qr) div 2; j++) { re[size(re) + 1] = std(qr[2 * j]); } } } setring P; re = imap(gnir, re); re=phi(re); option(set, op); setring(P0); list re=imap(P,re); return(re); } // Lines changed if ((facstdOption == "noFacstd") || (dim(i) == 0)) { if (algorithm == "GTZ") { re[1] = decomp(i, 2); } else { re[1] = minAssSL(i); } re = union(re); option(set, op); re=phi(re); setring(P0); list re=imap(P,re); return(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); // Debug dbprint(printlevel - voice, "Components returned by facstd", size(q), q); for(j = 1; j <= size(q); j++) { if(a == 0){attrib(q[j], "isSB", 1);} // Debug dbprint(printlevel - voice, "We compute the decomp of component", j); // Lines changed if (algorithm == "GTZ") { re[j] = decomp(q[j], 2); } else { re[j] = minAssSL(q[j]); } // Debug dbprint(printlevel - voice, "Number of components obtained for this component:", size(re[j]) div 2); dbprint(printlevel - voice, "re[j]:", re[j]); } re = union(re); re=phi(re); setring(P0); list re=imap(P,re); return(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]) div 2;j++) { if(deg(l[k][2*j][1])!=0) { i++; @erg[i]=l[k][2*j]; } } } list @wos; i=0; ideal i1,i2; while(isize(@erg)) { @wos[size(@wos)+1]=@erg[i]; } } } if(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(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } intvec op ; def P = basering; list eq; intvec w; int n,m; int g=size(i); int a=attrib(i,"isSB"); int homo=homog(i); if(size(#)!=0) { m=1; } if(((homo==1)||(a==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(attrib(basering,"global")!=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); //---Anfang Provisorium if((size(i)==2) && (w==2)) { option(redSB); ideal J=std(i); option(noredSB); if((size(J)==2)&&(deg(J[1])==1)) { ideal keep; poly f; int j; for(j=1;j<=nvars(basering);j++) { f=J[2]; while((f/var(j))*var(j)-f==0) { f=f/var(j); keep=keep,var(j); } J[2]=f; } ideal K=factorize(J[2],1); if(deg(K[1])==0){K=0;} K=K+std(keep); ideal L; list resu; for(j=1;j<=size(K);j++) { L=J[1],K[j]; resu[j]=L; } return(resu); } } //---Ende Provisorium 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); ideal I=subst(i,var(nvars(basering)),0); int j; for(j=1;j<=ncols(i);j++) { if(i[j]!=I[j]){break;} } if((j>ncols(i))&&(deg(p)==1)) { setring R; kill RH; kill gnir; string gnir="ring RH="+string(char(R))+",("+varstr(R)+"),dp;"; execute(gnir); ideal i=imap(R,i); ideal J; } else { i=i,p; } list pr; if(w==0) { pr=decomp(i); } if(w==1) { pr=prim_dec(i,1); pr=reconvList(pr); } if(w==2) { pr=minAssPrimes(i); } if(n0) { if((#[1]==1)||(#[1]==2)||(#[1]==3)) { @wr=#[1]; if(@wr==3){abspri=1;@wr=0;} if(size(#)>1) { seri=1; peek=#[2]; ser=#[3]; } } else { seri=1; peek=#[1]; ser=#[2]; } } if(abspri) { list absprimary,abskeep,absprimarytmp,abskeeptmp; } homo=homog(i); if(homo==1) { if(attrib(i,"isSB")!=1) { //ltras=mstd(i); tras=groebner(i); ltras=tras,tras; 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) && (!abspri)) { 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=ideal(0),ideal(0); if (abspri) { return(prepare_absprimdec(primary));} 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)&&(!abspri)) { ideal fried; @n=size(@j); for(@k=1;@k<=@n;@k++) { if(deg(lead(@j[@k]))==1) { fried[size(fried)+1]=@j[@k]; @j[@k]=0; } } if(size(fried)==nvars(basering)) { setring @P; primary[1]=i; primary[2]=i; if (abspri) { return(prepare_absprimdec(primary));} return(primary); } if(size(fried)>0) { string newva; string newma; poly f; for(@k=1;@k<=nvars(basering);@k++) { @n1=0; for(@n=1;@n<=size(fried);@n++) { if(leadmonom(fried[@n])==var(@k)) { @n1=1; break; } } if(@n1==0) { newva=newva+string(var(@k))+","; newma=newma+string(var(@k))+","; } else { newma=newma+string(0)+","; fried[@n]=fried[@n]/leadcoef(fried[@n]); f=fried[@n]-lead(fried[@n]); @j=subst(@j,var(@k),-f); } } newva[size(newva)]=")"; newma[size(newma)]=";"; execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;"); execute("map @kappa=gnir,"+newma); ideal @j= @kappa(@j); @j=std(@j); list pr=decomp(@j); setring gnir; list pr=imap(@deirf,pr); for(@k=1;@k<=size(pr);@k++) { @j=pr[@k]+fried; pr[@k]=@j; } setring @P; primary=imap(gnir,pr); if (abspri) { return(prepare_absprimdec(primary));} return(primary); } } //---------------------------------------------------------------- //j is the ring //---------------------------------------------------------------- if (dim(@j)==-1) { setring @P; primary=ideal(1),ideal(1); if (abspri) { return(prepare_absprimdec(primary));} return(primary); } //---------------------------------------------------------------- // the case of one variable //---------------------------------------------------------------- if(nvars(basering)==1) { list fac=factor(@j[1]); list gprimary; for(@k=1;@k<=size(fac[1]);@k++) { if(@wr==0) { gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } else { gprimary[2*@k-1]=ideal(fac[1][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } } setring @P; primary=fetch(gnir,gprimary); //HIER if (abspri) { return(prepare_absprimdec(primary));} return(primary); } //------------------------------------------------------------------ //the zero-dimensional case //------------------------------------------------------------------ if (dim(@j)==0) { op=option(get); option(redSB); list gprimary= zero_decomp(@j,ser,@wr); setring @P; primary=fetch(gnir,gprimary); if(size(ser)>0) { primary=cleanPrimary(primary); } //HIER if(abspri) { setring gnir; list primary=imap(@P,primary); list resu,tempo; string absotto; map sigma,invsigma; ideal II,jmap; nn=nvars(basering); for(ab=1;ab<=size(primary) div 2;ab++) { II=primary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])==vdim(II)) { absotto= absFactorize(primary[2*ab][1],77); tempo= primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); } else { invsigma=basering,maxideal(1); jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); absotto = absFactorize(II[1],77); II=var(nn); tempo= primary[2*ab-1],primary[2*ab],absotto,string(invsigma(II)); } resu[ab]=tempo; } primary=resu; setring @P; primary=imap(gnir,primary); } option(set,op); 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; ideal @Ptest=1; di=dim(jwork); keepdi=di; setring gnir; for(@m=1;@m<=size(indep);@m++) { isat=0; @n2=0; if((indep[@m][1]==varstr(basering))&&(@m==1)) //this is the good case, nothing to do, just to have the same notations //change the ring { 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)); if(@m==1) { @j=fetch(@P,i); } execute("ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),(" +indep[@m][2]+");"); execute("map phi=gnir,"+@va+";"); op=option(get); option(redSB); 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..] //--------------------------------------------------------------------- ideal @jj=lead(@j); //!! vorn vereinbaren execute(quotring); ideal @jj=imap(gnir1,@jj); @vv=clearSBNeu(@jj,fett); //!! vorn vereinbaren setring gnir1; @k=size(@j); for (lauf=1;lauf<=@k;lauf++) { if(@vv[lauf]==1) { @j[lauf]=0; } } @j=simplify(@j,2); setring quring; // @j considered in the quotientring ideal @j=imap(gnir1,@j); ideal ser=imap(gnir1,ser); kill gnir1; //j is a standardbasis in the quotientring but usually not minimal //here it becomes minimal attrib(@j,"isSB",1); //we need later ggt(h[1],...)=gh for saturation ideal @h; if(deg(@j[1])>0) { for(@n=1;@n<=size(@j);@n++) { @h[@n]=leadcoef(@j[@n]); } //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] op=option(get); option(redSB); list uprimary= zero_decomp(@j,ser,@wr); //HIER if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); } else { list uprimary; uprimary[1]=ideal(1); uprimary[2]=ideal(1); } //we need the intersection of the ideals in the list quprimary with the //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal //but fi polynomials, then the intersection of q with the polynomialring //is the saturation of the ideal generated by f1,...,fr with respect to //h which is the lcm of the leading coefficients of the fi considered in //in the quotientring: this is coded in saturn list saturn; ideal hpl; for(@n=1;@n<=size(uprimary);@n++) { uprimary[@n]=interred(uprimary[@n]); // temporary fix hpl=0; for(@n1=1;@n1<=size(uprimary[@n]);@n1++) { hpl=hpl,leadcoef(uprimary[@n][@n1]); } saturn[@n]=hpl; } //-------------------------------------------------------------------- //we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //back to the polynomialring //--------------------------------------------------------------------- setring gnir; collectprimary=imap(quring,uprimary); lsau=imap(quring,saturn); @h=imap(quring,@h); kill quring; @n2=size(quprimary); @n3=@n2; for(@n1=1;@n1<=size(collectprimary) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } if(size(@h)>0) { //--------------------------------------------------------------- //we change to @Phelp to have the ordering dp for saturation //--------------------------------------------------------------- setring @Phelp; @h=imap(gnir,@h); if(@wr!=1) { if(defined(@LL)){kill @LL;} list @LL=minSat(jwork,@h); @Ptest=intersect(@Ptest,@LL[1]); @q=@LL[2]; } else { fac=ideal(0); for(lauf=1;lauf<=ncols(@h);lauf++) { if(deg(@h[lauf])>0) { fac=fac+factorize(@h[lauf],1); } } fac=simplify(fac,6); @q=1; for(lauf=1;lauf<=size(fac);lauf++) { @q=@q*fac[lauf]; } } jwork=std(jwork,@q); keepdi=dim(jwork); if(keepdi0) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); if(@wr==0) { //HIER STATT DURCHSCHNITT SATURIEREN! ideal htest=@Ptest; } else { ideal htest=hquprimary[2]; for (@n1=2;@n1<=size(hquprimary) div 2;@n1++) { htest=intersect(htest,hquprimary[2*@n1]); } } if(size(ser)>0) { ser=intersect(htest,ser); } else { ser=htest; } setring gnir; ser=imap(@Phelp,ser); } if(size(reduce(ser,peek,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); //HIER if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); //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) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } if(@n2>=@n3+2) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(@wr==0) { ser=intersect(ser,hquprimary[2*@n-1]); } else { ser=intersect(ser,hquprimary[2*@n]); } } setring gnir; ser=imap(@Phelp,ser); } // } } //HIER if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { if (deg(quprimary[2*ab][1])!=0) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } } quprimary=resu; @wr=3; } if(size(reduce(ser,peek,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]; } } } } else { if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } quprimary=resu; } } //--------------------------------------------------------------------------- //back to the ring we started with //the final result: primary //--------------------------------------------------------------------------- setring @P; primary=imap(gnir,quprimary); if(!abspri) { primary=cleanPrimary(primary); } if (abspri && (typeof(primary[1][1])=="poly")) { return(prepare_absprimdec(primary));} 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); } /////////////////////////////////////////////////////////////////////////////// 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))); if((reduce(f,std(h))!=0)||(reduce(diff(f,var(i)),std(h))!=0)) { ERROR("FEHLER IN GCD"); } poly g1=lift(h,f)[1][1]; // f/h poly h1; while(h!=h1) { h1=h; h=gcd(h,diff(h,var(i))); } if(deg(h1)==0){return(list(g1,k));} //in characteristic 0 we return here k++; ideal ma=maxideal(1); ma[i]=var(i)^char(R); map phi=R,ma; ideal hh=h; //this is technical because preimage works only for ideals poly u=preimage(R,phi,hh)[1]; //h=u(x(i)^p) list g2=sep(u,i,k); //we consider u(t(1)^(p^-1),...,t(m)^(p^-1)) g1=powerCoeffs(g1,g2[2]-k+1); //to have g1 over the same field as g2[1] list g3=sep(g1*g2[1],i,g2[2]); return(g3); } example { "EXAMPLE:"; echo = 2; ring R=(5,t,s),(x,y,z),dp; poly f=(x^25-t*x^5+t)*(x^3+s); sep(f,1); } /////////////////////////////////////////////////////////////////////////////// proc zeroRad(ideal I,list #) "USAGE: zeroRad(I) , I a zero-dimensional ideal RETURN: the radical of I NOTE: Algorithm of Kemper EXAMPLE: example zeroRad; shows an example { if(homog(I)==1){return(maxideal(1));} //I needs to be a reduced standard basis def R=basering; int m=npars(R); int n=nvars(R); int p=char(R); int d=vdim(I); int i,k; list l; if(((p==0)||(p>d))&&(d==deg(I[1]))) { intvec e=leadexp(I[1]); for(i=1;i<=nvars(basering);i++) { if(e[i]!=0) break; } I[1]=sep(I[1],i)[1]; return(interred(I)); } intvec op=option(get); option(redSB); ideal F=finduni(I);//F[i] generates I intersected with K[var(i)] option(set,op); if(size(#)>0){I=#[1];} for(i=1;i<=n;i++) { l[i]=sep(F[i],i); F[i]=l[i][1]; if(l[i][2]>k){k=l[i][2];} //computation of the maximal k } if((k==0)||(m==0)) //the separable case { intvec save=option(get);option(redSB); I=interred(I+F);option(set,save);return(I); } 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); } /////////////////////////////////////////////////////////////////////////////// proc radicalEHV(ideal i) "USAGE: radicalEHV(i); i ideal. RETURN: ideal, the radical of i. NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos, which reduces the computation to the complete intersection case, by taking, in the general case, a generic linear combination of the input. Works only in characteristic 0 or p large. EXAMPLE: example radicalEHV; shows an example " { if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } if((char(basering)<100)&&(char(basering)!=0)) { "WARNING: The characteristic is too small, the result may be wrong"; } ideal J,I,I0,radI0,L,radI1,I2,radI2; int l,n; intvec op=option(get); matrix M; option(redSB); list m=mstd(i); I=m[2]; option(set,op); int cod=nvars(basering)-dim(m[1]); //-------------------complete intersection case:---------------------- if(cod==size(m[2])) { J=minor(jacob(I),cod); return(quotient(I,J)); } //-----first codim elements of I are a complete intersection:--------- for(l=1;l<=cod;l++) { I0[l]=I[l]; } n=dim(std(I0))+cod-nvars(basering); //-----last codim elements of I are a complete intersection:---------- if(n!=0) { for(l=1;l<=cod;l++) { I0[l]=I[size(I)-l+1]; } n=dim(std(I0))+cod-nvars(basering); } //-----taking a generic linear combination of the input:-------------- if(n!=0) { M=transpose(sparsetriag(size(m[2]),cod,95,1)); I0=ideal(M*transpose(I)); n=dim(std(I0))+cod-nvars(basering); } //-----taking a more generic linear combination of the input:--------- if(n!=0) { M=transpose(sparsetriag(size(m[2]),cod,0,100)); I0=ideal(M*transpose(I)); n=dim(std(I0))+cod-nvars(basering); } if(n==0) { J=minor(jacob(I0),cod); radI0=quotient(I0,J); L=quotient(radI0,I); radI1=quotient(radI0,L); if(size(reduce(radI1,m[1],1))==0) { return(I); } I2=sat(I,radI1)[1]; if(deg(I2[1])<=0) { return(radI1); } return(intersect(radI1,radicalEHV(I2))); } //---------------------general case------------------------------------- return(radical(I)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; ideal pr= radicalEHV(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc Ann(module M) "USAGE: Ann(M); M module RETURN: ideal, the annihilator of coker(M) NOTE: The output is the ideal of all elements a of the basering R such that a * R^m is contained in M (m=number of rows of M). EXAMPLE: example Ann; shows an example " { M=prune(M); //to obtain a small embedding ideal ann=quotient1(M,freemodule(nrows(M))); return(ann); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; module M = x2-y2,z3; Ann(M); M = [1,x2],[y,x]; Ann(M); qring Q=std(xy-1); module M=imap(r,M); Ann(M); } /////////////////////////////////////////////////////////////////////////////// //computes the equidimensional part of the ideal i of codimension e static proc int_ass_primary_e(ideal i, int e) { if(homog(i)!=1) { i=std(i); } list re=sres(i,0); //the resolution re=minres(re); //minimized resolution ideal ann=AnnExt_R(e,re); if(nvars(basering)-dim(std(ann))!=e) { return(ideal(1)); } return(ann); } /////////////////////////////////////////////////////////////////////////////// //computes the annihilator of Ext^n(R/i,R) with given resolution re //n is not necessarily the number of variables static proc AnnExt_R(int n,list re) { if(n1;ii--) { for(jj=1;jj0; n--) { m=coef(p,var(n)); if(m[1,1]!=1) { p=m[2,1]; break; } } if(deg(p)==0) { p=0; } return(p); } /////////////////////////////////////////////////////// // min_ass_prim_charsets // input: generators of an ideal PS and an integer cho // If cho=0, the given ordering of the variables is used. // Otherwise, the system tries to find an "optimal ordering", // which in some cases may considerably speed up the algorithm // output: the minimal associated primes of PS // algorithm: via characteriostic sets ////////////////////////////////////////////////////// static proc min_ass_prim_charsets (ideal PS, int cho) { if((cho<0) and (cho>1)) { ERROR(" must be 0 or 1"); } intvec saveopt=option(get); option(notWarnSB); list L; if(cho==0) { L=min_ass_prim_charsets0(PS); } else { L=min_ass_prim_charsets1(PS); } option(set,saveopt); return(L); } /////////////////////////////////////////////////////// // min_ass_prim_charsets0 // input: generators of an ideal PS // output: the minimal associated primes of PS // algorithm: via characteristic sets // the given ordering of the variables is used ////////////////////////////////////////////////////// static proc min_ass_prim_charsets0 (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--) { I=I+ini_mod(PHI[i][j]); } JS=std(PHI[i]); sizeJS=size(JS); for(j=size(I);j>0;j--) { II=0; sizeII=0; k=0; while(k<=sizeII) // successive saturation { op=option(get); option(returnSB); II=quotient(JS,I[j]); option(set,op); sizeII=size(II); if(sizeII==sizeJS) { for(k=1;k<=sizeII;k++) { if(leadexp(II[k])!=leadexp(JS[k])) break; } } JS=II; sizeJS=sizeII; } } PSI=insert(PSI,JS); } int sizePSI=size(PSI); // We eliminate redundant ideals for(i=1;i=1;i--) { if(size(PSI[i])==0) { PSI=delete(PSI,i); } } 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--) { I=I+ITPHI[i][j]; } JS=std(PHI[i]); sizeJS=size(JS); for(j=size(I);j>0;j--) { II=0; sizeII=0; k=0; while(k<=sizeII) // successive iteration { op=option(get); option(returnSB); II=quotient(JS,I[j]); option(set,op); //std // II=std(II); sizeII=size(II); if(sizeII==sizeJS) { for(k=1;k<=sizeII;k++) { if(leadexp(II[k])!=leadexp(JS[k])) break; } } JS=II; sizeJS=sizeII; } } PSI=insert(PSI,JS); } int sizePSI=size(PSI); // We eliminate redundant ideals for(i=1;i=1;i--) { if(size(PSI[i])==0) { PSI=delete(PSI,i); } } 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("ERROR: must be 0 or 1 or 2 or 3"); } ideal H=1; // The intersection of the primary components list U; // the leaves of the decomposition tree, i.e., // pairs consisting of a primary component of I // and the corresponding associated prime list W; // the non-leaf vertices in the decomposition tree. // every entry has 6 components: // 1- the vertex itself , i.e., a standard bais of the // given ideal I (type 1), or a standard basis of a // pseudo-primary component arising from // pseudo-primary decomposition (type 2), or a // standard basis of a remaining component arising from // pseudo-primary decomposition or extraction (type 3) // 2- the type of the vertex as indicated above // 3- the weighted_tree_depth of the vertex // 4- the tester of the vertex // 5- a standard basis of the associated prime // of a vertex of type 2, or 0 otherwise // 6- a list of pairs consisting of a standard // basis of a minimal associated prime ideal // of the father of the vertex and the // irreducible factors of the "minimal // divisor" of the seperator or extractor // corresponding to the prime ideal // as computed by the procedure minsat, // if the vertex is of type 3, or // the empty list otherwise ideal SI=std(I); if(SI[1]==1) // primdecSY(ideal(1)) { return(list()); } intvec save=option(get); option(notWarnSB); int ncolsSI=ncols(SI); int ncolsH=1; W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree int weighted_tree_depth; int i,j; int check; list V; // current vertex list VV; // new vertex list QQ; list WI; ideal Qi,SQ,SRest,fac; poly tester; while(1) { i=1; while(1) { while(i<=size(W)) // find vertex V of smallest weighted tree-depth { if (W[i][3]<=weighted_tree_depth) break; i++; } if (i<=size(W)) break; i=1; weighted_tree_depth++; } V=W[i]; W=delete(W,i); // delete V from W // now proceed by type of vertex V if (V[2]==2) // extraction needed { SQ,SRest,fac=extraction(V[1],V[5]); // standard basis of primary component, // standard basis of remaining component, // irreducible factors of // the "minimal divisor" of the extractor // as computed by the procedure minsat, check=0; for(j=1;j<=ncolsH;j++) { if (NF(H[j],SQ,1)!=0) // Q is not redundant { check=1; break; } } if(check==1) // Q is not redundant { QQ=list(); QQ[1]=list(SQ,V[5]); // primary component, associated prime, // i.e., standard bases thereof U=U+QQ; H=intersect(H,SQ); H=std(H); ncolsH=ncols(H); check=0; if(ncolsH==ncolsSI) { for(j=1;j<=ncolsSI;j++) { if(leadexp(H[j])!=leadexp(SI[j])) { check=1; break; } } } else { check=1; } if(check==0) // H==I => U is a primary decomposition { option(set,save); return(U); } } if (SRest[1]!=1) // the remaining component is not // the whole ring { if (rad_con(V[4],SRest)==0) // the new vertex is not the // root of a redundant subtree { VV[1]=SRest; // remaining component VV[2]=3; // pseudoprimdec_special VV[3]=V[3]+1; // weighted depth VV[4]=V[4]; // the tester did not change VV[5]=ideal(0); VV[6]=list(list(V[5],fac)); W=insert(W,VV,size(W)); } } } else { if (V[2]==3) // pseudo_prim_dec_special is needed { QQ,SRest=pseudo_prim_dec_special_charsets(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)); } } } } option(set,save); } ////////////////////////////////////////////////////////////////////////// // proc pseudo_prim_dec_charsets // input: Generators of an arbitrary ideal I, a standard basis SI of I, // and an integer choo // If choo=0, min_ass_prim_charsets with the given // ordering of the variables is used. // If choo=1, min_ass_prim_charsets with the "optimized" // ordering of the variables is used. // If choo=2, minAssPrimes from primdec.lib is used // If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used // output: a pseudo primary decomposition of I, i.e., a list // of pseudo primary components together with a standard basis of the // remaining component. Each pseudo primary component is // represented by a quadrupel: A standard basis of the component, // a standard basis of the corresponding associated prime, the // seperator of the component, and the irreducible factors of the // "minimal divisor" of the seperator as computed by the procedure minsat, // calls proc pseudo_prim_dec_i ////////////////////////////////////////////////////////////////////////// static proc pseudo_prim_dec_charsets (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--) { 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+2); // saturation SQi,f0,f,fac=minsat_ppd(SI,fac); I'=I',f; QP=SQi,L[i],f0,fac; // the quadrupel: // a standard basis of Q_i, // a standard basis of P_i, // sep_i, // irreducible factors of // the "minimal divisor" of the seperator // as computed by the procedure minsat, Q[i]=QP; } I'=std(I'); return (Q, I'); // I' = remaining component } //////////////////////////////////////////////////////////////// // proc extraction // input: A standard basis of a pseudo primary ideal I, and a standard // basis of the unique minimal associated prime P of I // output: an extraction of I, i.e., a standard basis of the primary // component Q of I with associated prime P, a standard basis of the // remaining component, and the irreducible factors of the // "minimal divisor" of the extractor as computed by the procedure minsat //////////////////////////////////////////////////////////////// static proc extraction (ideal SI, ideal SP) { 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) div 2;i++) { he=l[2*i-1],l[2*i]; re[i]=he; } return(re); } /////////////////////////////////////////////////////////////////////////////// static proc reconvList(list l) { int i; list re; for(i=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, list #) "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!) - For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring - For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument EXAMPLE: example primdecGTZ; shows an example " { if(size(#)>0) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=primdecGTZ(i); setring r; def li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k][2]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } 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 absPrimdecGTZ(ideal I, list #) "USAGE: absPrimdecGTZ(I); I ideal ASSUME: Ground field has characteristic 0. RETURN: a ring containing two lists: @code{absolute_primes}, the absolute prime components of I, and @code{primary_decomp}, the output of @code{primdecGTZ(I)}. The list absolute_primes has to be interpreted as follows: each entry describes a class of conjugated absolute primes, @format absolute_primes[i][1] the absolute prime component, absolute_primes[i][2] the number of conjugates. @end format The first entry of @code{absolute_primes[i][1]} is the minimal polynomial of a minimal finite field extension over which the absolute prime component is defined. For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring. For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument NOTE: Algorithm of Gianni/Trager/Zacharias combined with the @code{absFactorize} command. SEE ALSO: primdecGTZ; absFactorize EXAMPLE: example absPrimdecGTZ; shows an example " { if (char(basering) != 0) { ERROR("primdec.lib::absPrimdecGTZ is only implemented for "+ +"characteristic 0"); } if(size(#)>0) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithm automatically passes to the global case // hence prepare to go back to an appropriate new ring def r=basering; ideal max_of_r=maxideal(1); def s=changeord(list(list("dp",1:nvars(basering)))); setring s; def I=imap(r,I); def S=absPrimdecGTZ(I); setring S; ring r1=char(basering),var(nvars(r)+1),dp; def rS=r+r1; // move objects to appropriate ring and clean up setring rS; def max_of_r=imap(r,max_of_r); attrib(max_of_r,"isSB",1); def absolute_primes=imap(S,absolute_primes); def primary_decomp=imap(S,primary_decomp); if(!defined(keep_comp)) { ideal tempid; for(int k=size(absolute_primes);k>=1;k--) { tempid=absolute_primes[k][1]; tempid[1]=0; // ignore minimal polynomial if(size(reduce(lead(tempid),max_of_r))!=0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering absolute_primes=delete(absolute_primes,k); } } for(k=size(primary_decomp);k>=1;k--) { if(mindeg(std(lead(primary_decomp[k][2]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering primary_decomp=delete(primary_decomp,k); } } kill tempid; } export(primary_decomp); export(absolute_primes); return(rS); } if(minpoly!=0) { //return(algeDeco(i,0)); ERROR( "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal" ); } def R=basering; int n=nvars(R); list L=decomp(I,3); string newvar=L[1][3]; int k=find(newvar,",",find(newvar,",")+1); newvar=newvar[k+1..size(newvar)]; list lR=ringlist(R); int i,de,ii; intvec vv=1:n; //for(i=1;i<=n;i++){vv[i]=1;} list orst; orst[1]=list("dp",vv); orst[2]=list("dp",intvec(1)); orst[3]=list("C",0); lR[3]=orst; lR[2][n+1] = newvar; def Rz = ring(lR); setring Rz; list L=imap(R,L); list absolute_primes,primary_decomp; ideal I,M,N,K; M=maxideal(1); N=maxideal(1); poly p,q,f,g; map phi,psi; string tvar; for(i=1;i<=size(L);i++) { tvar=L[i][4]; ii=find(tvar,"+"); while(ii) { tvar=tvar[ii+1..size(tvar)]; ii=find(tvar,"+"); } for(ii=1;ii<=nvars(basering);ii++) { if(tvar==string(var(ii))) break; } I=L[i][2]; execute("K="+L[i][3]+";"); p=K[1]; q=K[2]; execute("f="+L[i][4]+";"); g=2*var(ii)-f; M[ii]=f; N[ii]=g; de=deg(p); psi=Rz,M; phi=Rz,N; I=phi(I),p,q; I=std(I); absolute_primes[i]=list(psi(I),de); primary_decomp[i]=list(L[i][1],L[i][2]); } export(primary_decomp); export(absolute_primes); setring R; dbprint( printlevel-voice+3," // 'absPrimdecGTZ' created a ring, in which two lists absolute_primes (the // absolute prime components) and primary_decomp (the primary and prime // components over the current basering) are stored. // To access the list of absolute prime components, type (if the name S was // assigned to the return value): setring S; absolute_primes; "); return(Rz); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; def S = absPrimdecGTZ(i); setring S; absolute_primes; } /////////////////////////////////////////////////////////////////////////////// proc primdecSY(ideal i, list #) "USAGE: primdecSY(I, c); I ideal, c int (optional) RETURN: a list pr of primary ideals and their associated primes: @format pr[i][1] the i-th primary component, pr[i][2] the i-th prime component. @end format NOTE: Algorithm of Shimoyama/Yokoyama. @format if c=0, the given ordering of the variables is used, if c=1, minAssChar tries to use an optimal ordering (default), if c=2, minAssGTZ is used, if c=3, minAssGTZ and facstd are used. @end format For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring. For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as third argument EXAMPLE: example primdecSY; shows an example " { if(size(#)>1) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=primdecSY(i); setring r; def li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k][2]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } i=simplify(i,2); if ((i[1]==0)||(i[1]==1)) { list L=list(ideal(i[1]),ideal(i[1])); return(list(L)); } if(minpoly!=0) { return(algeDeco(i,1)); } if (size(#)!=0) { 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[, l]); I ideal, l list (optional) @* Optional parameters in list l (can be entered in any order): @* 0, \"facstd\" -> uses facstd to first decompose the ideal (default) @* 1, \"noFacstd\" -> does not use facstd @* \"GTZ\" -> the original algorithm by Gianni, Trager and Zacharias is used @* \"SL\" -> GTZ algorithm with modificiations by Laplagne is used (default) RETURN: a list, the minimal associated prime ideals of I. NOTE: - Designed for characteristic 0, works also in char k > 0 based on an algorithm of Yokoyama - For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring - For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as second argument EXAMPLE: example minAssGTZ; shows an example " { if(size(#)>0) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=minAssGTZ(i); setring r; def li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } int j; string algorithm; string facstdOption; int useFac; // Set input parameters algorithm = "SL"; // Default: SL algorithm facstdOption = "facstd"; if(size(#) > 0) { int valid; for(j = 1; j <= size(#); j++) { valid = 0; if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) { if (#[j] == 1) {facstdOption = "noFacstd"; valid = 1;} // If #[j] == 1, facstd is not used. if (#[j] == 0) {facstdOption = "facstd"; valid = 1;} // If #[j] == 0, facstd is used. } if(typeof(#[j]) == "string") { if((#[j] == "GTZ") || (#[j] == "SL")) { algorithm = #[j]; valid = 1; } if((#[j] == "noFacstd") || (#[j] == "facstd")) { facstdOption = #[j]; valid = 1; } } if(valid == 0) { dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); } } } if(minpoly!=0) { return(algeDeco(i,2)); } list result = minAssPrimes(i, facstdOption, algorithm); return(result); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; list pr = minAssGTZ(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc minAssChar(ideal i, list #) "USAGE: minAssChar(I[,c]); i ideal, c int (optional). 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. @* For local orderings, the result is considered in the localization of the polynomial ring, not in the power series ring For local and mixed orderings, the decomposition in the corresponding global ring is returned if the string 'global' is specified as third argument EXAMPLE: example minAssChar; shows an example " { if(size(#)>1) { int keep_comp=1; } if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // decompose and go back list li=minAssChar(i); setring r; def li=imap(s,li); // clean up if(!defined(keep_comp)) { for(int k=size(li);k>=1;k--) { if(mindeg(std(lead(li[k]))[1])==0) { // 1 contained in ideal, i.e. component does not meet origin in local ordering li=delete(li,k); } } } return(li); } if (size(#)>0) { return(min_ass_prim_charsets(i,#[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 (with modifications by Laplagne) and Kemper is used. Works also in positive characteristic (Kempers algorithm). EXAMPLE: example equiRadical; shows an example " { if(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } return(radical(i, 1)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; ideal pr= equiRadical(i); pr; } /////////////////////////////////////////////////////////////////////////////// proc radical(ideal i, list #) "USAGE: radical(I[, l]); I ideal, l list (optional) @* Optional parameters in list l (can be entered in any order): @* 0, \"fullRad\" -> full radical is computed (default) @* 1, \"equiRad\" -> equiRadical is computed @* \"KL\" -> Krick/Logar algorithm is used @* \"SL\" -> modifications by Laplagne are used (default) @* \"facstd\" -> uses facstd to first decompose the ideal (default for non homogeneous ideals) @* \"noFacstd\" -> does not use facstd (default for homogeneous ideals) RETURN: ideal, the radical of I (or the equiradical if required in the input parameters) NOTE: A combination of the algorithms of Krick/Logar (with modifications by Laplagne) and Kemper is used. Works also in positive characteristic (Kempers algorithm). EXAMPLE: example radical; shows an example " { dbprint(printlevel - voice, "Radical, version 2006.05.08"); if(attrib(basering,"global")!=1) { // algorithms only work in global case! // pass to appropriate global ring def r=basering; def s=changeord(list(list("dp",1:nvars(basering)))); setring s; ideal i=imap(r,i); // compute radical and go back def j=radical(i); setring r; def j=imap(s,j); return(j); } if(size(i) == 0){return(ideal(0));} int j; def P0 = basering; list Pl=ringlist(P0); intvec dp_w; for(j=nvars(P0);j>0;j--) {dp_w[j]=1;} Pl[3]=list(list("dp",dp_w),list("C",0)); def @P=ring(Pl); setring @P; ideal i=imap(P0,i); int il; string algorithm; int useFac; // Set input parameters algorithm = "SL"; // Default: SL algorithm il = 0; // Default: Full radical (not only equiRadical) if (homog(i) == 1) { // Default: facStd is used, except if the ideal is homogeneous. useFac = 0; } else { useFac = 1; } if(size(#) > 0) { int valid; for(j = 1; j <= size(#); j++) { valid = 0; if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) { il = #[j]; // If il == 1, equiRadical is computed valid = 1; } if(typeof(#[j]) == "string") { if(#[j] == "KL") { algorithm = "KL"; valid = 1; } if(#[j] == "SL") { algorithm = "SL"; valid = 1; } if(#[j] == "noFacstd") { useFac = 0; valid = 1; } if(#[j] == "facstd") { useFac = 1; valid = 1; } if(#[j] == "equiRad") { il = 1; valid = 1; } if(#[j] == "fullRad") { il = 0; valid = 1; } } if(valid == 0) { dbprint(1, "Warning! The following input parameter was not recognized:", #[j]); } } } ideal rad = 1; intvec op = option(get); list qr = simplifyIdeal(i); map phi = @P, qr[2]; option(redSB); i = groebner(qr[1]); option(set, op); int di = dim(i); if(di == 0) { i = zeroRad(i, qr[1]); option(redSB); i=interred(phi(i)); option(set, op); setring(P0); i=imap(@P,i); return(i); } option(redSB); list pr; if(useFac == 1) { pr = facstd(i); } else { pr = i; } option(set, op); int s = size(pr); if(useFac == 1) { dbprint(printlevel - voice, "Number of components returned by facstd: ", s); } for(j = 1; j <= s; j++) { attrib(pr[s + 1 - j], "isSB", 1); if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il)) { // SL Debug messages dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]); dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j])); if(algorithm == "KL") { rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il)); } if(algorithm == "SL") { rad = intersect(rad, radicalSL(pr[s + 1 - j], il)); } } else { // SL Debug dbprint(printlevel-voice, "The radical of this component is not needed."); dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 1))); dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j])); dbprint(printlevel-voice, "il", il); } } rad=interred(phi(rad)); setring(P0); i=imap(@P,rad); return(i); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly p = z2+1; poly q = z3+2; ideal i = p*q^2,y-z2; ideal pr = radical(i); pr; } /////////////////////////////////////////////////////////////////////////////// // // Computes the radical of I using KL algorithm. // The only difference with the previous implementation of KL algorithm is // that now it uses block dp instead of lp ordering for the reduction to the // zerodimensional case. // The reduction step has been moved to the new routine radicalReduction, so that it can be // used also by radicalSL procedure. // static proc radicalKL(ideal I, ideal ser, list #) { // ideal I The ideal for which the radical is computed // ideal ser Used to reduce components already obtained // list # If #[1] = 1, equiradical is computed. // I needs to be a Groebner basis. if (attrib(I, "isSB") != 1) { I = groebner(I); } ideal rad; // The radical int allIndep = 1; // All max independent sets are used list result = radicalReduction(I, ser, allIndep, #); int done = result[3]; rad = result[1]; if (done == 0) { rad = intersect(rad, radicalKL(result[2], ideal(1), #)); } return(rad); } /////////////////////////////////////////////////////////////////////////////// // // Computes the radical of I via Laplagne algorithm, using zerodimensional radical in // the zero dimensional case. // For the reduction to the zerodimensional case, it uses the procedure // radical, with some modifications to avoid the recursion. // static proc radicalSL(ideal I, list #) // Input = I, ideal // #, list. If #[1] = 1, then computes only the equiradical. // Output = (P, primaryDec) where P = rad(I) and primaryDec is the list of the radicals // obtained in intermediate steps. { ideal rad = 1; ideal equiRad = 1; list primes; int k; // Counter int il; // If il = 1, only the equiradical is required. int iDim; // The dimension of I int stop = 0; // Checks if the radical has been obtained if (attrib(I, "isSB") != 1) { I = groebner(I); } iDim = dim(I); // Checks if only equiradical is required if (size(#) > 0) { il = #[1]; } while(stop == 0) { dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals."); primes = radicalSLIteration(I, rad); // A list of primes or intersections of primes, not included in P dbprint (printlevel - voice, "// Output of Iteration Step:"); dbprint (printlevel - voice, primes); if (size(primes) > 0) { dbprint (printlevel - voice, "// We intersect P with the ideal just obtained."); for(k = 1; k <= size(primes); k++) { rad = intersect(rad, primes[k]); if (il == 1) { if (attrib(primes[k], "isSB") != 1) { primes[k] = groebner(primes[k]); } if (iDim == dim(primes[k])) { equiRad = intersect(equiRad, primes[k]); } } } } else { stop = 1; } } if (il == 0) { return(rad); } else { return(equiRad); } } ////////////////////////////////////////////////////////////////////////// // Based on radicalKL. // It contains all of old version of proc radicalKL except the recursion call. // // Output: // #1 -> output ideal, the part of the radical that has been computed // #2 -> complementary ideal, the part of the ideal I whose radical remains to be computed // = (I, h) in KL algorithm // This is not used in the new algorithm. It is part of KL algorithm // #3 -> done, 1: output = radical, there is no need to continue // 0: radical = output \cap \sqrt{complementary ideal} // This is not used in the new algorithm. It is part of KL algorithm static proc radicalReduction(ideal I, ideal ser, int allIndep, list #) { // allMaximal 1 -> Indicates that the reduction to the zerodim case // must be done for all indep set of the leading terms ideal // 0 -> Otherwise // ideal ser Only for radicalKL. (Same as in radicalKL) // list # Only for radicalKL (If #[1] = 1, // only equiradical is required. // It is used to set the value of done.) 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); // Computes the dimension of I int homo = homog(I); // Finds out if I is homogeneous ideal rad = ideal(1); // The unit ideal ideal te = ser; if(size(#) > 0) { @wr = #[1]; } if(homo == 1) { for(@n = 1; @n <= nvars(basering); @n++) { @w[@n] = ord(var(@n)); } @hilb = hilb(I, 1, @w); } // SL 2006.04.11 1 Debug messages dbprint(printlevel-voice, "//Computes the radical of the ideal:", I); // SL 2006.04.11 2 Debug messages //--------------------------------------------------------------------------- //j is the ring //--------------------------------------------------------------------------- if (jdim==-1) { return(ideal(1), ideal(1), 1); } //--------------------------------------------------------------------------- //the zero-dimensional case //--------------------------------------------------------------------------- if (jdim==0) { return(zeroRad(I), ideal(1), 1); } //------------------------------------------------------------------------- //search for a maximal independent set indep,i.e. //look for subring such that the intersection with the ideal is zero //j intersected with K[var(indep[3]+1),...,var(nvar)] is zero, //indep[1] is the new varstring, indep[2] the string for the block-ordering //------------------------------------------------------------------------- // SL 2006-04-24 1 If allIndep = 0, then it only computes one maximal // independent set. // This looks better for the new algorithm but not for KL // algorithm list parameters = allIndep; indep = newMaxIndependSetDp(I, parameters); // SL 2006-04-24 2 for(@m = 1; @m <= size(indep); @m++) { if((indep[@m][1] == varstr(basering)) && (@m == 1)) //this is the good case, nothing to do, just to have the same notations //change the ring { 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) < jdim)) { setring @P; break; } for (lauf = 1; lauf <= size(@j); lauf++) { fett[lauf] = size(@j[lauf]); } //------------------------------------------------------------------------ // We have now the following situation: // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass // to this quotientring, j is there still a standardbasis, the // leading coefficients of the polynomials there (polynomials in // K[var(nnp+1),..,var(nva)]) are collected in the list h, // we need their LCM, gh, because of the following: // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..] // intersected with K[var(1),...,var(nva)] is (j:gh^n) // on the other hand j = ((j, gh^n) intersected with (j : gh^n)) //------------------------------------------------------------------------ // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..] // and the map phi:K[var(1),...,var(nva)] -----> // K(var(nnpr+1),..,var(nva))[..the rest..] //------------------------------------------------------------------------ quotring = prepareQuotientRingDp(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 LCM(h[1],...) = gh for saturation ideal @h; if(deg(@j[1]) > 0) { for(@n = 1; @n <= size(@j); @n++) { @h[@n] = leadcoef(@j[@n]); } op = option(get); option(redSB); @j = std(@j); //to obtain a reduced standardbasis option(set, op); // SL 1 Debug messages dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j))); ideal zero_rad = zeroRad(@j); dbprint(printlevel - voice, "zero_rad passed"); // SL 2 } else { ideal zero_rad = ideal(1); } // We need the intersection of the ideals in the list quprimary with the // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal // but fi polynomials, then the intersection of q with the polynomialring // is the saturation of the ideal generated by f1,...,fr with respect to // h which is the lcm of the leading coefficients of the fi considered in // the quotientring: this is coded in saturn zero_rad = std(zero_rad); ideal hpl; for(@n = 1; @n <= size(zero_rad); @n++) { hpl = hpl, leadcoef(zero_rad[@n]); } //------------------------------------------------------------------------ // We leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] // back to the polynomialring //------------------------------------------------------------------------ setring @P; collectrad = imap(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, 6); @q = 1; for(lauf = 1; lauf <= size(fac); lauf++) { @q = @q * fac[lauf]; } op = option(get); option(returnSB); option(redSB); I = quotient(I + ideal(@q), rad); attrib(I, "isSB", 1); option(set, op); } if((deg(rad[1]) > 0) && (deg(collectrad[1]) > 0)) { rad = intersect(rad, collectrad); te = intersect(te, collectrad); te = simplify(reduce(te, I, 1), 2); } else { if(deg(collectrad[1]) > 0) { rad = collectrad; te = intersect(te, collectrad); te = simplify(reduce(te, I, 1), 2); } } if((dim(I) < jdim)||(size(te) == 0)) { break; } if(homo==1) { @hilb = hilb(I, 1, @w); } } // SL 2006.04.11 1 Debug messages dbprint (printlevel-voice, "// Part of the Radical already computed:", rad); dbprint (printlevel-voice, "// Dimension:", dim(groebner(rad))); // SL 2006.04.11 2 Debug messages // SL 2006.04.21 1 New variable "done". // It tells if the radical is already computed or // if it still has to be computed the radical of the new ideal I int done; if(((@wr == 1) && (dim(I) 0) { allMaximal = #[1]; } else { allMaximal = 1; } int nMax; if (allMaximal == 1) { nMax = size(v); } else { nMax = 1; } for(n = 1; n <= nMax; n++) // SL 2006.04.21 2 { di = 0; var1 = ""; var2 = ""; for(k = 1; k <= size(v[n]); k++) { if(v[n][k] != 0) { di++; var2 = var2 + "var(" + string(k) + "), "; } else { var1 = var1 + "var(" + string(k) + "), "; } } if(di > 0) { var1 = var1 + var2; var1 = var1[1..size(var1) - 2]; // The "- 2" removes the trailer comma hilf[1] = var1; // SL 2006.21.04 1 The order is now block dp instead of lp hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")"; // SL 2006.21.04 2 hilf[3] = di; resu[n] = hilf; } else { resu[n] = varstr(basering), ordstr(basering), 0; } } return(resu); } example { "EXAMPLE:"; echo = 2; ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp; ideal i = ea - fbg, fa + be, ec - fdg, fc + de; i = std(i); list l = newMaxIndependSetDp(i); l; i = i, g; l = newMaxIndependSetDp(i); l; ring s = 0, (x, y, z), lp; ideal i = z, yx; list l = newMaxIndependSetDp(i); l; } /////////////////////////////////////////////////////////////////////////////// // based on prepareQuotientring // The order returned is now (C, dp) instead of (C, lp) static proc prepareQuotientRingDp (int nnp) "USAGE: prepareQuotientRingDp(nnp); nnp int RETURN: string = to define Kvar(nnp+1),...,var(nvars)[..rest ] NOTE: EXAMPLE: example prepareQuotientRingDp; 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); } // SL 2006-04-21 1 The order returned is now (C, dp) instead of (C, lp) quotring = quotring + "), (C, dp);"; // SL 2006-04-21 2 return(quotring); } example { "EXAMPLE:"; echo = 2; ring s1=(0,x),(a,b,c,d,e,f,g),lp; def @Q=basering; list l= prepareQuotientRingDp(3); l; execute(l[1]); execute(l[2]); basering; phi; setring @Q; } /////////////////////////////////////////////////////////////////////////////// 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(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } ideal j=std(i); int cod=nvars(basering)-dim(j); int e; list er; ideal ann; if(homog(i)==1) { 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(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } ideal j=groebner(i); int cod=nvars(basering)-dim(j); 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) div 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(attrib(basering,"global")!=1) { ERROR( "// Not implemented for this ordering, please change to global ordering." ); } def R=basering; poly q; int j,time; matrix m; list re; poly va=var(1); ideal J=groebner(I); ideal ba=kbase(J); int d=vdim(J); dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d)); //------ compute matrix of multiplication on R/I with generic element p ----- int e=nvars(basering); poly p=randomLast(100)[e]+random(-50,50); //the generic element matrix n[d][d]; time = timer; for(j=2;j<=e;j++) { va=va*var(j); } for(j=1;j<=d;j++) { q=reduce(p*ba[j],J); m=coeffs(q,ba,va); n[j,1..d]=m[1..d,1]; } dbprint(printlevel-voice+2, "// time for computing multiplication matrix (with generic element) : "+ string(timer-time)); //---------------- compute characteristic polynomial of matrix -------------- execute("ring P1=("+charstr(R)+"),T,dp;"); matrix n=imap(R,n); time = timer; poly charpol=det(n-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; } /////////////////////////////////////////////////////////////////////////////// static proc newDecompStep(ideal i, list #) "USAGE: newDecompStep(i); i ideal (for primary decomposition) newDecompStep(i,1); (for the associated primes of dimension of i) newDecompStep(i,2); (for the minimal associated primes) newDecompStep(i,3); (for the absolute primary decomposition (not tested!)) "oneIndep"; (for using only one max indep set) "intersect"; (returns alse the intersection of the components founded) RETURN: list = list of primary ideals and their associated primes (at even positions in the list) (resp. a list of the minimal associated primes) NOTE: Algorithm of Gianni/Trager/Zacharias EXAMPLE: example newDecompStep; shows an example " { intvec op,@vv; def @P = basering; list primary,indep,ltras; intvec @vh,isat,@w; int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn; ideal peek=i; ideal ser,tras; list data; list result; intvec @hilb; int isS=(attrib(i,"isSB")==1); // Debug dbprint(printlevel - voice, "newDecompStep, v2.0"); string indepOption = "allIndep"; string intersectOption = "noIntersect"; if(size(#)>0) { int count = 1; if(typeof(#[count]) == "string") { if ((#[count] == "oneIndep") or (#[count] == "allIndep")) { indepOption = #[count]; count++; } } if(typeof(#[count]) == "string") { if ((#[count] == "intersect") or (#[count] == "noIntersect")) { intersectOption = #[count]; count++; } } if((typeof(#[count]) == "int") or (typeof(#[count]) == "number")) { if ((#[count]==1)||(#[count]==2)||(#[count]==3)) { @wr=#[count]; if(@wr==3){abspri = 1; @wr = 0;} count++; } } if(size(#)>count) { seri=1; peek=#[count + 1]; ser=#[count + 2]; } } if(abspri) { list absprimary,abskeep,absprimarytmp,abskeeptmp; } homo=homog(i); if(homo==1) { if(attrib(i,"isSB")!=1) { //ltras=mstd(i); tras=groebner(i); ltras=tras,tras; 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[2]=maxideal(1); l[1]=maxideal(1); if (intersectOption == "intersect") { return(list(l, maxideal(1))); } else { return(l); } } if (intersectOption == "intersect") { return(list(primary, primary[1])); } else { return(primary); } } for(@n=1;@n<=nvars(basering);@n++) { @w[@n]=ord(var(@n)); } @hilb=hilb(tras,1,@w); intvec keephilb=@hilb; } //---------------------------------------------------------------- //i is the zero-ideal //---------------------------------------------------------------- if(size(i)==0) { primary=i,i; if (intersectOption == "intersect") { return(list(primary, i)); } else { return(primary); } } //---------------------------------------------------------------- //pass to the lexicographical ordering and compute a standardbasis //---------------------------------------------------------------- int lp=islp(); 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)&&(!abspri)) { ideal fried; @n=size(@j); for(@k=1;@k<=@n;@k++) { if(deg(lead(@j[@k]))==1) { fried[size(fried)+1]=@j[@k]; @j[@k]=0; } } if(size(fried)==nvars(basering)) { setring @P; primary[1]=i; primary[2]=i; if (intersectOption == "intersect") { return(list(primary, i)); } else { return(primary); } } if(size(fried)>0) { string newva; string newma; for(@k=1;@k<=nvars(basering);@k++) { @n1=0; for(@n=1;@n<=size(fried);@n++) { if(leadmonom(fried[@n])==var(@k)) { @n1=1; break; } } if(@n1==0) { newva=newva+string(var(@k))+","; newma=newma+string(var(@k))+","; } else { newma=newma+string(0)+","; } } newva[size(newva)]=")"; newma[size(newma)]=";"; execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;"); execute("map @kappa=gnir,"+newma); ideal @j= @kappa(@j); @j=simplify(@j, 2); attrib(@j,"isSB",1); result = newDecompStep(@j, indepOption, intersectOption, @wr); if (intersectOption == "intersect") { list pr = result[1]; ideal intersection = result[2]; } else { list pr = result; } setring gnir; list pr=imap(@deirf,pr); for(@k=1;@k<=size(pr);@k++) { @j=pr[@k]+fried; pr[@k]=@j; } if (intersectOption == "intersect") { ideal intersection = imap(@deirf, intersection); @j = intersection + fried; intersection = @j; } setring @P; if (intersectOption == "intersect") { return(list(imap(gnir,pr), imap(gnir,intersection))); } else { return(imap(gnir,pr)); } } } //---------------------------------------------------------------- //j is the ring //---------------------------------------------------------------- if (dim(@j)==-1) { setring @P; primary=ideal(1),ideal(1); if (intersectOption == "intersect") { return(list(primary, ideal(1))); } else { return(primary); } } //---------------------------------------------------------------- // the case of one variable //---------------------------------------------------------------- if(nvars(basering)==1) { list fac=factor(@j[1]); list gprimary; poly generator; ideal gIntersection; for(@k=1;@k<=size(fac[1]);@k++) { if(@wr==0) { gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } else { gprimary[2*@k-1]=ideal(fac[1][@k]); gprimary[2*@k]=ideal(fac[1][@k]); } if (intersectOption == "intersect") { generator = generator * fac[1][@k]; } } if (intersectOption == "intersect") { gIntersection = generator; } setring @P; primary=fetch(gnir,gprimary); if (intersectOption == "intersect") { ideal intersection = fetch(gnir,gIntersection); } //HIER if(abspri) { list resu,tempo; string absotto; for(ab=1;ab<=size(primary) div 2;ab++) { absotto= absFactorize(primary[2*ab][1],77); tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); resu[ab]=tempo; } primary=resu; intersection = 1; for(ab=1;ab<=size(primary);ab++) { intersection = intersect(intersection, primary[ab][2]); } } if (intersectOption == "intersect") { return(list(primary, intersection)); } else { return(primary); } } //------------------------------------------------------------------ //the zero-dimensional case //------------------------------------------------------------------ if (dim(@j)==0) { op=option(get); option(redSB); list gprimary= newZero_decomp(@j,ser,@wr); setring @P; primary=fetch(gnir,gprimary); if(size(ser)>0) { primary=cleanPrimary(primary); } //HIER if(abspri) { list resu,tempo; string absotto; for(ab=1;ab<=size(primary) div 2;ab++) { absotto= absFactorize(primary[2*ab][1],77); tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering))); resu[ab]=tempo; } primary=resu; } if (intersectOption == "intersect") { return(list(primary, fetch(gnir,@j))); } else { 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 = newMaxIndependSetLp(@j, indepOption); for(@m=1;@m<=size(allindep);@m++) { if(allindep[@m][3]==jdim) { di++; indep[di]=allindep[@m]; } else { lauf++; restindep[lauf]=allindep[@m]; } } } else { indep = newMaxIndependSetLp(@j, indepOption); } ideal jkeep=@j; if(ordstr(@P)[1]=="w") { 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; //Aenderung================ ideal @Ptest=1; //========================= di=dim(jwork); keepdi=di; ser = 1; setring gnir; for(@m=1; @m<=size(indep); @m++) { data[1] = indep[@m]; result = newReduction(@j, ser, @hilb, @w, jdim, abspri, @wr, data); quprimary = quprimary + result[1]; if(abspri) { absprimary = absprimary + result[2]; abskeep = abskeep + result[3]; } @h = result[5]; ser = result[4]; if(size(@h)>0) { //--------------------------------------------------------------- //we change to @Phelp to have the ordering dp for saturation //--------------------------------------------------------------- setring @Phelp; @h=imap(gnir,@h); //Aenderung================================== if(defined(@LL)){kill @LL;} list @LL=minSat(jwork,@h); @Ptest=intersect(@Ptest,@LL[1]); ser = intersect(ser, @LL[1]); //=========================================== if(@wr!=1) { //Aenderung================================== @q=@LL[2]; //=========================================== //@q=minSat(jwork,@h)[2]; } else { fac=ideal(0); for(lauf=1;lauf<=ncols(@h);lauf++) { if(deg(@h[lauf])>0) { fac=fac+factorize(@h[lauf],1); } } fac=simplify(fac,6); @q=1; for(lauf=1;lauf<=size(fac);lauf++) { @q=@q*fac[lauf]; } } jwork = std(jwork,@q); keepdi = dim(jwork); if(keepdi < di) { setring gnir; @j = imap(@Phelp, jwork); ser = imap(@Phelp, ser); break; } if(homo == 1) { @hilb = hilb(jwork, 1, @w); } setring gnir; ser = imap(@Phelp, ser); @j = imap(@Phelp, jwork); } } if((size(quprimary)==0)&&(@wr==1)) { @j=ideal(1); quprimary[1]=ideal(1); quprimary[2]=ideal(1); } if((size(quprimary)==0)) { keepdi = di - 1; quprimary[1]=ideal(1); quprimary[2]=ideal(1); } //--------------------------------------------------------------- //notice that j=sat(j,gh) intersected with (j,gh^n) //we finished with sat(j,gh) and have to start with (j,gh^n) //--------------------------------------------------------------- if((deg(@j[1])!=0)&&(@wr!=1)) { if(size(quprimary)>0) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); if(@wr==0) { //Aenderung==================================================== //HIER STATT DURCHSCHNITT SATURIEREN! ideal htest=@Ptest; /* ideal htest=hquprimary[1]; for (@n1=2;@n1<=size(hquprimary) div 2;@n1++) { htest=intersect(htest,hquprimary[2*@n1-1]); } */ //============================================================= } else { ideal htest=hquprimary[2]; for (@n1=2;@n1<=size(hquprimary) div 2;@n1++) { htest=intersect(htest,hquprimary[2*@n1]); } } if(size(ser)>0) { ser=intersect(htest,ser); } else { ser=htest; } setring gnir; ser=imap(@Phelp,ser); } if(size(reduce(ser,peek,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= newZero_decomp(@j,ser,@wr); //HIER if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); //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); //================NEU========================================= if(deg(quprimary[1][1])<=0){ @n2=0; } //============================================================ @n3=@n2; for(@n1=1;@n1<=size(collectprimary) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } if(@n2>=@n3+2) { setring @Phelp; ser=imap(gnir,ser); hquprimary=imap(gnir,quprimary); for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(@wr==0) { ser=intersect(ser,hquprimary[2*@n-1]); } else { ser=intersect(ser,hquprimary[2*@n]); } } setring gnir; ser=imap(@Phelp,ser); } // } } //HIER if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { if (deg(quprimary[2*ab][1])!=0) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } } quprimary=resu; @wr=3; } if(size(reduce(ser,peek,1))!=0) { if(@wr>0) { // The following line was dropped to avoid the recursion step: //htprimary=newDecompStep(@j,@wr,peek,ser); htprimary = list(); } else { // The following line was dropped to avoid the recursion step: //htprimary=newDecompStep(@j,peek,ser); htprimary = list(); } // here we collect now both results primary(sat(j,gh)) // and primary(j,gh^n) @n=size(quprimary); if (deg(quprimary[1][1])<=0) { @n=0; } for (@k=1;@k<=size(htprimary);@k++) { quprimary[@n+@k]=htprimary[@k]; } } } } else { if(abspri) { list resu,tempo; for(ab=1;ab<=size(quprimary) div 2;ab++) { tempo=quprimary[2*ab-1],quprimary[2*ab], absprimary[ab],abskeep[ab]; resu[ab]=tempo; } quprimary=resu; } } //--------------------------------------------------------------------------- //back to the ring we started with //the final result: primary //--------------------------------------------------------------------------- setring @P; primary=imap(gnir,quprimary); if (intersectOption == "intersect") { return(list(primary, imap(gnir, ser))); } else { return(primary); } } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; list pr= newDecompStep(i); pr; testPrimary( pr, i); } // This was part of proc decomp. // In proc newDecompStep, used for the computation of the minimal associated primes, // this part was separated as a soubrutine to make the code more clear. // Also, since the reduction is performed twice in proc newDecompStep, it should use both times this routine. // This is not yet implemented, since the reduction is not exactly the same and some changes should be made. static proc newReduction(ideal @j, ideal ser, intvec @hilb, intvec @w, int jdim, int abspri, int @wr, list data) { string @va; string quotring; intvec op; intvec @vv; def gnir = basering; ideal isat=0; int @n; int @n1 = 0; int @n2 = 0; int @n3 = 0; int homo = homog(@j); int lauf; int @k; list fett; int keepdi; list collectprimary; list lsau; list lnew; ideal @h; list indepInfo = data[1]; list quprimary = list(); //if(abspri) //{ int ab; list absprimarytmp,abskeeptmp; list absprimary, abskeep; //} // Debug dbprint(printlevel - voice, "newReduction, v2.0"); if((indepInfo[1]==varstr(basering))) // &&(@m==1) //this is the good case, nothing to do, just to have the same notations //change the ring { 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)); //Aenderung============== //if(@m==1) //{ // @j=fetch(@P,i); //} //======================= execute("ring gnir1 = ("+charstr(basering)+"),("+indepInfo[1]+"),(" +indepInfo[2]+");"); execute("map phi=gnir,"+@va+";"); op=option(get); option(redSB); if(homo==1) { ideal @j=std(phi(@j),@hilb,@w); } else { ideal @j=groebner(phi(@j)); } ideal ser=phi(ser); option(set,op); } if((deg(@j[1])==0)||(dim(@j)K(var(nnpr+1),..,var(nva))[..rest..] //------------------------------------------------------------------------ quotring=prepareQuotientring(nvars(basering)-indepInfo[3]); //--------------------------------------------------------------------- //we pass to the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //--------------------------------------------------------------------- ideal @jj=lead(@j); //!! vorn vereinbaren execute(quotring); ideal @jj=imap(gnir1,@jj); @vv=clearSBNeu(@jj,fett); //!! vorn vereinbaren setring gnir1; @k=size(@j); for (lauf=1;lauf<=@k;lauf++) { if(@vv[lauf]==1) { @j[lauf]=0; } } @j=simplify(@j,2); setring quring; // @j considered in the quotientring ideal @j=imap(gnir1,@j); ideal ser=imap(gnir1,ser); kill gnir1; //j is a standardbasis in the quotientring but usually not minimal //here it becomes minimal attrib(@j,"isSB",1); //we need later ggt(h[1],...)=gh for saturation ideal @h; if(deg(@j[1])>0) { for(@n=1;@n<=size(@j);@n++) { @h[@n]=leadcoef(@j[@n]); } //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] op=option(get); option(redSB); int zeroMinAss = @wr; if (@wr == 2) {zeroMinAss = 1;} list uprimary= newZero_decomp(@j, ser, zeroMinAss); //HIER if(abspri) { ideal II; ideal jmap; map sigma; nn=nvars(basering); map invsigma=basering,maxideal(1); for(ab=1;ab<=size(uprimary) div 2;ab++) { II=uprimary[2*ab]; attrib(II,"isSB",1); if(deg(II[1])!=vdim(II)) { jmap=randomLast(50); sigma=basering,jmap; jmap[nn]=2*var(nn)-jmap[nn]; invsigma=basering,jmap; II=groebner(sigma(II)); } absprimarytmp[ab]= absFactorize(II[1],77); II=var(nn); abskeeptmp[ab]=string(invsigma(II)); invsigma=basering,maxideal(1); } } option(set,op); } else { list uprimary; uprimary[1]=ideal(1); uprimary[2]=ideal(1); } //we need the intersection of the ideals in the list quprimary with the //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal //but fi polynomials, then the intersection of q with the polynomialring //is the saturation of the ideal generated by f1,...,fr with respect to //h which is the lcm of the leading coefficients of the fi considered in //in the quotientring: this is coded in saturn list saturn; ideal hpl; for(@n=1;@n<=size(uprimary);@n++) { uprimary[@n]=interred(uprimary[@n]); // temporary fix hpl=0; for(@n1=1;@n1<=size(uprimary[@n]);@n1++) { hpl=hpl,leadcoef(uprimary[@n][@n1]); } saturn[@n]=hpl; } //-------------------------------------------------------------------- //we leave the quotientring K(var(nnp+1),..,var(nva))[..the rest..] //back to the polynomialring //--------------------------------------------------------------------- setring gnir; collectprimary=imap(quring,uprimary); lsau=imap(quring,saturn); @h=imap(quring,@h); kill quring; @n2=size(quprimary); @n3=@n2; for(@n1=1;@n1<=size(collectprimary) div 2;@n1++) { if(deg(collectprimary[2*@n1][1])>0) { @n2++; quprimary[@n2]=collectprimary[2*@n1-1]; lnew[@n2]=lsau[2*@n1-1]; @n2++; lnew[@n2]=lsau[2*@n1]; quprimary[@n2]=collectprimary[2*@n1]; if(abspri) { absprimary[@n2 div 2]=absprimarytmp[@n1]; abskeep[@n2 div 2]=abskeeptmp[@n1]; } } } //here the intersection with the polynomialring //mentioned above is really computed for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++) { if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; quprimary[2*@n]=quprimary[2*@n-1]; } else { if(@wr==0) { quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1]; } quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; } } return(quprimary, absprimary, abskeep, ser, @h); } //////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // Based on minAssGTZ proc minAss(ideal i,list #) "USAGE: minAss(I[, l]); i ideal, l list (optional) of parameters, same as minAssGTZ RETURN: a list, the minimal associated prime ideals of I. NOTE: Designed for characteristic 0, works also in char k > 0 based on an algorithm of Yokoyama EXAMPLE: example minAss; shows an example " { return(minAssGTZ(i,#)); } example { "EXAMPLE:"; echo = 2; ring r = 0, (x, y, z), dp; poly p = z2 + 1; poly q = z3 + 2; ideal i = p * q^2, y - z2; list pr = minAss(i); pr; } /////////////////////////////////////////////////////////////////////////////// // // Computes the minimal associated primes of I via Laplagne algorithm, // using primary decomposition in the zero dimensional case. // For reduction to the zerodimensional case, it uses the procedure // decomp, with some modifications to avoid the recursion. // static proc minAssSL(ideal I) // Input = I, ideal // Output = primaryDec where primaryDec is the list of the minimal // associated primes and the primary components corresponding to these primes. { ideal P = 1; list pd = list(); int k; int stop = 0; list primaryDec = list(); while (stop == 0) { // Debug dbprint(printlevel - voice, "// We call minAssSLIteration to find new prime ideals!"); pd = minAssSLIteration(I, P); // Debug dbprint(printlevel - voice, "// Output of minAssSLIteration:"); dbprint(printlevel - voice, pd); if (size(pd[1]) > 0) { primaryDec = primaryDec + pd[1]; // Debug dbprint(printlevel - voice, "// We intersect the prime ideals obtained."); P = intersect(P, pd[2]); // Debug dbprint(printlevel - voice, "// Intersection finished."); } else { stop = 1; } } // Returns only the primary components, not the radical. return(primaryDec); } /////////////////////////////////////////////////////////////////////////////// // Given an ideal I and an ideal P (intersection of some minimal prime ideals // associated to I), it calculates new minimal prime ideals associated to I // which were not used to calculate P. // This version uses Primary Decomposition in the zerodimensional case. static proc minAssSLIteration(ideal I, ideal P); { int k = 1; int good = 0; list primaryDec = list(); // Debug dbprint (printlevel-voice, "// We search for an element in P - sqrt(I)."); while ((k <= size(P)) and (good == 0)) { good = 1 - rad_con(P[k], I); k++; } k--; if (good == 0) { // Debug dbprint (printlevel - voice, "// No element was found, P = sqrt(I)."); return (list(primaryDec, ideal(0))); } // Debug dbprint (printlevel - voice, "// We found h = ", P[k]); dbprint (printlevel - voice, "// We calculate the saturation of I with respect to the element just founded."); ideal J = sat(I, P[k])[1]; // Uses decomp from primdec, modified to avoid the recursion. // Debug dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via decomp."); primaryDec = newDecompStep(J, "oneIndep", "intersect", 2); // Debug dbprint(printlevel - voice, "// Proc decomp has found", size(primaryDec) div 2, "new primary components."); return(primaryDec); } /////////////////////////////////////////////////////////////////////////////////// // Based on maxIndependSet // Added list # as parameter // If the first element of # is 0, the output is only 1 max indep set. // If no list is specified or #[1] = 1, the output is all the max indep set of the // leading terms ideal. This is the original output of maxIndependSet proc newMaxIndependSetLp(ideal j, list #) "USAGE: newMaxIndependentSetLp(i); i ideal (returns all maximal independent sets of the corresponding leading terms ideal) newMaxIndependentSetLp(i, 0); i ideal (returns only one maximal independent set) RETURN: list = #1. new varstring with the maximal independent set at the end, #2. ordstring with the lp ordering, #3. the number of independent variables NOTE: EXAMPLE: example newMaxIndependentSetLp; shows an example " { int n, k, di; list resu, hilf; string var1, var2; list v = indepSet(j, 0); // SL 2006.04.21 1 Lines modified to use only one independent Set string indepOption; if (size(#) > 0) { indepOption = #[1]; } else { indepOption = "allIndep"; } int nMax; if (indepOption == "allIndep") { nMax = size(v); } else { nMax = 1; } for(n = 1; n <= nMax; n++) // SL 2006.04.21 2 { di = 0; var1 = ""; var2 = ""; for(k = 1; k <= size(v[n]); k++) { if(v[n][k] != 0) { di++; var2 = var2 + "var(" + string(k) + "), "; } else { var1 = var1 + "var(" + string(k) + "), "; } } if(di > 0) { var1 = var1 + var2; var1 = var1[1..size(var1) - 2]; // The "- 2" removes the trailer comma hilf[1] = var1; // SL 2006.21.04 1 The order is now block dp instead of lp //hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")"; // SL 2006.21.04 2 // For decomp, lp ordering is needed. Nothing is changed. hilf[2] = "lp"; hilf[3] = di; resu[n] = hilf; } else { resu[n] = varstr(basering), ordstr(basering), 0; } } return(resu); } example { "EXAMPLE:"; echo = 2; ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp; ideal i = ea - fbg, fa + be, ec - fdg, fc + de; i = std(i); list l = newMaxIndependSetLp(i); l; i = i, g; l = newMaxIndependSetLp(i); l; ring s = 0, (x, y, z), lp; ideal i = z, yx; list l = newMaxIndependSetLp(i); l; } /////////////////////////////////////////////////////////////////////////////// proc newZero_decomp (ideal j, ideal ser, int @wr, list #) "USAGE: newZero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1 (@wr=0 for primary decomposition, @wr=1 for computation of associated primes) if #[1] = "nest", then #[2] indicates the nest level (number of recursive calls) When the nest level is high it indicates that the computation is difficult, and different methods are applied. RETURN: list = list of primary ideals and their radicals (at even positions in the list) if the input is zero-dimensional and a standardbases with respect to lex-ordering If ser!=(0) and ser is contained in j or if j is not zero-dimen- sional then ideal(1),ideal(1) is returned NOTE: Algorithm of Gianni/Trager/Zacharias EXAMPLE: example newZero_decomp; shows an example " { def @P = basering; int uytrewq; int nva = nvars(basering); int @k,@s,@n,@k1,zz; list primary,lres0,lres1,act,@lh,@wh; map phi,psi,phi1,psi1; ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1; intvec @vh,@hilb; string @ri; poly @f; // Debug dbprint(printlevel - voice, "proc newZero_decomp"); if (dim(j)>0) { primary[1]=ideal(1); primary[2]=ideal(1); return(primary); } j=interred(j); attrib(j,"isSB",1); int nestLevel = 0; if (size(#) > 0) { if (typeof(#[1]) == "string") { if (#[1] == "nest") { nestLevel = #[2]; } # = list(); } } if(vdim(j)==deg(j[1])) { act=factor(j[1]); for(@k=1;@k<=size(act[1]);@k++) { @qh=j; if(@wr==0) { @qh[1]=act[1][@k]^act[2][@k]; } else { @qh[1]=act[1][@k]; } primary[2*@k-1]=interred(@qh); @qh=j; @qh[1]=act[1][@k]; primary[2*@k]=interred(@qh); attrib( primary[2*@k-1],"isSB",1); if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],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(@wr!=0) { primary[1]=std(j,act[1][1]); } if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1]))) { primary[2]=primary[1]; } else { primary[2]=primaryTest(primary[1],act[1][1]); } } if(size(#)==0) { primary=splitPrimary(primary,ser,@wr,act); } if((voice>=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 1){primary=extF(primary);} //test whether all ideals in the decomposition are primary and //in general position //if not after a random coordinate transformation of the last //variable the corresponding ideal is decomposed again. if((npars(basering)>0)&&(nestLevel > 1)) { poly randp; for(zz=1;zz0)&&(nestLevel > 1)) { jmap[size(jmap)]=randp; } for(@n=2;@n<=size(primary[2*@k-1]);@n++) { if(deg(lead(primary[2*@k-1][@n]))==1) { for(zz=1;zz<=nva;zz++) { if(lead(primary[2*@k-1][@n])/var(zz)!=0) { jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n] +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]); jmap2[zz]=primary[2*@k-1][@n]; @qht[@n]=var(zz); } } jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0); } } if(size(subst(jmap[nva],var(1),0)-var(nva))!=0) { // jmap[nva]=subst(jmap[nva],var(1),0); //hier geaendert +untersuchen!!!!!!!!!!!!!! } phi1=@P,jmap1; phi=@P,jmap; for(@n=1;@n<=nva;@n++) { jmap[@n]=-(jmap[@n]-2*var(@n)); } psi=@P,jmap; psi1=@P,jmap2; @qh=phi(@qht); //=================== the new part ============================ if (npars(basering)>1) { @qh=groebner(@qh,"par2var"); } else { @qh=groebner(@qh); } //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // ideal @qh=homog(imap(@P,@qht),@t); // // ideal @qh1=std(@qh); // @hilb=hilb(@qh1,1); // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // ideal @qh=homog(imap(@P,@qh),@t); // kill @Phelp; // @qh=std(@qh,@hilb); // @qh=subst(@qh,@t,1); // setring @P; // @qh=imap(@Phelp1,@qh); // kill @Phelp1; // @qh=clearSB(@qh); // attrib(@qh,"isSB",1); //============================================================= ser1=phi1(ser); @lh=newZero_decomp (@qh,phi(ser1),@wr, list("nest", nestLevel + 1)); kill lres0; list lres0; if(size(@lh)==2) { helpprim=@lh[2]; lres0[1]=primary[2*@k-1]; 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) { if (npars(basering)>1) { lres0[2]=groebner(lres0[2],"par2var"); } else { lres0[2]=groebner(lres0[2]); } } else { for(@n=1;@n<=size(lres0) div 2;@n++) { if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1) { lres0[2*@n-1]=groebner(lres0[2*@n-1]); lres0[2*@n]=lres0[2*@n-1]; attrib(lres0[2*@n],"isSB",1); } else { lres0[2*@n-1]=groebner(lres0[2*@n-1]); lres0[2*@n]=groebner(lres0[2*@n]); } } } primary=primary+lres0; //============================================================= // if(npars(@P)>0) // { // @ri= "ring @Phelp =" // +string(char(@P))+", // ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; // } // else // { // @ri= "ring @Phelp =" // +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; // } // execute(@ri); // list @lvec; // list @lr=imap(@P,lres0); // ideal @lr1; // // if(size(@lr)==2) // { // @lr[2]=homog(@lr[2],@t); // @lr1=std(@lr[2]); // @lvec[2]=hilb(@lr1,1); // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lvec[2*@n]=@lvec[2*@n-1]; // } // else // { // @lr[2*@n-1]=homog(@lr[2*@n-1],@t); // @lr1=std(@lr[2*@n-1]); // @lvec[2*@n-1]=hilb(@lr1,1); // @lr[2*@n]=homog(@lr[2*@n],@t); // @lr1=std(@lr[2*@n]); // @lvec[2*@n]=hilb(@lr1,1); // // } // } // } // @ri= "ring @Phelp1 =" // +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; // execute(@ri); // list @lr=imap(@Phelp,@lr); // // kill @Phelp; // if(size(@lr)==2) // { // @lr[2]=std(@lr[2],@lvec[2]); // @lr[2]=subst(@lr[2],@t,1); // // } // else // { // for(@n=1;@n<=size(@lr) div 2;@n++) // { // if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1) // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=@lr[2*@n-1]; // attrib(@lr[2*@n],"isSB",1); // } // else // { // @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]); // @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1); // @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]); // @lr[2*@n]=subst(@lr[2*@n],@t,1); // } // } // } // kill @lvec; // setring @P; // lres0=imap(@Phelp1,@lr); // kill @Phelp1; // for(@n=1;@n<=size(lres0);@n++) // { // lres0[@n]=clearSB(lres0[@n]); // attrib(lres0[@n],"isSB",1); // } // // primary[2*@k-1]=lres0[1]; // primary[2*@k]=lres0[2]; // @s=size(primary) div 2; // for(@n=1;@n<=size(lres0) div 2-1;@n++) // { // primary[2*@s+2*@n-1]=lres0[2*@n+1]; // primary[2*@s+2*@n]=lres0[2*@n+2]; // } // @k--; //============================================================= } } return(primary); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; poly p = z2+1; poly q = z4+2; ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4; i=std(i); list pr= newZero_decomp(i,ideal(0),0); pr; } /////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// /* //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); */