/////////////////// Deciding 3-jet ///////////////////////////////// LIB"primdec.lib"; LIB"sing.lib"; proc decompose(ideal I) { int i; I=jet(I,3); ideal J=slocus(I); list LI=primdecGTZ(I); list LJ=primdecGTZ(J); if((size(LI)==3)&&(size(LJ)==1)){return(6);} if((size(LI)==3)&&(size(LJ)==3)){return(1);} if((size(LI)==2)&&(size(LJ)==3)){return(2);} if((size(LI)==2)&&(size(LJ)==1)){return(5);} if((size(LI)==1)&&(size(LJ)==1)){return(7);} if(size(LI)==1&&size(LJ)==2) { for(i=1;i<=size(LJ);i++) { if(vdim(std(LJ[i][1]))==6){return(3);} if(vdim(std(LJ[i][1]))==4){return(4);} } } return(0); //this should not occur } ////////////////////// Transform 3-jet into //////////////// proc linearNF(poly f) { // computes the map phi to obtain jet(phi(f),3) in one of the following normal forms // x3+xyz,x2z+yz2,x3+xz2,x3+y3+xyz,x3+yz2 // it is assumed that jet(f,3) has one of this normal forms if(jet(f,3)==0){return(-1);} //modality too high ideal I=factorize(jet(f,3),1); //the factors of the 3-jet ideal K=reduce(maxideal(1),lead(std(I))); //the variables not in I if(size(K)==3) //3-jet is irreducible { I=jacob(jet(f,3)); int mu=milnor(jet(f,3)); //mu=8 is the case x^3+y^3+z^3+axyz if(mu==-1) { poly a=x3+yz2; intvec h=hilbPoly(I); mu=h[1]; //mu=1 x^3+y^3+xyz mu=2 x3+yz^2 if(mu==1){a=x3+y3+xyz;} I=findMap(a,jet(f,3)); if(I[1]==-77){return(-77);} matrix M=lift(I,maxideal(1)); matrix N[3][1]=var(1),var(2),var(3); ideal J=transpose(M)*N; map phi=basering,J; return(phi(f)); } return(-77); //the case x^3+y^3+z^3+axyz } if((size(I)==3)&&(size(std(I))<3)) //the case x(x^2+z^2) if x^2+z^2 splits { poly a=x3+xz2; I=findMap(a,jet(f,3)); if(I[1]==-77){return(-77);} matrix M=lift(I,maxideal(1)); matrix N[3][1]=var(1),var(2),var(3); ideal J=transpose(M)*N; map phi=basering,J; return(phi(f)); } ideal L=simplify(reduce(maxideal(1),std(K)),2); //the variables in I int n=size(L); int m=size(I); if((n==2)&&(m!=3)){return(-1);} //the case x^2y if((n==1)&&(m==1)){return(-1);} //the case x^3 if((n==1)&&(m==2))//3-jet has one linear factor and one quadratic { int l=size(minAssGTZ(I)); //l=2 x(x^2+yz) l=1 z(x^2+yz), x(x^2+z^2) poly s,t,g,h,h1,h2; int mu; if(l==1) { h=I[1]; if(deg(h)==1){h=I[2];} mu=milnor(h); //mu=1 z(x^2+yz) mu=-1 x(x^2+z^2) } if(leadmonom(I[1])==L[1]){h1=I[2];h2=I[1];} else{h2=I[2];h1=I[1];} s=jet(f,3)/(h1*h2); if((l==2)||(mu==-1)) // the cases x(x^2+yz), x(x^2+z^2) { //K=factorize(h1-s*h2^2,1);// h2 corresponds to x and K to y,z resp. z K=factorize(s*h1-h2^2,1);// h2 corresponds to x and K to y,z resp. z if(mu==-1){s=K[1];K[1]=var(2);K[2]=s;} // K corresponds to z I=simplify(jet(I,1)+K,2); } else { // h2 corresponds to z h=reduce(h1,std(ideal(h2))); // this corresponds to x^2 g=factorize(h,1)[1]; // this corresponds to x h=(h1-h)/h2; // this corresponds to y K=g,h; I=simplify(K+jet(I,1),2); } } if(size(I)!=3){return(-77);} // this should not occur matrix M=lift(I,maxideal(1)); matrix N[3][1]=var(1),var(2),var(3); ideal J=transpose(M)*N; map phi=basering,J; f=phi(f); if(lead(f)==-var(1)*var(2)*var(3)){f=subst(f,var(1),-var(1));} if(lead(f)==-var(1)^3){f=subst(f,var(1),-var(1));} if((lead(f)==var(1)^3)&&(lead(f[2])==-var(1)*var(2)*var(3))) { f=subst(f,var(2),-var(2)); } if(leadmonom(f)==var(1)^2*var(3)) { f=subst(f,var(3),1/leadcoef(f)*var(3)); f=subst(f,var(2),1/leadcoef(f[2])*var(2)); } return(f); } proc findMap(poly p, poly q) { // computes an ideal T such that for the map phi defined by T // phi(p)=q and phi is an automorphism, p and q are supposed to be // homogeneous polynomials in 3 variables and it is assumed that such a map exists. int j,l; def R=basering; int pp=char(R); ring S=(pp,a,b,c,d,e,f,g,h,i),(x,y,z),dp; ideal H; map phi=S,a*x+b*y+c*z,d*x+e*y+f*z,g*x+h*y+i*z; //Ansatz for the map poly p=imap(R,p); poly q=imap(R,q); p=phi(p); matrix M=coef(p-q,xyz); for(j=1;j<=ncols(M);j++){H[j]=M[2,j];} //H the ideal of the conditions for the map ring S1=pp,(a,b,c,d,e,f,g,h,i),dp; ideal H=imap(S,H); ideal K,N,Q; poly cc,dd; list L=minAssGTZ(H); //the associated prime ideals of H matrix E[3][3]=a,b,c,d,e,f,g,h,i; poly de=det(E); // we first test if one of the associated primes to H gives such a map directly for(j=1;j<=size(L);j++) { K=std(L[j]); dd=reduce(de,K); if((dd!=0)&&(vdim(K)==1)){break;} //we have already the map dd=0; } if(vdim(K)!=1) { for(j=1;j<=size(L);j++) { K=std(L[j]); Q=reduce(maxideal(1),K); for(l=1;l<=nvars(basering);l++) { if(deg(Q[l])!=1){Q[l]=0;} else{Q[l]=Q[l]+random(-100,100);} } K=std(K+Q); dd=reduce(de,K); if((dd!=0)&&(vdim(K)==1)){break;} //we have already the map dd=0; } } j=0; // all associated primes of H are not maximal, we will add a suitable condition // for one variable while((dd==0)&&(j1)){return(-77);} //this should never occur setring S; ring S2=pp,(a,b,c,d,e,f,g,h,i,x,y,z),dp; ideal N=imap(S1,N); ideal T=a*x+b*y+c*z,d*x+e*y+f*z,g*x+h*y+i*z; T=reduce(T,std(N)); setring R; ideal T=imap(S2,T); return(T); } //////////////////////////// Blowing-Ups/////////////////////////// proc bLowing1(poly f) { poly g,h,k; int a,b,c; map phi1=basering,x,xy,xz; g=phi1(f)/x3; a=milnor(g); int q=a+4; map phi2=basering,xy,y,yz; h=phi2(f)/y3; b=milnor(h); int r=b+4; map phi3=basering,xz,yz,z; k=phi3(f)/z3; c=milnor(k); int s=c+4; list L=q,r,s; return(L); } ////////////////////////////////////////////////////////////////////////////////// proc bLowing2(poly f) { poly g,h,k; int b,c; map phi2=basering,xy,y,yz; h=phi2(f)/y3; b=milnor(h); int r=b+4; map phi3=basering,xz,yz,z; k=phi3(f)/z3; c=milnor(k); int s=c+4; list L=r,s; return(L); } ///////////////////////////////////////////////////////////////////////////////// proc bLowing3(poly f) { poly g,h,k; int c; map phi3=basering,xz,yz,z; k=phi3(f)/z3; c=milnor(k); int s=c+4; list L=s; return(L); } ///////////////////////////////////////////////////////////////////////////////// proc BlowingQ1(poly f) { poly g,h; map phi=basering,xy,y,yz; g=phi(f)/y3; h=phi(g)/y2; ideal I=h; I=jet(I,2); ideal J=slocus(I); list LI=primdecGTZ(I); list LJ=primdecGTZ(J); if((size(LI)==2)&&(size(LJ)==1)){return(8);} if((size(LI)==1)&&(size(LJ)==1)){return(9);} } /////////////////////////////////////////////////////////////////////////////////// proc BlowingQ2(poly f) { poly g,h,k; map phi=basering,xy,y,yz; g=phi(f)/y3; h=phi(g)/y2; k=phi(h)/y2; return(k); } ///////////////////////////////////////////////////////////////////////////////// proc BlowingQ3(poly f) { poly g,h; map phi=basering,xy,y,yz; g=phi(f)/y3; ideal I=g; I=jet(I,2); ideal J=slocus(I); list LI=primdecGTZ(I); list LJ=primdecGTZ(J); if((size(LI)==1)&&(size(LJ)==1)){return(10);} if((size(LI)==2)&&(size(LJ)==1)){return(11);} } /////////////////////////////////////////////////////////////////////////////////// proc BlowingQ4(poly f) { poly g,h; map phi=basering,xy,y,yz; g=phi(f)/y3; h=phi(g)/y2; return(h); } //////////////////////////// Tripple (q,r,s)/////////////////////////// proc Tripple1(list L) { list M; if(size(L)==3) { if(L[1]<=L[2]&&L[2]<=L[3]) { M[1]=L[3]; M[2]=L[2]; M[3]=L[1]; return(M); } if(L[1]<=L[3]&&L[3]<=L[2]) { M[1]=L[2]; M[2]=L[3]; M[3]=L[1]; return(M); } if(L[2]<=L[1]&&L[1]<=L[3]) { M[1]=L[3]; M[2]=L[1]; M[3]=L[2]; return(M); } if(L[2]<=L[3]&&L[3]<=L[1]) { M[1]=L[1]; M[2]=L[3]; M[3]=L[2]; return(M); } if(L[3]<=L[1]&&L[1]<=L[2]) { M[1]=L[2]; M[2]=L[1]; M[3]=L[3]; return(M); } if(L[3]<=L[2]&&L[2]<=L[1]) { M[1]=L[1]; M[2]=L[2]; M[3]=L[3]; return(M); } } else{return(0);} } ///////////////////////////////////////////////////////////////////////////////////////// proc Pair1(list L) { list M; if(size(L)==2) { if(L[1]<=L[2]) { M[1]=L[1]; M[2]=L[2]; return(M); } if(L[2]<=L[1]) { M[1]=L[2]; M[2]=L[1]; return(M); } } else{return(0);} } /////////////////// Classifier /////////////////////////////////// proc classiFy(poly f) { int p=char(basering); poly g,h,k,l; int mu,q,m,n,u; list L,M; poly g=linearNF(f); ideal I=g; ideal J=f; int t=decompose(I); int u=decompose(J); if(t==1) { L=bLowing1(g); M=Tripple1(L); if(M[1]3"); } } if(t==2) { L=bLowing2(g); M=Pair1(L); if(M[2]3"); } } if(t==3) { L=bLowing3(g); if(L[1]3"); } } if(t==4) { mu=milnor(g); if(10<=mu&&mu<=12) { return("f is unimodal of type Q"+string(mu)+""); } if(mu==14||mu==15||19<=mu) { return("f is bimodal of type Q2,"+string(mu-14)+""); } if(mu==16) { q=BlowingQ1(g); if(q==8) { return("f is bimodal of type Q16"); } if(q==9) { return("f is bimodal of type Q2,2"); } } if(mu==17||mu==18) { h=BlowingQ2(g); m=milnor(h); if(m==0) { return("f is bimodal of type Q"+string(mu)+""); } if(m!=0) { return("f is bimodal of type Q2,"+string(mu-14)+""); } } } if(t==5) { mu=milnor(g); if(mu==11||mu==12) { return("f is unimodal of type S"+string(mu)+""); } if(mu==14) { return("f is bimodal of type S1,0"); } if(mu==15) { q=BlowingQ3(g); if(q==10) { return("f is bimodal of type S1,1"); } if(q==11) { return("f is bimodal of type S^#1,1"); } } if(mu==16) { q=BlowingQ3(g); k=BlowingQ4(g); m=milnor(k); if(q==10&&m==0) { return("f is bimodal of type S16"); } if(q==10&&m!=0) { return("f is bimodal of type S1,2"); } if(q==11) { return("f is bimodal of type S^#1,2"); } } if(mu==17) { q=BlowingQ3(g); k=BlowingQ4(g); m=milnor(k); if(q==10&&m==0) { return("f is bimodal of type S17"); } if(q==10&&m!=0) { return("f is bimodal of type S1,3"); } if(q==11) { return("f is bimodal of type S^#1,3"); } } if(16