////////////////////////////////////////////////////////////////////////////// version="$Id$"; category="Commutative Algebra"; info=" LIBRARY: reszeta.lib topological Zeta-function and some other applications of desingularization AUTHORS: A. Fruehbis-Krueger, anne@mathematik.uni-kl.de, @* G. Pfister, pfister@mathematik.uni-kl.de MAIN PROCEDURES: intersectionDiv(L) computes intersection form and genera of exceptional divisors (isolated singularities of surfaces) spectralNeg(L) computes negative spectral numbers (isolated hypersurface singularity) discrepancy(L) computes discrepancy of given resolution zetaDL(L,d) computes Denef-Loeser zeta function (hypersurface singularity of dimension 2) AUXILLIARY PROCEDURES: collectDiv(L[,iv]) identify exceptional divisors in different charts (embedded and non-embedded case) prepEmbDiv(L[,b]) prepare list of divisors (including components of strict transform, embedded case) abstractR(L) pass from embedded to non-embedded resolution "; LIB "resolve.lib"; LIB "solve.lib"; LIB "normal.lib"; /////////////////////////////////////////////////////////////////////////////// static proc spectral1(poly h,list re, list DL,intvec v, intvec n) "Internal procedure - no help and no example available " { //--- compute one spectral number //--- DL is output of prepEmbDiv int i; intvec w=computeH(h,re,DL); number gw=number(w[1]+v[1])/number(n[1]); for(i=2;i<=size(v);i++) { if(gw>number(w[i]+v[i])/number(n[i])) { gw=number(w[i]+v[i])/number(n[i]); } } return(gw-1); } /////////////////////////////////////////////////////////////////////////////// proc spectralNeg(list re,list #) "USAGE: spectralNeg(L); @* L = list of rings ASSUME: L is output of resolution of singularities RETURN: list of numbers, each a spectral number in (-1,0] EXAMPLE: example spectralNeg; shows an example " { //----------------------------------------------------------------------------- // Initialization and Sanity Checks //----------------------------------------------------------------------------- int i,j,l; number bound; list resu; if(size(#)>0) { //--- undocumented feature: //--- if # is not empty it computes numbers up to this bound, //--- not necessarily spectral numbers bound=number(#[1]); } //--- get list of embedded divisors list DL=prepEmbDiv(re,1); int k=1; ideal I,delI; number g; int m=nvars(basering); //--- prepare the multiplicities of exceptional divisors N and nu intvec v=computeV(re,DL); // nu intvec n=computeN(re,DL); // N //--------------------------------------------------------------------------- // start computation, first case separately, then loop //--------------------------------------------------------------------------- resu[1]=spectral1(1,re,DL,v,n); // first number, corresponding to // volume form itself if(resu[1]>=bound) { //--- exceeds bound ==> not a spectral number resu=delete(resu,1); return(resu); } delI=std(ideal(0)); while(k) { //--- now run through all monomial x volume form, degree by degree j++; k=0; I=maxideal(j); I=reduce(I,delI); for(i=1;i<=size(I);i++) { //--- all monomials in degree j g=spectral1(I[i],re,DL,v,n); if(g not a spectral number k=1; l=1; while(resu[l]K[y_1,...,y_n]/J defined by x_i ---> I_i. //--- Let basering=K[y_1,...,y_n], l=n-dim(basering/J), //--- I=, J= //--- For each subset v in {1,...,n} of l elements and //--- w in {1,...,r} of l elements //--- let K_v,w be the ideal generated by the n-l-minors of the matrix //--- (diff(I_i,y_j)+ //--- \sum_k diff(I_i,y_v[k])*diff(J_w[k],y_j))_{j not in v multiplied with //--- the determinant of (diff(J_w[i],y_v[j])) //--- the sum of all such ideals K_v,w plus J is returned. //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int n=nvars(basering); int i,j,k; intvec u,v,w,x; matrix MI[ncols(I)][n]=jacob(I); matrix N=unitmat(n); matrix L; ideal K=J; if(size(J)==0) { K=minor(MI,n); } //--------------------------------------------------------------------------- // Do calculation as described above. // separately for case size(J)=1 //--------------------------------------------------------------------------- if(size(J)==1) { matrix MJ[ncols(J)][n]=jacob(J); N=concat(N,transpose(MJ)); v=1..n; for(i=1;i<=n;i++) { L=transpose(permcol(N,i,n+1)); if(i==1){w=2..n;} if(i==n){w=1..n-1;} if((i!=1)&&(i!=n)){w=1..i-1,i+1..n;} L=submat(L,v,w); L=MI*L; K=K+minor(L,n-1)*MJ[1,i]; } } if(size(J)>1) { matrix MJ[ncols(J)][n]=jacob(J); matrix SMJ; N=concat(N,transpose(MJ)); ideal Jstd=std(J); int l=n-dim(Jstd); int r=ncols(J); list L1=indexSet(n,l); list L2=indexSet(r,l); for(i=1;i<=size(L1);i++) { for(j=1;j<=size(L2);j++) { for(k=1;k<=size(L1[i]);k++) { if(L1[i][k]){v[size(v)+1]=k;} } v=v[2..size(v)]; for(k=1;k<=size(L2[j]);k++) { if(L2[j][k]){w[size(w)+1]=k;} } w=w[2..size(w)]; SMJ=submat(MJ,w,v); L=N; for(k=1;k<=l;k++) { L=permcol(L,v[k],n+w[k]); } u=1..n; x=1..n; v=sort(v)[1]; for(k=l;k>=1;k--) { if(v[k]) { u=deleteInt(u,v[k],1); } } L=transpose(submat(L,u,x)); L=MI*L; K=K+minor(L,n-l)*det(SMJ); } } } return(K); } /////////////////////////////////////////////////////////////////////////////// static proc computeH(ideal h,list re,list DL) "Internal procedure - no help and no example available " { //--- additional procedure to computeV, allows //--- computation for polynomial x volume form //--- by computing the contribution of the polynomial h //--- Note: DL is output of prepEmbDiv //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- def R=basering; ideal II=h; list iden=DL; def T=re[2][1]; setring T; int i,k; intvec v; v[size(iden)]=0; if(deg(II[1])==0){return(v);} //---------------------------------------------------------------------------- // Run through all exceptional divisors //---------------------------------------------------------------------------- for(k=1;k<=size(iden);k++) { for(i=1;i<=size(iden[k]);i++) { if(defined(S)){kill S;} def S=re[2][iden[k][i][1]]; setring S; if((!v[k])&&(defined(EList))) { if(defined(JJ)){kill JJ;} if(defined(phi)){kill phi;} map phi=T,BO[5]; ideal JJ=phi(II); if(size(JJ)!=0) { v[k]=ordE(JJ,EList[iden[k][i][2]],BO[1]); } } setring R; } } return(v); } ////////////////////////////////////////////////////////////////////////////// static proc computeN(list re,list DL) "Internal procedure - no help and no example available " { //--- computes for every (Q-irred.) divisor E_i its multiplicity in f \circ pi //--- DL is output of prepEmbDiv //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- def R=basering; list iden=DL; def T=re[2][1]; setring T; ideal J=BO[2]; int i,k; intvec v; v[size(iden)]=0; //---------------------------------------------------------------------------- // Run through all exceptional divisors //---------------------------------------------------------------------------- for(k=1;k<=size(iden);k++) { for(i=1;i<=size(iden[k]);i++) { if(defined(S)){kill S;} def S=re[2][iden[k][i][1]]; setring S; if((!v[k])&&(defined(EList))) { if(defined(II)){kill II;} if(defined(phi)){kill phi;} map phi=T,BO[5]; ideal II=phi(J); if(size(II)!=0) { v[k]=ordE(II,EList[iden[k][i][2]],BO[1]); } } setring R; } } return(v); } example {"EXAMPLE:"; echo = 2; ring R=0,(x,y,z),dp; ideal I=(x-y)*(x-z)*(y-z)-z4; list re=resolve(I,1); intvec v=computeN(re); v; } ////////////////////////////////////////////////////////////////////////////// static proc countEijk(list re,list iden,intvec iv,list #) "Internal procedure - no help and no example available " { //--- count the number of points in the intersection of 3 exceptional //--- hyperplanes (of dimension 2) - one of them is allowed to be a component //--- of the strict transform //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,j,k,comPa,numPts,localCase; intvec ituple,jtuple,ktuple; list chList,tmpList; def R=basering; if(size(#)>0) { if(string(#[1])=="local") { localCase=1; } } //---------------------------------------------------------------------------- // Find common charts //---------------------------------------------------------------------------- for(i=1;i<=size(iden[iv[1]]);i++) { //--- find common charts - only for final charts if(defined(S)) {kill S;} def S=re[2][iden[iv[1]][i][1]]; setring S; if(!defined(EList)) { i++; setring R; continue; } setring R; kill ituple,jtuple,ktuple; intvec ituple=iden[iv[1]][i]; intvec jtuple=findInIVList(1,ituple[1],iden[iv[2]]); intvec ktuple=findInIVList(1,ituple[1],iden[iv[3]]); if((size(jtuple)!=1)&&(size(ktuple)!=1)) { //--- chList contains all information about the common charts, //--- each entry represents a chart and contains three intvecs from iden //--- one for each E_l kill tmpList; list tmpList=ituple,jtuple,ktuple; chList[size(chList)+1]=tmpList; i++; if(i<=size(iden[iv[1]])) { continue; } else { break; } } } if(size(chList)==0) { //--- no common chart !!! return(int(0)); } //---------------------------------------------------------------------------- // Count points in common charts //---------------------------------------------------------------------------- for(i=1;i<=size(chList);i++) { //--- run through all common charts if(defined(S)) { kill S;} def S=re[2][chList[i][1][1]]; setring S; //--- intersection in this chart if(defined(interId)){kill interId;} if(localCase==1) { //--- in this case we need to intersect with \pi^-1(0) ideal interId=EList[chList[i][1][2]]+EList[chList[i][2][2]] +EList[chList[i][3][2]]+BO[5]; } else { ideal interId=EList[chList[i][1][2]]+EList[chList[i][2][2]] +EList[chList[i][3][2]]; } interId=std(interId); if(defined(otherId)) {kill otherId;} ideal otherId=1; for(j=1;jncols(path))||(comPa>ncols(opath))) break; } comPa=int(leadcoef(path[1,comPa-1])); otherId=otherId+interId; otherId=intersect(otherId, fetchInTree(re,chList[j][1][1], comPa,chList[i][1][1],"interId",iden)); } otherId=std(otherId); //--- do not count each point more than once interId=sat(interId,otherId)[1]; export(interId); if(dim(interId)>0) { ERROR("CountEijk: intersection not zerodimensional"); } //--- add the remaining number of points to the total point count numPts numPts=numPts+vdim(interId); } return(numPts); } ////////////////////////////////////////////////////////////////////////////// static proc chiEij(list re, list iden, intvec iv) "Internal procedure - no help and no example available " { //!!! Copy of chiEij_local adjusted for non-local case //!!! changes must be made in both copies //--- compute the Euler characteristic of the intersection //--- curve of two exceptional hypersurfaces (of dimension 2) //--- one of which is allowed to be a component of the strict transform //--- using the formula chi(Eij)=2-2g(Eij) //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,j,k,chi,g; intvec ituple,jtuple,inters; def R=basering; //---------------------------------------------------------------------------- // Find a common chart in which they intersect //---------------------------------------------------------------------------- for(i=1;i<=size(iden[iv[1]]);i++) { //--- find a common chart in which they intersect: only for final charts if(defined(S)) {kill S;} def S=re[2][iden[iv[1]][i][1]]; setring S; if(!defined(EList)) { i++; setring R; continue; } setring R; kill ituple,jtuple; intvec ituple=iden[iv[1]][i]; intvec jtuple=findInIVList(1,ituple[1],iden[iv[2]]); if(size(jtuple)==1) { if(i1) { ERROR("genus_Eij: higher dimensional intersection"); } if(dim(interId)>=0) { //--- save the index of the current chart for future use export(interId); inters[size(inters)+1]=iden[iv[1]][i][1]; } BO[1]=std(BO[1]); if(((dim(interId)<=0)&&(dim(BO[1])>2))|| ((dim(interId)<0)&&(dim(BO[1])==2))) { if(i1) { ERROR("genus_Eij: higher dimensional intersection"); } if(dim(interId)>=0) { //--- save the index of the current chart for future use export(interId); inters[size(inters)+1]=iden[iv[1]][i][1]; } BO[1]=std(BO[1]); if(dim(interId)<0) { if(i2)) { //--- for sets of points the Euler characteristic is just //--- the number of points //--- fat points are impossible, since everything is smooth and n.c. chi=chi+vdim(interId); points=1; } else { if(points==1) { ERROR("components of intersection do not have same dimension"); } g=genus(interId); //--- chi is the Euler characteristic of the (disjoint !!!) union of the //--- considered components //--- remark: components are disjoint, because the E_i are normal crossing!!! chi=chi+(2-2*g); } } return(chi); } ////////////////////////////////////////////////////////////////////////////// static proc computeChiE(list re, list iden) "Internal procedure - no help and no example available " { //--- compute the Euler characteristic of the exceptional hypersurfaces //--- (of dimension 2), not considering the components of the strict //--- transform //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,j,k,m,thisE,otherE; def R=basering; intvec nulliv,chi_temp,kvec; nulliv[size(iden)]=0; list chi_E; for(i=1;i<=size(iden);i++) { chi_E[i]=list(); } //--------------------------------------------------------------------------- // Run through the list of charts and compute the Euler characteristic of // the new exceptional hypersurface and change the values for the old ones // according to the blow-up which has just been performed // For initialization reasons, treat the case of the first blow-up separately //--------------------------------------------------------------------------- for(i=2;i<=size(re[2]);i++) { //--- run through all charts if(defined(S)){kill S;} def S=re[2][i]; setring S; m=int(leadcoef(path[1,ncols(path)])); if(defined(Spa)){kill Spa;} def Spa=re[2][m]; if(size(BO[4])==1) { //--- just one exceptional divisor thisE=1; setring Spa; if(i==2) { //--- have not set the initial value of chi(E_1) yet if(dim(std(cent))==0) { //--- center was point ==> new except. div. is a P^2 list templist=3*vdim(std(BO[1]+cent)),nulliv; } else { //--- center was curve ==> new except. div. is curve x P^1 list templist=4-4*genus(BO[1]+cent),nulliv; } chi_E[1]=templist; kill templist; } setring S; i++; if(i new except. div. is a P^2 list templist=3*vdim(std(BO[1]+cent)),nulliv; } else { //--- center was curve ==> new except. div. is a C x P^1 list templist=4-4*genus(BO[1]+cent),nulliv; } chi_E[thisE]=templist; kill templist; } for(j=1;j<=size(BO[4]);j++) { //--- we are in the parent ring ==> thisE is not yet born //--- all the other E_i have already been initialized, but the chi //--- might change with the current blow-up at cent if(BO[6][j]==1) { //--- ignore empty sets j++; if(j<=size(BO[4])) { continue; } else { break; } } for(k=1;k<=size(iden);k++) { //--- find global index of BO[4][j] if(inIVList(intvec(m,j),iden[k])) { otherE=k; break; } } if(chi_E[otherE][2][thisE]==1) { //--- already considered this one j++; if(j<=size(BO[4])) { continue; } else { break; } } //--------------------------------------------------------------------------- // update chi according to the formula // chi(E_k^transf)=chi(E_k) - chi(C \cap E_k) + chi(E_k \cap E_new) // for convenience of implementation, we first compute // chi(E_k) - chi(C \cap E_k) // and afterwards add the last term chi(E_k \cap E_new) //--------------------------------------------------------------------------- ideal CinE=std(cent+BO[4][j]+BO[1]); // this is C \cap E_k if(dim(CinE)==1) { //--- center meets E_k in a curve chi_temp[otherE]=chi_E[otherE][1]-(2-2*genus(CinE)); } if(dim(CinE)==0) { //--- center meets E_k in points chi_temp[otherE]=chi_E[otherE][1]-vdim(std(CinE)); } kill CinE; setring S; //--- now we are back in the i-th ring in the list ideal CinE=std(BO[4][j]+BO[4][size(BO[4])]+BO[1]); // this is E_k \cap E_new if(dim(CinE)==1) { //--- if the two divisors meet, they meet in a curve chi_E[otherE][1]=chi_temp[otherE]+(2-2*genus(CinE)); chi_E[otherE][2][thisE]=1; // this blow-up of E_k is done } kill CinE; setring Spa; } } setring R; return(chi_E); } ////////////////////////////////////////////////////////////////////////////// static proc computeChiE_local(list re, list iden) "Internal procedure - no help and no example available " { //--- compute the Euler characteristic of the intersection of the //--- exceptional hypersurfaces with \pi^-1(0) which can be of //--- dimension 1 or 2 - not considering the components of the strict //--- transform //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,j,k,aa,m,n,thisE,otherE; def R=basering; intvec nulliv,chi_temp,kvec,dimEi,endiv; nulliv[size(iden)]=0; dimEi[size(iden)]=0; endiv[size(re[2])]=0; list chi_E; for(i=1;i<=size(iden);i++) { chi_E[i]=list(); } //--------------------------------------------------------------------------- // Run through the list of charts and compute the Euler characteristic of // the new exceptional hypersurface and change the values for the old ones // according to the blow-up which has just been performed // For initialization reasons, treat the case of the first blow-up separately //--------------------------------------------------------------------------- for(i=2;i<=size(re[2]);i++) { //--- run through all charts if(defined(S)){kill S;} def S=re[2][i]; setring S; if(defined(EList)) { endiv[i]=1; } m=int(leadcoef(path[1,ncols(path)])); if(defined(Spa)){kill Spa;} def Spa=re[2][m]; if(size(BO[4])==1) { //--- just one exceptional divisor thisE=1; setring Spa; if(i==2) { //--- have not set the initial value of chi(E_1) yet //--- in the local case, we need to know whether the center contains 0 if(size(reduce(cent,std(maxideal(1))))!=0) { //--- first center does not meet 0 list templist=0,nulliv; dimEi[1]=-1; } else { if(dim(std(cent))==0) { //--- center was point ==> new except. div. is a P^2 list templist=3*vdim(std(BO[1]+cent)),nulliv; dimEi[1]=2; } else { //--- center was curve ==> intersection of new exceptional divisor //--- with \pi^-1(0) is a curve, namely P^1 setring S; list templist=2,nulliv; dimEi[1]=1; } } chi_E[1]=templist; kill templist; } setring S; i++; if(i new except. div. is a P^2 list templist=3*vdim(std(BO[1]+cent)),nulliv; dimEi[thisE]=2; } else { //--- center was curve ==> new except. div. is a C x P^1 if(dim(std(cent+BO[5]))==1) { //--- whole curve is in \pi^-1(0) list templist=4-4*genus(BO[1]+cent),nulliv; dimEi[thisE]=2; } else { //--- curve meets \pi^-1(0) in points //--- in S, the intersection will be a curve!!! setring S; list templist=2-2*genus(BO[1]+BO[4][size(BO[4])]+BO[5]),nulliv; dimEi[thisE]=1; setring Spa; } } } if(defined(templist)) { chi_E[thisE]=templist; kill templist; } } for(j=1;j<=size(BO[4]);j++) { //--- we are in the parent ring ==> thisE is not yet born //--- all the other E_i have already been initialized, but the chi //--- might change with the current blow-up at cent if(BO[6][j]==1) { //--- ignore empty sets j++; if(j<=size(BO[4])) { continue; } else { break; } } for(k=1;k<=size(iden);k++) { //--- find global index of BO[4][j] if(inIVList(intvec(m,j),iden[k])) { otherE=k; break; } } if(dimEi[otherE]<=1) { //--- dimEi[otherE]==-1: center leading to this E does not meet \pi^-1(0) //--- dimEi[otherE]== 0: center leading to this E does not meet \pi^-1(0) //--- in any previously visited charts //--- maybe in some other branch later, but has nothing //--- to do with this center //--- dimEi[otherE]== 1: E \cap \pi^-1(0) is curve //--- ==> chi is birational invariant j++; if(j<=size(BO[4])) { continue; } break; } if(chi_E[otherE][2][thisE]==1) { //--- already considered this one j++; if(j<=size(BO[4])) { continue; } else { break; } } //--------------------------------------------------------------------------- // update chi according to the formula // chi(E_k^transf)=chi(E_k) - chi(C \cap E_k) + chi(E_k \cap E_new) // for convenience of implementation, we first compute // chi(E_k) - chi(C \cap E_k) // and afterwards add the last term chi(E_k \cap E_new) //--------------------------------------------------------------------------- ideal CinE=std(cent+BO[4][j]+BO[1]); // this is C \cap E_k if(dim(CinE)==1) { //--- center meets E_k in a curve chi_temp[otherE]=chi_E[otherE][1]-(2-2*genus(CinE)); } if(dim(CinE)==0) { //--- center meets E_k in points chi_temp[otherE]=chi_E[otherE][1]-vdim(std(CinE)); } kill CinE; setring S; //--- now we are back in the i-th ring in the list ideal CinE=std(BO[4][j]+BO[4][size(BO[4])]+BO[1]); // this is E_k \cap E_new if(dim(CinE)==1) { //--- if the two divisors meet, they meet in a curve chi_E[otherE][1]=chi_temp[otherE][1]+(2-2*genus(CinE)); chi_E[otherE][2][thisE]=1; // this blow-up of E_k is done } kill CinE; setring Spa; } } //--- we still need to clean-up the 1-dimensional E_i \cap \pi^-1(0) for(i=1;i<=size(dimEi);i++) { if(dimEi[i]!=1) { //--- not 1-dimensional ==> skip i++; if(i>size(dimEi)) break; continue; } if(defined(myCharts)) {kill myCharts;} intvec myCharts; chi_E[i]=0; for(j=1;j<=size(re[2]);j++) { if(endiv[j]==0) { //--- not an endChart ==> skip j++; if(j>size(re[2])) break; continue; } if(defined(mtuple)) {kill mtuple;} intvec mtuple=findInIVList(1,j,iden[i]); if(size(mtuple)==1) { //-- nothing to do with this Ei ==> skip j++; if(j>size(re[2])) break; continue; } myCharts[size(myCharts)+1]=j; if(defined(S)){kill S;} def S=re[2][j]; setring S; if(defined(interId)){kill interId;} //--- all components ideal interId=std(BO[4][mtuple[2]]+BO[5]); if(defined(myPts)){kill myPts;} ideal myPts=1; export(myPts); export(interId); if(defined(doneId)){kill doneId;} if(defined(donePts)){kill donePts;} ideal donePts=1; ideal doneId=1; for(k=2;k0) { chi_E[i][1]=chi_E[i][1]+(2-2*genus(pr[k])); } } myPts=sat(myPts,donePts)[1]; chi_E[i][1]=chi_E[i][1]-vdim(myPts); } } setring R; return(chi_E); } ////////////////////////////////////////////////////////////////////////////// static proc chi_ast(list re,list iden,list #) "Internal procedure - no help and no example available " { //--- compute the Euler characteristic of the Ei,Eij,Eijk and the //--- corresponding Ei^*,Eij^*,Eijk^* by preparing the input to the //--- specialized auxilliary procedures and then recombining the results //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,j,k,g; intvec tiv; list chi_ijk,chi_ij,chi_i,ast_ijk,ast_ij,ast_i,tmplist,g_ij,emptylist; list leererSchnitt; def R=basering; ring Rhelp=0,@t,dp; setring R; //---------------------------------------------------------------------------- // first compute the chi(Eij) and at the same time // check whether E_i \cap E_j is empty // the formula is // chi_ij=2-2*genus(E_i \cap E_j) //---------------------------------------------------------------------------- if(size(#)>0) { "Entering chi_ast"; } for(i=1;i<=size(iden)-1;i++) { for(j=i+1;j<=size(iden);j++) { if(defined(blub)){kill blub;} def blub=chiEij(re,iden,intvec(i,j)); if(typeof(blub)=="int") { tmplist=intvec(i,j),blub; } else { leererSchnitt[size(leererSchnitt)+1]=intvec(i,j); tmplist=intvec(i,j),0; } chi_ij[size(chi_ij)+1]=tmplist; } } if(size(#)>0) { "chi_ij computed"; } //----------------------------------------------------------------------------- // compute chi(Eijk)=chi^*(Eijk) by counting the points in the intersection // chi_ijk=#(E_i \cap E_j \cap E_k) // ast_ijk=chi_ijk //----------------------------------------------------------------------------- for(i=1;i<=size(iden)-2;i++) { for(j=i+1;j<=size(iden)-1;j++) { for(k=j+1;k<=size(iden);k++) { if(inIVList(intvec(i,j),leererSchnitt)) { tmplist=intvec(i,j,k),0; } else { tmplist=intvec(i,j,k),countEijk(re,iden,intvec(i,j,k)); } chi_ijk[size(chi_ijk)+1]=tmplist; } } } ast_ijk=chi_ijk; if(size(#)>0) { "chi_ijk computed"; } //---------------------------------------------------------------------------- // construct chi(Eij^*) by the formula // ast_ij=chi_ij - sum_ijk chi_ijk, // where k runs over all indices != i,j //---------------------------------------------------------------------------- for(i=1;i<=size(chi_ij);i++) { ast_ij[i]=chi_ij[i]; for(k=1;k<=size(chi_ijk);k++) { if(((chi_ijk[k][1][1]==chi_ij[i][1][1])|| (chi_ijk[k][1][2]==chi_ij[i][1][1]))&& ((chi_ijk[k][1][2]==chi_ij[i][1][2])|| (chi_ijk[k][1][3]==chi_ij[i][1][2]))) { ast_ij[i][2]=ast_ij[i][2]-chi_ijk[k][2]; } } } if(size(#)>0) { "ast_ij computed"; } //---------------------------------------------------------------------------- // construct ast_i according to the following formulae // ast_i=0 if E_i is (Q- resp. C-)component of the strict transform // chi_i=3*n if E_i originates from blowing up a Q-point, // which consists of n (different) C-points // chi_i=2-2g(C) if E_i originates from blowing up a (Q-)curve C // (chi_i=n*(2-2g(C_i))=2-2g(C), // where C=\cup C_i, C_i \cap C_j = \emptyset) // if E_i is not a component of the strict transform, then // ast_i=chi_i - sum_{j!=i} ast_ij //---------------------------------------------------------------------------- for(i=1;i<=size(iden);i++) { if(defined(S)) {kill S;} def S=re[2][iden[i][1][1]]; setring S; if(iden[i][1][2]>size(BO[4])) { i--; break; } } list idenE=iden; while(size(idenE)>i) { idenE=delete(idenE,size(idenE)); } list cl=computeChiE(re,idenE); for(i=1;i<=size(idenE);i++) { chi_i[i]=list(intvec(i),cl[i][1]); } if(size(#)>0) { "chi_i computed"; } for(i=1;i<=size(idenE);i++) { ast_i[i]=chi_i[i]; for(j=1;j<=size(ast_ij);j++) { if((ast_ij[j][1][1]==i)||(ast_ij[j][1][2]==i)) { ast_i[i][2]=ast_i[i][2]-chi_ij[j][2]; } } for(j=1;j<=size(ast_ijk);j++) { if((ast_ijk[j][1][1]==i)||(ast_ijk[j][1][2]==i) ||(ast_ijk[j][1][3]==i)) { ast_i[i][2]=ast_i[i][2]+chi_ijk[j][2]; } } } for(i=size(idenE)+1;i<=size(iden);i++) { ast_i[i]=list(intvec(i),0); } //--- results are in ast_i, ast_ij and ast_ijk //--- all are of the form intvec(indices),int(value) list result=ast_i,ast_ij,ast_ijk; return(result); } ////////////////////////////////////////////////////////////////////////////// static proc chi_ast_local(list re,list iden,list #) "Internal procedure - no help and no example available " { //--- compute the Euler characteristic of the Ei,Eij,Eijk and the //--- corresponding Ei^*,Eij^*,Eijk^* by preparing the input to the //--- specialized auxilliary procedures and then recombining the results //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,j,k,g; intvec tiv; list chi_ijk,chi_ij,chi_i,ast_ijk,ast_ij,ast_i,tmplist,g_ij,emptylist; list leererSchnitt; def R=basering; ring Rhelp=0,@t,dp; setring R; //---------------------------------------------------------------------------- // first compute // if E_i \cap E_j \cap \pi^-1(0) is a curve: // chi(Eij) and at the same time // check whether E_i \cap E_j is empty // the formula is // chi_ij=2-2*genus(E_i \cap E_j) // otherwise (points): // chi(E_ij) by counting the points //---------------------------------------------------------------------------- if(size(#)>0) { "Entering chi_ast_local"; } for(i=1;i<=size(iden)-1;i++) { for(j=i+1;j<=size(iden);j++) { if(defined(blub)){kill blub;} def blub=chiEij_local(re,iden,intvec(i,j)); if(typeof(blub)=="int") { tmplist=intvec(i,j),blub; } else { leererSchnitt[size(leererSchnitt)+1]=intvec(i,j); tmplist=intvec(i,j),0; } chi_ij[size(chi_ij)+1]=tmplist; } } if(size(#)>0) { "chi_ij computed"; } //----------------------------------------------------------------------------- // compute chi(Eijk)=chi^*(Eijk) by counting the points in the intersection // chi_ijk=#(E_i \cap E_j \cap E_k \cap \pi^-1(0)) // ast_ijk=chi_ijk //----------------------------------------------------------------------------- for(i=1;i<=size(iden)-2;i++) { for(j=i+1;j<=size(iden)-1;j++) { for(k=j+1;k<=size(iden);k++) { if(inIVList(intvec(i,j),leererSchnitt)) { tmplist=intvec(i,j,k),0; } else { tmplist=intvec(i,j,k),countEijk(re,iden,intvec(i,j,k),"local"); } chi_ijk[size(chi_ijk)+1]=tmplist; } } } ast_ijk=chi_ijk; if(size(#)>0) { "chi_ijk computed"; } //---------------------------------------------------------------------------- // construct chi(Eij^*) by the formula // ast_ij=chi_ij - sum_ijk chi_ijk, // where k runs over all indices != i,j //---------------------------------------------------------------------------- for(i=1;i<=size(chi_ij);i++) { ast_ij[i]=chi_ij[i]; for(k=1;k<=size(chi_ijk);k++) { if(((chi_ijk[k][1][1]==chi_ij[i][1][1])|| (chi_ijk[k][1][2]==chi_ij[i][1][1]))&& ((chi_ijk[k][1][2]==chi_ij[i][1][2])|| (chi_ijk[k][1][3]==chi_ij[i][1][2]))) { ast_ij[i][2]=ast_ij[i][2]-chi_ijk[k][2]; } } } if(size(#)>0) { "ast_ij computed"; } //---------------------------------------------------------------------------- // construct ast_i according to the following formulae // ast_i=0 if E_i is (Q- resp. C-)component of the strict transform // if E_i \cap \pi^-1(0) is of dimension 2: // chi_i=3*n if E_i originates from blowing up a Q-point, // which consists of n (different) C-points // chi_i=2-2g(C) if E_i originates from blowing up a (Q-)curve C // (chi_i=n*(2-2g(C_i))=2-2g(C), // where C=\cup C_i, C_i \cap C_j = \emptyset) // if E_i \cap \pi^-1(0) is a curve: // use the formula chi_i=2-2*genus(E_i \cap \pi^-1(0)) // // for E_i not a component of the strict transform we have // ast_i=chi_i - sum_{j!=i} ast_ij //---------------------------------------------------------------------------- for(i=1;i<=size(iden);i++) { if(defined(S)) {kill S;} def S=re[2][iden[i][1][1]]; setring S; if(iden[i][1][2]>size(BO[4])) { i--; break; } } list idenE=iden; while(size(idenE)>i) { idenE=delete(idenE,size(idenE)); } list cl=computeChiE_local(re,idenE); for(i=1;i<=size(cl);i++) { if(size(cl[i])==0) { cl[i][1]=0; } } for(i=1;i<=size(idenE);i++) { chi_i[i]=list(intvec(i),cl[i][1]); } if(size(#)>0) { "chi_i computed"; } for(i=1;i<=size(idenE);i++) { ast_i[i]=chi_i[i]; for(j=1;j<=size(ast_ij);j++) { if((ast_ij[j][1][1]==i)||(ast_ij[j][1][2]==i)) { ast_i[i][2]=ast_i[i][2]-chi_ij[j][2]; } } for(j=1;j<=size(ast_ijk);j++) { if((ast_ijk[j][1][1]==i)||(ast_ijk[j][1][2]==i) ||(ast_ijk[j][1][3]==i)) { ast_i[i][2]=ast_i[i][2]+chi_ijk[j][2]; } } } for(i=size(idenE)+1;i<=size(iden);i++) { ast_i[i]=list(intvec(i),0); } //--- results are in ast_i, ast_ij and ast_ijk //--- all are of the form intvec(indices),int(value) //"End of chi_ast_local"; //~; list result=ast_i,ast_ij,ast_ijk; return(result); } ////////////////////////////////////////////////////////////////////////////// proc discrepancy(list re) "USAGE: discrepancy(L); @* L = list of rings ASSUME: L is the output of resolution of singularities RETRUN: discrepancies of the given resolution" { //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- def R=basering; int i,j; list iden=prepEmbDiv(re); //--- identify the E_i intvec Vvec=computeV(re,iden); //--- nu intvec Nvec=computeN(re,iden); //--- N intvec Avec; //--- only look at exceptional divisors, not at strict transform for(i=1;i<=size(iden);i++) { if(defined(S)) {kill S;} def S=re[2][iden[i][1][1]]; setring S; if(iden[i][1][2]>size(BO[4])) { i--; break; } } j=i; //--- discrepancies are a_i=nu_i-N_i for(i=1;i<=j;i++) { Avec[i]=Vvec[i]-Nvec[i]-1; } return(Avec); } example {"EXAMPLE:"; echo = 2; ring R=0,(x,y,z),dp; ideal I=x2+y2+z3; list re=resolve(I); discrepancy(re); } ////////////////////////////////////////////////////////////////////////////// proc zetaDL(list re,int d,list #) "USAGE: zetaDL(L,d[,s1][,s2][,a]); L = list of rings; d = integer; s1,s2 = string; a = integer ASSUME: L is the output of resolution of singularities COMPUTE: local Denef-Loeser zeta function, if string s1 is present and has the value 'local'; global Denef-Loeser zeta function otherwise if string s1 or s2 has the value "A", additionally the characteristic polynomial of the monodromy is computed RETURN: list l if a is not present: l[1]: string specifying the top. zeta function l[2]: string specifying characteristic polynomial of monodromy, if "A" was specified if a is present: l[1]: string specifying the top. zeta function l[2]: list ast, ast[1]=chi(Ei^*) ast[2]=chi(Eij^*) ast[3]=chi(Eijk^*) l[3]: intvec nu of multiplicites as needed in computation of zeta function l[4]: intvec N of multiplicities as needed in compuation of zeta function l[5]: string specifying characteristic polynomial of monodromy, if "A" was specified EXAMPLE: example zetaDL; shows an example " { //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- def R=basering; int show_all,i; if(size(#)>0) { if((typeof(#[1])=="int")||(size(#)>2)) { show_all=1; } if(typeof(#[1])=="string") { if((#[1]=="local")||(#[1]=="lokal")) { // ERROR("Local case not implemented yet"); "Local Case: Assuming that no (!) charts were dropped"; "during calculation of the resolution (option \"A\")"; int localComp=1; if(size(#)>1) { if(#[2]=="A") { int aCampoFormula=1; } } } else { if(#[1]=="A") { int aCampoFormula=1; } "Computing global zeta function"; } } } //---------------------------------------------------------------------------- // Identify the embedded divisors and chi(Ei^*), chi(Eij^*) and chi(Eijk^*) // as well as the integer vector V(=nu) and N //---------------------------------------------------------------------------- list iden=prepEmbDiv(re); //--- identify the E_i //!!! TIMING: E8 takes 520 sec ==> needs speed up if(!defined(localComp)) { list ast_list=chi_ast(re,iden); //--- compute chi(E^*) } else { list ast_list=chi_ast_local(re,iden); } intvec Vvec=computeV(re,iden); //--- nu intvec Nvec=computeN(re,iden); //--- N //---------------------------------------------------------------------------- // Build a new ring with one parameter s // and compute Zeta_top^(d) in its ground field //---------------------------------------------------------------------------- ring Qs=(0,s),x,dp; number zetaTop=0; number enum,denom; denom=1; for(i=1;i<=size(Nvec);i++) { denom=denom*(Vvec[i]+s*Nvec[i]); } //--- factors for which index set J consists of one element //--- (do something only if d divides N_j) for(i=1;i<=size(ast_list[1]);i++) { if((((Nvec[ast_list[1][i][1][1]]/d)*d)-Nvec[ast_list[1][i][1][1]]==0)&& (ast_list[1][i][2]!=0)) { enum=enum+ast_list[1][i][2]*(denom/(Vvec[ast_list[1][i][1][1]]+s*Nvec[ast_list[1][i][1][1]])); } } //--- factors for which index set J consists of two elements //--- (do something only if d divides both N_i and N_j) //!!! TIMING: E8 takes 690 sec and has 703 elements //!!! ==> need to implement a smarter way to do this //!!! e.g. build up enumerator and denominator separately, thus not //!!! searching for common factors in each step for(i=1;i<=size(ast_list[2]);i++) { if((((Nvec[ast_list[2][i][1][1]]/d)*d)-Nvec[ast_list[2][i][1][1]]==0)&& (((Nvec[ast_list[2][i][1][2]]/d)*d)-Nvec[ast_list[2][i][1][2]]==0)&& (ast_list[2][i][2]!=0)) { enum=enum+ast_list[2][i][2]*(denom/((Vvec[ast_list[2][i][1][1]]+s*Nvec[ast_list[2][i][1][1]])*(Vvec[ast_list[2][i][1][2]]+s*Nvec[ast_list[2][i][1][2]]))); } } //--- factors for which index set J consists of three elements //--- (do something only if d divides N_i, N_j and N_k) //!!! TIMING: E8 takes 490 sec and has 8436 elements //!!! ==> same kind of improvements as in the previous case needed for(i=1;i<=size(ast_list[3]);i++) { if((((Nvec[ast_list[3][i][1][1]]/d)*d)-Nvec[ast_list[3][i][1][1]]==0)&& (((Nvec[ast_list[3][i][1][2]]/d)*d)-Nvec[ast_list[3][i][1][2]]==0)&& (((Nvec[ast_list[3][i][1][3]]/d)*d)-Nvec[ast_list[3][i][1][3]]==0)&& (ast_list[3][i][2]!=0)) { enum=enum+ast_list[3][i][2]*(denom/((Vvec[ast_list[3][i][1][1]]+s*Nvec[ast_list[3][i][1][1]])*(Vvec[ast_list[3][i][1][2]]+s*Nvec[ast_list[3][i][1][2]])*(Vvec[ast_list[3][i][1][3]]+s*Nvec[ast_list[3][i][1][3]]))); } } zetaTop=enum/denom; zetaTop=numerator(zetaTop)/denominator(zetaTop); string zetaStr=string(zetaTop); if(show_all) { list result=zetaStr,ast_list[1],ast_list[2],ast_list[3],Vvec,Nvec; } else { list result=zetaStr; } //--- compute characteristic polynomial of the monodromy //--- by the A'Campo formula if(defined(aCampoFormula)) { poly charP=1; for(i=1;i<=size(ast_list[1]);i++) { charP=charP*((s^Nvec[i]-1)^ast_list[1][i][2]); } string charPStr=string(charP/(s-1)); result[size(result)+1]=charPStr; } setring R; return(result); } example {"EXAMPLE:"; echo = 2; ring R=0,(x,y,z),dp; ideal I=x2+y2+z3; list re=resolve(I,"K"); zetaDL(re,1); I=(xz+y2)*(xz+y2+x2)+z5; list L=resolve(I,"K"); zetaDL(L,1); //===== expected zeta function ========= // (20s^2+130s+87)/((1+s)*(3+4s)*(29+40s)) //====================================== } ////////////////////////////////////////////////////////////////////////////// proc abstractR(list re) "USAGE: abstractR(L); @* L = list of rings ASSUME: L is output of resolution of singularities NOTE: currently only implemented for isolated surface singularities RETURN: list l l[1]: intvec, where l[1][i]=1 if the corresponding ring is a final chart of non-embedded resolution l[1][i]=0 otherwise l[2]: intvec, where l[2][i]=1 if the corresponding ring does not occur in the non-embedded resolution l[2][i]=0 otherwise l[3]: list L EXAMPLE: example abstractR; shows an example " { //--------------------------------------------------------------------------- // Initialization and sanity checks //--------------------------------------------------------------------------- def R=basering; //---Test whether we are in the irreducible surface case def S=re[2][1]; setring S; BO[2]=BO[2]+BO[1]; if(dim(std(BO[2]))!=2) { ERROR("NOT A SURFACE"); } if(dim(std(slocus(BO[2])))>0) { ERROR("NOT AN ISOLATED SINGULARITY"); } setring R; int i,j,k,l,i0; intvec deleted; intvec endiv; endiv[size(re[2])]=0; deleted[size(re[2])]=0; //----------------------------------------------------------------------------- // run through all rings, only consider final charts // for each final chart follow the list of charts leading up to it until // we encounter a chart which is not finished in the non-embedded case //----------------------------------------------------------------------------- for(i=1;i<=size(re[2]);i++) { if(defined(S)){kill S;} def S=re[2][i]; setring S; if(size(reduce(cent,std(BO[2]+BO[1])))!=0) { //--- only consider endrings i++; continue; } i0=i; for(j=ncols(path);j>=2;j--) { //--- walk backwards through history if(j==2) { endiv[i0]=1; break; } k=int(leadcoef(path[1,j])); if((deleted[k]==1)||(endiv[k]==1)) { deleted[i0]=1; break; } if(defined(SPa)){kill SPa;} def SPa=re[2][k]; setring SPa; l=int(leadcoef(path[1,ncols(path)])); if(defined(SPa2)){kill SPa2;} def SPa2=re[2][l]; setring SPa2; if((deleted[l]==1)||(endiv[l]==1)) { //--- parent was already treated via different final chart //--- we may safely inherit the data deleted[i0]=1; setring S; i0=k; j--; continue; } setring SPa; //!!! Idea of Improvement: //!!! BESSER: rueckwaerts gehend nur testen ob glatt //!!! danach vorwaerts bis zum ersten Mal abstractNC //!!! ACHTUNG: rueckweg unterwegs notieren - wir haben nur vergangenheit! if((deg(std(slocus(BO[2]))[1])!=0)||(!abstractNC(BO))) { //--- not finished in the non-embedded case endiv[i0]=1; break; } //--- unnecessary chart in non-embedded case setring S; deleted[i0]=1; i0=k; } } //----------------------------------------------------------------------------- // Clean up the intvec deleted and return the result //----------------------------------------------------------------------------- setring R; for(i=1;i<=size(endiv);i++) { if(endiv[i]==1) { if(defined(S)) {kill S;} def S=re[2][i]; setring S; for(j=3;j0) { //--- we are not interested in the strict transform of X return(iden); } //---------------------------------------------------------------------------- // Run through all final charts and collect and identify all components of // the strict transform //---------------------------------------------------------------------------- //--- first final chart - to be used for initialization def S=re[2][iden0[size(iden0)][1][1]]; setring S; ncomps=size(EList)-size(BO[4]); if((ncomps==1)&&(deg(std(EList[size(EList)])[1])==0)) { ncomps=0; } offset=size(BO[4]); for(i=1;i<=ncomps;i++) { //--- add components of strict transform tmpList[1]=intvec(iden0[size(iden0)][1][1],size(BO[4])+i); iden[size(iden)+1]=tmpList; } //--- now run through the other final charts for(i=2;i<=size(iden0[size(iden0)]);i++) { if(defined(S2)){kill S2;} def S2=re[2][iden0[size(iden0)][i][1]]; setring S2; //--- determine common parent of this ring and re[2][iden0[size(iden0)][1][1]] if(defined(opath)){kill opath;} def opath=imap(S,path); j=1; while(opath[1,j]==path[1,j]) { j++; if((j>ncols(path))||(j>ncols(opath))) break; } if(defined(li1)){kill li1;} list li1; //--- fetch the components we have considered in //--- re[2][iden0[size(iden0)][1][1]] //--- via the resolution tree for(k=1;k<=ncomps;k++) { if(defined(id1)){kill id1;} string tempstr="EList["+string(eval(k+offset))+"]"; ideal id1=fetchInTree(re,iden0[size(iden0)][1][1], int(leadcoef(path[1,j-1])), iden0[size(iden0)][i][1],tempstr,iden0,1); kill tempstr; li1[k]=id1; kill id1; } //--- do the comparison for(k=size(BO[4])+1;k<=size(EList);k++) { //--- only components of the strict transform are interesting if((size(BO[4])+1==size(EList))&&(deg(std(EList[size(EList)])[1])==0)) { break; } found=0; for(j=1;j<=size(li1);j++) { if((size(reduce(li1[j],std(EList[k])))==0)&& (size(reduce(EList[k],std(li1[j])))==0)) { //--- found a match li1[j]=ideal(1); iden[size(iden0)-1+j][size(iden[size(iden0)-1+j])+1]= intvec(iden0[size(iden0)][i][1],k); found=1; break; } } if(!found) { //--- no match yet, maybe there are entries not corresponding to the //--- initialization of the list -- collected in list repair if(!defined(repair)) { //--- no entries in repair, we add the very first one list repair; repair[1]=list(intvec(iden0[size(iden0)][i][1],k)); } else { //--- compare against repair, and add the item appropriately //--- steps of comparison as before for(c=1;c<=size(repair);c++) { for(d=1;d<=size(repair[c]);d++) { if(defined(opath)) {kill opath;} def opath=imap(re[2][repair[c][d][1]],path); b=0; while(path[1,b+1]==opath[1,b+1]) { b++; if((b>ncols(path)-1)||(b>ncols(opath)-1)) break; } b=int(leadcoef(path[1,b])); string tempstr="EList["+string(eval(repair[c][d][2])) +"]"; if(defined(id1)){kill id1;} ideal id1=fetchInTree(re,repair[c][d][1],b, iden0[size(iden0)][i][1],tempstr,iden0,1); kill tempstr; if((size(reduce(EList[k],std(id1)))==0)&& (size(reduce(id1,std(EList[k])))==0)) { repair[c][size(repair[c])+1]=intvec(iden0[size(iden0)][i][1],k); break; } } if(d<=size(repair[c])) { break; } } if(c>size(repair)) { repair[size(repair)+1]=list(intvec(iden0[size(iden0)][i][1],k)); } } } } } if(defined(repair)) { //--- there were further components, add them for(c=1;c<=size(repair);c++) { iden[size(iden)+1]=repair[c]; } kill repair; } //--- up to now only Q-irred components - not C-irred components !!! return(iden); } example {"EXAMPLE:"; echo = 2; ring R=0,(x,y,z),dp; ideal I=x2+y2+z11; list L=resolve(I); prepEmbDiv(L); } /////////////////////////////////////////////////////////////////////////////// static proc decompEinX(list BO) "Internal procedure - no help and no example available " { //--- decomposition of exceptional divisor, non-embedded resolution. //--- even a single exceptional divisor may be Q-reducible when considered //--- as divisor on the strict transform //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,j,k,de,contact; intmat interMat; list dcE,tmpList,prList,sa,nullList; string mpol,compList; def R=basering; ideal I; //---------------------------------------------------------------------------- // pass to divisors on V(J) and throw away components already present as // previous exceptional divisors //---------------------------------------------------------------------------- for(i=1;i<=size(BO[4]);i++) { I=BO[4][i]+BO[2]; for(j=i+1;j<=size(BO[4]);j++) { sa=sat(I,BO[4][j]+BO[2]); if(sa[2]) { I=sa[1]; } } //!!! Practical improvement - not yet implemented: //!!!hier den Input besser aufbereiten (cf. J. Wahl's example) //!!!I[1]=x(2)^15*y(2)^9+3*x(2)^10*y(2)^6+3*x(2)^5*y(2)^3+x(2)+1; //!!!I[2]=x(2)^8*y(2)^6+y(0); //!!!heuristisch die Ordnung so waehlen, dass y(0) im Prinzip eliminiert //!!!wird. //----------------------------------------------------------------------------- // 1) decompose exceptional divisor (over Q) // 2) check whether there are C-reducible Q-components // 3) if necessary, find appropriate field extension of Q to decompose // 4) in each chart collect information in list dcE and export it //----------------------------------------------------------------------------- prList=primdecGTZ(I); for(j=1;j<=size(prList);j++) { tmpList=grad(prList[j][2]); de=tmpList[1]; interMat=tmpList[2]; mpol=tmpList[3]; compList=tmpList[4]; nullList=tmpList[5]; contact=Kontakt(prList[j][1],BO[2]); tmpList=prList[j][2],de,contact,interMat,mpol,compList,nullList; prList[j]=tmpList; } dcE[i]=prList; } return(dcE); } ////////////////////////////////////////////////////////////////////////////// static proc getMinpoly(poly p) "Internal procedure - no help and no example available " { //---assume that p is a polynomial in 2 variables and irreducible //---over Q. Computes an irreducible polynomial mp in one variable //---over Q such that p splits completely over the splitting field of mp //---returns mp as a string //---use a variant of the algorithm of S. Gao def R=basering; int i,j,k,a,b,m,n; intvec v; string mp="poly p=t-1;"; list Li=string(1); list re=mp,Li,1; //---check which variables occur in p for(i=1;i<=nvars(basering);i++) { if(p!=subst(p,var(i),0)){v[size(v)+1]=i;} } //---the polynomial is constant if(size(v)==1){return(re);} //---the polynomial depends only on one variable or is homogeneous //---in 2 variables if((size(v)==2)||((size(v)==3)&&(homog(p)))) { if((size(v)==3)&&(homog(p))) { p=subst(p,var(v[3]),1); } ring Rhelp=0,var(v[2]),dp; poly p=imap(R,p); ring Shelp=0,t,dp; poly p=fetch(Rhelp,p); int de=deg(p); p=simplifyMinpoly(p); Li=getNumZeros(p); short=0; mp="poly p="+string(p)+";"; re=mp,Li,de; setring R; return(re); } v=v[2..size(v)]; if(size(v)>2){ERROR("getMinpoly:input depends on more then 2 variables");} //---the general case, the polynomial is considered as polynomial in x an y now ring T=0,(x,y),lp; ideal M,N; M[nvars(R)]=0; N[nvars(R)]=0; M[v[1]]=x; N[v[1]]=y; M[v[2]]=y; N[v[2]]=x; map phi=R,M; map psi=R,N; poly p=phi(p); poly q=psi(p); ring Thelp=(0,x),y,dp; poly p=imap(T,p); poly q=imap(T,q); n=deg(p); //---the degree with respect to y m=deg(q); //---the degree with respect to x setring T; ring A=0,(u(1..m*(n+1)),v(1..(m+1)*n),x,y,t),dp; poly f=imap(T,p); poly g; poly h; for(i=0;i<=m-1;i++) { for(j=0;j<=n;j++) { g=g+u(i*(n+1)+j+1)*x^i*y^j; } } for(i=0;i<=m;i++) { for(j=0;j<=n-1;j++) { h=h+v(i*n+j+1)*x^i*y^j; } } poly L=f*(diff(g,y)-diff(h,x))+h*diff(f,x)-g*diff(f,y); //---according to the theory f is absolutely irreducible if and only if //---L(g,h)=0 has no non-trivial solution g,h //---(g=diff(f,x),h=diff(f,y) is always a solution) //---therefore we compute a vector space basis of G //---G={g in Q[x,y],deg_x(g)) //---according to the theory f=product over all c in C of the //---gcd(f,g-c*diff(f,x)) //---let g_1,...,g_a be a basis of G and write //---g*g_i=sum a_ij*g_j*diff(f,x) mod f //---let B=(a_ij) and ch=det(t*unitmat(a)-B) the characteristic //---polynomial then the number of distinct irreducible factors //---of gcd(f,g-c*diff(f,x)) in C[x,y] is equal to the multiplicity //---of c as a root of ch. //---in our special situation (f is irreducible over Q) ch should //---be irreducible and the different roots of ch lead to the //---factors of f, i.e. ch is the minpoly we are looking for poly fh=homog(f,t); //---homogenization is used to obtain a constant matrix using lift ideal Gh=homog(G,t); int dh,df; df=deg(fh); for(i=1;i<=a;i++) { if(deg(Gh[i])>dh){dh=deg(Gh[i]);} } for(i=1;i<=a;i++) { Gh[i]=t^(dh-deg(Gh[i]))*Gh[i]; } ideal GF=simplify(diff(fh,x),1)*Gh,fh; poly ch; matrix LI; matrix B[a][a]; matrix E=unitmat(a); poly gran; ideal fac; for(i=1;i<=a;i++) { LI=lift(GF,t^(df-1-dh)*Gh[i]*Gh); B=LI[1..a,1..a]; ch=det(t*E-B); //---irreducibility test fac=factorize(ch,1); if(deg(fac[1])==a) { ch=simplifyMinpoly(ch); Li=getNumZeros(ch); int de=deg(ch); short=0; mp="poly p="+string(ch)+";"; re=mp,Li,de; setring R; return(re); } } ERROR("getMinpoly:not found:please send the example to the authors"); } ////////////////////////////////////////////////////////////////////////////// static proc getNumZeros(poly p) "Internal procedure - no help and no example available " { //--- compute numerically (!!!) the zeros of the minimal polynomial def R=basering; ring S=0,t,dp; poly p=imap(R,p); def L=laguerre_solve(p,30); //!!! practical improvement: //!!! testen ob die Nullstellen signifikant verschieden sind //!!! und im Notfall Genauigkeit erhoehen list re; int i; for(i=1;i<=size(L);i++) { re[i]=string(L[i]); } setring R; return(re); } ////////////////////////////////////////////////////////////////////////////// static proc simplifyMinpoly(poly p) "Internal procedure - no help and no example available " { //--- describe field extension in a simple way p=cleardenom(p); int n=int(leadcoef(p)); int d=deg(p); int i,k; int re=1; number s=1; list L=primefactors(n); for(i=1;i<=size(L[1]);i++) { k=L[2][i] mod d; s=1/number((L[1][i])^(L[2][i]/d)); if(!k){p=subst(p,t,s*t);} } p=cleardenom(p); n=int(leadcoef(subst(p,t,0))); L=primefactors(n); for(i=1;i<=size(L[1]);i++) { k=L[2][i] mod d; s=(L[1][i])^(L[2][i]/d); if(!k){p=subst(p,t,s*t);} } p=cleardenom(p); return(p); } /////////////////////////////////////////////////////////////////////////////// static proc grad(ideal I) "Internal procedure - no help and no example available " { //--- computes the number of components over C //--- for a prime ideal of height 1 over Q def R=basering; int n=nvars(basering); string mp="poly p=t-1;"; string str=string(1); list zeroList=string(1); int i,j,k,l,d,e,c,mi; ideal Istd=std(I); intmat interMat; d=dim(Istd); if(d==-1){return(list(0,0,mp,str,zeroList));} if(d!=1){ERROR("ideal is not one-dimensional");} ideal Sloc=std(slocus(I)); if(deg(Sloc[1])>0) { //---This is only to test that in case of singularities we have only //---one singular point which is a normal crossing //---consider the different singular points ideal M; list pr=minAssGTZ(Sloc); if(size(pr)>1){ERROR("grad:more then one singular point");} for(l=1;l<=size(pr);l++) { M=std(pr[l]); d=vdim(M); if(d!=1) { //---now we have to extend the field if(defined(S)){kill S;} ring S=0,x(1..n),lp; ideal M=fetch(R,M); ideal I=fetch(R,I); ideal jmap; map phi=S,maxideal(1);; ideal Mstd=std(M); //---M has to be in general position with respect to lp, i.e. //---vdim(M)=deg(M[1]) poly p=Mstd[1]; e=vdim(Mstd); while(e!=deg(p)) { jmap=randomLast(100); phi=S,jmap; Mstd=std(phi(M)); p=Mstd[1]; } I=phi(I); kill phi; //---now it is in general position an M[1] defines the field extension //---Q[x]/M over Q ring Shelp=0,t,dp; ideal helpmap; helpmap[n]=t; map psi=S,helpmap; poly p=psi(p); ring T=(0,t),x(1..n),lp; poly p=imap(Shelp,p); //---we are now in the polynomial ring over the field Q[x]/M minpoly=leadcoef(p); ideal M=imap(S,Mstd); M=M,var(n)-t; ideal I=fetch(S,I); } //---we construct a map phi which maps M to maxideal(1) option(redSB); ideal Mstd=-simplify(std(M),1); option(noredSB); for(i=1;i<=n;i++) { Mstd=subst(Mstd,var(i),-var(i)); M[n-i+1]=Mstd[i]; } M=M[1..n]; //---go to the localization with respect to if(d!=1) { ring Tloc=(0,t),x(1..n),ds; poly p=imap(Shelp,p); minpoly=leadcoef(p); ideal M=fetch(T,M); map phi=T,M; } else { ring Tloc=0,x(1..n),ds; ideal M=fetch(R,M); map phi=R,M; } ideal I=phi(I); ideal Istd=std(I); mi=mi+milnor(Istd); if(mi>l) { ERROR("grad:divisor is really singular"); } setring R; } } intvec ind=indepSet(Istd,1)[1]; for(i=1;i<=n;i++){if(ind[i]) break;} //---the i-th variable is the independent one ring Shelp=0,x(1..n),dp; ideal I=fetch(R,I); if(defined(S)){kill S;} if(i==1){ring S=(0,x(1)),x(2..n),lp;} if(i==n){ring S=(0,x(n)),x(1..n-1),lp;} if((i!=1)&&(i!=n)){ring S=(0,x(i)),(x(1..i-1),x(i+1..n)),lp;} //---I is zero-dimensional now ideal I=imap(Shelp,I); ideal Istd=std(I); ideal jmap; map phi; poly p=Istd[1]; e=vdim(Istd); if(e==1) { setring R; str=string(I); list resi=1,interMat,mp,str,zeroList; return(resi); } //---move I to general position with respect to lp if(e!=deg(p)) { jmap=randomLast(5); phi=S,jmap; Istd=std(phi(I)); p=Istd[1]; } while(e!=deg(p)) { jmap=randomLast(100); phi=S,jmap; Istd=std(phi(I)); p=Istd[1]; } setring Shelp; poly p=imap(S,p); list Q=getMinpoly(p); int de=Q[3]; mp=Q[1]; //!!!diese Stelle effizienter machen //!!!minAssGTZ vermeiden durch direkte Betrachtung von //!!!p und mp und evtl. Quotientenbildung //!!!bisher nicht zeitkritisch string Tesr="ring Tes=(0,t),("+varstr(R)+"),dp;"; execute(Tesr); execute(mp); minpoly=leadcoef(p); ideal I=fetch(R,I); list pr=minAssGTZ(I); ideal allgEbene=randomLast(100)[nvars(basering)]; int minpts=vdim(std(I+allgEbene)); ideal tempi; j=1; for(i=1;i<=size(pr);i++) { tempi=std(pr[i]+allgEbene); if(vdim(tempi), ai in a field //---extension of Q, containing I+K and an integer a //---such that in the localization of the polynomial ring with //---respect to M the ideal I is not contained in K+M^a+1 but in M^a in def R=basering; int n=nvars(basering); int i,j,k,d,e; ideal J=std(I+K); if(dim(J)==-1){return(0);} ideal W; //---choice of the maximal ideal M for(i=1;i<=n;i++) { W=std(J,var(i)); d=dim(W); if(d==0) break; } i=1;k=2; while((d)&&(i if(d!=1) { ring Tloc=(0,t),x(1..n),ds; poly p=imap(Shelp,p); minpoly=leadcoef(p); ideal M=fetch(T,M); map phi=T,M; } else { ring Tloc=0,x(1..n),ds; ideal M=fetch(R,M); map phi=R,M; } ideal K=phi(K); ideal I=phi(I); //---compute the order of I in (Q[x]/M)[[x]]/K k=1;d=0; while(!d) { k++; d=size(reduce(I,std(maxideal(k)+K))); } setring R; return(k-1); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; ideal I=x4+z4+1; ideal K=x+y2+z2; Kontakt(I,K); } ////////////////////////////////////////////////////////////////////////////// static proc abstractNC(list BO) "Internal procedure - no help and no example available " { //--- check normal crossing property //--- used for passing from embedded to non-embedded resolution //---------------------------------------------------------------------------- // Initialization //---------------------------------------------------------------------------- int i,k,j,flag; list L; ideal J; if(dim(std(cent))>0){return(1);} //---------------------------------------------------------------------------- // check each exceptional divisor on V(J) //---------------------------------------------------------------------------- for(i=1;i<=size(BO[4]);i++) { if(dim(std(BO[2]+BO[4][i]))>0) { //--- really something to do J=radical(BO[4][i]+BO[2]); if(deg(std(slocus(J))[1])!=0) { if(!nodes(J)) { //--- really singular, not only nodes ==> not normal crossing return(0); } } for(k=1;k<=size(L);k++) { //--- run through previously considered divisors //--- we do not want to bother with the same one twice if((size(reduce(J,std(L[k])))==0)&&(size(reduce(L[k],std(J)))==0)) { //--- already considered this one flag=1;break; } //--- drop previously considered exceptional divisors from the current one J=sat(J,L[k])[1]; if(deg(std(J)[1])==0) { //--- nothing remaining flag=1;break; } } if(flag==0) { //--- add exceptional divisor to the list L[size(L)+1]=J; } flag=0; } } //--------------------------------------------------------------------------- // check intersection properties between different exceptional divisors //--------------------------------------------------------------------------- for(k=1;k not normal crossing return(0); } for(j=i+1;j<=size(L);j++) { if(deg(std(L[i]+L[j]+L[k])[1])>0) { //--- three divisors meet simultaneously ==> not normal crossing return(0); } } } } //--- we reached this point ==> normal crossing return(1); } ////////////////////////////////////////////////////////////////////////////// static proc nodes(ideal J) "Internal procedure - no help and no example available " { //--- check whether at most nodes occur as singularities ideal K=std(slocus(J)); if(deg(K[1])==0){return(1);} if(dim(K)>0){return(0);} if(vdim(K)!=vdim(std(radical(K)))){return(0);} return(1); } ////////////////////////////////////////////////////////////////////////////// proc intersectionDiv(list re) "USAGE: intersectionDiv(L); @* L = list of rings ASSUME: L is output of resolution of singularities (only case of isolated surface singularities) COMPUTE: intersection matrix and genera of the exceptional divisors (considered as curves on the strict transform) RETURN: list l, where l[1]: intersection matrix of exceptional divisors l[2]: intvec, genera of exceptional divisors l[3]: divisorList, encoding the identification of the divisors EXAMPLE: example intersectionDiv; shows an example " { //---------------------------------------------------------------------------- //--- Computes in case of surface singularities (non-embedded resolution): //--- the intersection of the divisors (on the surface) //--- assuming that re=resolve(J) //---------------------------------------------------------------------------- def R=basering; //---Test whether we are in the irreducible surface case def S=re[2][1]; setring S; BO[2]=BO[2]+BO[1]; // make sure we are living in the smooth W if(dim(std(BO[2]))!=2) { ERROR("The given original object is not a surface"); } if(dim(std(slocus(BO[2])))>0) { ERROR("The given original object has non-isolated singularities."); } setring R; //---------------------------------------------------------------------------- // Compute a non-embedded resolution from the given embedded one by // dropping redundant trailing blow-ups //---------------------------------------------------------------------------- list resu,tmpiden,templist; intvec divcomp; int i,j,k,offset1,offset2,a,b,c,d,q,found; //--- compute non-embedded resolution list abst=abstractR(re); intvec endiv=abst[1]; intvec deleted=abst[2]; //--- identify the divisors in the various final charts list iden=collectDiv(re,deleted)[2]; // list of final divisors list iden0=iden; // backup copy of iden for later use iden=delete(iden,size(iden)); // drop list of endRings from iden //--------------------------------------------------------------------------- // In iden, only the final charts should be listed, whereas iden0 contains // everything. //--------------------------------------------------------------------------- for(i=1;i<=size(iden);i++) { k=size(iden[i]); tmpiden=iden[i]; for(j=k;j>0;j--) { if(!endiv[iden[i][j][1]]) { //---not a final chart tmpiden=delete(tmpiden,j); } } if(size(tmpiden)==0) { //--- oops, this divisor does not appear in final charts iden=delete(iden,i); continue; } else { iden[i]=tmpiden; } } //--------------------------------------------------------------------------- // Even though the exceptional divisors were irreducible in the embedded // case, they may very well have become reducible after intersection with // the strict transform of the original object. // ===> compute a decomposition for each divisor in each of the final charts // and change the entries of iden accordingly // In particular, it is important to keep track of the identification of the // components of the divisors in each of the charts //--------------------------------------------------------------------------- int n=size(iden); for(i=1;i<=size(re[2]);i++) { if(endiv[i]) { def SN=re[2][i]; setring SN; if(defined(dcE)){kill dcE;} list dcE=decompEinX(BO); // decomposition of exceptional divisors export(dcE); setring R; kill SN; } } if(defined(tmpiden)){kill tmpiden;} list tmpiden=iden; for(i=1;i<=size(iden);i++) { for(j=size(iden[i]);j>0;j--) { def SN=re[2][iden[i][j][1]]; setring SN; if(size(dcE[iden[i][j][2]])==1) { if(dcE[iden[i][j][2]][1][2]==0) { tmpiden[i]=delete(tmpiden[i],j); } } setring R; kill SN; } } for(i=size(tmpiden);i>0;i--) { if(size(tmpiden[i])==0) { tmpiden=delete(tmpiden,i); } } iden=tmpiden; kill tmpiden; list tmpiden; //--- change entries of iden accordingly for(i=1;i<=size(iden);i++) { //--- first set up new entries in iden if necessary - using the first chart //--- in which we see the respective exceptional divisor if(defined(S)){kill S;} def S=re[2][iden[i][1][1]]; //--- considering first entry for i-th divisor setring S; a=size(dcE[iden[i][1][2]]); for(j=1;j<=a;j++) { //--- reducible - add to the list considering each component as an exceptional //--- divisor in its own right list tl; tl[1]=intvec(iden[i][1][1],iden[i][1][2],j); tmpiden[size(tmpiden)+1]=tl; kill tl; } //--- now identify the components in the other charts w.r.t. the ones in the //--- first chart which have already been added to the list for(j=2;j<=size(iden[i]);j++) { //--- considering remaining entries for the same original divisor if(defined(S2)){kill S2;} def S2=re[2][iden[i][j][1]]; setring S2; //--- determine common parent of this ring and re[2][iden[i][1][1]] if(defined(opath)){kill opath;} def opath=imap(S,path); b=1; while(opath[1,b]==path[1,b]) { b++; if((b>ncols(path))||(b>ncols(opath))) break; } if(defined(li1)){kill li1;} list li1; //--- fetch the components we have considered in re[2][iden[i][1][1]] //--- via the resolution tree for(k=1;k<=a;k++) { string tempstr="dcE["+string(eval(iden[i][1][2]))+"]["+string(k)+"][1]"; if(defined(id1)){kill id1;} ideal id1=fetchInTree(re,iden[i][1][1],int(leadcoef(path[1,b-1])), iden[i][j][1],tempstr,iden0,1); kill tempstr; li1[k]=radical(id1); // for comparison only the geometric // object matters kill id1; } //--- compare the components we have fetched with the components in the //--- current ring for(k=1;k<=size(dcE[iden[i][j][2]]);k++) { found=0; for(b=1;b<=size(li1);b++) { if((size(reduce(li1[b],std(dcE[iden[i][j][2]][k][1])))==0)&& (size(reduce(dcE[iden[i][j][2]][k][1],std(li1[b]+BO[2])))==0)) { li1[b]=ideal(1); tmpiden[size(tmpiden)-a+b][size(tmpiden[size(tmpiden)-a+b])+1]= intvec(iden[i][j][1],iden[i][j][2],k); found=1; break; } } if(!found) { if(!defined(repair)) { list repair; repair[1]=list(intvec(iden[i][j][1],iden[i][j][2],k)); } else { for(c=1;c<=size(repair);c++) { for(d=1;d<=size(repair[c]);d++) { if(defined(opath)) {kill opath;} def opath=imap(re[2][repair[c][d][1]],path); q=0; while(path[1,q+1]==opath[1,q+1]) { q++; if((q>ncols(path)-1)||(q>ncols(opath)-1)) break; } q=int(leadcoef(path[1,q])); string tempstr="dcE["+string(eval(repair[c][d][2]))+"]["+string(eval(repair[c][d][3]))+"][1]"; if(defined(id1)){kill id1;} ideal id1=fetchInTree(re,repair[c][d][1],q, iden[i][j][1],tempstr,iden0,1); kill tempstr; //!!! sind die nicht schon radical? id1=radical(id1); // for comparison // only the geometric // object matters if((size(reduce(dcE[iden[i][j][2]][k][1],std(id1+BO[2])))==0)&& (size(reduce(id1+BO[2],std(dcE[iden[i][j][2]][k][1])))==0)) { repair[c][size(repair[c])+1]=intvec(iden[i][j][1],iden[i][j][2],k); break; } } if(d<=size(repair[c])) { break; } } if(c>size(repair)) { repair[size(repair)+1]=list(intvec(iden[i][j][1],iden[i][j][2],k)); } } } } } if(defined(repair)) { for(c=1;c<=size(repair);c++) { tmpiden[size(tmpiden)+1]=repair[c]; } kill repair; } } setring R; for(i=size(tmpiden);i>0;i--) { if(size(tmpiden[i])==0) { tmpiden=delete(tmpiden,i); continue; } } iden=tmpiden; // store the modified divisor list kill tmpiden; // and clean up temporary objects //--------------------------------------------------------------------------- // Now we have decomposed everything into irreducible components over Q, // but over C there might still be some reducible ones left: // Determine the number of components over C. //--------------------------------------------------------------------------- n=0; for(i=1;i<=size(iden);i++) { if(defined(S)) {kill S;} def S=re[2][iden[i][1][1]]; setring S; divcomp[i]=ncols(dcE[iden[i][1][2]][iden[i][1][3]][4]); // number of components of the Q-irreducible curve dcE[iden[i][1][2]] n=n+divcomp[i]; setring R; } //--------------------------------------------------------------------------- // set up the entries Inters[i,j] , i!=j, in the intersection matrix: // we have to compute the intersection of the exceptional divisors (over C) // i.e. we have to work in over appropriate algebraic extension of Q. // (1) plug the intersection matrices of the components of the same Q-irred. // divisor into the correct position in the intersection matrix // (2) for comparison of Ei,k and Ej,l move to a chart where both divisors // are present, fetch the components from the very first chart containing // the respective divisor and then compare by using intersComp // (4) put the result into the correct position in the integer matrix Inters //--------------------------------------------------------------------------- //--- some initialization int comPai,comPaj; intvec v,w; intmat Inters[n][n]; //--- run through all Q-irreducible exceptional divisors for(i=1;i<=size(iden);i++) { if(divcomp[i]>1) { //--- (1) put the intersection matrix for Ei,k with Ei,l into the correct place for(k=1;k<=size(iden[i]);k++) { if(defined(tempmat)){kill tempmat;} intmat tempmat=imap(re[2][iden[i][k][1]],dcE)[iden[i][k][2]][iden[i][k][3]][4]; if(size(ideal(tempmat))!=0) { Inters[i+offset1..(i+offset1+divcomp[i]-1), i+offset1..(i+offset1+divcomp[i]-1)]= tempmat[1..nrows(tempmat),1..ncols(tempmat)]; break; } kill tempmat; } } offset2=offset1+divcomp[i]-1; //--- set up the components over C of the i-th exceptional divisor if(defined(S)){kill S;} def S=re[2][iden[i][1][1]]; setring S; if(defined(idlisti)) {kill idlisti;} list idlisti; idlisti[1]=dcE[iden[i][1][2]][iden[i][1][3]][6]; export(idlisti); setring R; //--- run through the remaining exceptional divisors and check whether they //--- have a chart in common with the i-th divisor for(j=i+1;j<=size(iden);j++) { kill templist; list templist; for(k=1;k<=size(iden[i]);k++) { intvec tiv2=findInIVList(1,iden[i][k][1],iden[j]); if(size(tiv2)!=1) { //--- tiv2[1] is a common chart for the divisors i and j tiv2[4..6]=iden[i][k]; templist[size(templist)+1]=tiv2; } kill tiv2; } if(size(templist)==0) { //--- the two (Q-irred) divisors do not appear in any chart simultaneously offset2=offset2+divcomp[j]-1; j++; continue; } for(k=1;k<=size(templist);k++) { if(defined(S)) {kill S;} //--- set up the components over C of the j-th exceptional divisor def S=re[2][iden[j][1][1]]; setring S; if(defined(idlistj)) {kill idlistj;} list idlistj; idlistj[1]=dcE[iden[j][1][2]][iden[j][1][3]][6]; export(idlistj); if(defined(opath)){kill opath;} def opath=imap(re[2][templist[k][1]],path); comPaj=1; while(opath[1,comPaj]==path[1,comPaj]) { comPaj++; if((comPaj>ncols(opath))||(comPaj>ncols(path))) break; } comPaj=int(leadcoef(path[1,comPaj-1])); setring R; kill S; def S=re[2][iden[i][1][1]]; setring S; if(defined(opath)){kill opath;} def opath=imap(re[2][templist[k][1]],path); comPai=1; while(opath[1,comPai]==path[1,comPai]) { comPai++; if((comPai>ncols(opath))||(comPai>ncols(path))) break; } comPai=int(leadcoef(opath[1,comPai-1])); setring R; kill S; def S=re[2][templist[k][1]]; setring S; if(defined(il)) {kill il;} if(defined(jl)) {kill jl;} if(defined(str1)) {kill str1;} if(defined(str2)) {kill str2;} string str1="idlisti"; string str2="idlistj"; attrib(str1,"algext",imap(re[2][iden[i][1][1]],dcE)[iden[i][1][2]][iden[i][1][3]][5]); attrib(str2,"algext",imap(re[2][iden[j][1][1]],dcE)[iden[j][1][2]][iden[j][1][3]][5]); list il=fetchInTree(re,iden[i][1][1],comPai, templist[k][1],str1,iden0,1); list jl=fetchInTree(re,iden[j][1][1],comPaj, templist[k][1],str2,iden0,1); list nulli=imap(re[2][iden[i][1][1]],dcE)[iden[i][1][2]][iden[i][1][3]][7]; list nullj=imap(re[2][iden[j][1][1]],dcE)[iden[j][1][2]][iden[j][1][3]][7]; string mpi=imap(re[2][iden[i][1][1]],dcE)[iden[i][1][2]][iden[i][1][3]][5]; string mpj=imap(re[2][iden[j][1][1]],dcE)[iden[j][1][2]][iden[j][1][3]][5]; if(defined(tintMat)){kill tintMat;} intmat tintMat=intersComp(il[1],mpi,nulli,jl[1],mpj,nullj); kill mpi; kill mpj; kill nulli; kill nullj; for(a=1;a<=divcomp[i];a++) { for(b=1;b<=divcomp[j];b++) { if(tintMat[a,b]!=0) { Inters[i+offset1+a-1,j+offset2+b-1]=tintMat[a,b]; Inters[j+offset2+b-1,i+offset1+a-1]=tintMat[a,b]; } } } } offset2=offset2+divcomp[j]-1; } offset1=offset1+divcomp[i]-1; } Inters=addSelfInter(re,Inters,iden,iden0,endiv); intvec GenusIden; list tl_genus; a=1; for(i=1;i<=size(iden);i++) { tl_genus=genus_E(re,iden0,iden[i][1]); for(j=1;j<=tl_genus[2];j++) { GenusIden[a]=tl_genus[1]; a++; } } list retlist=Inters,GenusIden,iden,divcomp; return(retlist); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp(3); ideal J=x(3)^5+x(2)^4+x(1)^3+x(1)*x(2)*x(3); list re=resolve(J); list di=intersectionDiv(re); di; } ////////////////////////////////////////////////////////////////////////////// static proc intersComp(string str1, string mp1, list null1, string str2, string mp2, list null2) "Internal procedure - no help and no example available " { //--- format of input //--- str1 : ideal (over field extension 1) //--- mp1 : minpoly of field extension 1 //--- null1: numerical zeros of minpoly //--- str2 : ideal (over field extension 2) //--- mp2 : minpoly of field extension 2 //--- null2: numerical zeros of minpoly //--- determine intersection matrix of the C-components defined by the input //--------------------------------------------------------------------------- // Initialization //--------------------------------------------------------------------------- int ii,jj,same; def R=basering; intmat InterMat[size(null1)][size(null2)]; ring ringst=0,(t,s),dp; //--------------------------------------------------------------------------- // Add new variables s and t and compare the minpolys and ideals // to find out whether they are identical //--------------------------------------------------------------------------- def S=R+ringst; setring S; if((mp1==mp2)&&(str1==str2)) { same=1; } //--- define first Q-component/C-components, substitute t by s string tempstr="ideal id1="+str1+";"; execute(tempstr); execute(mp1); id1=subst(id1,t,s); poly q=subst(p,t,s); kill p; //--- define second Q-component/C-components tempstr="ideal id2="+str2+";"; execute(tempstr); execute(mp2); //--- do the intersection ideal interId=id1+id2+ideal(p)+ideal(q); if(same) { interId=quotient(interId,t-s); } interId=std(interId); //--- refine the comparison by passing to each of the numerical zeros //--- of the two minpolys ideal stid=nselect(interId,1..nvars(R)); ring compl_st=complex,(s,t),dp; def stid=imap(S,stid); ideal tempid,tempid2; for(ii=1;ii<=size(null1);ii++) { tempstr="number numi="+null1[ii]+";"; execute(tempstr); tempid=subst(stid,s,numi); kill numi; for(jj=1;jj<=size(null2);jj++) { tempstr="number numj="+null2[jj]+";"; execute(tempstr); tempid2=subst(tempid,t,numj); kill numj; if(size(tempid2)==0) { InterMat[ii,jj]=1; } } } //--- sanity check; as both Q-components were Q-irreducible, //--- summation over all entries of a single row must lead to the same //--- result, no matter which row is chosen //--- dito for the columns int cou,cou1; for(ii=1;ii<=ncols(InterMat);ii++) { cou=0; for(jj=1;jj<=nrows(InterMat);jj++) { cou=cou+InterMat[jj,ii]; } if(ii==1){cou1=cou;} if(cou1!=cou){ERROR("intersComp:matrix has wrong entries");} } for(ii=1;ii<=nrows(InterMat);ii++) { cou=0; for(jj=1;jj<=ncols(InterMat);jj++) { cou=cou+InterMat[ii,jj]; } if(ii==1){cou1=cou;} if(cou1!=cou){ERROR("intersComp:matrix has wrong entries");} } return(InterMat); } ///////////////////////////////////////////////////////////////////////////// static proc addSelfInter(list re,intmat Inters,list iden,list iden0,intvec endiv) "Internal procedure - no help and no example available " { //--------------------------------------------------------------------------- // Initialization //--------------------------------------------------------------------------- def R=basering; int i,j,k,l,a,b; int n=size(iden); intvec v,w; list satlist; def T=re[2][1]; setring T; poly p; p=var(1); //any linear form will do, //but this one is most convenient ideal F=ideal(p); //---------------------------------------------------------------------------- // lift linear form to every end ring, determine the multiplicity of // the exceptional divisors and store it in Flist //---------------------------------------------------------------------------- list templist; intvec tiv; for(i=1;i<=size(endiv);i++) { if(endiv[i]==1) { kill v; intvec v; a=0; if(defined(S)) {kill S;} def S=re[2][i]; setring S; map resi=T,BO[5]; ideal F=resi(F)+BO[2]; ideal Ftemp=F; list Flist; if(defined(satlist)){kill satlist;} list satlist; for(a=1;a<=size(dcE);a++) { for(b=1;b<=size(dcE[a]);b++) { Ftemp=sat(Ftemp,dcE[a][b][1])[1]; } } F=sat(F,Ftemp)[1]; Flist[1]=Ftemp; Ftemp=1; list pr=primdecGTZ(F); v[size(pr)]=0; for(j=1;j<=size(pr);j++) { for(a=1;a<=size(dcE);a++) { if(j==1) { kill tiv; intvec tiv; tiv[size(dcE[a])]=0; templist[a]=tiv; if(v[j]==1) { a++; continue; } } if(dcE[a][1][2]==0) { a++; continue; } for(b=1;b<=size(dcE[a]);b++) { if((size(reduce(dcE[a][b][1],std(pr[j][2])))==0)&& (size(reduce(pr[j][2],std(dcE[a][b][1])))==0)) { templist[a][b]=Vielfachheit(pr[j][1],pr[j][2]); v[j]=1; break; } } if((v[j]==1)&&(j>1)) break; } } kill v; intvec v; Flist[2]=templist; } } //----------------------------------------------------------------------------- // Now set up all the data: // 1. run through all exceptional divisors in iden and determine the // coefficients c_i of the divisor of F. ===> civ // 2. determine the intersection locus of F^bar and the Ei and from this data // the F^bar.Ei . ===> intF //----------------------------------------------------------------------------- intvec civ; intvec intF; intF[ncols(Inters)]=0; int offset,comPa,ncomp,vd; for(i=1;i<=size(iden);i++) { ncomp=0; for(j=1;j<=size(iden[i]);j++) { if(defined(S)) {kill S;} def S=re[2][iden[i][j][1]]; setring S; if((size(civ)0) { civ[i+offset+k]=civ[i+k]; } } } if(defined(interId)) {kill interId;} ideal interId=dcE[iden[i][j][2]][iden[i][j][3]][1]+Flist[1]; if(defined(interList)) {kill interList;} list interList; interList[1]=string(interId); interList[2]=ideal(0); export(interList); if(defined(doneId)) {kill doneId;} if(defined(tempId)) {kill tempId;} ideal doneId=ideal(1); if(defined(dl)) {kill dl;} list dl; for(k=1;kncols(path))||(comPa>ncols(opath))) break; } comPa=int(leadcoef(path[1,comPa-1])); if(defined(str)) {kill str;} string str="interList"; attrib(str,"algext","poly p=t-1;"); dl=fetchInTree(re,iden[i][k][1],comPa,iden[i][j][1],str,iden0,1); if(defined(tempId)){kill tempId;} str="ideal tempId="+dl[1]+";"; execute(str); doneId=intersect(doneId,tempId); str="interId="+interList[1]+";"; execute(str); interId=sat(interId,doneId)[1]; interList[1]=string(interId); } interId=std(interId); if(dim(interId)>0) { "oops, intersection not a set of points"; ~; } vd=vdim(interId); if(vd>0) { for(k=i+offset;k<=i+offset+ncomp-1;k++) { intF[k]=intF[k]+(vd/ncomp); } } } offset=size(civ)-i-1; } if(defined(tiv)){kill tiv;} intvec tiv=civ[2..size(civ)]; civ=tiv; kill tiv; //----------------------------------------------------------------------------- // Using the F_total= sum c_i Ei + F^bar, the intersection matrix Inters and // the f^bar.Ei, determine the selfintersection numbers of the Ei from the // equation F_total.Ei=0 and store it in the diagonal of Inters. //----------------------------------------------------------------------------- intvec diag=Inters*civ+intF; for(i=1;i<=size(diag);i++) { Inters[i,i]=-diag[i]/civ[i]; } return(Inters); } ////////////////////////////////////////////////////////////////////////////// static proc invSort(list re, list #) "Internal procedure - no help and no example available " { int i,j,k,markier,EZeiger,offset; intvec v,e; intvec deleted; if(size(#)>0) { deleted=#[1]; } else { deleted[size(re[2])]=0; } list LE,HI; def R=basering; //---------------------------------------------------------------------------- // Go through all rings //---------------------------------------------------------------------------- for(i=1;i<=size(re[2]);i++) { if(deleted[i]){i++;continue} def S=re[2][i]; setring S; //---------------------------------------------------------------------------- // Determine Invariant //---------------------------------------------------------------------------- if((size(BO[3])==size(BO[9]))||(size(BO[3])==size(BO[9])+1)) { if(defined(merk2)){kill merk2;} intvec merk2; EZeiger=0; for(j=1;j<=size(BO[9]);j++) { offset=0; if(BO[7][j]==-1) { BO[7][j]=size(BO[4])-EZeiger; } for(k=EZeiger+1;(k<=EZeiger+BO[7][j])&&(k<=size(BO[4]));k++) { if(BO[6][k]==2) { offset++; } } EZeiger=EZeiger+BO[7][1]; merk2[3*j-2]=BO[3][j]; merk2[3*j-1]=BO[9][j]-offset; if(size(invSat[2])>j) { merk2[3*j]=-invSat[2][j]; } else { if(jsize(BO[9]))) { merk2[size(merk2)+1]=BO[3][size(BO[3])]; } if((size(merk2)%3)==0) { intvec tintvec=merk2[1..size(merk2)-1]; merk2=tintvec; kill tintvec; } } else { ERROR("This situation should not occur, please send the example to the authors."); } //---------------------------------------------------------------------------- // Save invariant describing current center as an object in this ring // We also store information on the intersection with the center and the // exceptional divisors //---------------------------------------------------------------------------- cent=std(cent); kill e; intvec e; for(j=1;j<=size(BO[4]);j++) { if(size(reduce(BO[4][j],std(cent+BO[1])))==0) { e[j]=1; } else { e[j]=0; } } if(size(ideal(merk2))==0) { markier=1; } if((size(merk2)%3==0)&&(merk2[size(merk2)]==0)) { intvec blabla=merk2[1..size(merk2)-1]; merk2=blabla; kill blabla; } if(defined(invCenter)){kill invCenter;} list invCenter=cent,merk2,e; export invCenter; //---------------------------------------------------------------------------- // Insert it into correct place in the list //---------------------------------------------------------------------------- if(i==1) { if(!markier) { HI=intvec(merk2[1]+1),intvec(1); } else { HI=intvec(778),intvec(1); // some really large integer // will be changed at the end!!! } LE[1]=HI; i++; setring R; kill S; continue; } if(markier==1) { if(i==2) { HI=intvec(777),intvec(2); // same really large integer-1 LE[2]=HI; i++; setring R; kill S; continue; } else { if(ncols(path)==2) { LE[2][2][size(LE[2][2])+1]=i; i++; setring R; kill S; continue; } else { markier=0; } } } j=1; def SOld=re[2][int(leadcoef(path[1,ncols(path)]))]; setring SOld; merk2=invCenter[2]; setring S; kill SOld; while(merk2size(LE)) break; } HI=merk2,intvec(i); if(j<=size(LE)) { if(merk2>LE[j][1]) { LE=insert(LE,HI,j-1); } else { while((merk2==LE[j][1])&&(size(merk2)size(LE)) break; } if(j<=size(LE)) { if((merk2!=LE[j][1])||(size(merk2)!=size(LE[j][1]))) { LE=insert(LE,HI,j-1); } else { LE[j][2][size(LE[j][2])+1]=i; } } else { LE[size(LE)+1]=HI; } } } else { LE[size(LE)+1]=HI; } setring R; kill S; } if((LE[1][1]==intvec(778)) && (size(LE)>2)) { LE[1][1]=intvec(LE[3][1][1]+2); // by now we know what 'sufficiently LE[2][1]=intvec(LE[3][1][1]+1); // large' is } return(LE); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp(3); ideal J=x(1)^3-x(1)*x(2)^3+x(3)^2; list re=resolve(J,1); list di=invSort(re); di; } ///////////////////////////////////////////////////////////////////////////// static proc addToRE(intvec v,int x,list RE) "Internal procedure - no help and no example available " { //--- auxilliary procedure for collectDiv, //--- inserting an entry at the correct place int i=1; while(i<=size(RE)) { if(v==RE[i][1]) { RE[i][2][size(RE[i][2])+1]=x; return(RE); } if(v>RE[i][1]) { list templist=v,intvec(x); RE=insert(RE,templist,i-1); return(RE); } i++; } list templist=v,intvec(x); RE=insert(RE,templist,size(RE)); return(RE); } //////////////////////////////////////////////////////////////////////////// proc collectDiv(list re,list #) "USAGE: collectDiv(L); @* L = list of rings ASSUME: L is output of resolution of singularities COMPUTE: list representing the identification of the exceptional divisors in the various charts RETURN: list l, where l[1]: intmat, entry k in position i,j implies BO[4][j] of chart i is divisor k (if k!=0) if k==0, no divisor corresponding to i,j l[2]: list ll, where each entry of ll is a list of intvecs entry i,j in list ll[k] implies BO[4][j] of chart i is divisor k l[3]: list L EXAMPLE: example collectDiv; shows an example " { //------------------------------------------------------------------------ // Initialization //------------------------------------------------------------------------ int i,j,k,l,m,maxk,maxj,mPa,oPa,interC,pa,ignoreL,iTotal; int mLast,oLast=1,1; intvec deleted; //--- sort the rings by the invariant which controlled the last of the //--- exceptional divisors if(size(#)>0) { deleted=#[1]; } else { deleted[size(re[2])]=0; } list LE=invSort(re,deleted); list LEtotal=LE; intmat M[size(re[2])][size(re[2])]; intvec invar,tempiv; def R=basering; list divList; list RE,SE; intvec myEi,otherEi,tempe; int co=2; while(size(LE)>0) { //------------------------------------------------------------------------ // Run through the sorted list LE whose entries are lists containing // the invariant and the numbers of all rings corresponding to it //------------------------------------------------------------------------ for(i=co;i<=size(LE);i++) { //--- i==1 in first iteration: //--- the original ring which did not arise from a blow-up //--- hence there are no exceptional divisors to be identified there ; //------------------------------------------------------------------------ // For each fixed value of the invariant, run through all corresponding // rings //------------------------------------------------------------------------ for(l=1;l<=size(LE[i][2]);l++) { if(defined(S)){kill S;} def S=re[2][LE[i][2][l]]; setring S; if(size(BO[4])>maxj){maxj=size(BO[4]);} //--- all exceptional divisors, except the last one, were previously //--- identified - hence we can simply inherit the data from the parent ring for(j=1;j0) { k=int(leadcoef(path[1,ncols(path)])); k=M[k,j]; if(k==0) { RE=addToRE(LE[i][1],LE[i][2][l],RE); ignoreL=1; break; } M[LE[i][2][l],j]=k; tempiv=LE[i][2][l],j; divList[k][size(divList[k])+1]=tempiv; } } if(ignoreL){ignoreL=0;l++;continue;} //---------------------------------------------------------------------------- // In the remaining part of the procedure, the identification of the last // exceptional divisor takes place. // Step 1: check whether there is a previously considered ring with the // same parent; if this is the case, we can again inherit the data // Step 1':check whether the parent had a stored center which it then used // in this case, we are dealing with an additional component of this // divisor: store it in the integer otherComp // Step 2: if no appropriate ring was found in step 1, we check whether // there is a previously considered ring, in the parent of which // the center intersects the same exceptional divisors as the center // in our parent. // if such a ring does not exist: new exceptional divisor // if it exists: see below //---------------------------------------------------------------------------- if(path[1,ncols(path)-1]==0) { //--- current ring originated from very first blow-up //--- hence exceptional divisor is the first one M[LE[i][2][l],1]=1; if(size(divList)>0) { divList[1][size(divList[1])+1]=intvec(LE[i][2][l],j); } else { divList[1]=list(intvec(LE[i][2][l],j)); } l++; continue; } if(l==1) { list TE=addToRE(LE[i][1],1,SE); if(size(TE)!=size(SE)) { //--- new value of invariant hence new exceptional divisor SE=TE; divList[size(divList)+1]=list(intvec(LE[i][2][l],j)); M[LE[i][2][l],j]=size(divList); } kill TE; } for(k=1;k<=size(LEtotal);k++) { if(LE[i][1]==LEtotal[k][1]) { iTotal=k; break; } } //--- Step 1 k=1; while(LEtotal[iTotal][2][k]size(LEtotal[iTotal][2])) {break;} } if(ignoreL){ignoreL=0;l++;continue;} //--- Step 1', if necessary if(M[LE[i][2][l],j]==0) { int savedCent; def SPa1=re[2][int(leadcoef(path[1,ncols(path)]))]; // parent ring setring SPa1; if(size(BO)>9) { if(size(BO[10])>0) { savedCent=1; } } if(!savedCent) { def SPa2=re[2][int(leadcoef(path[1,ncols(path)]))]; map lMa=SPa2,lastMap; // map leading from grandparent to parent list transBO=lMa(BO); // actually we only need BO[10], but this is an // object not a name list tempsat; if(size(transBO)>9) { //--- there were saved centers while((k<=size(transBO[10])) & (savedCent==0)) { tempsat=sat(transBO[10][k][1],BO[4][size(BO[4])]); if(deg(tempsat[1][1])!=0) { //--- saved center can be seen in this affine chart if((size(reduce(tempsat[1],std(cent)))==0) && (size(reduce(cent,tempsat[1]))==0)) { //--- this was the saved center which was used savedCent=1; } } k++; } } kill lMa; // clean up temporary objects kill tempsat; kill transBO; } setring S; // back to the ring which we want to consider if(savedCent==1) { vector otherComp; otherComp[M[int(leadcoef(path[1,ncols(path)])),size(BO[4])-1]] =1; } kill savedCent; if (defined(SPa2)){kill SPa2;} kill SPa1; } //--- Step 2, if necessary if(M[LE[i][2][l],j]==0) { //--- we are not done after step 1 and 2 pa=int(leadcoef(path[1,ncols(path)])); // parent ring tempe=imap(re[2][pa],invCenter)[3]; // intersection there kill myEi; intvec myEi; for(k=1;k<=size(tempe);k++) { if(tempe[k]==1) { //--- center meets this exceptional divisor myEi[size(myEi)+1]=M[pa,k]; mLast=k; } } //--- ring in which the last divisor we meet is new-born mPa=int(leadcoef(path[1,mLast+2])); k=1; while(LEtotal[iTotal][2][k]size(LEtotal[iTotal][2])) { break; } else { continue; } } //---------------------------------------------------------------------------- // Current situation: // 1. the last exceptional divisor could not be identified by simply // considering its parent // 2. it could not be proved to be a new one by considering its intersections // with previous exceptional divisors //---------------------------------------------------------------------------- if(defined(bool1)) { kill bool1;} int bool1= compareE(re,LE[i][2][l],LEtotal[iTotal][2][k],divList); if(bool1) { //--- found some non-empty intersection if(bool1==1) { //--- it is really the same exceptional divisor m=size(imap(re[2][LEtotal[iTotal][2][k]],BO)[4]); m=M[LEtotal[iTotal][2][k],m]; if(m==0) { RE=addToRE(LE[i][1],LE[i][2][l],RE); ignoreL=1; break; } M[LE[i][2][l],j]=m; tempiv=LE[i][2][l],j; divList[m][size(divList[m])+1]=tempiv; break; } else { m=size(imap(re[2][LEtotal[iTotal][2][k]],BO)[4]); m=M[LEtotal[iTotal][2][k],m]; if(m!=0) { otherComp[m]=1; } } } k++; if(k>size(LEtotal[iTotal][2])) { break; } } if(ignoreL){ignoreL=0;l++;continue;} if( M[LE[i][2][l],j]==0) { divList[size(divList)+1]=list(intvec(LE[i][2][l],j)); M[LE[i][2][l],j]=size(divList); } } setring R; kill S; } } LE=RE; co=1; kill RE; list RE; } //---------------------------------------------------------------------------- // Add the strict transform to the list of divisors at the last place // and clean up M //---------------------------------------------------------------------------- //--- add strict transform for(i=1;i<=size(re[2]);i++) { if(defined(S)){kill S;} def S=re[2][i]; setring S; if(size(reduce(cent,std(BO[2])))==0) { tempiv=i,0; RE[size(RE)+1]=tempiv; } setring R; } divList[size(divList)+1]=RE; //--- drop trailing zero-columns of M intvec iv0; iv0[nrows(M)]=0; for(i=ncols(M);i>0;i--) { if(intvec(M[1..nrows(M),i])!=iv0) break; } intmat N[nrows(M)][i]; for(i=1;i<=ncols(N);i++) { N[1..nrows(M),i]=M[1..nrows(M),i]; } kill M; intmat M=N; list retlist=cleanUpDiv(re,M,divList); return(retlist); } example {"EXAMPLE:"; echo = 2; ring R=0,(x,y,z),dp; ideal I=xyz+x4+y4+z4; //we really need to blow up curves even if the generic point of //the curve the total transform is n.c. //this occurs here in r[2][5] list re=resolve(I); list di=collectDiv(re); di[1]; di[2]; } ////////////////////////////////////////////////////////////////////////////// static proc cleanUpDiv(list re,intmat M,list divList) "Internal procedure - no help and no example available " { //--- It may occur that two different entries of invSort coincide on the //--- first part up to the last entry of the shorter one. In this case //--- exceptional divisors may appear in both entries of the invSort-list. //--- To correct this, we now compare the final collection of Divisors //--- for coinciding ones. int i,j,k,a,oPa,mPa,comPa,mdim,odim; def R=basering; for(i=1;i<=size(divList)-2;i++) { if(defined(Sm)){kill Sm;} def Sm=re[2][divList[i][1][1]]; setring Sm; mPa=int(leadcoef(path[1,ncols(path)])); if(defined(SmPa)){kill SmPa;} def SmPa=re[2][mPa]; setring SmPa; mdim=dim(std(BO[1]+cent)); setring Sm; if(mPa==1) { //--- very first divisor originates exactly from the first blow-up //--- there cannot be any mistake here i++; continue; } for(j=i+1;j<=size(divList)-1;j++) { setring Sm; for(k=1;k<=size(divList[j]);k++) { if(size(findInIVList(1,divList[j][k][1],divList[i]))>1) { //--- same divisor cannot appear twice in the same chart k=-1; break; } } if(k==-1) { j++; if(j>size(divList)-1) break; continue; } if(defined(opath)){kill opath;} def opath=imap(re[2][divList[j][1][1]],path); oPa=int(leadcoef(opath[1,ncols(opath)])); if(defined(SoPa)){kill SoPa;} def SoPa=re[2][oPa]; setring SoPa; odim=dim(std(BO[1]+cent)); setring Sm; if(mdim!=odim) { //--- different dimension ==> cannot be same center j++; if(j>size(divList)-1) break; continue; } comPa=1; while(path[1,comPa]==opath[1,comPa]) { comPa++; if((comPa>ncols(path))||(comPa>ncols(opath))) break; } comPa=int(leadcoef(path[1,comPa-1])); if(defined(SPa)){kill SPa;} def SPa=re[2][mPa]; setring SPa; if(defined(tempIdE)){kill tempIdE;} ideal tempIdE=fetchInTree(re,oPa,comPa,mPa,"cent",divList); if((size(reduce(cent,std(tempIdE)))!=0)|| (size(reduce(tempIdE,std(cent)))!=0)) { //--- it is not the same divisor! j++; if(j>size(divList)) { break; } else { continue; } } for(k=1;k<=size(divList[j]);k++) { //--- append the entries of the j-th divisor (which is actually also the i-th) //--- to the i-th divisor divList[i][size(divList[i])+1]=divList[j][k]; } divList=delete(divList,j); //kill obsolete entry from the list for(k=1;k<=nrows(M);k++) { for(a=1;a<=ncols(M);a++) { if(M[k,a]==j) { //--- j-th divisor is actually the i-th one M[k,a]=i; } if(M[k,a]>j) { //--- index j was deleted from the list ==> all subsequent indices dropped by //--- one M[k,a]=M[k,a]-1; } } } j--; //do not forget to consider new j-th entry } } setring R; list retlist=M,divList; return(retlist); } ///////////////////////////////////////////////////////////////////////////// static proc findTrans(ideal Z, ideal E, list notE, list #) "Internal procedure - no help and no example available " { //---Auxilliary procedure for fetchInTree! //---Assume E prime ideal, Z+E eqidimensional, //---ht(E)+r=ht(Z+E). Compute P= in Z+E, and polynomial f, //---such that radical(Z+E)=radical((E+P):f) int i,j,d,e; ideal Estd=std(E); //!!! alternative to subsequent line: //!!! ideal Zstd=std(radical(Z+E)); ideal Zstd=std(Z+E); ideal J=1; if(size(#)>0) { J=#[1]; } if(deg(Zstd[1])==0){return(list(ideal(1),poly(1)));} for(i=1;i<=size(notE);i++) { notE[i]=std(notE[i]); } ideal Zred=simplify(reduce(Z,Estd),2); if(size(Zred)==0){Z,Estd;~;ERROR("Z is contained in E");} ideal P,Q,Qstd; Q=Estd; attrib(Q,"isSB",1); d=dim(Estd); e=dim(Zstd); for(i=1;i<=size(Zred);i++) { Qstd=std(Q,Zred[i]); if(dim(Qstd)20) { ~; ERROR("findTrans:Hier ist was faul"); } } list resu=P,f; return(resu); } ///////////////////////////////////////////////////////////////////////////// static proc compareE(list L, int m, int o, list DivL) "Internal procedure - no help and no example available " { //---------------------------------------------------------------------------- // We want to compare the divisors BO[4][size(BO[4])] of the rings // L[2][m] and L[2][o]. // In the initialization step, we collect all necessary data from those // those rings. In particular, we determine at what point (in the resolution // history) the branches for L[2][m] and L[2][o] were separated, denoting // the corresponding ring indices by mPa, oPa and comPa. //---------------------------------------------------------------------------- def R=basering; int i,j,k,len; //-- find direct parents and branching point in resolution history matrix tpm=imap(L[2][m],path); matrix tpo=imap(L[2][o],path); int m1,o1=int(leadcoef(tpm[1,ncols(tpm)])), int(leadcoef(tpo[1,ncols(tpo)])); while((i they cannot give rise //--- to the same divisor return(0); } def T=L[2][o1]; setring T; int dimCo1=dim(std(cent+BO[1])); def S=L[2][m1]; setring S; int dimCm1=dim(std(cent+BO[1])); if(dimCm1!=dimCo1) { //--- centers do not have same dimension ==> they cannot give rise //--- to the same divisor return(0); } //---------------------------------------------------------------------------- // fetch the center via the tree for comparison //---------------------------------------------------------------------------- if(defined(invLocus0)) {kill invLocus0;} ideal invLocus0=fetchInTree(L,o1,comPa,m1,"cent",DivL); // blow down from L[2][o1] to L[2][comPa] and then up to L[2][m1] if(deg(std(invLocus0+invCenter[1]+BO[1])[1])!=0) { setring R; return(int(1)); } if(size(BO)>9) { for(i=1;i<=size(BO[10]);i++) { if(deg(std(invLocus0+BO[10][i][1]+BO[1])[1])!=0) { if(dim(std(BO[10][i][1]+BO[1])) > dim(std(invLocus0+BO[10][i][1]+BO[1]))) { ERROR("Internal Error: Please send this example to the authors."); } setring R; return(int(2)); } } } setring R; return(int(0)); //---------------------------------------------------------------------------- // Return-Values: // TRUE (=1) if the exceptional divisors coincide, // TRUE (=2) if the exceptional divisors originate from different // components of the same center // FALSE (=0) otherwise //---------------------------------------------------------------------------- } ////////////////////////////////////////////////////////////////////////////// proc fetchInTree(list L, int o1, int comPa, int m1, string idname, list DivL, list #); "Internal procedure - no help and no example available " { //---------------------------------------------------------------------------- // Initialization and Sanity Checks //---------------------------------------------------------------------------- int i,j,k,m,branchPos,inJ,exception; string algext; //--- we need to be in L[2][m1] def R=basering; ideal test_for_the_same_ring=-77; def Sm1=L[2][m1]; setring Sm1; if(!defined(test_for_the_same_ring)) { //--- we are not in L[2][m1] ERROR("basering has to coincide with L[2][m1]"); } else { //--- we are in L[2][m1] kill test_for_the_same_ring; } //--- non-embedded case? if(size(#)>0) { inJ=1; } //--- do parameter values make sense? if(comPa<1) { ERROR("Common Parent should at least be the first ring!"); } //--- do we need to pass to an algebraic field extension of Q? if(typeof(attrib(idname,"algext"))=="string") { algext=attrib(idname,"algext"); } //--- check wheter comPa is in the history of m1 //--- same test for o1 can be done later on (on the fly) if(m1==comPa) { j=1; i=ncols(path)+1; } else { for(i=1;i<=ncols(path);i++) { if(int(leadcoef(path[1,i]))==comPa) { //--- comPa occurs in the history j=1; break; } } } branchPos=i; if(j==0) { ERROR("L[2][comPa] not in history of L[2][m1]!"); } //---------------------------------------------------------------------------- // Blow down ideal "idname" from L[2][o1] to L[2][comPa], where the latter // is assumed to be the common parent of L[2][o1] and L[2][m1] //---------------------------------------------------------------------------- if(size(algext)>0) { //--- size(algext)>0: case of algebraic extension of base field if(defined(tstr)){kill tstr;} string tstr="ring So1=(0,t),("+varstr(L[2][o1])+"),("+ordstr(L[2][o1])+");"; execute(tstr); setring So1; execute(algext); minpoly=leadcoef(p); if(defined(id1)) { kill id1; } if(defined(id2)) { kill id2; } if(defined(idlist)) { kill idlist; } execute("int bool2=defined("+idname+");"); if(bool2==0) { execute("list ttlist=imap(L[2][o1],"+idname+");"); } else { execute("list ttlist="+idname+";"); } kill bool2; def BO=imap(L[2][o1],BO); def path=imap(L[2][o1],path); def lastMap=imap(L[2][o1],lastMap); ideal id2=1; if(defined(notE)){kill notE;} list notE; intvec nE; list idlist; for(i=1;i<=size(ttlist);i++) { if((i==size(ttlist))&&(typeof(ttlist[i])!="string")) break; execute("ideal tid="+ttlist[i]+";"); idlist[i]=list(tid,ideal(1),nE); kill tid; } } else { //--- size(algext)==0: no algebraic extension of base needed def So1=L[2][o1]; setring So1; if(defined(id1)) { kill id1; } if(defined(id2)) { kill id2; } if(defined(idlist)) { kill idlist; } execute("ideal id1="+idname+";"); if(deg(std(id1)[1])==0) { //--- problems with findTrans if id1 is empty set //!!! todo: also correct in if branch!!! setring R; return(ideal(1)); } // id1=radical(id1); ideal id2=1; list idlist; if(defined(notE)){kill notE;} list notE; intvec nE; idlist[1]=list(id1,id2,nE); } if(defined(tli)){kill tli;} list tli; if(defined(id1)) { kill id1; } if(defined(id2)) { kill id2; } ideal id1; ideal id2; if(defined(Etemp)){kill Etemp;} ideal Etemp; for(m=1;m<=size(idlist);m++) { //!!! Duplicate Block!!! All changes also needed below!!! //!!! no subprocedure due to large data overhead!!! //--- run through all ideals to be fetched id1=idlist[m][1]; id2=idlist[m][2]; nE=idlist[m][3]; for(i=branchPos-1;i<=size(BO[4]);i++) { //--- run through all relevant exceptional divisors if(size(reduce(BO[4][i],std(id1+BO[1])))==0) { //--- V(id1) is contained in except. div. i in this chart if(size(reduce(id1,std(BO[4][i])))!=0) { //--- V(id1) does not equal except. div. i of this chart Etemp=BO[4][i]; if(npars(basering)>0) { //--- we are in an algebraic extension of the base field if(defined(prtemp)){kill prtemp;} list prtemp=minAssGTZ(BO[4][i]); // C-comp. of except. div. j=1; if(size(prtemp)>1) { //--- more than 1 component Etemp=ideal(1); for(j=1;j<=size(prtemp);j++) { //--- find correct component if(size(reduce(prtemp[j],std(id1)))==0) { Etemp=prtemp[j]; break; } } if(deg(std(Etemp)[1])==0) { ERROR("fetchInTree:something wrong in field extension"); } } prtemp=delete(prtemp,j); // remove this comp. from list while(size(prtemp)>1) { //--- collect all the others into prtemp[1] prtemp[1]=intersect(prtemp[1],prtemp[size(prtemp)]); prtemp=delete(prtemp,size(prtemp)); } } //--- determine tli[1] and tli[2] such that //--- V(id1) \cap D(id2) = V(tli[1]) \cap D(tli[2]) \cap BO[4][i] //--- inside V(BO[1]) (and if necessary inside V(BO[1]+BO[2])) if(inJ) { tli=findTrans(id1+BO[2]+BO[1],Etemp,notE,BO[2]); } else { tli=findTrans(id1+BO[1],Etemp,notE); } if(npars(basering)>0) { //--- in algebraic extension: make sure we stay outside the other components if(size(prtemp)>0) { for(j=1;j<=ncols(prtemp[1]);j++) { //--- find the (univariate) generator of prtemp[1] which is the remaining //--- factor from the factorization over the extension field if(size(reduce(prtemp[1][j],std(id1)))>0) { tli[2]=tli[2]*prtemp[1][j]; } } } } } else { //--- V(id1) equals except. div. i of this chart tli[1]=ideal(0); tli[2]=ideal(1); } id1=tli[1]; id2=id2*tli[2]; notE[size(notE)+1]=BO[4][i]; for(j=1;j<=size(DivL);j++) { if(inIVList(intvec(o1,i),DivL[j])) { nE[size(nE)+1]=j; break; } } if(size(nE)1) { while(int(leadcoef(path[1,ncols(path)]))>=comPa) { if((int(leadcoef(path[1,ncols(path)]))>comPa)&& (int(leadcoef(path[1,ncols(path)-1]))0) { if(defined(T0)){kill T0;} def T0=L[2][int(leadcoef(path[1,ncols(path)]))]; if(defined(tstr)){kill tstr;} string tstr="ring T=(0,t),(" +varstr(L[2][int(leadcoef(path[1,ncols(path)]))])+"),(" +ordstr(L[2][int(leadcoef(path[1,ncols(path)]))])+");"; execute(tstr); setring T; execute(algext); minpoly=leadcoef(p); kill tstr; def BO=imap(T0,BO); if(!defined(und_jetzt_raus)) { def path=imap(T0,path); def lastMap=imap(T0,lastMap); } if(defined(idlist)){kill idlist;} list idlist=list(list(ideal(1),ideal(1))); } else { def T=L[2][int(leadcoef(path[1,ncols(path)]))]; setring T; if(defined(id1)) { kill id1; } if(defined(id2)) { kill id2; } if(defined(idlist)){kill idlist;} list idlist=list(list(ideal(1),ideal(1))); } setring S; if(defined(phi)) { kill phi; } map phi=T,lastMap; //--- now do the actual blowing down ... for(m=1;m<=size(idlist);m++) { //--- ... for each entry of idlist separately if(defined(id1)){kill id1;} if(defined(id2)){kill id2;} ideal id1=idlist[m][1]+BO[1]; ideal id2=idlist[m][2]; nE=idlist[m][3]; if(defined(debug_fetchInTree)>0) { "Blowing down entry",m,"of idlist:"; setring S; "Abbildung:";phi; "before preimage"; id1; id2; } setring T; ideal id1=preimage(S,phi,id1); ideal id2=preimage(S,phi,id2); if(defined(debug_fetchInTree)>0) { "after preimage"; id1; id2; } if(size(id2)==0) { //--- preimage of (principal ideal) id2 was zero, i.e. //--- generator of previous id2 not in image setring S; //--- it might just be one offending factor ==> factorize ideal id2factors=factorize(id2[1])[1]; int zzz=size(id2factors); ideal curfactor; setring T; id2=ideal(1); ideal curfactor; for(int mm=1;mm<=zzz;mm++) { //--- blow down each factor separately setring S; curfactor=id2factors[mm]; setring T; curfactor=preimage(S,phi,curfactor); if(size(curfactor)>0) { id2[1]=id2[1]*curfactor[1]; } } kill curfactor; setring S; kill curfactor; kill id2factors; setring T; kill mm; kill zzz; if(defined(debug_fetchInTree)>0) { "corrected id2:"; id2; } } idlist[m]=list(id1,id2,nE); kill id1,id2; setring S; } setring T; //--- after blowing down we might again be sitting inside a relevant //--- exceptional divisor for(m=1;m<=size(idlist);m++) { //!!! Duplicate Block!!! All changes also needed above!!! //!!! no subprocedure due to large data overhead!!! //--- run through all ideals to be fetched if(defined(id1)) {kill id1;} if(defined(id2)) {kill id2;} if(defined(notE)) {kill notE;} if(defined(notE)) {kill notE;} list notE; ideal id1=idlist[m][1]; ideal id2=idlist[m][2]; nE=idlist[m][3]; for(i=branchPos-1;i<=size(BO[4]);i++) { //--- run through all relevant exceptional divisors if(size(reduce(BO[4][i],std(id1)))==0) { //--- V(id1) is contained in except. div. i in this chart if(size(reduce(id1,std(BO[4][i])))!=0) { //--- V(id1) does not equal except. div. i of this chart if(defined(Etemp)) {kill Etemp;} ideal Etemp=BO[4][i]; if(npars(basering)>0) { //--- we are in an algebraic extension of the base field if(defined(prtemp)){kill prtemp;} list prtemp=minAssGTZ(BO[4][i]); // C-comp.except.div. if(size(prtemp)>1) { //--- more than 1 component Etemp=ideal(1); for(j=1;j<=size(prtemp);j++) { //--- find correct component if(size(reduce(prtemp[j],std(id1)))==0) { Etemp=prtemp[j]; break; } } if(deg(std(Etemp)[1])==0) { ERROR("fetchInTree:something wrong in field extension"); } } prtemp=delete(prtemp,j); // remove this comp. from list while(size(prtemp)>1) { //--- collect all the others into prtemp[1] prtemp[1]=intersect(prtemp[1],prtemp[size(prtemp)]); prtemp=delete(prtemp,size(prtemp)); } } if(defined(tli)) {kill tli;} //--- determine tli[1] and tli[2] such that //--- V(id1) \cap D(id2) = V(tli[1]) \cap D(tli[2]) \cap BO[4][i] //--- inside V(BO[1]) (and if necessary inside V(BO[1]+BO[2])) if(inJ) { def tli=findTrans(id1+BO[2]+BO[1],Etemp,notE,BO[2]); } else { def tli=findTrans(id1+BO[1],Etemp,notE); } if(npars(basering)>0) { //--- in algebraic extension: make sure we stay outside the other components if(size(prtemp)>0) { for(j=1;j<=ncols(prtemp[1]);j++) { //--- find the (univariate) generator of prtemp[1] which is the remaining //--- factor from the factorization over the extension field if(size(reduce(prtemp[1][j],std(id1)))>0) { tli[2]=tli[2]*prtemp[1][j]; } } } } } else { tli[1]=ideal(0); tli[2]=ideal(1); } id1=tli[1]; id2=id2*tli[2]; notE[size(notE)+1]=BO[4][i]; for(j=1;j<=size(DivL);j++) { if(inIVList(intvec(o1,i),DivL[j])) { nE[size(nE)+1]=j; break; } } if(size(nE)0) { "idlist after current blow down step:"; idlist; } } if(defined(debug_fetchInTree)>0) { "Blowing down ended"; } //---------------------------------------------------------------------------- // Blow up ideal id1 from L[2][comPa] to L[2][m1]. To this end, first // determine the path to follow and save it in path_togo. //---------------------------------------------------------------------------- if(m1==comPa) { //--- no further blow ups needed if(size(algext)==0) { //--- no field extension ==> we are done return(idlist[1][1]); } else { //--- field extension ==> we need to encode the result list retlist; for(m=1;m<=size(idlist);m++) { retlist[m]=string(idlist[m][1]); } return(retlist); } } //--- we need to blow up if(defined(path_m1)) { kill path_m1; } matrix path_m1=imap(Sm1,path); intvec path_togo; for(i=1;i<=ncols(path_m1);i++) { if(path_m1[1,i]>=comPa) { path_togo=path_togo,int(leadcoef(path_m1[1,i])); } } path_togo=path_togo[2..size(path_togo)],m1; i=1; while(i0) { //--- in an algebraic extension of the base field if(defined(T0)){kill T0;} def T0=L[2][path_togo[i+1]]; if(defined(tstr)){kill tstr;} string tstr="ring T=(0,t),(" +varstr(T0)+"),(" +ordstr(T0)+");"; execute(tstr); setring T; execute(algext); minpoly=leadcoef(p); kill tstr; def path=imap(T0,path); def BO=imap(T0,BO); def lastMap=imap(T0,lastMap); if(defined(phi)){kill phi;} map phi=S,lastMap; list idlist=phi(idlist); if(defined(debug_fetchInTree)>0) { "in blowing up (algebraic extension case):"; phi; idlist; } } else { def T=L[2][path_togo[i+1]]; setring T; if(defined(phi)) { kill phi; } map phi=S,lastMap; if(defined(idlist)) {kill idlist;} list idlist=phi(idlist); idlist[1][1]=radical(idlist[1][1]); idlist[1][2]=radical(idlist[1][2]); if(defined(debug_fetchInTree)>0) { "in blowing up (case without field extension):"; phi; idlist; } } for(m=1;m<=size(idlist);m++) { //--- get rid of new exceptional divisor idlist[m][1]=sat(idlist[m][1]+BO[1],BO[4][size(BO[4])])[1]; idlist[m][2]=sat(idlist[m][2],BO[4][size(BO[4])])[1]; } if(defined(debug_fetchInTree)>0) { "after saturation:"; idlist; } if((size(algext)==0)&&(deg(std(idlist[1][1])[1])==0)) { //--- strict transform empty in this chart, it will stay empty till the end setring Sm1; return(ideal(1)); } kill S; i++; } if(defined(debug_fetchInTree)>0) { "End of blowing up steps"; } //--------------------------------------------------------------------------- // prepare results for returning them //--------------------------------------------------------------------------- ideal E,bla; intvec kv; list retlist; for(m=1;m<=size(idlist);m++) { for(j=2;j<=size(idlist[m][3]);j++) { kv=findInIVList(1,path_togo[size(path_togo)],DivL[idlist[m][3][j]]); if(kv!=intvec(0)) { E=E+BO[4][kv[2]]; } } bla=quotient(idlist[m][1]+E,idlist[m][2]); retlist[m]=string(bla); } if(size(algext)==0) { return(bla); } return(retlist); } ///////////////////////////////////////////////////////////////////////////// static proc findInIVList(int pos, int val, list ivl) "Internal procedure - no help and no example available " { //--- find entry with value val at position pos in list of intvecs //--- and return the corresponding entry int i; for(i=1;i<=size(ivl);i++) { if(ivl[i][pos]==val) { return(ivl[i]); } } return(intvec(0)); } ///////////////////////////////////////////////////////////////////////////// //static proc inIVList(intvec iv, list li) "Internal procedure - no help and no example available " { //--- if intvec iv is contained in list li return 1, 0 otherwise int i; int s=size(iv); for(i=1;i<=size(li);i++) { if(typeof(li[i])!="intvec"){ERROR("Not integer vector in the list");} if(s==size(li[i])) { if(iv==li[i]){return(1);} } } return(0); } ////////////////////////////////////////////////////////////////////////////// static proc Vielfachheit(ideal J,ideal I) "Internal procedure - no help and no example available " { //--- auxilliary procedure for addSelfInter //--- compute multiplicity, suitable for the special situation there int d=1; int vd; int c; poly p; ideal Ip,Jp; while((d>0)||(!vd)) { p=randomLast(100)[nvars(basering)]; Ip=std(I+ideal(p)); c++; if(c>20){ERROR("Vielfachheit: Dimension is wrong");} d=dim(Ip); vd=vdim(Ip); } Jp=std(J+ideal(p)); return(vdim(Jp)/vdim(Ip)); } ///////////////////////////////////////////////////////////////////////////// static proc genus_E(list re, list iden0, intvec Eindex) "Internal procedure - no help and no example available " { int i,ge,gel,num; def R=basering; ring Rhelp=0,@t,dp; def S=re[2][Eindex[1]]; setring S; def Sh=S+Rhelp; //---------------------------------------------------------------------------- //--- The Q-component X is reducible over C, decomposes into s=num components //--- X_i, we assume they have n.c. //--- s*g(X_i)=g(X)+s-1. //---------------------------------------------------------------------------- if(defined(I2)){kill I2;} ideal I2=dcE[Eindex[2]][Eindex[3]][1]; num=ncols(dcE[Eindex[2]][Eindex[3]][4]); setring Sh; if(defined(I2)){kill I2;} ideal I2=imap(S,I2); I2=homog(I2,@t); ge=genus(I2); gel=(ge+(num-1))/num; if(gel*num-ge-num+1!=0){ERROR("genus_E: not divisible by num");} setring R; return(gel,num); }