//////////////////////////////////////////////////////////////////////////// version="$Id$"; category="Resolution of singularities"; info=" LIBRARY: resbinomial.lib Combinatorial algorithm of resolution of singularities of binomial ideals in arbitrary characteristic. Binomial resolution algorithm of Blanco AUTHORS: R. Blanco, mariarocio.blanco@uclm.es, @* G. Pfister, pfister@mathematik.uni-kl.de PROCEDURES: BINresol(J); computes a E-resolution of singularities of (J) (THE SECOND PART IS NOT IMPLEMENTED YET) Eresol(J); computes a E-resolution of singularities of (J) in char 0 determinecenter(L1,L2,c,n,Y,a,mb,flag,control3); computes the next blowing-up center Blowupcenter(L1,id,m,L2,c,n,h); makes the blowing-up Nonhyp(Coef,expJ,sJ,n,flag,sums); computes the ideal generated by the non hyperbolic generators of expJ identifyvar(); identifies status of variables Edatalist(Coef,Exp,k,n,flag); gives the E-order of each term in Exp EOrdlist(Coef,Exp,k,n,flag); computes the E-order of an ideal (giving in the language of lists) maxEord(Coef,Exp,k,n,flag); computes de maximum E-order of an ideal given by Coef and Exp ECoef(Coef,expP,sP,V,auxc,n,flag); Computes a simplified version of the E-Coeff ideal. The E-orders are correct, but tranformations of coefficients of the generators and powers of binomials cannot be computed easily in terms of lists. elimrep(L); removes repeated terms from a list Emaxcont(Coef,Exp,k,n,flag); computes a list of hypersurfaces of E-maximal contact cleanunit(mon,n,flag); clean the units in a monomial mon resfunction(t,auxinv,nchart,n); composes the E-resolution function calculateI(Coef,J,c,n,Y,a,b,D); computes the order of the non monomial part of an ideal J Maxord(L,n); computes the maximum exponent of an exceptional monomial ideal Gamma(L,c,n); computes the Gamma function for an exceptional monomial ideal given by L convertdata(C,L,n,flag); computes the ideal corresponding to C,L lcmofall(nchart,mobile); computes the lcm of the denominators of the E-orders for all the charts computemcm(Eolist); computes the lcm of the denominators of the E-orders for one chart constructH(Hhist,n,flag); construct the list of exceptional divisors accumulated at this chart constructblwup(blwhist,n,chy,flag); construct the ideal defining the map K[W] --> K[Wi], which gives the composition map of all the blowing up leading to this chart constructlastblwup(blwhist,n,chy,flag); construct the ideal defining the last blowup leading to this chart genoutput(chart,mobile,nchart,nsons,n,q,p); generates the output for visualization salida(idchart,chart,mobile,numson,previousa,n,q); generates the output for one chart iniD(n); creates a list of lists of zeros of size n sumlist(L1,L2); sums two lists component to component reslist(L1,L2); subtracts two lists component to component multiplylist(L,a); multiplies a list by a number, component to component dividelist(L1,L2); divides two lists component to component createlist(L1,L2); creates a list of lists of two elements "; // inidata(K,k); verifies input data, a binomial ideal K of k generators // data(K,k,n); transforms data on lists of length n // list0(n); creates a list of zeros of size n LIB "general.lib"; LIB "qhmoduli.lib"; LIB "inout.lib"; LIB "poly.lib"; LIB "resolve.lib"; LIB "reszeta.lib"; LIB "resgraph.lib"; //////////////////////////////////////////////////////////////////////////// static proc inidata(ideal K,int k) "USAGE: inidata(K,k); K any ideal, k integer (!=0) COMPUTE: Verifies the input data RETURN: flag indicating if the ideal is binomial or not EXAMPLE: example inidata; shows an example " { int i; for (i=1;i<=k; i++) { if (size(K[i])>2){return(0);} } return(1); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; ideal J1=x(1)^4*x(2)^2, x(1)^2+x(3)^3; inidata(J1,2); ideal J2=x(1)^4*x(2)^2, x(1)^2+x(2)^3+x(3)^5; inidata(J2,2); } ///////////////////////////////////////////////////////////////////////////////// proc changeoriginalvar() "USAGE: changeoriginalvar(); COMPUTE: Change the name of the variables to x(1...n), only necessary at the beginning RETURN: the new ring with the suitable names EXAMPLE: example changeoriginalvar; shows an example " { int i,n,cont; n=nvars(basering); cont=0; def r=basering; // check the name of the variables for (i=1;i<=n; i++){if (varstr(i)[1]=="x" or varstr(i)[1]=="y"){cont=cont+1;}} // change them if there exists some variable different from x(i) or y(i) if (cont!=n or n<=2){ // making the change def Rnew=changevar ("x()"); setring Rnew; // print("INVERTIBLE VARIABLES NOT CONSIDERED AT THE BEGINNING"); return(Rnew,1); } else{ // print("INVERTIBLE VARIABLES ALREADY CONSIDERED AT THE BEGINNING"); return(r,0); } } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; changeoriginalvar(); ring r = 0,(x,y,z,w),dp; changeoriginalvar(); } ///////////////////////////////////////////////////////////////////////////////// proc identifyvar() "USAGE: identifyvar(); COMPUTE: Asign 0 to variables x and 1 to variables y, only necessary at the beginning RETURN: list, say l, of size the dimension of the basering l[i] is: 0 if the i-th variable is x(i), 1 if the i-th variable is y(i) EXAMPLE: example identifyvar; shows an example " { int i,n; list flaglist; n=nvars(basering); flaglist=list0(n); for (i=1;i<=n; i++){if (varstr(i)[1]=="y"){flaglist[i]=1;}} return(flaglist); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; identifyvar(); } ///////////////////////////////////////////////////////////////////////////////// proc data(ideal K,int k,int n) "USAGE: data(K,k,n); K any ideal, k integer (!=0), n integer (!=0) COMPUTE: Construcs a list with the coefficients and exponents of one ideal RETURN: lists of coefficients and exponents of K EXAMPLE: example data; shows an example " {int i,j,lon; number aa; intvec cc; list bb,dd,aux,ddaux,Coef,Exp; for (i=1;i<=k; i++) { lon=size(K[i]); // binomial if (lon==2){aa=leadcoef(K[i][1]); bb=aa; Coef[i]=bb; // coefficients cc=leadexp(K[i][1]); // exponents // cc is an intvec, transform cc in dd, a list of lists dd=cc[1..n]; aux[1]=dd; // the same for the second term aa=leadcoef(K[i][2]); bb=aa; Coef[i]=Coef[i] + bb; // all the coefficients of i-th generator of K cc=leadexp(K[i][2]); dd=cc[1..n]; aux[2]=dd; Exp[i]=aux;} // monomial if (lon==1){aux=list(); aa=leadcoef(K[i][1]); bb=aa; Coef[i]=bb; cc=leadexp(K[i][1]); dd=cc[1..n]; aux[1]=dd; Exp[i]=aux;} } //end for return(Coef,Exp); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3; data(J,2,3); } ////////////////////////////////////////////////////// proc Edatalist(list Coef,list Exp,int k,int n,list flaglist) "USAGE: Edatalist(Coef,Exp,k,n,flaglist); Coef,Exp,flaglist lists, k,n, integers Exp is a list of lists of exponents, k=size(Exp) COMPUTE: computes a list with the E-order of each term RETURN: a list with the E-order of each term EXAMPLE: example Edatalist; shows an example " {int i,j,lon,mm; list dd,ss,sums; number aux,aux1,aux2; for (i=1;i<=k;i++){lon=size(Coef[i]); if (lon==1) { for (j=1;j<=n;j++){if (flaglist[j]==0){aux=aux+Exp[i][1][j];}} ss=aux; aux=0;} // monomial else { for (j=1;j<=n;j++){if (flaglist[j]==0){ aux1=aux1+Exp[i][1][j]; aux2=aux2+Exp[i][2][j];}} ss=aux1,aux2; aux1=0; aux2=0; } // binomial sums[i]=ss;} return(sums); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; list flag=identifyvar(); ideal J=x(1)^3*x(3)-y(2)*y(4)^2,x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5); list L=data(J,3,8); list EL=Edatalist(L[1],L[2],3,8,flag); EL; // E-order of each term ring r = 2,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; list flag=identifyvar(); ideal J=x(1)^3*x(3)-y(2)*y(4)^2,x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5); list L=data(J,3,8); list EL=Edatalist(L[1],L[2],3,8,flag); EL; // E-order of each term IN CHAR 2, COMPUTATIONS NEED TO BE DONE IN CHAR 0 ring r = 0,(x(1..3)),dp; list flag=identifyvar(); ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3; list L=data(J,2,3); list EL=Edatalist(L[1],L[2],2,3,flag); EL; // E-order of each term } /////////////////////////////////////////////////////////////////////////////////// proc EOrdlist(list Coef,list Exp,int k,int n,list flaglist) "USAGE: EOrdlist(Coef,Exp,k,n,flaglist); Coef,Exp,flaglist lists, k,n, integers Exp is a list of lists of exponents, k=size(Exp) COMPUTE: computes de E-order of an ideal given by a list (Coef,Exp) and extra information RETURN: maximal E-order, and its position=number of generator and term EXAMPLE: example EOrdlist; shows an example " {int i,can,canpost,lon; number canmin; list sums; sums=Edatalist(Coef,Exp,k,n,flaglist); canmin=sums[1][1]; // inicializating, works also with a monomial for (i=1;i<=k; i++){lon=size(sums[i]); // this is 2 for binomial and 1 for monomial generators if (sums[i][1]<=canmin and Coef[i][1]!=0){canmin=sums[i][1]; can=i; canpost=1;} // if the generator is a binomial we check the second term if (lon==2) {if (sums[i][2]0){ for (i=1;i<=k; i++){lon=size(sums[i]); if (lon==2){if (sums[i][1]==maxEo) // variables of the first term {for (j=1;j<=n; j++){if(Exp[i][1][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}}} if (sums[i][2]==maxEo) // variables of the second term {for (j=1;j<=n; j++){if(Exp[i][2][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}}}} else {if (sums[i][1]==maxEo) {for (j=1;j<=n; j++){if(Exp[i][1][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}}}} }} else {maxvar=list();} // eliminating repeated terms maxvar=elimrep(maxvar); // It is necessary to check if flag[j]==0 in order to avoid the selection of y variables return(maxEo,maxvar); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; list flag=identifyvar(); ideal J=x(1)^3*x(3)-y(2)*y(4)^2,x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5),x(7)^4*y(8)^2; list L=data(J,4,8); list hyp=Emaxcont(L[1],L[2],4,8,flag); hyp[1]; // max E-order=0 hyp[2]; // There are no hypersurfaces of E-maximal contact ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; list flag=identifyvar(); ideal J=x(1)^3*x(3)-y(2)*y(4)^2*x(3),x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5),x(7)^4*y(8)^2; list L=data(J,4,8); list hyp=Emaxcont(L[1],L[2],4,8,flag); hyp[1]; // the E-order is 1 hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of E-maximal contact } /////////////////////////////////////////////////////// proc cleanunit(list mon,int n,list flaglist) "USAGE: cleanunit(mon,n,flaglist); mon, flaglist lists, n integer COMPUTE: We clean (or forget) the units in a monomial, given by "y" variables RETURN: The list defining the monomial ideal already cleaned EXAMPLE: example cleanunit; shows an example " { int i; for (i=1;i<=n;i++){if (flaglist[i]==1){mon[i]=0;}} // coef[1]=coef[1]*y(i)^mon[i]; IS NOT ALLOWED because mon[i] can be a number // therefore, the coefficients remain constant return(mon); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1),y(2),x(3),y(4)),dp; list flag=identifyvar(); ideal J=x(1)^3*y(2)*x(3)^5*y(4)^8; list L=data(J,1,4); L[2][1][1]; // list of exponents of the monomial J list M=cleanunit(L[2][1][1],4,flag); M; // new list without units } ////////////////////////////////////////////////////// // Classification of the ideal E-Coeff_V(P): // ccase=1, E-Coeff_V(P)=0 // 2,3 Bold regular case // 4 P=1 monomial case (detected before) // 0 Otherwise proc ECoef(list Coef,list expP,int sP,int V,number auxc,int n,list flaglist) "USAGE: ECoef(Coef,expP,sP,V,auxc,n,flaglist); Coef, expP, flaglist lists, sP, V, n integers, auxc number COMPUTE: The ideal E-Coeff_V(P), where V is a permissible hypersurface which belongs to the center RETURN: list of exponents, list of coefficients and classification of the ideal E-Coeff_V(P) EXAMPLE: example ECoef; shows an example " { int i,j,k,l,numg,ccase,cont2,cont3,val; number aa; list Eco,newcoef,auxexp,newL,rs,rs2,aux,aux2,aux3,aux4,L; auxexp=expP; l=1; for (i=1;i<=sP;i++) {rs[i]=size(Coef[i]); if (rs[i]==2){ // binomials if (auxexp[i][1][V]!=auxexp[i][2][V]) // no common factors for the variable in V {for (j=1;j<=2;j++){if (auxexp[i][j][V]1){ccase=3;} // CASE P=x^{\alpha},x^{\beta}, IN FACT, BOLD REGULAR if (cont2==val and cont3==0){ccase=4;} // P=1, then I=1 monomial case // Case BOLD REGULAR P=x^{\alpha}*(1-\mu y^{\delta}) // IT IS NON NECESSARY TO CHECK IT, Eco=empty list, already done! L=maxEord(newcoef,Eco,numg,n,flaglist); // L[1] is the E-order of Eco if (L[1]==0){ccase=2; print("E-order zero!");} // BOLD REGULAR CASE // we leave it to check the computations } // close else return(Eco,newcoef,ccase); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp; list flag=identifyvar(); ideal P=x(1)^2*x(3)^5-x(5)^7*y(4),x(6)^3*y(2)^5-x(7)^5,x(5)^3*x(6)-y(4)^3*x(1)^5; list L=data(P,3,7); list L2=ECoef(L[1],L[2],3,1,3,7,flag); L2[1]; // exponents of the E-Coefficient ideal respect to x(1) L2[2]; // its coefficients L2[3]; // classify the type of ideal obtained ring r = 0,(x(1),y(2),x(3),y(4)),dp; list flag=identifyvar(); ideal J=x(1)^3*(1-2*y(2)*y(4)^2); // Bold regular case list L=data(J,1,4); list L2=ECoef(L[1],L[2],1,1,3,4,flag); L2; ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp; list flag=identifyvar(); ideal J=x(1)^3-x(3)^2*y(4)^2,x(1)*x(7)*y(2)-x(6)^3*x(5)*y(4)^3,x(5)^3-x(5)^3*y(2)^2; list L=data(J,3,7); list L2=ECoef(L[1],L[2],3,1,2,7,flag); L2; ring r = 3,(x(1),y(2),x(3),y(4),x(5..7)),dp; list flag=identifyvar(); ideal J=x(1)^3-x(3)^2*y(4)^2,x(1)*x(7)*y(2)-x(6)^3*x(5)*y(4)^3,x(5)^3-x(5)^3*y(2)^2; list L=data(J,3,7); list L2=ECoef(L[1],L[2],3,1,2,7,flag); L2; // THE COMPUTATIONS ARE NOT CORRECT IN CHARACTERISTIC p>0 // because numbers are treated as 0 in assignments } //////////////////////////////////////////////////////////////////////////// // The intvec a indicates the previous center // Hhist = intvec of exceptional divisors of the parent chart proc determinecenter(list Coef,list expJ,number c,int n,int Y,intvec a,list listmb,list flag,int control3,intvec Hhist) "USAGE: determinecenter(Coef,expJ,c,n,Y,a,listmb,flag,control3,Hhist); Coef, expJ, listmb, flag lists, c number, n, Y, control3 integers, a, Hhist intvec COMPUTE: next center of blowing up and related information, see example RETURN: several lists defining the center and related information EXAMPLE: example determinecenter; shows an example " {int i,j,rstep,l,mm,cont,cont1,cont2,cont3,a4,sI,sP,V,V2,ccase,b,Mindx,tip,mval; number auxc,a1,a2,ex,maxEo,aux; list D,H,auxJ; // lists of D_n,D_n-1,...,D_1; H_n,H_n-1,...,H_1; J_n,J_n-1,...,J_1 list oldOlist,oldC,oldt,oldD,oldH,allH; // information of the previous step list Olist,C,t,Dstar,center,expI,expP,newJ,maxset; list maxvar,auxlist,aux3,auxD,auxolist,auxdiv,auxaux,L,rs,auxgamma,auxg2,aux1; // auxiliary lists list auxinvlist,newcoef,EL,Ecoaux,Hplus,transH,Hsum,auxset,sumnewH; // auxiliary lists list auxcoefI,auxcent,center2; intvec oldinfobo7,infobo7; int infaux,leh,leh2,leh3; tip=listmb[1]; // It is not used in this procedure, it is used to compute the lcm of the denominators oldOlist=listmb[2]; oldC=listmb[3]; oldt=listmb[4]; // t= resolution function oldD=listmb[5]; oldH=listmb[6]; allH=listmb[7]; oldinfobo7=listmb[8]; // auxiliary intvec, it is used to define BO[7] // inicializating lists Olist=list(); C=list(); auxinvlist=list(); auxJ[1]=expJ; rstep=n; // we are in dimension rstep auxc=c; cont=1; if (Y==0) {D=iniD(n); H=iniD(n); infobo7=-1;} // first center, inicializate previous information if (Y!=0 and rstep==n) // In dimension n, D'_n is always of this form { auxdiv=list0(n); Dstar[1]=oldD[1]; b=size(a); for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}} Dstar[1][Y]=aux; aux=0; auxdiv[Y]=oldOlist[1]-oldC[1]; D[1]=sumlist(Dstar[1],auxdiv);} // list defining D_n // computing strict transforms of the exceptional divisors H if (Y!=0){transH=iniD(n); for (i=1;i<=size(oldH);i++){transH[i]=oldH[i]; transH[i][Y]=0;} // Note: size(oldH)<=n allH[Y]=1;} // transform of |H|=H_nU...UH_1 // We put here size(oldH) instead of n because maybe we have not // calculated all the dimensions in the previous step // STARTING THE LOOP while (rstep>=1) { if (Y!=0 and rstep!=n) // transformation law of D_i for ia4){cont=0; oldt[n-rstep+1]=0; } // VERIFICAR!!!! if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1] and control3==0){H[n-rstep+1]=transH[n-rstep+1]; // we fill now the value for BO[7] if (oldinfobo7[n-rstep+1]==-1){leh=size(Hhist); infobo7[n-rstep+1]=Hhist[leh];} // suitable index !!! else{ infaux=oldinfobo7[n-rstep+1]; infobo7[n-rstep+1]=infaux;} // the same as the previous step } else { if (rstep=auxc) {expP=expI; Mindx=0;} // The coefficients also remain constant else {ex=maxEo/(auxc-maxEo); auxlist=list(); Mindx=1; auxlist[1]=multiplylist(D[n-rstep+1],ex); // weighted monomial part: D[n-rstep+1]^ex; expP=insert(expI,auxlist); // P=I+D[n-rstep+1]^ex; auxcoefI=Coef; Coef=insert(Coef,list(1));} // Adding the coefficient for M // NOTE: IT IS NECESSARY TO ADD COEFFICIENT 1 TO THE MONOMIAL PART M // E-ord(P_i)=E-ord(I_i) so to compute the E-order of P_i we can compute E-ord(I_i) // Calculate variables of E-maximal contact, ALWAYS WITH RESPECT TO THE IDEAL I !! sP=size(expP); // Can be different from size(expI) if (Mindx==1){ maxvar=Emaxcont(auxcoefI,expI,sI,n,flag);} else{ maxvar=Emaxcont(Coef,expP,sP,n,flag);} auxc=maxvar[1]; // E-order of P, critical value for the next step, ALSO VALID auxc=maxEo; if (auxc!=maxEo){print("ERROR, the E-order is not well computed");} maxset=maxvar[2]; // center=center + maxset; // HACER DESPUES Y A?ADIR SOLO V!!!!!! // Cleaning the center: eliminating repeated variables // center=elimrep(center); // if (rstep==1) {break;} // Induction finished, is not necessary to compute the rest // Calculate Hplus=set of non permissible hypersurfaces // RESET Hplus if c has dropped or we have eliminated hyperbolic generators // ES NECESARIO PONER CONDICION DE SI INVARIANTE BAJA O NO??? SI BAJA HPLUS NO SE USA... if (Y==0 or c gives (max E-order,sums) proc Nonhyp(list Coef,list expJ,int sJ,int n,list flaglist,list sums) "USAGE: Nonhyp(Coef,expJ,sJ,n,flaglist,sums); Coef, expJ, flaglist, sums lists, sJ, n integers COMPUTE: The "ideal" generated by the non hyperbolic generators of J RETURN: lists with the following information newcoef,newJ: coefficients and exponents of the non hyperbolic generators totalhyp,totalgen: coefficients and exponents of the hyperbolic generators flaglist: new list saying status of variables NOTE: the basering r is supposed to be a polynomial ring K[x,y], in fact, we work in a localization of K[x,y], of type K[x,y]_y with y invertible variables. EXAMPLE: example Nonhyp; shows an example " { int i,j,k,h,lon,lon2,cont; number eordcontrol; list genhyp,listgen,listid,posnumJ,newJ,newcoef,hypcoef,hyp,aux1,aux2,aux3,aux,midlist; list totalhyp,totalgen; eordcontrol=0; while (eordcontrol==0 and sJ!=0) { // Give a positional number/flag to each generator of expJ for (i=1;i<=sJ; i++){listgen=expJ[i]; listid=i; posnumJ[i]=listgen+listid; } // Select the non hyperbolic and hyperbolic generators for (j=1;j<=sJ; j++){lon=size(Coef[j]); if (lon==1){ // IS NOT NECESSARY TO CHECK IF THERE EXIST A MONOMIAL WITH ONLY UNITS, ALREADY DONE!! aux1=aux1+posnumJ[j]; aux3=list(); aux3[1]=expJ[j]; newJ=newJ+aux3; aux3[1]=Coef[j]; newcoef=newcoef+aux3; } else{ // CHECKING BINOMIALS, ONE TERM WITH E-ORDER ZERO GIVES HYPERBOLIC EQ if (sums[j][1]==0 or sums[j][2]==0){aux2=aux2+posnumJ[j]; aux3=list(); aux3[1]=expJ[j]; genhyp=genhyp+aux3; aux3[1]=Coef[j]; hypcoef=hypcoef+aux3; if (sums[j][1]==0){aux3[1]=expJ[j][2]; hyp=hyp+aux3;} if (sums[j][2]==0){aux3[1]=expJ[j][1]; hyp=hyp+aux3;} } else {aux1=aux1+posnumJ[j]; aux3=list(); aux3[1]=expJ[j]; newJ=newJ+aux3; aux3[1]=Coef[j]; newcoef=newcoef+aux3;} } } // NOTE: aux1 and aux2 are no needed right now! // Identify new y variables, that is, x variables in the monomials contained in hyp h=size(hyp); for (k=1;k<=h; k++){ for(i=1;i<=n; i++){ if (hyp[k][i]!=0 and flaglist[i]==0) {flaglist[i]=1;}}} // To replace x by y IT IS NECESSARY TO CHANGE THE BASERING!!! We change only the list flaglist // CHECK IF THE IDEAL IS ALREADY GENERATED BY MONOMIALS, in this case // WE HAVE FINISHED THE E-RESOLUTION PART, J GENERATED BY MONOMIALS AND HYPERBOLIC EQS cont=0; lon2=size(newJ); for (j=1;j<=lon2; j++){if (size(newJ[j])==1){cont=cont+1;}} if (cont==lon2){newcoef=list(); newJ=list(); totalgen=totalgen+genhyp; totalhyp=totalhyp+hypcoef; break;} // CHECK IF THERE ARE MORE HYPERBOLIC EQUATIONS AFTER UPDATE THE FLAG LIST // CHECK THE MAXIMAL E-ORDER AGAIN if (lon2==0){ // we are in the previous case, newJ=empty list, save values and exit totalgen=totalgen+genhyp; totalhyp=totalhyp+hypcoef; break; } midlist=maxEord(newcoef,newJ,lon2,n,flaglist); eordcontrol=midlist[1]; if (eordcontrol==0){ // new input for the loop Coef=newcoef; expJ=newJ; sJ=lon2; sums=midlist[2]; // flaglist is already updated totalgen=totalgen+genhyp; totalhyp=totalhyp+hypcoef; hypcoef=list(); genhyp=list(); newJ=list(); newcoef=list(); } else{ // If the process is already finished we save the values and exit totalgen=totalgen+genhyp; totalhyp=totalhyp+hypcoef; } } // closing while return(newcoef,newJ,totalhyp,totalgen,flaglist); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp; list flag=identifyvar(); // List giving flag=1 to invertible variables: y(2),y(4) ideal J=x(1)^3-x(3)^2*y(4)^2,x(1)*x(7)*y(2)-x(6)^3*x(5)*y(4)^3,1-x(5)^2*y(2)^2; list L=data(J,3,7); list L2=maxEord(L[1],L[2],3,7,flag); L2[1]; // Maximum E-order list New=Nonhyp(L[1],L[2],3,7,flag,L2[2]); New[1]; // Coefficients of the non hyperbolic part New[2]; // Exponents of the non hyperbolic part New[3]; // Coefficients of the hyperbolic part New[4]; // New hyperbolic equations New[5]; // New list giving flag=1 to invertible variables: y(2),y(4),y(5) ring r = 0,(x(1..4)),dp; list flag=identifyvar(); ideal J=1-x(1)^5*x(2)^2*x(3)^5, x(1)^2*x(3)^3+x(1)^4*x(4)^6; list L=data(J,2,4); list L2=maxEord(L[1],L[2],2,4,flag); L2[1]; // Maximum E-order list New=Nonhyp(L[1],L[2],2,4,flag,L2[2]); New; } ////////////////////////////////////////////////////////////// proc calculateI(list Coef,list expJ,number c,int n,int Y,intvec a,number oldordI,list oldD) "USAGE: calculateI(Coef,expJ,c,n,Y,a,b,D); Coef, expJ, D lists, c, b numbers, n,Y integers, a intvec RETURN: ideal I, non monomial part of J EXAMPLE: example calculateI; shows an example " { int i,cont1,b,j; number EordI,aux; list D,L,expI; list auxdiv,Dstar,aux1,rs; // WE NEED THE MONOMIAL PART, BUT ONLY IN DIMENSION n auxdiv=list0(n); auxdiv[Y]=oldordI-c; Dstar[1]=oldD[1]; b=size(a); for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}} Dstar[1][Y]=aux; aux=0; D[1]=sumlist(Dstar[1],auxdiv); cont1=0; for (i=1;i<=n;i++) {if (D[1][i]==0) {cont1=cont1+1;}} // if it fails write listO(n)[i] if (cont1==n) {expI=expJ;} else { for (i=1;i<=size(expJ);i++) {rs[i]=size(Coef[i]); if (rs[i]==2){ aux1=list(); aux1[1]=reslist(expJ[i][1],D[1]); aux1[2]=reslist(expJ[i][2],D[1]); expI[i]=aux1;} // binomial else {aux1=list(); aux1[1]=reslist(expJ[i][1],D[1]); expI[i]=aux1;}} // monomial } return(expI); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; list flag=identifyvar(); ideal J=x(1)^4*x(2)^2, x(3)^3; list Lmb=1,list(0,0,0),list(0,0,0),list(3),iniD(3),iniD(3),list(0,0,0),-1; list L=data(J,2,3); list LL=determinecenter(L[1],L[2],3,3,0,0,Lmb,flag,0,-1); // Calculate the center module auxpath=[0,-1]; list infochart=0,0,0,L[2],L[1],flag,0,list(0,0,0),auxpath,list(),list(); list L3=Blowupcenter(LL[1],1,1,infochart,3,3,0); // blowing-up and looking to the x(3) chart calculateI(L3[2][1][5],L3[2][1][4],3,3,3,L3[2][1][3],3,iniD(3)); // (I_3) // looking to the x(1) chart calculateI(L3[2][2][5],L3[2][2][4],3,3,1,L3[2][2][3],3,iniD(3)); // (I_3) } ////////////////////////////////////////////////////////////////////////////////////// // // // E-RESOLUTION: Eresol(J) subroutine computing the E-resolution of J, char 0 // // // ////////////////////////////////////////////////////////////////////////////////////// proc Eresol(ideal J) "USAGE: Eresol(J); J ideal RETURN: The E-resolution of singularities of J in terms of the affine charts, see example EXAMPLE: example Eresol; shows an example " {int i,n,k,idchart,nchart,parent,Y,oldid,tnum,s,cont,control,control2,control3,cont2,val,rs2,l,cont3,tip; intvec a,Hhist; number c,EordJ,EordI,oldordI; list L,LL,oldD,t,auxL,finalchart,chart,auxchart,newL,auxp,auxfchart,L2; list Coef,expJ,expI,sons,oldOlist,oldC,oldt,oldH,allH,auxordJ,auxordI,auxmb,mobile,invariant; list step,nsons,auxinv,extraL,totalinv,auxsum; string empstring; module auxpath; // ADDED LATER list flag,newflag,blwhist,hipercoef,hiperexp,hipercoefson,hiperexpson; intvec infobo7; export finalchart; // export nsons; // export tnum; // export nchart; // export step; export invariant; export auxinv; export mobile; n=nvars(basering); flag=identifyvar(); k=size(J); // Checking input data if (inidata(J,k)==0){return("This library only works for binomial ideals.");} idchart=1; nchart=1; parent=0; step=0; control=0; control2=0; control3=0; // Translate the input ideal to a list auxL=data(J,k,n); // data gives (Coef,Exp) // THEREAFTER WE WORK ALL THE TIME WITH LISTS L=maxEord(auxL[1],auxL[2],k,n,flag); // gives (max E-ord, sums) EordJ=L[1]; // before the first blow up I=J EordI=EordJ; // main loop AT EACH CHART WE MUST INICIALIZATE ALL THE VALUES AND // CONSTRUCT THE FIRST CHART chart[1] BEFORE THE LOOP // at the first step, before the blow up, there are no exceptional divisors, Y=0 Y=0; expJ=auxL[2]; Coef=auxL[1]; Hhist=0; blwhist=list0(n); auxpath=[0,-1]; hipercoef=list(); // this is for the first chart hiperexp=list(); auxp=parent,Y,a,expJ,Coef,flag,Hhist,blwhist,auxpath,hipercoef,hiperexp; chart[1]=auxp; // information of the first chart tip=1; oldOlist=list0(n); oldC=list0(n); oldC[1]=EordJ; // non necessary here c=EordJ; // the value c is given by the previous step oldt=list0(n); oldD=iniD(n); oldH=iniD(n); allH=list0(n); for (i=1;i<=n;i++){infobo7[i]=-1;} auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; mobile[1]=auxmb; // mobile corresponding to the first chart auxinv[1]=list(0); // NOTE: oldC[1] is the value c to classify the chart in one of the next cases // HERE BEGIN THE LOOP while (idchart<=nchart) // WE PROCEED WHILE THERE EXIST UNSOLVED CHARTS { if (idchart!=1) // WE ARE NOT IN THE FIRST CHART, INICIALIZATE ALL THE VALUES { parent=chart[idchart][1]; Y=chart[idchart][2]; a=chart[idchart][3]; expJ=chart[idchart][4]; Coef=chart[idchart][5]; flag=chart[idchart][6]; Hhist=chart[idchart][7]; // it is not necessary for the computations blwhist=chart[idchart][8]; auxpath=chart[idchart][9]; hipercoef=chart[idchart][10]; hiperexp=chart[idchart][11]; k=size(Coef); // IT IS NECESSARY TO COMPUTE IT BECAUSE IT DECREASES IF THERE ARE HYPERBOLIC EQS auxordJ=maxEord(Coef,expJ,k,n,flag); EordJ=auxordJ[1]; if (control==0){c=mobile[parent+1][3][1];} // we keep c from the last step else {c=EordJ; control=0; } // we reset the value of c if (control2==1){c=EordJ; control2=0; control3=1;} // we reset the value of c // NOTE: oldC[1] is the value c to classify the chart in one of the next cases } // The E-order must be computed here oldid=idchart; if (EordJ<0) {print("ERROR in J in chart"); print(idchart); ~; break;} //------------------------------------------------------------- // CASE J=1, if we reset c, can happen Eord=c=0 // or if there are hyperbolic equations at the beginning!!! A?ADIR!!!! // if (EordJ==0){auxfchart[1]=chart[idchart]; // WE HAVE FINISHED // finalchart=finalchart+auxfchart; // empstring="#"; print("reset c and Eord=c=0"); print(idchart); // invariant[idchart]=empstring; // auxinv[idchart]=list(0); // nsons[idchart]=0; // idchart=idchart+1;} //---------------------------------------------------------------------- if (EordJ>=c and EordJ!=0) // subroutine: E-RESOLUTION OF PAIRS { if (parent>0) { LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,chart[parent][7]); } else { LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,Hhist); } // determinecenter gives (center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7) // save current values, before the blow up oldOlist=LL[3]; tip=computemcm(oldOlist); oldC=LL[4]; oldt=LL[5]; oldD=LL[6]; oldH=LL[7]; allH=LL[8]; auxinv[idchart]=LL[9]; infobo7=LL[10]; auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; mobile[idchart+1]=auxmb; invariant[idchart]=oldt; newL=Blowupcenter(LL[1],idchart,nchart,chart[idchart],c,n,step[idchart]); // Blowupcenter gives (info,auxchart,nchart,auxstep,num) // IMPORTANT: ADD THE NEW CHARTS AFTER EACH BLOW UP, IN ORDER TO KEEP THEM CORRECTLY step=step+newL[4]; nsons[idchart]=newL[5]; chart=chart+newL[2]; finalchart=finalchart+newL[1]; // new input for the loop idchart=idchart+1; nchart=newL[3]; control3=0; } // END OF CASE EordJ>=c //--------------------------------------------------------------------- else{ // compute EordI=max E-order(I) expI=calculateI(Coef,expJ,c,n,Y,a,mobile[parent+1][2][1],mobile[parent+1][5]); k=size(expJ); // probably non necessary auxordI=maxEord(Coef,expI,k,n,flag); EordI=auxordI[1]; auxsum=auxordI[2]; // CASE EordI>0 DROP c AND CONTINUE if (EordI>0){idchart=idchart; // keep the chart and back to the main loop while, dropping the value of c control=1;} else{ // EordI=0, so check if I=1 or not cont2=0; // If cont2=val then all the entries of expI are zero!! val=0; for (i=1;i<=k;i++) {rs2=size(Coef[i]); if (rs2==1){if (auxsum[i][1]==0){cont2=val; break;} // THERE EXIST A MONOMIAL WITH ONLY UNITS val=val+n; // monomials for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;}} } else{val=val+(2*n); // binomials for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;} if (expI[i][2][l]==0) {cont2=cont2+1;}} } } // CASE EordI==0 AND I=1 THIS CHART IS DONE, FINISH // NOTE: THIS CASE IS NOT MONOMIAL BECAUSE E-Sing(J,c) is empty if (cont2==val){auxfchart[1]=chart[idchart]; finalchart=finalchart+auxfchart; empstring="#"; invariant[idchart]=empstring; auxinv[idchart]=list(0); nsons[idchart]=0; // information for the mobile tip=1; oldOlist=list(0); oldC=list(0); oldt=list(0); oldD=list(0); oldH=list(0); allH=list(0); // the value of the parent + the new one infobo7=-1; auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; mobile[idchart+1]=auxmb; idchart=idchart+1;} else{ // CASE EordI==0 AND I!=1 --> HYPERBOLIC EQUATIONS // COMPUTE THE IDEAL OF NON HYPERBOLIC ELEMENTS extraL=Nonhyp(Coef,expI,k,n,flag,auxordI[2]); // gives (newcoef,newI,hypcoef,genhyp,flaglist) // CHECK IF ALL THE VARIABLES ARE ALREADY INVERTIBLE newflag=extraL[5]; chart[idchart][6]=extraL[5]; // update the status of variables cont3=0; for (i=1;i<=n;i++){if (newflag[i]==1){cont3=cont3+1;}} if (cont3==n){ // ALL THE VARIABLES ARE INVERTIBLE auxfchart[1]=chart[idchart]; finalchart=finalchart+auxfchart; empstring="@"; invariant[idchart]=empstring; auxinv[idchart]=list(0); nsons[idchart]=0; // information for the mobile tip=1; oldOlist=list(0); oldC=list(0); oldt=list(0); oldD=list(0); oldH=list(0); allH=list(0); infobo7=-1; auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; mobile[idchart+1]=auxmb; idchart=idchart+1;} else{ // OTHERWISE, CONTINUE CHEKING IF newI=0 or not Coef=extraL[1]; expI=extraL[2]; hipercoefson=extraL[3]; // Information about hyperbolic generators hiperexpson=extraL[4]; k=size(expI); if (k==0){auxfchart[1]=chart[idchart]; // WE HAVE FINISHED finalchart=finalchart+auxfchart; empstring="#"; // no more non-hyperbolic generators in this chart invariant[idchart]=empstring; auxinv[idchart]=list(0); nsons[idchart]=0; // information for the mobile tip=1; oldOlist=list(0); oldC=list(0); oldt=list(0); oldD=list(0); oldH=list(0); allH=list(0); infobo7=-1; auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; mobile[idchart+1]=auxmb; idchart=idchart+1;} else{ // CONTINUE WITH THE IDEAL OF NON HYPERBOLIC EQS chart[idchart][4]=expI; // new input ideal and coefficients chart[idchart][5]=Coef; chart[idchart][10]=hipercoef+hipercoefson; chart[idchart][11]=hiperexp+hiperexpson; idchart=idchart; control2=1; // it is necessary to reset the value of c control3=1; // and the previous exceptional divisors } // PROBABLY IT IS NEC MORE INFORMATION !!! } // closing else otherwise } // closing else case I!=1 } // closing else for EordI=0 if (EordI<0) {print("ERROR in chart"); print(idchart); ~; break;} //----------------------- guardar de momento-------- // if (EordI==0) {auxfchart[1]=chart[idchart]; // finalchart=finalchart+auxfchart; // L2=Gamma(expJ,c,n); // HAY QUE APLICARLO AL M NO AL J // invariant[idchart]=L2[2]; // auxinv[idchart]=list(0); // nsons[idchart]=0; // idchart=idchart+1;} //------------------------------------------------ } // END ELSE //--------------------------------------------------- } // END LOOP WHILE tnum=step[nchart]; totalinv=resfunction(invariant,auxinv,nchart,n); return(chart,finalchart,invariant,nchart,step,nsons,auxinv,mobile,totalinv); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..2)),dp; ideal J=x(1)^2-x(2)^3; list L=Eresol(J); "Please press return after each break point to see the next element of the output list"; L[1][1]; // information of the first chart, L[1] list of charts ~; L[2]; // list of charts with information about sons ~; L[3]; // invariant, "#" means solved chart ~; L[4]; // number of charts, 7 in this example ~; L[5]; // height corresponding to each chart ~; L[6]; // number of sons ~; L[7]; // auxiliary invariant ~; L[8]; // H exceptional divisors and more information ~; L[9]; // complete resolution function "Second example, write L[i] to see the i-th component of the list"; ring r = 0,(x(1..3)),dp; ideal J=x(1)^2*x(2),x(3)^3; // SOLVED! list L=Eresol(J); L[4]; // 16 charts L[9]; // complete resolution function ~; "Third example, write L[i] to see the i-th component of the list"; ring r = 0,(x(1..2)),dp; ideal J=x(1)^3-x(1)*x(2)^3; list L=Eresol(J); L[4]; // 8 charts, rational exponents L[9]; // complete resolution function ~; } ////////////////////////////////////////////////////////////////////////////////////// proc resfunction(list invariant, list auxinv, int nchart,int n) "USAGE: resfunction(invariant,auxinv,nchart,n); invariant, auxinv lists, nchart, n integers COMPUTE: Patch the resolution function RETURN: The complete resolution function EXAMPLE: example resfunction; shows an example " { int i,j,l,k; list patchfun,aux; for (i=1;i<=nchart;i++){patchfun[i]=invariant[i];} for (i=1;i<=nchart;i++){if (auxinv[i][1]!=0 and size(auxinv[i])==3){l=size(invariant[i]); for (j=1;j<=l;j++){ if (invariant[i][j]==0){aux=auxinv[i]; patchfun[i][j]=aux; if (l0){list Lring=ringlist(basering); Lring[1]=0; // def r=basering; def Rnew=ring(Lring); setring Rnew; ideal chy=maxideal(1); map fRnew=r,chy; ideal J=fRnew(J); // E-RESOLUTION, Computations in char 0 list L=Eresol(J); // STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL // not implemented yet, CHAR p !!!! // STEP 3: DO THE E-RESOLUTION AGAIN (char 0 again) // generating output in char p int q=lcmofall(L[4],L[8]); // lcm of the denominators list B=genoutput(L[1],L[8],L[4],L[6],n,q,p); // generate output needed for visualization // setring r; // Back to the basering // ideal chy=maxideal(1); // map fr=Rnew,chy; // list L=fr(L); // list B=fr(B); } else{ // E-RESOLUTION list L=Eresol(J); // STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL // not implemented yet // STEP 3: DO THE E-RESOLUTION AGAIN // generating output int q=lcmofall(L[4],L[8]); list B=genoutput(L[1],L[8],L[4],L[6],n,q,p); } return(B); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..2)),dp; ideal J=x(1)^2-x(2)^3; list B=BINresol(J); B[1]; // list of final charts B[2]; // list of all charts ring r = 2,(x(1..3)),dp; ideal J=x(1)^2-x(2)^2*x(3)^2; list B=BINresol(J); B[2]; // list of all charts } /////////////////////////////////////////////////////// proc Maxord(list L,int n) "USAGE: Maxord(L,n); L list, n integer COMPUTE: Find the maximal entry of a list, input is a list defining a monomial RETURN: maximum entry of a list and its position EXAMPLE: example Maxord; shows an example " {int i,can; number canmax; list aux; canmax=1; can=1; for (i=1;i<=n;i++) { if (L[i]>=canmax and i>=can) {canmax=L[i]; can=i;}} return(canmax,can); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; ideal J=x(1)^2*x(2)*x(3)^5; list L=data(J,1,3); L[2]; // list of exponents Maxord(L[2][1][1],3); } /////////////////////////////////////////////////////// proc Gamma(list expM,number c,int n) "USAGE: Gamma(L,c,n); L list, c number, n integer COMPUTE: The Gamma function, resolution function corresponding to the monomial case RETURN: lists of maximum exponents in L, value of Gamma function, center of blow up EXAMPLE: example Gamma; shows an example " {int i,j,k,l,cont,can; intvec upla; number canmax; list expM2,gamma,L,aux,maxlist,center,aux2; i=1; cont=0; expM2=expM; while (cont==0 and i<=n) { L=Maxord(expM2,n); aux=L[1]; maxlist=maxlist + aux; can=L[2]; if (i==1) {upla=can; center=can;} else {upla=upla,can; aux2=can; center=center+aux2;} canmax=sum(maxlist); if (canmax>=c) {gamma[1]=-i; gamma[2]=canmax/c; gamma[3]=upla; cont=1;} else {expM2[can]=0;} i=i+1; } return(maxlist,gamma,center); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..5)),dp; ideal J=x(1)^2*x(2)*x(3)^5*x(4)^2*x(5)^3; list L=data(J,1,5); list G=Gamma(L[2][1][1],9,5); // critical value c=9 G[1]; // maximum exponents in the ideal G[2]; // maximal value of Gamma function G[3]; // center given by Gamma } /////////////////////////////////////////////////////// proc convertdata(list C,list L, int n, list flaglist) "USAGE: convertdata(C,L,n,flaglist); C, L, flaglist lists, n integer COMPUTE: Compute the ideal corresponding to the given lists C,L RETURN: an ideal whose coefficients are given by C, exponents given by L EXAMPLE: example convertdata; shows an example " {int i,j,k,a,b,lon; poly aux,aux1,aux2,aux3,f; ideal J; aux=poly(0); aux1=poly(1); aux2=poly(0); aux3=poly(1); k=size(L); for (i=1;i<=k;i++){lon=size(C[i]); if (lon==1){ // variables in the monomial for (j=1;j<=n;j++){a=int(poly(L[i][1][j])); if (a!=0){ if (flaglist[j]==0){aux=poly(x(j)^a); aux1=aux1*aux;} else {aux=poly(y(j)^a); aux1=aux1*aux;} } } if (C[i][1]!=0){aux1=C[i][1]*aux1;} // we add the coefficient else {aux1=0;} J[i]=aux1; aux1=poly(1); } else{ // variables in the binomial for (j=1;j<=n;j++){a=int(poly(L[i][1][j])); b=int(poly(L[i][2][j])); if (a!=0){ if (flaglist[j]==0){aux=poly(x(j)^a); aux1=aux1*aux;} else {aux=poly(y(j)^a); aux1=aux1*aux;} } if (b!=0){ if (flaglist[j]==0){aux2=poly(x(j)^b); aux3=aux3*aux2;} else {aux2=poly(y(j)^b); aux3=aux3*aux2;} } } // we add the coefficients if (C[i][1]!=0){aux1=C[i][1]*aux1;} else {aux1=0;} if (C[i][2]!=0){aux3=C[i][2]*aux3;} else {aux3=0;} f=aux1+aux3; J[i]=f; aux1=poly(1); aux3=poly(1); } } return(J); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..4),y(5)),dp; list M=identifyvar(); ideal J=x(1)^2*y(5)^2-x(2)^2*x(3)^2,6*x(4)^2; list L=data(J,2,5); L[1]; // Coefficients L[2]; // Exponents ideal J2=convertdata(L[1],L[2],5,M); J2; } ///////////////////////////////////////////////////////////////////////////// proc lcmofall(int nchart,list mobile) "USAGE: lcmofall(nchart,mobile); nchart integer, mobile list of lists COMPUTE: Compute the lcm of the denominators of the E-orders of all the charts RETURN: an integer given the lcm NOTE: CALL BEFORE salida EXAMPLE: example lcmofall; shows an example " { int i,m,tip,mcmall; intvec numall; for (i=2;i<=nchart+1;i++){ tip=mobile[i][1]; if (tip!=1){numall=numall,tip;} } m=size(numall); if (m==1){mcmall=1;} else{ if (numall[1]==0){numall=numall[2..m];} mcmall=lcm(numall);} return(mcmall); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..2)),dp; ideal J=x(1)^3-x(1)*x(2)^3; list L=Eresol(J); L[4]; // 8 charts, rational exponents L[8][2][2]; // E-orders at the first chart lcmofall(8,L[8]); } ///////////////////////////////////////////////////////////////////////////// proc salida(int idchart,list chart,list mobile,int numson,intvec previousa,int n,int q,int p) "USAGE: salida(idchart,chart,mobile,numson,previousa,n,q,p); idchart, numson, n, q, p integers, chart, mobile, lists, previousa intvec COMPUTE: CONVERT THE OUTPUT OF A CHART IN A RING, WHERE DEFINE A BASIC OBJECT (BO) RETURN: the ring corresponding to the chart EXAMPLE: example salida; shows an example " { int l,i,m,aux,parent,m4,j; intvec Hhist,EOhist,aux7,aux9; list expJ,Coef,BO,blwhist,Eolist,hipercoef,hiperexp; list flag; // chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp // mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2]; // we need to define the suitable ring at this chart list Lring=ringlist(basering); def RR2=basering; flag=chart[6]; string newl; for (l=1;l<=n; l++){if (flag[l]==1){newl=string(l); Lring[2][l]="y("+newl+")";} } def RRnew=ring(Lring); setring RRnew; ideal chy=maxideal(1); map fRnew=RR2,chy; list chart=fRnew(chart); list mobile2=fRnew(mobile); flag=chart[6]; // we need to convert expJ and Coef to an ideal expJ=chart[4]; Coef=chart[5]; Hhist=chart[7]; blwhist=chart[8]; // now the ideal will be correctly defined in the ring Rnew ideal J2=convertdata(Coef,expJ,n,flag); // Computations in RRnew //------------------------------------------------------------------------------ // START TO CREATE THE BO corresponding to this chart BO=createBO(J2); // MODIFY BO WITH THE INFORMATION OF THE CHART // BO[1] an ideal, say W_i, defining the ambient space of the i-th chart of the blowing up // If there are hyperbolic equations, we put them here hipercoef=chart[10]; hiperexp=chart[11]; if (size(hipercoef)!=0){ ideal ambJ=convertdata(hipercoef,hiperexp,n,flag); BO[1]=ambJ; } // BO[2] an ideal defining the controlled transform BO[2]=J2; // BO[3] intvec, tupla containing the maximal E-order of BO[2] if (numson==0){BO[3]=1;} // we write 1 if the chart is a final chart else{ Eolist=mobile2[2]; // otherwise, convert the list of E-orders in an intvec m=size(Eolist); aux=int(Eolist[1]*q); EOhist=aux; if (m>1){for (i=2;i<=m;i++){aux=int(Eolist[i]*q); EOhist=EOhist,aux;}} BO[3]=EOhist; } // BO[4] the list of exceptional divisors given by Hhist BO[4]=constructH(Hhist,n,flag); // BO[5] an ideal defining the map K[W] ----> K[Wi] given by blwhist BO[5]=constructblwup(blwhist,n,chy,flag); // BO[6] an intvec, BO[6][j]=1 indicates that =1, i.e. the // strict transform does not meet the j-th exceptional divisor m4=size(BO[4]); ideal auxydeal; ideal Jint; for (j=1;j<=m4;j++){ auxydeal=BO[4][j]+J2; Jint=std(auxydeal); if (size(Jint)==1 and Jint[1]==1){BO[6][j]=1;} else{BO[6][j]=0;} } // BO[7] intvec, the index of the first blown-up object in the resolution process // leading to this object for which the value of b was BO[3] // the subsequent ones are the indices for the Coeff-Objects // of BO[2] used when determining the center // index of last element of H^- in H if (numson!=0){BO[7]=mobile2[8];} // it is always -1 at the final charts // BO[8] a matrix indicating that BO[4][i] meets BO[4][j] by BO[8][i,j]=1 for i < j if (m4>0){ matrix aux8[m4][m4]; BO[8]=aux8; ideal auxydeal2; ideal Jint2; for (i=1;i<=m4;i++){ for (j=i+1;j<=m4;j++){ auxydeal2=BO[4][i]+BO[4][j]; Jint2=std(auxydeal2); if (size(Jint2)==1 and Jint2[1]==1){BO[8][i,j]=0;} else{ for (l=1;l0 if (p>0){ list Lring; Lring=ringlist(RRnew); Lring[1]=p; def auxRnew=ring(Lring); kill Lring; setring auxRnew; ideal chy=maxideal(1); map frnew=RRnew,chy; def BO=frnew(BO); // def chart=frr(chart); def invSat=frnew(invSat); def lastMap=frnew(lastMap); def cent=frnew(cent); def path=frnew(path); } // export everything needed export BO; export(invSat); export lastMap; export path; export cent; if (p==0){return(RRnew);} else{ return(auxRnew);} } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..2)),dp; ideal J=x(1)^2-x(2)^3; list L=Eresol(J); list B=salida(5,L[1][5],L[8][6],2,L[1][3][3],2,1,0); // chart 5 def RR=B[1]; setring RR; BO; "press return to see next example"; ~; ring r = 0,(x(1..2)),dp; ideal J=x(1)^2-x(2)^3; list L=Eresol(J); list B=salida(7,L[1][7],L[8][8],0,L[1][5][3],2,1,0); // chart 7 def RR=B[1]; setring RR; BO; showBO(BO); "press return to see next example"; ~; ring r = 0,(x(1..2)),dp; ideal J=x(1)^3-x(1)*x(2)^3; list L=Eresol(J); // 8 charts, rational exponents list B=salida(1,L[1][1],L[8][2],2,0,2,2,0); // CHART 1 def RR=B[1]; setring RR; BO; } ///////////////////////////////////////////////////////////////////////////// // CONVERT THE OUTPUT OF Eresol IN A LIST OF RINGS, WHERE A BASIC OBJECT (BO) IS DEFINED // IN ORDER TO INTEGRATE THIS LIBRARY INSIDE THE LIBRARY resolve.lib proc genoutput(list chart,list mobile,int nchart,list nsons,int n,int q, int p) "USAGE: genoutput(chart,mobile,nchart,nsons,n,q,p); chart, mobile, nsons lists, nchart, n,q, p integers RETURN: two lists, the first one gives the rings corresponding to the final charts, the second one is the list of all rings corresponding to the affine charts of the resolution process EXAMPLE: example genoutput; shows an example " { int idchart,parent; list auxlist,solvedrings,totalringlist,previousa; list auxlistenp,solvedringsenp,totalringenp; // chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp // mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2]; idchart=1; // first loop, construct list previousa while (idchart<=nchart) { if (idchart==1){previousa[1]=chart[2][3];} else { // if there are no sons, the next center is nothing if (nsons[idchart]==0){previousa[idchart]=0;} // always fill the parent parent=chart[idchart][1]; previousa[parent]=chart[idchart][3]; } idchart=idchart+1; } // HERE BEGIN THE LOOP idchart=1; while (idchart<=nchart) { def auxexit=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q,p); if (p>0) { // we need the computations in char 0 too def auxexitenp=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q,0); } else{def auxexitenp=auxexit;} // we add the ring to the list of all rings auxlist[1]=auxexit; totalringlist=totalringlist+auxlist; auxlistenp[1]=auxexitenp; totalringenp=totalringenp+auxlistenp; // if the chart has no sons, add it to the list of final charts if (nsons[idchart]==0) { solvedrings=solvedrings+auxlist; solvedringsenp=solvedringsenp+auxlistenp; } auxlist=list(); auxlistenp=list(); kill auxexit; kill auxexitenp; idchart=idchart+1; } // EXIT WHILE return(solvedrings,totalringlist,solvedringsenp,totalringenp); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..2)),dp; ideal J=x(1)^3-x(1)*x(2)^3; list L=Eresol(J); // 8 charts, rational exponents list B=genoutput(L[1],L[8],L[4],L[6],2,2,0); // generates the output presentTree(B); list iden0=collectDiv(B); ResTree(B,iden0[1]); // generates the resolution tree // Use presentTree(B); to see the final charts // To see the tree type in another shell // dot -Tjpg ResTree.dot -o ResTree.jpg // /usr/bin/X11/xv ResTree.jpg } ///////////////////////////////////////////////////////////////////// proc computemcm(list Eolist) "USAGE: computemcm(Eolist); Eolist list RETURN: an integer, the least common multiple of the denominators of the E-orders NOTE: Make the same as lcmofall but for one chart. NECESSARY BECAUSE THE E-ORDERS ARE OF TYPE NUMBER!! EXAMPLE: example computemcm; shows an example " { int m,i,aux,mcmchart; intvec num; m=size(Eolist); if (m==1){mcmchart=int(denominator(Eolist[1])); return(mcmchart);} if (m>1) { num=int(denominator(Eolist[1])); for (i=2;i<=m;i++) {aux=int(denominator(Eolist[i])); num=num,aux; } } mcmchart=lcm(num); return(mcmchart); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..2)),dp; ideal J=x(1)^3-x(1)*x(2)^3; list L=Eresol(J); // 8 charts, rational exponents L[8][2][2]; // maximal E-order at the first chart computemcm(L[8][2][2]); } ///////////////////////////////////////////////////////////////////// proc constructH(intvec Hhist,int n,list flag) "USAGE: constructH(Hhist,n,flag); Hhist intvec, n integer, flag list RETURN: the list of exceptional divisors accumulated at this chart EXAMPLE: example constructH; shows an example " { int i,j,m,l; list exceplist; ideal aux; m=size(Hhist); if (Hhist[1]==0 and m>1) { Hhist=Hhist[2..m]; m=m-1; for (i=1;i<=m;i++) { l=Hhist[i]; if (flag[l]==0){aux=ideal(poly(x(l))); } else {aux=ideal(poly(y(l))); } exceplist[i]=aux; } // eliminate repeated variables for (i=1;i<=m;i++) { for (j=1;j<=m;j++) { if (Hhist[i]==Hhist[j] and i!=j) { if (ij){exceplist[j]=ideal(1);} } } } } else {exceplist=list();} // else {exceplist=list(ideal(0));} // IF IT FAILS USE THIS return(exceplist); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; list flag=identifyvar(); ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; list L=Eresol(J); // 7 charts // history of the exceptional divisors at the 7-th chart L[1][7][7]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts constructH(L[1][7][7],3,flag); } ///////////////////////////////////////////////////////////////////// proc constructblwup(list blwhist,int n,ideal chy,list flag) "USAGE: constructblwup(blwhist,n,chy,flag); blwhist, flag lists, n integer, chy ideal RETURN: the ideal defining the map K[W] --> K[Wi], which gives the composition map of all the blowing up leading to this chart NOTE: NECESSARY START WITH COLUMNS EXAMPLE: example constructblwup; shows an example " { int i,j,m,m2; poly aux2; m=size(blwhist[1]); for (j=1;j<=m;j++) { for (i=1;i<=n;i++) { m2=blwhist[i][j]; // If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not if (m2!=0) { if (flag[m2]==0){aux2=poly(x(m2));} else {aux2=poly(y(m2));} // And then substitute this variable for the corresponding product in the whole ideal if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);} else {chy=subst(chy,y(i),y(i)*aux2);} } } } return(chy); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; list flag=identifyvar(); ideal chy=maxideal(1); ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; list L=Eresol(J); // 7 charts // history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts constructblwup(L[1][7][8],3,chy,flag); } ///////////////////////////////////////////////////////////////////// proc constructlastblwup(list blwhist,int n,ideal chy,list flag) "USAGE: constructlastblwup(blwhist,n,chy,flag); blwhist, flag lists, n integer, chy ideal RETURN: the ideal defining the last blow up NOTE: NECESSARY START WITH COLUMNS EXAMPLE: example constructlastblwup; shows an example " { int i,j,m,m2; poly aux2; m=size(blwhist[1]); if (m>0) { for (i=1;i<=n;i++){ m2=blwhist[i][m]; // If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not if (m2!=0) { if (flag[m2]==0){aux2=poly(x(m2));} else {aux2=poly(y(m2));} // And then substitute this variable for the corresponding product in the whole ideal if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);} else {chy=subst(chy,y(i),y(i)*aux2);} } } } return(chy); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; list flag=identifyvar(); ideal chy=maxideal(1); ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; list L=Eresol(J); // 7 charts // history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts constructlastblwup(L[1][7][8],3,chy,flag); } ///////////////////////////////////////////////////////////////////// proc tradtoideal(intvec a,ideal J2,list flag) "USAGE: tradtoideal(a,J2,flag); a intvec, J2 ideal, flag list COMPUTE: traslate to an ideal the intvec defining the center RETURN: the ideal of the center, given by the intvec a, or J2 if a=0 EXAMPLE: example tradtoideal; shows an example " { int i,m; ideal acenter,aux2; if (a==0) {acenter=J2;} else { m=size(a); for (i=1;i<=m;i++) { if (flag[a[i]]==0){aux2=poly(x(a[i]));} else {aux2=poly(y(a[i]));} acenter=acenter+aux2; } } return(acenter); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; list flag=identifyvar(); ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; intvec a=1,3; // first center of blowing up tradtoideal(a,J,flag); } ////////////////////////////////////////////////////////////////////////////////////// // OPERATIONS WITH LISTS ////////////////////////////////////////////////////////////////////////////////////// proc iniD(int n) "USAGE: iniD(n); n integer RETURN: list of lists of zeros of size n EXAMPLE: example iniD; shows an example " {int i,j; list D,auxD; for (j=1;j<=n; j++) {auxD[j]=0;} for (i=1;i<=n; i++) {D[i]=auxD;} return(D); } example {"EXAMPLE:"; echo = 2; iniD(3); } ///////////////////////////////////////////////////////// proc sumlist(list L1,list L2) "USAGE: sumlist(L1,L2); L1,L2 lists, (size(L1)==size(L2)) RETURN: a list, sum of L1 and L2 EXAMPLE: example sumlist; shows an example " { int i,k; list sumL; k=size(L1); if (size(L2)!=k) {return("ERROR en sumlist, lists must have the same size");} for (i=1;i<=k;i++) {sumL[i]=L1[i]+L2[i];} return(sumL); } example {"EXAMPLE:"; echo = 2; list L1=1,2,3; list L2=5,9,7; sumlist(L1,L2); } /////////////////////////////////////////////////////// proc reslist(list L1,list L2) "USAGE: reslist(L1,L2); L1,L2 lists, (size(L1)==size(L2)) RETURN: a list, subtraction of L1 and L2 EXAMPLE: example reslist; shows an example " { int i,k; list resL; k=size(L1); if (size(L2)!=k) {return("ERROR en reslist, lists must have the same size");} for (i=1;i<=k;i++) {resL[i]=L1[i]-L2[i];} return(resL); } example {"EXAMPLE:"; echo = 2; list L1=1,2,3; list L2=5,9,7; reslist(L1,L2); } ////////////////////////////////////////////////////// proc multiplylist(list L,number a) "USAGE: multiplylist(L,a); L list, a number RETURN: list of elements of type number, multiplication of L times a EXAMPLE: example multiplylist; shows an example " {int i,k; list newL,bb; number b; k=size(L); for (i=1;i<=k;i++) {b=L[i]*a; bb=b; newL=newL+bb;} return(newL); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; list L=1,2,3; multiplylist(L,1/5); } /////////////////////////////////////////////////////// proc dividelist(list L1,list L2) "USAGE: dividelist(L1,L2); L1,L2 lists RETURN: list of elements of type number, division of L1 by L2 EXAMPLE: example dividelist; shows an example " {int i,k,k1,k2; list LL,bb; number a1,a2,b; k1=size(L1); k2=size(L2); if (k2!=k1) {print("ERROR en dividelist, lists must have the same size");} if (k1<=k2) {k=k1;} else {k=k2;} for (i=1;i<=k;i++) {a1=L1[i]; a2=L2[i]; b=a1/a2; bb=b; LL=LL+bb;} return(LL); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x(1..3)),dp; list L1=1,2,3; list L2=5,9,7; dividelist(L1,L2); } /////////////////////////////////////////////////////// proc createlist(list L1,list L2) "USAGE: createlist(L1,L2); L1,L2 lists, (size(L1)==size(L2)) RETURN: list of lists of two elements, the first one of L1 and the second of L2 EXAMPLE: example createlist; shows an example " {int i,k; list L,aux; k=size(L1); if (size(L2)!=k) {return("ERROR en createlist, lists must have the same size");} L=list0(k); for (i=1;i<=k;i++) {if (L1[i]!=0) {aux=L1[i],L2[i]; L[i]=aux;} else {L=delete(L,i);}} return(L); } example {"EXAMPLE:"; echo = 2; list L1=1,2,3; list L2=5,9,7; createlist(L1,L2); } /////////////////////////////////////////////////////// static proc list0(int n) "USAGE: list0(n); n integer RETURN: list of n zeros EXAMPLE: example list0; shows an example " {int i; list L0; for (i=1;i<=n;i++) {L0[i]=0;} return(L0); } example {"EXAMPLE:"; echo = 2; list0(4); } ////////////////////////////////////////////////////////////