[7c0817] | 1 | // rocio, last modified 23.05.10 |
---|
| 2 | //////////////////////////////////////////////////////////////////////////// |
---|
| 3 | version="$Id: resbin.lib$"; |
---|
[f9b430] | 4 | category="Commutative algebra"; |
---|
[7c0817] | 5 | info=" |
---|
| 6 | LIBRARY: resbin.lib Combinatorial algorithm of resolution of singularities |
---|
| 7 | of binomial ideals in arbitrary characteristic. |
---|
| 8 | Binomial resolution algorithm of Blanco |
---|
| 9 | |
---|
| 10 | AUTHORS: R. Blanco, mariarocio.blanco@uclm.es, |
---|
| 11 | @* G. Pfister, pfister@mathematik.uni-kl.de |
---|
| 12 | |
---|
[f9b430] | 13 | SEE ALSO: resolve_lib |
---|
| 14 | |
---|
[7c0817] | 15 | PROCEDURES: |
---|
| 16 | BINresol(J); computes a E-resolution of singularities of (J) (THE SECOND PART IS NOT IMPLEMENTED YET) |
---|
| 17 | Eresol(J); computes a E-resolution of singularities of (J) in char 0 |
---|
| 18 | determinecenter(L1,L2,c,n,Y,a,mb,flag,control3); computes the next blowing-up center |
---|
| 19 | Blowupcenter(L1,id,m,L2,c,n,h); makes the blowing-up |
---|
| 20 | Nonhyp(Coef,expJ,sJ,n,flag,sums); computes the ideal generated by the non hyperbolic generators of expJ |
---|
| 21 | inidata(K,k); verifies input data, a binomial ideal K of k generators |
---|
| 22 | identifyvar(); identifies status of variables |
---|
| 23 | data(K,k,n); transforms data on lists of lenght n |
---|
| 24 | Edatalist(Coef,Exp,k,n,flag); gives the E-order of each term in Exp |
---|
| 25 | EOrdlist(Coef,Exp,k,n,flag); computes the E-order of an ideal (giving in the language of lists) |
---|
| 26 | maxEord(Coef,Exp,k,n,flag); computes de maximum E-order of an ideal given by Coef and Exp |
---|
| 27 | ECoef(Coef,expP,sP,V,auxc,n,flag); Computes a simplified version of the E-Coeff ideal. The E-orders are correct, |
---|
| 28 | but tranformations of coefficients of the generators and powers of binomials |
---|
| 29 | cannot be computed easily in terms of lists. |
---|
| 30 | elimrep(L); removes repeated terms from a list |
---|
| 31 | Emaxcont(Coef,Exp,k,n,flag); computes a list of hypersurfaces of E-maximal contact |
---|
| 32 | cleanunit(mon,n,flag); clean the units in a monomial mon |
---|
| 33 | resfunction(t,auxinv,nchart,n); composes the E-resolution function |
---|
| 34 | calculateI(Coef,J,c,n,Y,a,b,D); computes the order of the non monomial part of an ideal J |
---|
| 35 | Maxord(L,n); computes the maximum exponent of an exceptional monomial ideal |
---|
| 36 | Gamma(L,c,n); computes the Gamma function for an exceptional monomial ideal given by L |
---|
| 37 | convertdata(C,L,n,flag); computes the ideal corresponding to C,L |
---|
| 38 | tradblwup(blwhist,n,Y,a,num); composes the blowing up at this chart |
---|
| 39 | lcmofall(nchart,mobile); computes the lcm of the denominators of the E-orders for all the charts |
---|
| 40 | computemcm(Eolist); computes the lcm of the denominators of the E-orders for one chart |
---|
| 41 | constructH(Hhist,n,flag); construct the list of exceptional divisors accumulated at this chart |
---|
| 42 | constructblwup(blwhist,n,chy,flag); construct the ideal defining the map K[W] --> K[Wi], |
---|
| 43 | which gives the composition map of all the blowing up leading to this chart |
---|
| 44 | constructlastblwup(blwhist,n,chy,flag); construct the ideal defining the last blowup leading to this chart |
---|
| 45 | genoutput(chart,mobile,nchart,nsons,n,q,p); generates the output for visualization |
---|
| 46 | salida(idchart,chart,mobile,numson,previousa,n,q); generates the output for one chart |
---|
| 47 | iniD(n); creates a list of lists of zeros of size n |
---|
| 48 | sumlist(L1,L2); sums two lists component to component |
---|
| 49 | reslist(L1,L2); subtracts two lists component to component |
---|
| 50 | multiplylist(L,a); multiplies a list by a number, component to component |
---|
| 51 | dividelist(L1,L2); divides two lists component to component |
---|
| 52 | createlist(L1,L2); creates a list of lists of two elements |
---|
| 53 | list0(n); creates a list of zeros of size n |
---|
| 54 | |
---|
| 55 | "; |
---|
| 56 | |
---|
| 57 | LIB "general.lib"; |
---|
| 58 | LIB "qhmoduli.lib"; |
---|
| 59 | LIB "inout.lib"; |
---|
| 60 | LIB "poly.lib"; |
---|
| 61 | LIB "resolve.lib"; |
---|
| 62 | LIB "reszeta.lib"; |
---|
| 63 | LIB "resgraph.lib"; |
---|
| 64 | //////////////////////////////////////////////////////////////////////////// |
---|
| 65 | |
---|
| 66 | proc inidata(ideal K,int k) |
---|
| 67 | "USAGE: inidata(K,k); K any ideal, k integer (!=0) |
---|
| 68 | COMPUTE: Verifies the input data |
---|
| 69 | RETURN: flag indicating if the ideal is binomial or not |
---|
| 70 | EXAMPLE: example inidata; shows an example |
---|
| 71 | " |
---|
| 72 | { |
---|
| 73 | int i; |
---|
| 74 | for (i=1;i<=k; i++) |
---|
| 75 | { if (size(K[i])>2){return(0);} |
---|
| 76 | } |
---|
| 77 | return(1); |
---|
| 78 | } |
---|
| 79 | example |
---|
| 80 | {"EXAMPLE:"; echo = 2; |
---|
| 81 | ring r = 0,(x(1..3)),dp; |
---|
| 82 | ideal J1=x(1)^4*x(2)^2, x(1)^2+x(3)^3; |
---|
| 83 | inidata(J1,2); |
---|
| 84 | |
---|
| 85 | ideal J2=x(1)^4*x(2)^2, x(1)^2+x(2)^3+x(3)^5; |
---|
| 86 | inidata(J2,2); |
---|
| 87 | } |
---|
| 88 | ///////////////////////////////////////////////////////////////////////////////// |
---|
| 89 | |
---|
| 90 | proc changeoriginalvar() |
---|
| 91 | "USAGE: changeoriginalvar(); |
---|
| 92 | COMPUTE: Change the name of the variables to x(1...n), only necessary at the beginning |
---|
| 93 | RETURN: the new ring with the suitable names |
---|
| 94 | EXAMPLE: example changeoriginalvar; shows an example |
---|
| 95 | " |
---|
| 96 | { |
---|
| 97 | int i,n,cont; |
---|
| 98 | |
---|
| 99 | n=nvars(basering); |
---|
| 100 | cont=0; |
---|
| 101 | def r=basering; |
---|
| 102 | |
---|
| 103 | // check the name of the variables |
---|
| 104 | |
---|
| 105 | for (i=1;i<=n; i++){if (varstr(i)[1]=="x" or varstr(i)[1]=="y"){cont=cont+1;}} |
---|
| 106 | |
---|
| 107 | // change them if there exists some variable different from x(i) or y(i) |
---|
| 108 | |
---|
| 109 | if (cont!=n or n<=2){ |
---|
| 110 | |
---|
| 111 | // defining the new variables in Lring2[2] |
---|
| 112 | |
---|
| 113 | list Lring,Lring2; |
---|
| 114 | Lring=ringlist(basering); |
---|
| 115 | |
---|
| 116 | ring raux=0,(x(1..n)),dp; |
---|
| 117 | setring r; |
---|
| 118 | Lring2=ringlist(raux); |
---|
| 119 | |
---|
| 120 | // making the change |
---|
| 121 | |
---|
| 122 | for (i=1;i<=n; i++){ Lring[2][i]=Lring2[2][i];} |
---|
| 123 | |
---|
| 124 | def Rnew=ring(Lring); |
---|
| 125 | setring Rnew; |
---|
| 126 | // print("INVERTIBLE VARIABLES NOT CONSIDERED AT THE BEGINNING"); |
---|
| 127 | return(Rnew,1); |
---|
| 128 | } |
---|
| 129 | else{ // print("INVERTIBLE VARIABLES ALREADY CONSIDERED AT THE BEGINNING"); |
---|
| 130 | return(r,0); |
---|
| 131 | } |
---|
| 132 | } |
---|
| 133 | example |
---|
| 134 | {"EXAMPLE:"; echo = 2; |
---|
| 135 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 136 | changeoriginalvar(); |
---|
| 137 | |
---|
| 138 | ring r = 0,(x,y,z,w),dp; |
---|
| 139 | changeoriginalvar(); |
---|
| 140 | } |
---|
| 141 | |
---|
| 142 | ///////////////////////////////////////////////////////////////////////////////// |
---|
| 143 | |
---|
| 144 | proc identifyvar() |
---|
| 145 | "USAGE: identifyvar(); |
---|
| 146 | COMPUTE: Asign 0 to variables x and 1 to variables y, only necessary at the beginning |
---|
| 147 | RETURN: list, say l, of size the dimension of the basering |
---|
| 148 | l[i] is: 0 if the i-th variable is x(i), |
---|
| 149 | 1 if the i-th variable is y(i) |
---|
| 150 | EXAMPLE: example identifyvar; shows an example |
---|
| 151 | " |
---|
| 152 | { |
---|
| 153 | int i,n; |
---|
| 154 | list flaglist; |
---|
| 155 | |
---|
| 156 | n=nvars(basering); |
---|
| 157 | flaglist=list0(n); |
---|
| 158 | |
---|
| 159 | for (i=1;i<=n; i++){if (varstr(i)[1]=="y"){flaglist[i]=1;}} |
---|
| 160 | |
---|
| 161 | return(flaglist); |
---|
| 162 | } |
---|
| 163 | example |
---|
| 164 | {"EXAMPLE:"; echo = 2; |
---|
| 165 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 166 | identifyvar(); |
---|
| 167 | } |
---|
| 168 | ///////////////////////////////////////////////////////////////////////////////// |
---|
| 169 | |
---|
| 170 | proc data(ideal K,int k,int n) |
---|
| 171 | "USAGE: data(K,k,n); K any ideal, k integer (!=0), n integer (!=0) |
---|
| 172 | COMPUTE: Construcs a list with the coefficients and exponents of one ideal |
---|
| 173 | RETURN: lists of coefficients and exponents of K |
---|
| 174 | EXAMPLE: example data; shows an example |
---|
| 175 | " |
---|
| 176 | {int i,j,lon; |
---|
| 177 | number aa; |
---|
| 178 | intvec cc; |
---|
| 179 | list bb,dd,aux,ddaux,Coef,Exp; |
---|
| 180 | |
---|
| 181 | for (i=1;i<=k; i++) |
---|
| 182 | { lon=size(K[i]); |
---|
| 183 | |
---|
| 184 | // binomial |
---|
| 185 | if (lon==2){aa=leadcoef(K[i][1]); |
---|
| 186 | bb=aa; |
---|
| 187 | Coef[i]=bb; // coefficients |
---|
| 188 | cc=leadexp(K[i][1]); // exponents |
---|
| 189 | |
---|
| 190 | // cc is an intvec, transform cc in dd, a list of lists |
---|
| 191 | dd=cc[1..n]; |
---|
| 192 | aux[1]=dd; |
---|
| 193 | // the same for the second term |
---|
| 194 | |
---|
| 195 | aa=leadcoef(K[i][2]); |
---|
| 196 | bb=aa; |
---|
| 197 | Coef[i]=Coef[i] + bb; // all the coefficients of i-th generator of K |
---|
| 198 | cc=leadexp(K[i][2]); |
---|
| 199 | |
---|
| 200 | dd=cc[1..n]; |
---|
| 201 | aux[2]=dd; |
---|
| 202 | Exp[i]=aux;} |
---|
| 203 | |
---|
| 204 | // monomial |
---|
| 205 | if (lon==1){aux=list(); |
---|
| 206 | aa=leadcoef(K[i][1]); |
---|
| 207 | bb=aa; |
---|
| 208 | Coef[i]=bb; |
---|
| 209 | cc=leadexp(K[i][1]); |
---|
| 210 | dd=cc[1..n]; |
---|
| 211 | aux[1]=dd; |
---|
| 212 | Exp[i]=aux;} |
---|
| 213 | } //end for |
---|
| 214 | return(Coef,Exp); |
---|
| 215 | } |
---|
| 216 | example |
---|
| 217 | {"EXAMPLE:"; echo = 2; |
---|
| 218 | ring r = 0,(x(1..3)),dp; |
---|
| 219 | ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3; |
---|
| 220 | data(J,2,3); |
---|
| 221 | } |
---|
| 222 | ////////////////////////////////////////////////////// |
---|
| 223 | |
---|
| 224 | proc Edatalist(list Coef,list Exp,int k,int n,list flaglist) |
---|
| 225 | "USAGE: Edatalist(Coef,Exp,k,n,flaglist); |
---|
| 226 | Coef,Exp,flaglist lists, k,n, integers |
---|
| 227 | Exp is a list of lists of exponents, k=size(Exp) |
---|
| 228 | COMPUTE: computes a list with the E-order of each term |
---|
| 229 | RETURN: a list with the E-order of each term |
---|
| 230 | EXAMPLE: example Edatalist; shows an example |
---|
| 231 | " |
---|
| 232 | {int i,j,lon,mm; |
---|
| 233 | list dd,ss,sums; |
---|
| 234 | number aux,aux1,aux2; |
---|
| 235 | |
---|
| 236 | for (i=1;i<=k;i++){lon=size(Coef[i]); |
---|
| 237 | if (lon==1) { for (j=1;j<=n;j++){if (flaglist[j]==0){aux=aux+Exp[i][1][j];}} |
---|
| 238 | ss=aux; aux=0;} // monomial |
---|
| 239 | else { for (j=1;j<=n;j++){if (flaglist[j]==0){ aux1=aux1+Exp[i][1][j]; |
---|
| 240 | aux2=aux2+Exp[i][2][j];}} |
---|
| 241 | ss=aux1,aux2; aux1=0; aux2=0; } // binomial |
---|
| 242 | sums[i]=ss;} |
---|
| 243 | return(sums); |
---|
| 244 | } |
---|
| 245 | example |
---|
| 246 | {"EXAMPLE:"; echo = 2; |
---|
| 247 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 248 | list flag=identifyvar(); |
---|
| 249 | 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); |
---|
| 250 | list L=data(J,3,8); |
---|
| 251 | list EL=Edatalist(L[1],L[2],3,8,flag); |
---|
| 252 | EL; // E-order of each term |
---|
| 253 | |
---|
| 254 | |
---|
| 255 | ring r = 2,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 256 | list flag=identifyvar(); |
---|
| 257 | 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); |
---|
| 258 | list L=data(J,3,8); |
---|
| 259 | list EL=Edatalist(L[1],L[2],3,8,flag); |
---|
| 260 | EL; // E-order of each term IN CHAR 2, COMPUTATIONS NEED TO BE DONE IN CHAR 0 |
---|
| 261 | |
---|
| 262 | |
---|
| 263 | ring r = 0,(x(1..3)),dp; |
---|
| 264 | list flag=identifyvar(); |
---|
| 265 | ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3; |
---|
| 266 | list L=data(J,2,3); |
---|
| 267 | list EL=Edatalist(L[1],L[2],2,3,flag); |
---|
| 268 | EL; // E-order of each term |
---|
| 269 | } |
---|
| 270 | /////////////////////////////////////////////////////////////////////////////////// |
---|
| 271 | |
---|
| 272 | proc EOrdlist(list Coef,list Exp,int k,int n,list flaglist) |
---|
| 273 | "USAGE: EOrdlist(Coef,Exp,k,n,flaglist); |
---|
| 274 | Coef,Exp,flaglist lists, k,n, integers |
---|
| 275 | Exp is a list of lists of exponents, k=size(Exp) |
---|
| 276 | COMPUTE: computes de E-order of an ideal given by a list (Coef,Exp) and extra information |
---|
| 277 | RETURN: maximal E-order, and its position=number of generator and term |
---|
| 278 | EXAMPLE: example EOrdlist; shows an example |
---|
| 279 | " |
---|
| 280 | {int i,can,canpost,lon; |
---|
| 281 | number canmin; |
---|
| 282 | list sums; |
---|
| 283 | |
---|
| 284 | sums=Edatalist(Coef,Exp,k,n,flaglist); |
---|
| 285 | |
---|
| 286 | canmin=sums[1][1]; // inicializating, works also with a monomial |
---|
| 287 | for (i=1;i<=k; i++){lon=size(sums[i]); // this is 2 for binomial and 1 for monomial generators |
---|
| 288 | if (sums[i][1]<=canmin and Coef[i][1]!=0){canmin=sums[i][1]; |
---|
| 289 | can=i; canpost=1;} |
---|
| 290 | |
---|
| 291 | // if the generator is a binomial we check the second term |
---|
| 292 | |
---|
| 293 | if (lon==2) {if (sums[i][2]<canmin and Coef[i][2]!=0){canmin=sums[i][2]; |
---|
| 294 | can=i; canpost=2;}} |
---|
| 295 | } |
---|
| 296 | return(canmin,can,canpost); |
---|
| 297 | } |
---|
| 298 | example |
---|
| 299 | {"EXAMPLE:"; echo = 2; |
---|
| 300 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 301 | list flag=identifyvar(); |
---|
| 302 | 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; |
---|
| 303 | list L=data(J,4,8); |
---|
| 304 | list Eo=EOrdlist(L[1],L[2],4,8,flag); |
---|
| 305 | Eo[1]; // E-order |
---|
| 306 | Eo[2]; // generator giving the E-order |
---|
| 307 | Eo[3]; // term giving the E-order |
---|
| 308 | } |
---|
| 309 | |
---|
| 310 | ////////////////////////////////////////////////////// |
---|
| 311 | |
---|
| 312 | proc maxEord(list Coef,list Exp,int k,int n,list flaglist) |
---|
| 313 | "USAGE: maxEord(Coef,Exp,k,n,flaglist); |
---|
| 314 | Coef,Exp,flaglist lists, k,n, integers |
---|
| 315 | Exp is a list of lists of exponents, k=size(Exp) |
---|
| 316 | RETURN: computes de maximal E-order of an ideal given by Coef,Exp |
---|
| 317 | EXAMPLE: example maxEord; shows an example |
---|
| 318 | " |
---|
| 319 | { |
---|
| 320 | int i,lon; |
---|
| 321 | number canmin; // THE ASSIGNMENT IS NOT OK BECAUSE IT IS OF TYPE NUMBER |
---|
| 322 | list sums; |
---|
| 323 | |
---|
| 324 | sums=Edatalist(Coef,Exp,k,n,flaglist); |
---|
| 325 | |
---|
| 326 | canmin=sums[1][1]; // inicializating, works also with a monomial |
---|
| 327 | for (i=1;i<=k; i++){lon=size(sums[i]); // this is 2 for binomial and 1 for monomial generators |
---|
| 328 | if (sums[i][1]<=canmin and Coef[i][1]!=0){canmin=sums[i][1];} |
---|
| 329 | |
---|
| 330 | // if the generator is a binomial we check the second term |
---|
| 331 | |
---|
| 332 | if (lon==2) {if (sums[i][2]<canmin and Coef[i][2]!=0){canmin=sums[i][2];}} |
---|
| 333 | } |
---|
| 334 | return(canmin,sums); |
---|
| 335 | } |
---|
| 336 | example |
---|
| 337 | {"EXAMPLE:"; echo = 2; |
---|
| 338 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 339 | list flag=identifyvar(); |
---|
| 340 | 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; |
---|
| 341 | list L=data(J,4,8); |
---|
| 342 | list M=maxEord(L[1],L[2],4,8,flag); |
---|
| 343 | M[1]; // E-order |
---|
| 344 | } |
---|
| 345 | ////////////////////////////////////////////////////// |
---|
| 346 | |
---|
| 347 | proc elimrep(list maxvar) |
---|
| 348 | "USAGE: elimrep(L); L is a list |
---|
| 349 | COMPUTE: Eliminate repeated terms from a list |
---|
| 350 | RETURN: the same list without repeated terms |
---|
| 351 | EXAMPLE: example elimrep; shows an example |
---|
| 352 | " |
---|
| 353 | { |
---|
| 354 | int i,j; |
---|
| 355 | list aux2; |
---|
| 356 | |
---|
| 357 | aux2=maxvar; |
---|
| 358 | for (i=1;i<=size(aux2); i++) |
---|
| 359 | { for (j=i+1;j<=size(aux2); j++){if (aux2[i]==aux2[j] and i!=j){aux2=delete(aux2,j);}} |
---|
| 360 | } |
---|
| 361 | maxvar=aux2; |
---|
| 362 | return(maxvar); |
---|
| 363 | } |
---|
| 364 | example |
---|
| 365 | {"EXAMPLE:"; echo = 2; |
---|
| 366 | ring r = 0,(x(1..3)),dp; |
---|
| 367 | list L=4,5,2,5,7,8,6,3,2; |
---|
| 368 | elimrep(L); |
---|
| 369 | } |
---|
| 370 | ////////////////////////////////////////////////////// |
---|
| 371 | |
---|
| 372 | proc Emaxcont(list Coef,list Exp,int k,int n,list flag) |
---|
| 373 | "USAGE: Emaxcont(Coef,Exp,k,n,flag); |
---|
| 374 | Coef,Exp,flag lists, k,n, integers |
---|
| 375 | Exp is a list of lists of exponents, k=size(Exp) |
---|
| 376 | COMPUTE: Identify ALL the variables of E-maximal contact |
---|
| 377 | RETURN: a list with the indexes of the variables of E-maximal contact |
---|
| 378 | EXAMPLE: example Emaxcont; shows an example |
---|
| 379 | " |
---|
| 380 | { |
---|
| 381 | int i,j,lon; |
---|
| 382 | number maxEo; |
---|
| 383 | list L,sums,bx,maxvar; |
---|
| 384 | |
---|
| 385 | L=maxEord(Coef,Exp,k,n,flag); |
---|
| 386 | |
---|
| 387 | maxEo=L[1]; |
---|
| 388 | sums=L[2]; |
---|
| 389 | |
---|
| 390 | if (maxEo>0){ |
---|
| 391 | |
---|
| 392 | for (i=1;i<=k; i++){lon=size(sums[i]); |
---|
| 393 | if (lon==2){if (sums[i][1]==maxEo) // variables of the first term |
---|
| 394 | {for (j=1;j<=n; j++){if(Exp[i][1][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}}} |
---|
| 395 | |
---|
| 396 | if (sums[i][2]==maxEo) // variables of the second term |
---|
| 397 | {for (j=1;j<=n; j++){if(Exp[i][2][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}}}} |
---|
| 398 | else {if (sums[i][1]==maxEo) |
---|
| 399 | {for (j=1;j<=n; j++){if(Exp[i][1][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}}}} |
---|
| 400 | |
---|
| 401 | }} |
---|
| 402 | else {maxvar=list();} |
---|
| 403 | |
---|
| 404 | // eliminating repeated terms |
---|
| 405 | maxvar=elimrep(maxvar); |
---|
| 406 | |
---|
| 407 | // It is necessary to check if flag[j]==0 in order to avoid the selection of y variables |
---|
| 408 | |
---|
| 409 | return(maxEo,maxvar); |
---|
| 410 | } |
---|
| 411 | example |
---|
| 412 | {"EXAMPLE:"; echo = 2; |
---|
| 413 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 414 | list flag=identifyvar(); |
---|
| 415 | 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; |
---|
| 416 | list L=data(J,4,8); |
---|
| 417 | list hyp=Emaxcont(L[1],L[2],4,8,flag); |
---|
| 418 | hyp[1]; // max E-order=0 |
---|
| 419 | hyp[2]; // There are no hypersurfaces of E-maximal contact |
---|
| 420 | |
---|
| 421 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; |
---|
| 422 | list flag=identifyvar(); |
---|
| 423 | 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; |
---|
| 424 | list L=data(J,4,8); |
---|
| 425 | list hyp=Emaxcont(L[1],L[2],4,8,flag); |
---|
| 426 | hyp[1]; // the E-order is 1 |
---|
| 427 | hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of E-maximal contact |
---|
| 428 | |
---|
| 429 | } |
---|
| 430 | /////////////////////////////////////////////////////// |
---|
| 431 | |
---|
| 432 | proc cleanunit(list mon,int n,list flaglist) |
---|
| 433 | "USAGE: cleanunit(mon,n,flaglist); |
---|
| 434 | mon, flaglist lists, n integer |
---|
| 435 | COMPUTE: We clean (or forget) the units in a monomial, given by "y" variables |
---|
| 436 | RETURN: The list defining the monomial ideal already cleaned |
---|
| 437 | EXAMPLE: example cleanunit; shows an example |
---|
| 438 | " |
---|
| 439 | { |
---|
| 440 | int i; |
---|
| 441 | |
---|
| 442 | for (i=1;i<=n;i++){if (flaglist[i]==1){mon[i]=0;}} |
---|
| 443 | |
---|
| 444 | // coef[1]=coef[1]*y(i)^mon[i]; IS NOT ALLOWED because mon[i] can be a number |
---|
| 445 | // therefore, the coefficients remain constant |
---|
| 446 | |
---|
| 447 | return(mon); |
---|
| 448 | } |
---|
| 449 | example |
---|
| 450 | {"EXAMPLE:"; echo = 2; |
---|
| 451 | ring r = 0,(x(1),y(2),x(3),y(4)),dp; |
---|
| 452 | list flag=identifyvar(); |
---|
| 453 | ideal J=x(1)^3*y(2)*x(3)^5*y(4)^8; |
---|
| 454 | list L=data(J,1,4); |
---|
| 455 | L[2][1][1]; // list of exponents of the monomial J |
---|
| 456 | list M=cleanunit(L[2][1][1],4,flag); |
---|
| 457 | M; // new list without units |
---|
| 458 | } |
---|
| 459 | ////////////////////////////////////////////////////// |
---|
| 460 | // Classification of the ideal E-Coeff_V(P): |
---|
| 461 | // ccase=1, E-Coeff_V(P)=0 |
---|
| 462 | // 2,3 Bold regular case |
---|
| 463 | // 4 P=1 monomial case (detected before) |
---|
| 464 | // 0 Otherwise |
---|
| 465 | |
---|
| 466 | proc ECoef(list Coef,list expP,int sP,int V,number auxc,int n,list flaglist) |
---|
| 467 | "USAGE: ECoef(Coef,expP,sP,V,auxc,n,flaglist); |
---|
| 468 | Coef, expP, flaglist lists, sP, V, n integers, auxc number |
---|
| 469 | COMPUTE: The ideal E-Coeff_V(P), where V is a permissible hypersurface which belongs to the center |
---|
| 470 | RETURN: list of exponents, list of coefficients and classification of the ideal E-Coeff_V(P) |
---|
| 471 | EXAMPLE: example ECoef; shows an example |
---|
| 472 | " |
---|
| 473 | { |
---|
| 474 | int i,j,k,l,numg,ccase,cont2,cont3,val; |
---|
| 475 | number aa; |
---|
| 476 | list Eco,newcoef,auxexp,newL,rs,rs2,aux,aux2,aux3,aux4,L; |
---|
| 477 | |
---|
| 478 | auxexp=expP; |
---|
| 479 | |
---|
| 480 | l=1; |
---|
| 481 | for (i=1;i<=sP;i++) |
---|
| 482 | {rs[i]=size(Coef[i]); |
---|
| 483 | if (rs[i]==2){ // binomials |
---|
| 484 | if (auxexp[i][1][V]!=auxexp[i][2][V]) // no common factors for the variable in V |
---|
| 485 | |
---|
| 486 | {for (j=1;j<=2;j++){if (auxexp[i][j][V]<auxc){aa=auxc/(auxc-auxexp[i][j][V]); |
---|
| 487 | auxexp[i][j][V]=0; |
---|
| 488 | aux4[1]=multiplylist(auxexp[i][j],aa); |
---|
| 489 | Eco[l]=aux4; |
---|
| 490 | // newcoef[l]=Coef[i][j]^aa; IT IS NO ALLOWED!!! |
---|
| 491 | newcoef[l]=Coef[i][j]; // we leave it constant |
---|
| 492 | l=l+1;}}} |
---|
| 493 | |
---|
| 494 | else // common factors for the variable in V, of zero in both terms |
---|
| 495 | |
---|
| 496 | {if (auxexp[i][1][V]<auxc){aa=auxc/(auxc-auxexp[i][1][V]); |
---|
| 497 | auxexp[i][1][V]=0; auxexp[i][2][V]=0; |
---|
| 498 | |
---|
| 499 | // this generator is a power of a binomial |
---|
| 500 | // one possibility is Eco[l]=auxexp[i]; we leave it constant and add some extra number aa, or |
---|
| 501 | // define a binomial again. The E-order coincides!!! |
---|
| 502 | |
---|
| 503 | aux=multiplylist(auxexp[i][1],aa); |
---|
| 504 | aux2=multiplylist(auxexp[i][2],aa); |
---|
| 505 | aux3[1]=aux; |
---|
| 506 | aux3[2]=aux2; |
---|
| 507 | Eco[l]=aux3; |
---|
| 508 | newcoef[l]=Coef[i]; |
---|
| 509 | l=l+1;}} |
---|
| 510 | } |
---|
| 511 | |
---|
| 512 | else // monomials |
---|
| 513 | {if (auxexp[i][1][V]<auxc){aa=auxc/(auxc-auxexp[i][1][V]); |
---|
| 514 | auxexp[i][1][V]=0; |
---|
| 515 | aux4=list(); |
---|
| 516 | aux4[1]=multiplylist(auxexp[i][1],aa); |
---|
| 517 | Eco[l]=aux4; |
---|
| 518 | newcoef[l]=Coef[i]; |
---|
| 519 | l=l+1;}} |
---|
| 520 | } |
---|
| 521 | |
---|
| 522 | // cleaning units from the monomial generators of Eco |
---|
| 523 | // If there are hyperbolic equations in Eco, such that Eco=1, we detect it later, computing the E-order |
---|
| 524 | |
---|
| 525 | numg=size(Eco); |
---|
| 526 | for (k=1;k<=numg;k++){ if (size(newcoef[k])==1){Eco[k][1]=cleanunit(Eco[k][1],n,flaglist);}} |
---|
| 527 | |
---|
| 528 | // checking Eco |
---|
| 529 | |
---|
| 530 | ccase=0; |
---|
| 531 | cont2=0; |
---|
| 532 | cont3=0; |
---|
| 533 | val=0; |
---|
| 534 | |
---|
| 535 | // CASE Eco=0: If Eco=empty list then as ideal Eco=0 |
---|
| 536 | |
---|
| 537 | if (numg==0){ccase=1;} |
---|
| 538 | else |
---|
| 539 | { |
---|
| 540 | for (i=1;i<=numg;i++) {rs2[i]=size(newcoef[i]); |
---|
| 541 | if (rs2[i]==1){val=val+n; // monomials |
---|
| 542 | for (l=1;l<=n; l++) {if (Eco[i][1][l]==0) {cont2=cont2+1;}} |
---|
| 543 | } |
---|
| 544 | else{val=val+(2*n); // binomials |
---|
| 545 | for (l=1;l<=n; l++) {if (Eco[i][1][l]==0) {cont2=cont2+1;} |
---|
| 546 | if (Eco[i][2][l]==0) {cont2=cont2+1;}} |
---|
| 547 | } |
---|
| 548 | } |
---|
| 549 | |
---|
| 550 | // If cont2=val then all the entries of Eco are zero!! As ideal Eco=1 |
---|
| 551 | |
---|
| 552 | for (i=1;i<=sP;i++){if (rs[i]==2){ // binomials |
---|
| 553 | for (l=1;l<=n;l++) {if (expP[i][1][l]!=0) {cont3=cont3+1;} |
---|
| 554 | if (expP[i][2][l]!=0) {cont3=cont3+1;}} |
---|
| 555 | } |
---|
| 556 | else{ // monomials |
---|
| 557 | for (l=1;l<=n;l++) {if (expP[i][1][l]!=0) {cont3=cont3+1;}} |
---|
| 558 | } |
---|
| 559 | } |
---|
| 560 | |
---|
| 561 | // If cont3=0 all the entries of expP are zero!! As ideal P=1 this is detected before |
---|
| 562 | // If cont3=1 then P is bold regular |
---|
| 563 | |
---|
| 564 | |
---|
| 565 | // CASE Eco=1 |
---|
| 566 | |
---|
| 567 | if (cont2==val and cont3==1){ccase=2;} // BOLD REGULAR CASE |
---|
| 568 | if (cont2==val and cont3>1){ccase=3;} // CASE P=x^{\alpha},x^{\beta}, IN FACT, BOLD REGULAR |
---|
| 569 | if (cont2==val and cont3==0){ccase=4;} // P=1, then I=1 monomial case |
---|
| 570 | |
---|
| 571 | // Case BOLD REGULAR P=x^{\alpha}*(1-\mu y^{\delta}) |
---|
| 572 | // IT IS NON NECESSARY TO CHECK IT, Eco=empty list, already done! |
---|
| 573 | |
---|
| 574 | L=maxEord(newcoef,Eco,numg,n,flaglist); // L[1] is the E-order of Eco |
---|
| 575 | if (L[1]==0){ccase=2; print("E-order zero!");} // BOLD REGULAR CASE |
---|
| 576 | |
---|
| 577 | // we leave it to check the computations |
---|
| 578 | |
---|
| 579 | } // close else |
---|
| 580 | |
---|
| 581 | return(Eco,newcoef,ccase); |
---|
| 582 | } |
---|
| 583 | example |
---|
| 584 | {"EXAMPLE:"; echo = 2; |
---|
| 585 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp; |
---|
| 586 | list flag=identifyvar(); |
---|
| 587 | 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; |
---|
| 588 | list L=data(P,3,7); |
---|
| 589 | list L2=ECoef(L[1],L[2],3,1,3,7,flag); |
---|
| 590 | L2[1]; // exponents of the E-Coefficient ideal respect to x(1) |
---|
| 591 | L2[2]; // its coefficients |
---|
| 592 | L2[3]; // classify the type of ideal obtained |
---|
| 593 | |
---|
| 594 | ring r = 0,(x(1),y(2),x(3),y(4)),dp; |
---|
| 595 | list flag=identifyvar(); |
---|
| 596 | ideal J=x(1)^3*(1-2*y(2)*y(4)^2); // Bold regular case |
---|
| 597 | list L=data(J,1,4); |
---|
| 598 | list L2=ECoef(L[1],L[2],1,1,3,4,flag); |
---|
| 599 | L2; |
---|
| 600 | |
---|
| 601 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp; |
---|
| 602 | list flag=identifyvar(); |
---|
| 603 | 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; |
---|
| 604 | list L=data(J,3,7); |
---|
| 605 | list L2=ECoef(L[1],L[2],3,1,2,7,flag); |
---|
| 606 | L2; |
---|
| 607 | |
---|
| 608 | ring r = 3,(x(1),y(2),x(3),y(4),x(5..7)),dp; |
---|
| 609 | list flag=identifyvar(); |
---|
| 610 | 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; |
---|
| 611 | list L=data(J,3,7); |
---|
| 612 | list L2=ECoef(L[1],L[2],3,1,2,7,flag); |
---|
| 613 | L2; // THE COMPUTATIONS ARE NOT CORRECT IN CHARACTERISTIC p>0 |
---|
| 614 | // because numbers are treated as 0 in assignments |
---|
| 615 | |
---|
| 616 | } |
---|
| 617 | //////////////////////////////////////////////////////////////////////////// |
---|
| 618 | // The intvec a indicates the previous center |
---|
| 619 | // Hhist = intvec of exceptional divisors of the parent chart |
---|
| 620 | |
---|
| 621 | proc determinecenter(list Coef,list expJ,number c,int n,int Y,intvec a,list listmb,list flag,int control3,intvec Hhist) |
---|
| 622 | "USAGE: determinecenter(Coef,expJ,c,n,Y,a,listmb,flag,control3,Hhist); |
---|
| 623 | Coef, expJ, listmb, flag lists, c number, n, Y, control3 integers, a, Hhist intvec |
---|
| 624 | COMPUTE: next center of blowing up and related information, see example |
---|
| 625 | RETURN: several lists defining the center and related information |
---|
| 626 | EXAMPLE: example determinecenter; shows an example |
---|
| 627 | " |
---|
| 628 | {int i,j,rstep,l,mm,cont,cont1,cont2,cont3,a4,sI,sP,V,V2,ccase,b,Mindx,tip,mval; |
---|
| 629 | number auxc,a1,a2,ex,maxEo,aux; |
---|
| 630 | |
---|
| 631 | 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 |
---|
| 632 | |
---|
| 633 | list oldOlist,oldC,oldt,oldD,oldH,allH; // information of the previous step |
---|
| 634 | |
---|
| 635 | list Olist,C,t,Dstar,center,expI,expP,newJ,maxset; |
---|
| 636 | |
---|
| 637 | list maxvar,auxlist,aux3,auxD,auxolist,auxdiv,auxaux,L,rs,auxgamma,auxg2,aux1; // auxiliary lists |
---|
| 638 | list auxinvlist,newcoef,EL,Ecoaux,Hplus,transH,Hsum,auxset,sumnewH; // auxiliary lists |
---|
| 639 | list auxcoefI,auxcent,center2; |
---|
| 640 | |
---|
| 641 | intvec oldinfobo7,infobo7; |
---|
| 642 | int infaux,leh,leh2,leh3; |
---|
| 643 | |
---|
| 644 | tip=listmb[1]; // It is not used in this procedure, it is used to compute the lcm of the denominators |
---|
| 645 | oldOlist=listmb[2]; |
---|
| 646 | oldC=listmb[3]; |
---|
| 647 | oldt=listmb[4]; // t= resolution function |
---|
| 648 | oldD=listmb[5]; |
---|
| 649 | |
---|
| 650 | oldH=listmb[6]; |
---|
| 651 | allH=listmb[7]; |
---|
| 652 | |
---|
| 653 | oldinfobo7=listmb[8]; // auxiliary intvec, it is used to define BO[7] |
---|
| 654 | |
---|
| 655 | // inicializating lists |
---|
| 656 | Olist=list(); |
---|
| 657 | C=list(); |
---|
| 658 | auxinvlist=list(); |
---|
| 659 | |
---|
| 660 | auxJ[1]=expJ; |
---|
| 661 | rstep=n; // we are in dimension rstep |
---|
| 662 | auxc=c; |
---|
| 663 | cont=1; |
---|
| 664 | |
---|
| 665 | if (Y==0) {D=iniD(n); H=iniD(n); infobo7=-1;} // first center, inicializate previous information |
---|
| 666 | |
---|
| 667 | if (Y!=0 and rstep==n) // In dimension n, D'_n is always of this form |
---|
| 668 | { auxdiv=list0(n); |
---|
| 669 | Dstar[1]=oldD[1]; |
---|
| 670 | |
---|
| 671 | b=size(a); |
---|
| 672 | for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}} |
---|
| 673 | Dstar[1][Y]=aux; |
---|
| 674 | aux=0; |
---|
| 675 | |
---|
| 676 | auxdiv[Y]=oldOlist[1]-oldC[1]; |
---|
| 677 | D[1]=sumlist(Dstar[1],auxdiv);} // list defining D_n |
---|
| 678 | |
---|
| 679 | // computing strict transforms of the exceptional divisors H |
---|
| 680 | |
---|
| 681 | if (Y!=0){transH=iniD(n); |
---|
| 682 | for (i=1;i<=size(oldH);i++){transH[i]=oldH[i]; transH[i][Y]=0;} // Note: size(oldH)<=n |
---|
| 683 | allH[Y]=1;} // transform of |H|=H_nU...UH_1 |
---|
| 684 | |
---|
| 685 | // We put here size(oldH) instead of n because maybe we have not |
---|
| 686 | // calculated all the dimensions in the previous step |
---|
| 687 | |
---|
| 688 | // STARTING THE LOOP |
---|
| 689 | |
---|
| 690 | while (rstep>=1) |
---|
| 691 | { |
---|
| 692 | if (Y!=0 and rstep!=n) // transformation law of D_i for i<n |
---|
| 693 | { |
---|
| 694 | if (cont!=0) // the resolution function did not drop in higher dimensions |
---|
| 695 | { |
---|
| 696 | if (oldt[n-rstep]==a1/a2 and c==oldC[1] and control3==0) |
---|
| 697 | {auxD=list0(n); |
---|
| 698 | auxD[Y]=oldOlist[n-rstep+1]-oldC[n-rstep+1]; |
---|
| 699 | Dstar[n-rstep+1]=oldD[n-rstep+1]; |
---|
| 700 | |
---|
| 701 | for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[n-rstep+1][i];}}} |
---|
| 702 | Dstar[n-rstep+1][Y]=aux; |
---|
| 703 | aux=0; |
---|
| 704 | |
---|
| 705 | D[n-rstep+1]=sumlist(Dstar[n-rstep+1],auxD); |
---|
| 706 | |
---|
| 707 | } |
---|
| 708 | else |
---|
| 709 | {cont=0; |
---|
| 710 | for (j=n-rstep+1;j<=n; j++){D[j]=list0(n);} |
---|
| 711 | } |
---|
| 712 | } |
---|
| 713 | } |
---|
| 714 | |
---|
| 715 | // Factorizing J=M*I |
---|
| 716 | |
---|
| 717 | cont1=0; |
---|
| 718 | for (i=1;i<=n;i++) {if (D[n-rstep+1][i]==0) {cont1=cont1+1;}} // if it fails write: listO(n)[i] |
---|
| 719 | |
---|
| 720 | if (cont1==n) {expI=expJ;} // D[n-rstep+1]=0 (is a list of zeros) |
---|
| 721 | else { |
---|
| 722 | for (i=1;i<=size(expJ);i++) |
---|
| 723 | {rs[i]=size(Coef[i]); |
---|
| 724 | if (rs[i]==2){ aux1=list(); |
---|
| 725 | aux1[1]=reslist(expJ[i][1],D[n-rstep+1]); |
---|
| 726 | aux1[2]=reslist(expJ[i][2],D[n-rstep+1]); |
---|
| 727 | expI[i]=aux1;} // binomial |
---|
| 728 | else {aux1=list(); |
---|
| 729 | aux1[1]=reslist(expJ[i][1],D[n-rstep+1]); |
---|
| 730 | expI[i]=aux1;}} // monomial |
---|
| 731 | } |
---|
| 732 | |
---|
| 733 | // NOTE: coeficients of I = coeficients of J, because I and J differ in a monomial |
---|
| 734 | |
---|
| 735 | // Detecting errors, negative exponents in expI |
---|
| 736 | |
---|
| 737 | sI=size(expI); |
---|
| 738 | |
---|
| 739 | for (i=1;i<=sI;i++) |
---|
| 740 | {rs[i]=size(Coef[i]); |
---|
| 741 | if (rs[i]==2){for (j=1;j<=2;j++){for (l=1;l<=n; l++) |
---|
| 742 | {if (expI[i][j][l]<0) {print("ERROR, the BINOMIAL ideal I has negative components"); |
---|
| 743 | // print("M ideal"); print(D[n-rstep+1]); print(expI); print("dimension"); print(rstep); |
---|
| 744 | // print("previous chart"); print(size(finalchart)); ~; |
---|
| 745 | }}}} |
---|
| 746 | else {for (l=1;l<=n; l++) |
---|
| 747 | {if (expI[i][1][l]<0) {print("ERROR, the MONOMIAL ideal I has negative components"); |
---|
| 748 | // print("M ideal"); print(D[n-rstep+1]); print(expI); print("dimension"); print(rstep); |
---|
| 749 | // print("previous chart"); print(size(finalchart)); ~; |
---|
| 750 | }}} |
---|
| 751 | } |
---|
| 752 | |
---|
| 753 | // Compute the maximal E-order of I |
---|
| 754 | |
---|
| 755 | L=maxEord(Coef,expI,sI,n,flag); |
---|
| 756 | maxEo=L[1]; // E-order of I |
---|
| 757 | |
---|
| 758 | // Inicializating information |
---|
| 759 | |
---|
| 760 | auxolist=maxEo; |
---|
| 761 | a1=maxEo; |
---|
| 762 | a2=auxc; |
---|
| 763 | Olist=Olist+auxolist; // list of new maximal E-orders o_n,o_{n-1},...o_1 |
---|
| 764 | aux3=auxc; |
---|
| 765 | C=C+aux3; // list of new critical values c=c_{n+1},c_{n},...c_2 |
---|
| 766 | |
---|
| 767 | // It is necessary to check if the first coordinate of the invariant has dropped or not |
---|
| 768 | // NOTE: By construction, the first coordinate is always 1 !! |
---|
| 769 | // It has dropped is equivalent to: CURRENT C<c of the previous step |
---|
| 770 | |
---|
| 771 | // Calculate new H, this is done for every dimension |
---|
| 772 | |
---|
| 773 | if (Y!=0){a4=size(oldt); |
---|
| 774 | if (n-rstep+1>a4){cont=0; oldt[n-rstep+1]=0; } // VERIFICAR!!!! |
---|
| 775 | |
---|
| 776 | 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]; |
---|
| 777 | |
---|
| 778 | // we fill now the value for BO[7] |
---|
| 779 | if (oldinfobo7[n-rstep+1]==-1){leh=size(Hhist); |
---|
| 780 | infobo7[n-rstep+1]=Hhist[leh];} // suitable index !!! |
---|
| 781 | else{ infaux=oldinfobo7[n-rstep+1]; |
---|
| 782 | infobo7[n-rstep+1]=infaux;} // the same as the previous step |
---|
| 783 | |
---|
| 784 | } |
---|
| 785 | else { |
---|
| 786 | if (rstep<n) {sumnewH=list0(n); |
---|
| 787 | for (i=1;i<n-rstep+1;i++){sumnewH=sumlist(sumnewH,H[i]);} |
---|
| 788 | H[n-rstep+1]=reslist(allH,sumnewH);} |
---|
| 789 | else {H[n-rstep+1]=allH;} |
---|
| 790 | |
---|
| 791 | // we fill the value for BO[7] too, we complete it at the end if necessary |
---|
| 792 | infobo7[n-rstep+1]=-1; |
---|
| 793 | } |
---|
| 794 | } |
---|
| 795 | |
---|
| 796 | // It is necessary to detect the monomial case AFTER inicializate the information |
---|
| 797 | // OTHERWISE WE WILL HAVE EMPTY COMPONENTS IN THE RESOLUTION FUNCTION |
---|
| 798 | |
---|
| 799 | // If maxEo=0 but maxo!=0 MONOMIAL CASE (because E-Sing(J,c) still !=emptyset) |
---|
| 800 | // If maxEo=0 and maxo=0 then I=1, (real) monomial case, the same case for us |
---|
| 801 | // NOTE THAT IT DOESN'T MATTER IF THERE IS A p-TH POWER OF A HYPERBOLIC EQ, THE E-ORDER IS ZERO ANYWAY |
---|
| 802 | |
---|
| 803 | if (maxEo==0){auxgamma=Gamma(D[n-rstep+1],auxc,n); // Gamma gives (maxlist,gamma,center) |
---|
| 804 | auxg2=auxgamma[3]; |
---|
| 805 | center=center+auxg2; |
---|
| 806 | center=elimrep(center); |
---|
| 807 | auxinvlist=auxgamma[2]; |
---|
| 808 | |
---|
| 809 | // print("gamma"); print(auxg2); |
---|
| 810 | |
---|
| 811 | break;} |
---|
| 812 | |
---|
| 813 | // Calculate P // P=I+M^{o/(c-o)} with weight o |
---|
| 814 | |
---|
| 815 | if (maxEo>=auxc) {expP=expI; Mindx=0;} // The coefficients also remain constant |
---|
| 816 | else {ex=maxEo/(auxc-maxEo); |
---|
| 817 | auxlist=list(); |
---|
| 818 | Mindx=1; |
---|
| 819 | auxlist[1]=multiplylist(D[n-rstep+1],ex); // weighted monomial part: D[n-rstep+1]^ex; |
---|
| 820 | expP=insert(expI,auxlist); // P=I+D[n-rstep+1]^ex; |
---|
| 821 | auxcoefI=Coef; |
---|
| 822 | Coef=insert(Coef,list(1));} // Adding the coefficient for M |
---|
| 823 | |
---|
| 824 | // NOTE: IT IS NECESSARY TO ADD COEFFICIENT 1 TO THE MONOMIAL PART M |
---|
| 825 | // E-ord(P_i)=E-ord(I_i) so to compute the E-order of P_i we can compute E-ord(I_i) |
---|
| 826 | |
---|
| 827 | // Calculate variables of E-maximal contact, ALWAYS WITH RESPECT TO THE IDEAL I !! |
---|
| 828 | |
---|
| 829 | sP=size(expP); // Can be different from size(expI) |
---|
| 830 | |
---|
| 831 | if (Mindx==1){ maxvar=Emaxcont(auxcoefI,expI,sI,n,flag);} |
---|
| 832 | else{ maxvar=Emaxcont(Coef,expP,sP,n,flag);} |
---|
| 833 | |
---|
| 834 | auxc=maxvar[1]; // E-order of P, critical value for the next step, ALSO VALID auxc=maxEo; |
---|
| 835 | if (auxc!=maxEo){print("ERROR, the E-order is not well computed");} |
---|
| 836 | |
---|
| 837 | maxset=maxvar[2]; |
---|
| 838 | |
---|
| 839 | // center=center + maxset; // HACER DESPUES Y A?ADIR SOLO V!!!!!! |
---|
| 840 | // Cleaning the center: eliminating repeated variables |
---|
| 841 | // center=elimrep(center); |
---|
| 842 | |
---|
| 843 | // if (rstep==1) {break;} // Induction finished, is not necessary to compute the rest |
---|
| 844 | |
---|
| 845 | // Calculate Hplus=set of non permissible hypersurfaces |
---|
| 846 | // RESET Hplus if c has dropped or we have eliminated hyperbolic generators |
---|
| 847 | |
---|
| 848 | // ES NECESARIO PONER CONDICION DE SI INVARIANTE BAJA O NO??? SI BAJA HPLUS NO SE USA... |
---|
| 849 | |
---|
| 850 | if (Y==0 or c<oldC[1] or control3==1) {Hplus=list0(n);} |
---|
| 851 | else {Hsum=list0(n); |
---|
| 852 | Hplus=allH; |
---|
| 853 | for (i=1;i<=n-rstep+1;i++){Hsum=sumlist(Hsum,H[i]);} |
---|
| 854 | Hplus=reslist(Hplus,Hsum); // CHEQUEAR QUE NO SALEN -1'S |
---|
| 855 | } |
---|
| 856 | |
---|
| 857 | // Taking into account variables of maxset outside of Hplus (so inside Hminus) |
---|
| 858 | |
---|
| 859 | if (Y==0 or c<oldC[1] or control3==1){V=maxset[1]; // Hplus=0 so any variable is permissible |
---|
| 860 | maxset=delete(maxset,1);} // eliminating this variable V from maxset |
---|
| 861 | else{ |
---|
| 862 | // If the invariant remains constant V comes from the previous step |
---|
| 863 | |
---|
| 864 | if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1]){ |
---|
| 865 | if (Mindx==1){ |
---|
| 866 | //----------------------------USING HPLUS---------------------------------------- |
---|
| 867 | // REMIND THAT IN THIS CASE maxset=HYPERSURFACES OF E-MAXIMAL CONTACT FOR I, INSTEAD OF P |
---|
| 868 | |
---|
| 869 | V2=a[n-rstep+1]; // V can be different from the variable coming from the previous step |
---|
| 870 | // check that V2 belongs to maxset |
---|
| 871 | |
---|
| 872 | for (i=1;i<=size(maxset);i++){ |
---|
| 873 | if (V2==maxset[i]){mval=1; break;} |
---|
| 874 | else{mval=0;} |
---|
| 875 | } |
---|
| 876 | |
---|
| 877 | if (Hplus[V2]==0 and mval==1){V=V2;} // V2 is permissible |
---|
| 878 | else{ |
---|
| 879 | cont2=1; |
---|
| 880 | cont3=1; |
---|
| 881 | auxset=maxset; |
---|
| 882 | while (cont2!=0){mm=auxset[1]; |
---|
| 883 | if (Hplus[mm]!=0) {auxset=delete(auxset,1); cont3=cont3+1;} |
---|
| 884 | // eliminating non permissible variables from maxset |
---|
| 885 | else {cont2=0;}} |
---|
| 886 | V=maxset[cont3]; // first permissible variable |
---|
| 887 | maxset=delete(maxset,cont3); |
---|
| 888 | } |
---|
| 889 | } |
---|
| 890 | |
---|
| 891 | //------------------------------------------------------------------------------- |
---|
| 892 | else{ V=a[n-rstep+1];} |
---|
| 893 | } |
---|
| 894 | else {V=maxset[1]; // Hplus=0 so any variable is permissible |
---|
| 895 | maxset=delete(maxset,1); |
---|
| 896 | } |
---|
| 897 | |
---|
| 898 | } |
---|
| 899 | |
---|
| 900 | |
---|
| 901 | // if (V!=V2 and V2!=0){print(a); print(rstep); print(V); print(V2); print("num cartas"); print(size(finalchart)); ~;} |
---|
| 902 | |
---|
| 903 | V2=0; |
---|
| 904 | |
---|
| 905 | // Adding the new hypersurface of E-maximal contact to the center |
---|
| 906 | |
---|
| 907 | auxcent[1]=V; |
---|
| 908 | |
---|
| 909 | center=center + auxcent; // print("num cartas"); print(size(finalchart)); print(center); if (size(finalchart)==2){~~;} |
---|
| 910 | |
---|
| 911 | auxcent=list(); |
---|
| 912 | |
---|
| 913 | // Cleaning the center: eliminating repeated variables CREO QUE NO HACE FALTA |
---|
| 914 | |
---|
| 915 | center2=elimrep(center); // print(center2); print("-----------"); |
---|
| 916 | |
---|
| 917 | // if (size(center2)!=size(center)){print("MAL");} |
---|
| 918 | |
---|
| 919 | // for (i=1;i<=size(center);i++){if (center2[i]!=center[i]){print("cambia");}} |
---|
| 920 | |
---|
| 921 | |
---|
| 922 | if (rstep==1) {break;} // Induction finished, is not necessary to compute the rest |
---|
| 923 | |
---|
| 924 | |
---|
| 925 | // Calculate Eco=E-Coeff_V(P) where V is a permissible hypersurface which belongs to the center |
---|
| 926 | // Eco can have rational exponents |
---|
| 927 | |
---|
| 928 | Ecoaux=ECoef(Coef,expP,sP,V,auxc,n,flag); |
---|
| 929 | |
---|
| 930 | // SPECIAL CASES: BOLD REGULAR CASE |
---|
| 931 | //-------------------------------------------------------------------- |
---|
| 932 | |
---|
| 933 | if (Ecoaux[3]==1){ // Eco=EMPTY LIST, Eco=0 AS IDEAL |
---|
| 934 | aux1[1]=list0(n); |
---|
| 935 | newJ[1]=aux1; // monomial with zero entries, newJ=1 as ideal |
---|
| 936 | newcoef[1]=list(1); // the new coefficient is only 1 |
---|
| 937 | auxaux=list(); |
---|
| 938 | auxaux[1]=newJ; |
---|
| 939 | auxJ=auxJ+auxaux; // auxJ list of ideals J_i |
---|
| 940 | auxinvlist=1; |
---|
| 941 | break;} |
---|
| 942 | |
---|
| 943 | //----------------------------------------------------------- |
---|
| 944 | // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS |
---|
| 945 | |
---|
| 946 | if (Ecoaux[3]==2 or Ecoaux[3]==3){ // Eco=0 LIST, Eco=1 AS IDEAL |
---|
| 947 | aux1[1]=list0(n); |
---|
| 948 | newJ[1]=aux1; |
---|
| 949 | newcoef[1]=list(1); // print("Strange case happens"); ~; |
---|
| 950 | auxaux=list(); |
---|
| 951 | auxaux[1]=newJ; |
---|
| 952 | auxJ=auxJ + auxaux; // auxJ list of ideals J_i |
---|
| 953 | auxinvlist=1; |
---|
| 954 | break;} |
---|
| 955 | //----------------------------------------------------------- |
---|
| 956 | // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS |
---|
| 957 | |
---|
| 958 | // P=1 THIS CANNOT HAPPEN SINCE P=1 IFF I=1 (or I is equivalent to 1) |
---|
| 959 | // and this is the monomial case, already checked |
---|
| 960 | |
---|
| 961 | if (Ecoaux[3]==4){print("ERROR in ECoef"); break;} |
---|
| 962 | //----------------------------------------------------------- |
---|
| 963 | |
---|
| 964 | // If we are here Ecoaux[3]=0, then continue |
---|
| 965 | |
---|
| 966 | // Filling the list of "ideals", auxJ=J_n,J_{n-1},... |
---|
| 967 | |
---|
| 968 | newJ=Ecoaux[1]; |
---|
| 969 | newcoef=Ecoaux[2]; |
---|
| 970 | |
---|
| 971 | auxJ=insert(auxJ,newJ,n-rstep+1); // newJ is inserted after n-rstep+1 position, so in position n-rstep+2 |
---|
| 972 | |
---|
| 973 | // New input for the loop, if we are here newJ is different from 0 |
---|
| 974 | |
---|
| 975 | expJ=newJ; |
---|
| 976 | Coef=newcoef; |
---|
| 977 | |
---|
| 978 | newJ=list(); |
---|
| 979 | expI=list(); |
---|
| 980 | expP=list(); |
---|
| 981 | rstep=rstep-1; // print(size(auxJ)); |
---|
| 982 | } |
---|
| 983 | |
---|
| 984 | // EXIT LOOP "while" |
---|
| 985 | // we do NOT construct the center as an ideal because WE USE LISTS |
---|
| 986 | |
---|
| 987 | t=dividelist(Olist,C); // resolution function t |
---|
| 988 | |
---|
| 989 | // Complete the intvec infobo7 if necessary |
---|
| 990 | |
---|
| 991 | if (control3==1){infobo7=-1;} // We reset the value after clean hyperbolic equations |
---|
| 992 | leh2=size(Olist); |
---|
| 993 | leh3=size(infobo7); |
---|
| 994 | if (leh3<leh2){for (j=leh3+1;j<=leh2; j++){infobo7[j]=-1;}} |
---|
| 995 | |
---|
| 996 | // Auxiliary list to complete the resolution function in special cases |
---|
| 997 | if (size(auxinvlist)==0) {auxinvlist[1]=0;} |
---|
| 998 | |
---|
| 999 | return(center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7); |
---|
| 1000 | } |
---|
| 1001 | example |
---|
| 1002 | {"EXAMPLE:"; echo = 2; |
---|
| 1003 | ring r = 0,(x(1..4)),dp; |
---|
| 1004 | list flag=identifyvar(); |
---|
| 1005 | ideal J=x(1)^2-x(2)^2*x(3)^5, x(1)*x(3)^3+x(4)^6; |
---|
| 1006 | list Lmb=1,list0(4),list0(4),list0(4),iniD(4),iniD(4),list0(4),-1; |
---|
| 1007 | list L=data(J,2,4); |
---|
| 1008 | list LL=determinecenter(L[1],L[2],2,4,0,0,Lmb,flag,0,-1); // Compute the first center |
---|
| 1009 | LL[1]; // index of variables in the center |
---|
| 1010 | LL[2]; // exponents of ideals J_4,J_3,J_2,J_1 |
---|
| 1011 | LL[3]; // list of orders of J_4,J_3,J_2,J_1 |
---|
| 1012 | LL[4]; // list of critical values |
---|
| 1013 | LL[5]; // components of the resolution function t |
---|
| 1014 | LL[6]; // list of D_4,D_3,D_2,D_1 |
---|
| 1015 | LL[7]; // list of H_4,H_3,H_2,H_1 (exceptional divisors) |
---|
| 1016 | LL[8]; // list of all exceptional divisors acumulated |
---|
| 1017 | LL[9]; // auxiliary invariant |
---|
| 1018 | LL[10]; // intvec pointing out the last step where the function t has dropped |
---|
| 1019 | |
---|
| 1020 | ring r= 0,(x(1..4)),dp; |
---|
| 1021 | list flag=identifyvar(); |
---|
| 1022 | ideal J=x(1)^3-x(2)^2*x(3)^5, x(1)*x(3)^3+x(4)^5; |
---|
| 1023 | list Lmb=2,list0(4),list0(4),list0(4),iniD(4),iniD(4),list0(4),-1; |
---|
| 1024 | list L2=data(J,2,4); |
---|
| 1025 | list L3=determinecenter(L2[1],L2[2],2,4,0,0,Lmb,flag,0,-1); // Example with rational exponents in E-Coeff |
---|
| 1026 | L3[1]; // index of variables in the center |
---|
| 1027 | L3[2]; // exponents of ideals J_4,J_3,J_2,J_1 |
---|
| 1028 | L3[3]; // list of orders of J_4,J_3,J_2,J_1 |
---|
| 1029 | L3[4]; // list of critical values |
---|
| 1030 | L3[5]; // components of the resolution function |
---|
| 1031 | } |
---|
| 1032 | //////////////////////////////////////////////////////// |
---|
| 1033 | // idchart= identity number of the current chart |
---|
| 1034 | // infochart=chart[idchart] information related to the chart to blow up |
---|
| 1035 | // infochart= int parent,int Y,intvec a,list expJ,list Coef, list flag, // NEEDED FOR THE RESOLUTION |
---|
| 1036 | // intvec Hhist, list blwhist, module path, list hipercoef, list hiperexp // NEEDED FOR THE OUTPUT |
---|
| 1037 | |
---|
| 1038 | // NOTE: IT IS NOT NECESSARY TAKE INTO ACCOUNT "y" VARIABLES BECAUSE THE CENTER IS ALREADY GIVEN |
---|
| 1039 | |
---|
| 1040 | proc Blowupcenter(list center,int idchart,int nchart,list infochart,number c,int n,int currentstep) |
---|
| 1041 | "USAGE: Blowupcenter(center,id,nchart,infochart,c,n,cstep); |
---|
| 1042 | center, infochart lists, id, nchart, n, cstep integers, c number |
---|
| 1043 | COMPUTE: The blowing up at the chart IDCHART along the given center |
---|
| 1044 | RETURN: new affine charts and related information, see example |
---|
| 1045 | EXAMPLE: example Blowupcenter; shows an example |
---|
| 1046 | " |
---|
| 1047 | {int num,i,j,k,l,parent,Y,lon,m,m2; |
---|
| 1048 | intvec a,Hhist,auxHhist; |
---|
| 1049 | number auxsum, auxsum2; |
---|
| 1050 | list sons,aux,expJ,blexpJ,blD; |
---|
| 1051 | list auxstep,Coef; |
---|
| 1052 | list auxchart,auxchart1,info,flaglist; |
---|
| 1053 | list auxblwhist,blwhist,hipercoef,hiperexp; |
---|
| 1054 | module auxpath,auxp2; |
---|
| 1055 | |
---|
| 1056 | parent=idchart; |
---|
| 1057 | num=size(center); |
---|
| 1058 | |
---|
| 1059 | // Transform to intvec the list of variables defining the center |
---|
| 1060 | a=center[1]; |
---|
| 1061 | for (i=2;i<=num;i++){a=a,center[i];} |
---|
| 1062 | |
---|
| 1063 | expJ=infochart[4]; |
---|
| 1064 | Coef=infochart[5]; |
---|
| 1065 | flaglist=infochart[6]; |
---|
| 1066 | Hhist=infochart[7]; |
---|
| 1067 | blwhist=infochart[8]; |
---|
| 1068 | auxpath=infochart[9]; |
---|
| 1069 | hipercoef=infochart[10]; |
---|
| 1070 | hiperexp=infochart[11]; |
---|
| 1071 | |
---|
| 1072 | l=size(expJ); |
---|
| 1073 | |
---|
| 1074 | // input for the loop |
---|
| 1075 | blexpJ=expJ; |
---|
| 1076 | |
---|
| 1077 | // making the blowing up in the i-th chart |
---|
| 1078 | for (i=1;i<=num;i++) |
---|
| 1079 | { |
---|
| 1080 | // we assign the current number of charts +1 to the i-th chart |
---|
| 1081 | idchart=nchart+1; |
---|
| 1082 | nchart=nchart+1; |
---|
| 1083 | aux=idchart; |
---|
| 1084 | sons=sons+aux; |
---|
| 1085 | |
---|
| 1086 | auxstep[i]=currentstep+1; |
---|
| 1087 | |
---|
| 1088 | Y=center[i]; |
---|
| 1089 | |
---|
| 1090 | // The blowing up |
---|
| 1091 | |
---|
| 1092 | for (j=1;j<=l;j++){lon=size(Coef[j]); |
---|
| 1093 | if (lon==1){for (m=1;m<=n;m++){for (m2=1;m2<=num;m2++){ |
---|
| 1094 | if (m==center[m2]){auxsum=auxsum+ expJ[j][1][m];}}} |
---|
| 1095 | blexpJ[j][1][Y]=auxsum-c; |
---|
| 1096 | auxsum=0;} // monomial |
---|
| 1097 | else {for (m=1;m<=n;m++){for (m2=1;m2<=num;m2++){ |
---|
| 1098 | if (m==center[m2]){auxsum=auxsum+expJ[j][1][m]; |
---|
| 1099 | auxsum2=auxsum2+expJ[j][2][m];}}} |
---|
| 1100 | blexpJ[j][1][Y]=auxsum-c; |
---|
| 1101 | blexpJ[j][2][Y]=auxsum2-c; |
---|
| 1102 | auxsum=0; auxsum2=0;} // binomial |
---|
| 1103 | } |
---|
| 1104 | |
---|
| 1105 | |
---|
| 1106 | auxHhist=Hhist,Y; // history of the exceptional divisors in this chart |
---|
| 1107 | auxblwhist=tradblwup(blwhist,n,Y,a,num); // history of the blow ups in this chart |
---|
| 1108 | |
---|
| 1109 | auxp2=auxpath,[parent,i]; |
---|
| 1110 | |
---|
| 1111 | auxchart1=parent,Y,a,blexpJ,Coef,flaglist,auxHhist,auxblwhist,auxp2,hipercoef,hiperexp; |
---|
| 1112 | |
---|
| 1113 | // Coef, flaglist are not modified after the blowing-up, the hyperbolic information is the same as in the parent chart |
---|
| 1114 | |
---|
| 1115 | auxchart[i]=auxchart1; |
---|
| 1116 | |
---|
| 1117 | // Inicializating the exponents of J for the next chart |
---|
| 1118 | |
---|
| 1119 | blexpJ=expJ; |
---|
| 1120 | } |
---|
| 1121 | // end of the loop |
---|
| 1122 | |
---|
| 1123 | // we add its sons to the current chart |
---|
| 1124 | infochart=infochart+sons; |
---|
| 1125 | info[1]=infochart; |
---|
| 1126 | |
---|
| 1127 | return(info,auxchart,nchart,auxstep,num); |
---|
| 1128 | } |
---|
| 1129 | example |
---|
| 1130 | {"EXAMPLE:"; echo = 2; |
---|
| 1131 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp; |
---|
| 1132 | list flag=identifyvar(); |
---|
| 1133 | 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; |
---|
| 1134 | list Lmb=2,list0(7),list0(7),list0(7),iniD(7),iniD(7),list0(7),-1; |
---|
| 1135 | list L=data(J,3,7); |
---|
| 1136 | list L2=determinecenter(L[1],L[2],2,7,0,0,Lmb,flag,0,-1); // Computing the center |
---|
| 1137 | module auxpath=[0,-1]; |
---|
| 1138 | list infochart=0,0,0,L[2],L[1],flag,0,list0(7),auxpath,list(),list(); |
---|
| 1139 | list L3=Blowupcenter(L2[1],1,1,infochart,2,7,0); |
---|
| 1140 | L3[1]; // current chart (parent,Y,center,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp) with sons: [12],...,[16] |
---|
| 1141 | L3[2][1]; // information of its first son, write L3[2][2],...,L3[2][5] to see the other sons |
---|
| 1142 | L3[3]; // current number of charts |
---|
| 1143 | L3[4]; // step/level associated to each son |
---|
| 1144 | L3[5]; // number of variables in the center |
---|
| 1145 | } |
---|
| 1146 | ////////////////////////////////////////////////////////////// |
---|
| 1147 | |
---|
| 1148 | proc tradblwup(list blwhist,int n,int Y,intvec a,int num) |
---|
| 1149 | "Internal procedure - no help and no example available |
---|
| 1150 | " |
---|
| 1151 | { |
---|
| 1152 | int i,j,blwnew; |
---|
| 1153 | intvec aux,aux2; |
---|
| 1154 | |
---|
| 1155 | for (j=1;j<=n;j++){ |
---|
| 1156 | for (i=1;i<=num;i++){ |
---|
| 1157 | if (j==a[i] and a[i]!=Y){blwnew=Y; break;} |
---|
| 1158 | else {blwnew=0;} |
---|
| 1159 | } |
---|
| 1160 | aux=blwhist[j]; |
---|
| 1161 | aux2=aux,blwnew; |
---|
| 1162 | blwhist[j]=aux2; |
---|
| 1163 | } |
---|
| 1164 | return(blwhist); |
---|
| 1165 | } |
---|
| 1166 | ////////////////////////////////////////////////////////////// |
---|
| 1167 | // It is called only when Eord(J)=0, and J!=1 it is checked inside |
---|
| 1168 | // SO IT IS CALLED AFTER: maxEord(Coef,expJ,sJ,n,flaglist); --> gives (max E-order,sums) |
---|
| 1169 | |
---|
| 1170 | proc Nonhyp(list Coef,list expJ,int sJ,int n,list flaglist,list sums) |
---|
| 1171 | "USAGE: Nonhyp(Coef,expJ,sJ,n,flaglist,sums); |
---|
| 1172 | Coef, expJ, flaglist, sums lists, sJ, n integers |
---|
| 1173 | COMPUTE: The "ideal" generated by the non hyperbolic generators of J |
---|
| 1174 | RETURN: lists with the following information |
---|
| 1175 | newcoef,newJ: coefficients and exponents of the non hyperbolic generators |
---|
| 1176 | totalhyp,totalgen: coefficients and exponents of the hyperbolic generators |
---|
| 1177 | flaglist: new list saying status of variables |
---|
| 1178 | NOTE: the basering r is supposed to be a polynomial ring K[x,y], |
---|
| 1179 | in fact, we work in a localization of K[x,y], of type K[x,y]_y with y invertible variables. |
---|
| 1180 | EXAMPLE: example Nonhyp; shows an example |
---|
| 1181 | " |
---|
| 1182 | { |
---|
| 1183 | int i,j,k,h,lon,lon2,cont; |
---|
| 1184 | number eordcontrol; |
---|
| 1185 | list genhyp,listgen,listid,posnumJ,newJ,newcoef,hypcoef,hyp,aux1,aux2,aux3,aux,midlist; |
---|
| 1186 | list totalhyp,totalgen; |
---|
| 1187 | |
---|
| 1188 | eordcontrol=0; |
---|
| 1189 | |
---|
| 1190 | while (eordcontrol==0 and sJ!=0) |
---|
| 1191 | { |
---|
| 1192 | |
---|
| 1193 | // Give a positional number/flag to each generator of expJ |
---|
| 1194 | |
---|
| 1195 | for (i=1;i<=sJ; i++){listgen=expJ[i]; listid=i; posnumJ[i]=listgen+listid; } |
---|
| 1196 | |
---|
| 1197 | // Select the non hyperbolic and hyperbolic generators |
---|
| 1198 | |
---|
| 1199 | for (j=1;j<=sJ; j++){lon=size(Coef[j]); |
---|
| 1200 | if (lon==1){ |
---|
| 1201 | |
---|
| 1202 | // IS NOT NECESSARY TO CHECK IF THERE EXIST A MONOMIAL WITH ONLY UNITS, ALREADY DONE!! |
---|
| 1203 | |
---|
| 1204 | aux1=aux1+posnumJ[j]; |
---|
| 1205 | aux3=list(); |
---|
| 1206 | aux3[1]=expJ[j]; |
---|
| 1207 | newJ=newJ+aux3; |
---|
| 1208 | aux3[1]=Coef[j]; |
---|
| 1209 | newcoef=newcoef+aux3; |
---|
| 1210 | } |
---|
| 1211 | |
---|
| 1212 | else{ // CHECKING BINOMIALS, ONE TERM WITH E-ORDER ZERO GIVES HYPERBOLIC EQ |
---|
| 1213 | |
---|
| 1214 | if (sums[j][1]==0 or sums[j][2]==0){aux2=aux2+posnumJ[j]; |
---|
| 1215 | aux3=list(); |
---|
| 1216 | aux3[1]=expJ[j]; |
---|
| 1217 | genhyp=genhyp+aux3; |
---|
| 1218 | aux3[1]=Coef[j]; |
---|
| 1219 | hypcoef=hypcoef+aux3; |
---|
| 1220 | if (sums[j][1]==0){aux3[1]=expJ[j][2]; hyp=hyp+aux3;} |
---|
| 1221 | if (sums[j][2]==0){aux3[1]=expJ[j][1]; hyp=hyp+aux3;} |
---|
| 1222 | } |
---|
| 1223 | else {aux1=aux1+posnumJ[j]; |
---|
| 1224 | aux3=list(); |
---|
| 1225 | aux3[1]=expJ[j]; |
---|
| 1226 | newJ=newJ+aux3; |
---|
| 1227 | aux3[1]=Coef[j]; |
---|
| 1228 | newcoef=newcoef+aux3;} |
---|
| 1229 | |
---|
| 1230 | } |
---|
| 1231 | } |
---|
| 1232 | |
---|
| 1233 | // NOTE: aux1 and aux2 are no needed right now! |
---|
| 1234 | |
---|
| 1235 | // Identify new y variables, that is, x variables in the monomials contained in hyp |
---|
| 1236 | |
---|
| 1237 | h=size(hyp); |
---|
| 1238 | |
---|
| 1239 | for (k=1;k<=h; k++){ for(i=1;i<=n; i++){ if (hyp[k][i]!=0 and flaglist[i]==0) {flaglist[i]=1;}}} |
---|
| 1240 | |
---|
| 1241 | // To replace x by y IT IS NECESSARY TO CHANGE THE BASERING!!! We change only the list flaglist |
---|
| 1242 | |
---|
| 1243 | // CHECK IF THE IDEAL IS ALREADY GENERATED BY MONOMIALS, in this case |
---|
| 1244 | // WE HAVE FINISHED THE E-RESOLUTION PART, J GENERATED BY MONOMIALS AND HYPERBOLIC EQS |
---|
| 1245 | |
---|
| 1246 | cont=0; |
---|
| 1247 | lon2=size(newJ); |
---|
| 1248 | for (j=1;j<=lon2; j++){if (size(newJ[j])==1){cont=cont+1;}} |
---|
| 1249 | |
---|
| 1250 | if (cont==lon2){newcoef=list(); |
---|
| 1251 | newJ=list(); |
---|
| 1252 | totalgen=totalgen+genhyp; |
---|
| 1253 | totalhyp=totalhyp+hypcoef; |
---|
| 1254 | break;} |
---|
| 1255 | |
---|
| 1256 | // CHECK IF THERE ARE MORE HYPERBOLIC EQUATIONS AFTER UPDATE THE FLAG LIST |
---|
| 1257 | // CHECK THE MAXIMAL E-ORDER AGAIN |
---|
| 1258 | |
---|
| 1259 | if (lon2==0){ // we are in the previous case, newJ=empty list, save values and exit |
---|
| 1260 | |
---|
| 1261 | totalgen=totalgen+genhyp; |
---|
| 1262 | totalhyp=totalhyp+hypcoef; |
---|
| 1263 | break; |
---|
| 1264 | } |
---|
| 1265 | |
---|
| 1266 | midlist=maxEord(newcoef,newJ,lon2,n,flaglist); |
---|
| 1267 | |
---|
| 1268 | eordcontrol=midlist[1]; |
---|
| 1269 | |
---|
| 1270 | if (eordcontrol==0){ // new input for the loop |
---|
| 1271 | Coef=newcoef; |
---|
| 1272 | expJ=newJ; |
---|
| 1273 | sJ=lon2; |
---|
| 1274 | sums=midlist[2]; // flaglist is already updated |
---|
| 1275 | |
---|
| 1276 | totalgen=totalgen+genhyp; |
---|
| 1277 | totalhyp=totalhyp+hypcoef; |
---|
| 1278 | |
---|
| 1279 | hypcoef=list(); |
---|
| 1280 | genhyp=list(); |
---|
| 1281 | |
---|
| 1282 | newJ=list(); |
---|
| 1283 | newcoef=list(); |
---|
| 1284 | } |
---|
| 1285 | else{ // If the process is already finished we save the values and exit |
---|
| 1286 | |
---|
| 1287 | totalgen=totalgen+genhyp; |
---|
| 1288 | totalhyp=totalhyp+hypcoef; |
---|
| 1289 | } |
---|
| 1290 | |
---|
| 1291 | } // closing while |
---|
| 1292 | |
---|
| 1293 | return(newcoef,newJ,totalhyp,totalgen,flaglist); |
---|
| 1294 | } |
---|
| 1295 | example |
---|
| 1296 | {"EXAMPLE:"; echo = 2; |
---|
| 1297 | ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp; |
---|
| 1298 | list flag=identifyvar(); // List giving flag=1 to invertible variables: y(2),y(4) |
---|
| 1299 | 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; |
---|
| 1300 | list L=data(J,3,7); |
---|
| 1301 | list L2=maxEord(L[1],L[2],3,7,flag); |
---|
| 1302 | L2[1]; // Maximum E-order |
---|
| 1303 | list New=Nonhyp(L[1],L[2],3,7,flag,L2[2]); |
---|
| 1304 | New[1]; // Coefficients of the non hyperbolic part |
---|
| 1305 | New[2]; // Exponents of the non hyperbolic part |
---|
| 1306 | New[3]; // Coefficients of the hyperbolic part |
---|
| 1307 | New[4]; // New hyperbolic equations |
---|
| 1308 | New[5]; // New list giving flag=1 to invertible variables: y(2),y(4),y(5) |
---|
| 1309 | |
---|
| 1310 | ring r = 0,(x(1..4)),dp; |
---|
| 1311 | list flag=identifyvar(); |
---|
| 1312 | ideal J=1-x(1)^5*x(2)^2*x(3)^5, x(1)^2*x(3)^3+x(1)^4*x(4)^6; |
---|
| 1313 | list L=data(J,2,4); |
---|
| 1314 | list L2=maxEord(L[1],L[2],2,4,flag); |
---|
| 1315 | L2[1]; // Maximum E-order |
---|
| 1316 | list New=Nonhyp(L[1],L[2],2,4,flag,L2[2]); |
---|
| 1317 | New; |
---|
| 1318 | |
---|
| 1319 | } |
---|
| 1320 | ////////////////////////////////////////////////////////////// |
---|
| 1321 | |
---|
| 1322 | proc calculateI(list Coef,list expJ,number c,int n,int Y,intvec a,number oldordI,list oldD) |
---|
| 1323 | "USAGE: calculateI(Coef,expJ,c,n,Y,a,b,D); |
---|
| 1324 | Coef, expJ, D lists, c, b numbers, n,Y integers, a intvec |
---|
| 1325 | RETURN: ideal I, non monomial part of J |
---|
| 1326 | EXAMPLE: example calculateI; shows an example |
---|
| 1327 | " |
---|
| 1328 | { |
---|
| 1329 | int i,cont1,b,j; |
---|
| 1330 | number EordI,aux; |
---|
| 1331 | list D,L,expI; |
---|
| 1332 | list auxdiv,Dstar,aux1,rs; |
---|
| 1333 | |
---|
| 1334 | // WE NEED THE MONOMIAL PART, BUT ONLY IN DIMENSION n |
---|
| 1335 | |
---|
| 1336 | auxdiv=list0(n); |
---|
| 1337 | auxdiv[Y]=oldordI-c; |
---|
| 1338 | Dstar[1]=oldD[1]; |
---|
| 1339 | |
---|
| 1340 | b=size(a); |
---|
| 1341 | for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}} |
---|
| 1342 | Dstar[1][Y]=aux; |
---|
| 1343 | aux=0; |
---|
| 1344 | |
---|
| 1345 | D[1]=sumlist(Dstar[1],auxdiv); |
---|
| 1346 | |
---|
| 1347 | cont1=0; |
---|
| 1348 | for (i=1;i<=n;i++) {if (D[1][i]==0) {cont1=cont1+1;}} // if it fails write listO(n)[i] |
---|
| 1349 | |
---|
| 1350 | if (cont1==n) {expI=expJ;} |
---|
| 1351 | else { |
---|
| 1352 | for (i=1;i<=size(expJ);i++) |
---|
| 1353 | {rs[i]=size(Coef[i]); |
---|
| 1354 | if (rs[i]==2){ aux1=list(); |
---|
| 1355 | aux1[1]=reslist(expJ[i][1],D[1]); |
---|
| 1356 | aux1[2]=reslist(expJ[i][2],D[1]); |
---|
| 1357 | expI[i]=aux1;} // binomial |
---|
| 1358 | else {aux1=list(); |
---|
| 1359 | aux1[1]=reslist(expJ[i][1],D[1]); |
---|
| 1360 | expI[i]=aux1;}} // monomial |
---|
| 1361 | } |
---|
| 1362 | |
---|
| 1363 | return(expI); |
---|
| 1364 | } |
---|
| 1365 | example |
---|
| 1366 | {"EXAMPLE:"; echo = 2; |
---|
| 1367 | ring r = 0,(x(1..3)),dp; |
---|
| 1368 | list flag=identifyvar(); |
---|
| 1369 | ideal J=x(1)^4*x(2)^2, x(3)^3; |
---|
| 1370 | list Lmb=1,list0(3),list0(3),list0(3),iniD(3),iniD(3),list0(3),-1; |
---|
| 1371 | list L=data(J,2,3); |
---|
| 1372 | list LL=determinecenter(L[1],L[2],3,3,0,0,Lmb,flag,0,-1); // Calculate the center |
---|
| 1373 | module auxpath=[0,-1]; |
---|
| 1374 | list infochart=0,0,0,L[2],L[1],flag,0,list0(3),auxpath,list(),list(); |
---|
| 1375 | list L3=Blowupcenter(LL[1],1,1,infochart,3,3,0); // blowing-up and looking to the x(3) chart |
---|
| 1376 | calculateI(L3[2][1][5],L3[2][1][4],3,3,3,L3[2][1][3],3,iniD(3)); // (I_3) |
---|
| 1377 | // looking to the x(1) chart |
---|
| 1378 | calculateI(L3[2][2][5],L3[2][2][4],3,3,1,L3[2][2][3],3,iniD(3)); // (I_3) |
---|
| 1379 | } |
---|
| 1380 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1381 | // // |
---|
| 1382 | // E-RESOLUTION: Eresol(J) subroutine computing the E-resolution of J, char 0 // |
---|
| 1383 | // // |
---|
| 1384 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1385 | |
---|
| 1386 | proc Eresol(ideal J) |
---|
| 1387 | "USAGE: Eresol(J); J ideal |
---|
| 1388 | RETURN: The E-resolution of singularities of J in terms of the affine charts, see example |
---|
| 1389 | EXAMPLE: example Eresol; shows an example |
---|
| 1390 | " |
---|
| 1391 | {int i,n,k,idchart,nchart,parent,Y,oldid,tnum,s,cont,control,control2,control3,cont2,val,rs2,l,cont3,tip; |
---|
| 1392 | intvec a,Hhist; |
---|
| 1393 | number c,EordJ,EordI,oldordI; |
---|
| 1394 | list L,LL,oldD,t,auxL,finalchart,chart,auxchart,newL,auxp,auxfchart,L2; |
---|
| 1395 | list Coef,expJ,expI,sons,oldOlist,oldC,oldt,oldH,allH,auxordJ,auxordI,auxmb,mobile,invariant; |
---|
| 1396 | list step,nsons,auxinv,extraL,totalinv,auxsum; |
---|
| 1397 | string empstring; |
---|
| 1398 | module auxpath; |
---|
| 1399 | |
---|
| 1400 | // ADDED LATER |
---|
| 1401 | |
---|
| 1402 | list flag,newflag,blwhist,hipercoef,hiperexp,hipercoefson,hiperexpson; |
---|
| 1403 | intvec infobo7; |
---|
| 1404 | |
---|
| 1405 | export finalchart; |
---|
| 1406 | // export nsons; |
---|
| 1407 | // export tnum; |
---|
| 1408 | // export nchart; |
---|
| 1409 | // export step; |
---|
| 1410 | export invariant; |
---|
| 1411 | export auxinv; |
---|
| 1412 | export mobile; |
---|
| 1413 | |
---|
| 1414 | n=nvars(basering); |
---|
| 1415 | flag=identifyvar(); |
---|
| 1416 | |
---|
| 1417 | k=size(J); |
---|
| 1418 | // Checking input data |
---|
| 1419 | if (inidata(J,k)==0){return("This library only works for binomial ideals.");} |
---|
| 1420 | |
---|
| 1421 | idchart=1; |
---|
| 1422 | nchart=1; |
---|
| 1423 | parent=0; |
---|
| 1424 | step=0; |
---|
| 1425 | control=0; |
---|
| 1426 | control2=0; |
---|
| 1427 | control3=0; |
---|
| 1428 | |
---|
| 1429 | // Translate the input ideal to a list |
---|
| 1430 | auxL=data(J,k,n); // data gives (Coef,Exp) |
---|
| 1431 | |
---|
| 1432 | // THEREAFTER WE WORK ALL THE TIME WITH LISTS |
---|
| 1433 | |
---|
| 1434 | L=maxEord(auxL[1],auxL[2],k,n,flag); // gives (max E-ord, sums) |
---|
| 1435 | EordJ=L[1]; |
---|
| 1436 | |
---|
| 1437 | // before the first blow up I=J |
---|
| 1438 | EordI=EordJ; |
---|
| 1439 | |
---|
| 1440 | // main loop AT EACH CHART WE MUST INICIALIZATE ALL THE VALUES AND |
---|
| 1441 | // CONSTRUCT THE FIRST CHART chart[1] BEFORE THE LOOP |
---|
| 1442 | |
---|
| 1443 | // at the first step, before the blow up, there are no exceptional divisors, Y=0 |
---|
| 1444 | Y=0; |
---|
| 1445 | expJ=auxL[2]; |
---|
| 1446 | Coef=auxL[1]; |
---|
| 1447 | Hhist=0; |
---|
| 1448 | blwhist=list0(n); |
---|
| 1449 | auxpath=[0,-1]; |
---|
| 1450 | hipercoef=list(); // this is for the first chart |
---|
| 1451 | hiperexp=list(); |
---|
| 1452 | auxp=parent,Y,a,expJ,Coef,flag,Hhist,blwhist,auxpath,hipercoef,hiperexp; |
---|
| 1453 | chart[1]=auxp; // information of the first chart |
---|
| 1454 | |
---|
| 1455 | tip=1; |
---|
| 1456 | oldOlist=list0(n); |
---|
| 1457 | oldC=list0(n); |
---|
| 1458 | oldC[1]=EordJ; // non necessary here |
---|
| 1459 | c=EordJ; // the value c is given by the previous step |
---|
| 1460 | oldt=list0(n); |
---|
| 1461 | oldD=iniD(n); |
---|
| 1462 | oldH=iniD(n); |
---|
| 1463 | allH=list0(n); |
---|
| 1464 | |
---|
| 1465 | for (i=1;i<=n;i++){infobo7[i]=-1;} |
---|
| 1466 | |
---|
| 1467 | auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; |
---|
| 1468 | mobile[1]=auxmb; // mobile corresponding to the first chart |
---|
| 1469 | auxinv[1]=list(0); |
---|
| 1470 | |
---|
| 1471 | // NOTE: oldC[1] is the value c to classify the chart in one of the next cases |
---|
| 1472 | |
---|
| 1473 | // HERE BEGIN THE LOOP |
---|
| 1474 | |
---|
| 1475 | while (idchart<=nchart) // WE PROCEED WHILE THERE EXIST UNSOLVED CHARTS |
---|
| 1476 | { |
---|
| 1477 | if (idchart!=1) // WE ARE NOT IN THE FIRST CHART, INICIALIZATE ALL THE VALUES |
---|
| 1478 | { |
---|
| 1479 | |
---|
| 1480 | parent=chart[idchart][1]; |
---|
| 1481 | Y=chart[idchart][2]; |
---|
| 1482 | a=chart[idchart][3]; |
---|
| 1483 | expJ=chart[idchart][4]; |
---|
| 1484 | Coef=chart[idchart][5]; |
---|
| 1485 | flag=chart[idchart][6]; |
---|
| 1486 | Hhist=chart[idchart][7]; // it is not necessary for the computations |
---|
| 1487 | blwhist=chart[idchart][8]; |
---|
| 1488 | auxpath=chart[idchart][9]; |
---|
| 1489 | hipercoef=chart[idchart][10]; |
---|
| 1490 | hiperexp=chart[idchart][11]; |
---|
| 1491 | |
---|
| 1492 | k=size(Coef); // IT IS NECESSARY TO COMPUTE IT BECAUSE IT DECREASES IF THERE ARE HYPERBOLIC EQS |
---|
| 1493 | |
---|
| 1494 | auxordJ=maxEord(Coef,expJ,k,n,flag); |
---|
| 1495 | EordJ=auxordJ[1]; |
---|
| 1496 | |
---|
| 1497 | if (control==0){c=mobile[parent+1][3][1];} // we keep c from the last step |
---|
| 1498 | else {c=EordJ; control=0; } // we reset the value of c |
---|
| 1499 | |
---|
| 1500 | if (control2==1){c=EordJ; control2=0; control3=1;} // we reset the value of c |
---|
| 1501 | |
---|
| 1502 | // NOTE: oldC[1] is the value c to classify the chart in one of the next cases |
---|
| 1503 | |
---|
| 1504 | } |
---|
| 1505 | |
---|
| 1506 | // The E-order must be computed here |
---|
| 1507 | |
---|
| 1508 | oldid=idchart; |
---|
| 1509 | |
---|
| 1510 | if (EordJ<0) {print("ERROR in J in chart"); print(idchart); ~; break;} |
---|
| 1511 | |
---|
| 1512 | |
---|
| 1513 | //------------------------------------------------------------- |
---|
| 1514 | // CASE J=1, if we reset c, can happen Eord=c=0 |
---|
| 1515 | |
---|
| 1516 | // or if there are hyperbolic equations at the beginning!!! A?ADIR!!!! |
---|
| 1517 | |
---|
| 1518 | // if (EordJ==0){auxfchart[1]=chart[idchart]; // WE HAVE FINISHED |
---|
| 1519 | // finalchart=finalchart+auxfchart; |
---|
| 1520 | // empstring="#"; print("reset c and Eord=c=0"); print(idchart); |
---|
| 1521 | // invariant[idchart]=empstring; |
---|
| 1522 | // auxinv[idchart]=list(0); |
---|
| 1523 | // nsons[idchart]=0; |
---|
| 1524 | // idchart=idchart+1;} |
---|
| 1525 | |
---|
| 1526 | |
---|
| 1527 | //---------------------------------------------------------------------- |
---|
| 1528 | if (EordJ>=c and EordJ!=0) // subroutine: E-RESOLUTION OF PAIRS |
---|
| 1529 | { |
---|
| 1530 | if (parent>0) |
---|
| 1531 | { LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,chart[parent][7]); } |
---|
| 1532 | else { LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,Hhist); } |
---|
| 1533 | |
---|
| 1534 | // determinecenter gives (center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7) |
---|
| 1535 | |
---|
| 1536 | // save current values, before the blow up |
---|
| 1537 | oldOlist=LL[3]; |
---|
| 1538 | tip=computemcm(oldOlist); |
---|
| 1539 | oldC=LL[4]; |
---|
| 1540 | oldt=LL[5]; |
---|
| 1541 | oldD=LL[6]; |
---|
| 1542 | oldH=LL[7]; |
---|
| 1543 | allH=LL[8]; |
---|
| 1544 | auxinv[idchart]=LL[9]; |
---|
| 1545 | infobo7=LL[10]; |
---|
| 1546 | |
---|
| 1547 | auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; |
---|
| 1548 | mobile[idchart+1]=auxmb; |
---|
| 1549 | invariant[idchart]=oldt; |
---|
| 1550 | |
---|
| 1551 | newL=Blowupcenter(LL[1],idchart,nchart,chart[idchart],c,n,step[idchart]); |
---|
| 1552 | |
---|
| 1553 | // Blowupcenter gives (info,auxchart,nchart,auxstep,num) |
---|
| 1554 | |
---|
| 1555 | // IMPORTANT: ADD THE NEW CHARTS AFTER EACH BLOW UP, IN ORDER TO KEEP THEM CORRECTLY |
---|
| 1556 | |
---|
| 1557 | step=step+newL[4]; |
---|
| 1558 | nsons[idchart]=newL[5]; |
---|
| 1559 | |
---|
| 1560 | chart=chart+newL[2]; |
---|
| 1561 | finalchart=finalchart+newL[1]; |
---|
| 1562 | |
---|
| 1563 | // new input for the loop |
---|
| 1564 | |
---|
| 1565 | idchart=idchart+1; |
---|
| 1566 | nchart=newL[3]; |
---|
| 1567 | |
---|
| 1568 | control3=0; |
---|
| 1569 | |
---|
| 1570 | } // END OF CASE EordJ>=c |
---|
| 1571 | //--------------------------------------------------------------------- |
---|
| 1572 | |
---|
| 1573 | else{ |
---|
| 1574 | |
---|
| 1575 | // compute EordI=max E-order(I) |
---|
| 1576 | |
---|
| 1577 | expI=calculateI(Coef,expJ,c,n,Y,a,mobile[parent+1][2][1],mobile[parent+1][5]); |
---|
| 1578 | k=size(expJ); // probably non necessary |
---|
| 1579 | auxordI=maxEord(Coef,expI,k,n,flag); |
---|
| 1580 | EordI=auxordI[1]; |
---|
| 1581 | auxsum=auxordI[2]; |
---|
| 1582 | |
---|
| 1583 | // CASE EordI>0 DROP c AND CONTINUE |
---|
| 1584 | |
---|
| 1585 | if (EordI>0){idchart=idchart; // keep the chart and back to the main loop while, dropping the value of c |
---|
| 1586 | control=1;} |
---|
| 1587 | else{ // EordI=0, so check if I=1 or not |
---|
| 1588 | |
---|
| 1589 | cont2=0; // If cont2=val then all the entries of expI are zero!! |
---|
| 1590 | val=0; |
---|
| 1591 | |
---|
| 1592 | for (i=1;i<=k;i++) {rs2=size(Coef[i]); |
---|
| 1593 | if (rs2==1){if (auxsum[i][1]==0){cont2=val; break;} // THERE EXIST A MONOMIAL WITH ONLY UNITS |
---|
| 1594 | |
---|
| 1595 | val=val+n; // monomials |
---|
| 1596 | for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;}} |
---|
| 1597 | } |
---|
| 1598 | else{val=val+(2*n); // binomials |
---|
| 1599 | for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;} |
---|
| 1600 | if (expI[i][2][l]==0) {cont2=cont2+1;}} |
---|
| 1601 | } |
---|
| 1602 | } |
---|
| 1603 | |
---|
| 1604 | |
---|
| 1605 | // CASE EordI==0 AND I=1 THIS CHART IS DONE, FINISH |
---|
| 1606 | |
---|
| 1607 | // NOTE: THIS CASE IS NOT MONOMIAL BECAUSE E-Sing(J,c) is empty |
---|
| 1608 | |
---|
| 1609 | if (cont2==val){auxfchart[1]=chart[idchart]; |
---|
| 1610 | finalchart=finalchart+auxfchart; |
---|
| 1611 | empstring="#"; |
---|
| 1612 | invariant[idchart]=empstring; |
---|
| 1613 | auxinv[idchart]=list(0); |
---|
| 1614 | nsons[idchart]=0; |
---|
| 1615 | |
---|
| 1616 | // information for the mobile |
---|
| 1617 | tip=1; |
---|
| 1618 | oldOlist=list(0); |
---|
| 1619 | oldC=list(0); |
---|
| 1620 | oldt=list(0); |
---|
| 1621 | oldD=list(0); |
---|
| 1622 | oldH=list(0); |
---|
| 1623 | allH=list(0); // the value of the parent + the new one |
---|
| 1624 | infobo7=-1; |
---|
| 1625 | |
---|
| 1626 | auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; |
---|
| 1627 | mobile[idchart+1]=auxmb; |
---|
| 1628 | |
---|
| 1629 | idchart=idchart+1;} |
---|
| 1630 | |
---|
| 1631 | else{ // CASE EordI==0 AND I!=1 --> HYPERBOLIC EQUATIONS |
---|
| 1632 | |
---|
| 1633 | // COMPUTE THE IDEAL OF NON HYPERBOLIC ELEMENTS |
---|
| 1634 | |
---|
| 1635 | extraL=Nonhyp(Coef,expI,k,n,flag,auxordI[2]); // gives (newcoef,newI,hypcoef,genhyp,flaglist) |
---|
| 1636 | |
---|
| 1637 | // CHECK IF ALL THE VARIABLES ARE ALREADY INVERTIBLE |
---|
| 1638 | |
---|
| 1639 | newflag=extraL[5]; |
---|
| 1640 | chart[idchart][6]=extraL[5]; // update the status of variables |
---|
| 1641 | |
---|
| 1642 | cont3=0; |
---|
| 1643 | for (i=1;i<=n;i++){if (newflag[i]==1){cont3=cont3+1;}} |
---|
| 1644 | |
---|
| 1645 | if (cont3==n){ // ALL THE VARIABLES ARE INVERTIBLE |
---|
| 1646 | auxfchart[1]=chart[idchart]; |
---|
| 1647 | finalchart=finalchart+auxfchart; |
---|
| 1648 | empstring="@"; |
---|
| 1649 | invariant[idchart]=empstring; |
---|
| 1650 | auxinv[idchart]=list(0); |
---|
| 1651 | nsons[idchart]=0; |
---|
| 1652 | |
---|
| 1653 | // information for the mobile |
---|
| 1654 | tip=1; |
---|
| 1655 | oldOlist=list(0); |
---|
| 1656 | oldC=list(0); |
---|
| 1657 | oldt=list(0); |
---|
| 1658 | oldD=list(0); |
---|
| 1659 | oldH=list(0); |
---|
| 1660 | allH=list(0); |
---|
| 1661 | infobo7=-1; |
---|
| 1662 | |
---|
| 1663 | auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; |
---|
| 1664 | mobile[idchart+1]=auxmb; |
---|
| 1665 | |
---|
| 1666 | idchart=idchart+1;} |
---|
| 1667 | else{ // OTHERWISE, CONTINUE CHEKING IF newI=0 or not |
---|
| 1668 | |
---|
| 1669 | Coef=extraL[1]; |
---|
| 1670 | expI=extraL[2]; |
---|
| 1671 | |
---|
| 1672 | hipercoefson=extraL[3]; // Information about hyperbolic generators |
---|
| 1673 | hiperexpson=extraL[4]; |
---|
| 1674 | |
---|
| 1675 | k=size(expI); |
---|
| 1676 | |
---|
| 1677 | if (k==0){auxfchart[1]=chart[idchart]; // WE HAVE FINISHED |
---|
| 1678 | finalchart=finalchart+auxfchart; |
---|
| 1679 | empstring="#"; // no more non-hyperbolic generators in this chart |
---|
| 1680 | invariant[idchart]=empstring; |
---|
| 1681 | auxinv[idchart]=list(0); |
---|
| 1682 | nsons[idchart]=0; |
---|
| 1683 | |
---|
| 1684 | // information for the mobile |
---|
| 1685 | tip=1; |
---|
| 1686 | oldOlist=list(0); |
---|
| 1687 | oldC=list(0); |
---|
| 1688 | oldt=list(0); |
---|
| 1689 | oldD=list(0); |
---|
| 1690 | oldH=list(0); |
---|
| 1691 | allH=list(0); |
---|
| 1692 | infobo7=-1; |
---|
| 1693 | |
---|
| 1694 | auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; |
---|
| 1695 | mobile[idchart+1]=auxmb; |
---|
| 1696 | |
---|
| 1697 | idchart=idchart+1;} |
---|
| 1698 | |
---|
| 1699 | else{ // CONTINUE WITH THE IDEAL OF NON HYPERBOLIC EQS |
---|
| 1700 | |
---|
| 1701 | chart[idchart][4]=expI; // new input ideal and coefficients |
---|
| 1702 | chart[idchart][5]=Coef; |
---|
| 1703 | chart[idchart][10]=hipercoef+hipercoefson; |
---|
| 1704 | chart[idchart][11]=hiperexp+hiperexpson; |
---|
| 1705 | |
---|
| 1706 | idchart=idchart; |
---|
| 1707 | control2=1; // it is necessary to reset the value of c |
---|
| 1708 | control3=1; // and the previous exceptional divisors |
---|
| 1709 | } |
---|
| 1710 | |
---|
| 1711 | // PROBABLY IT IS NEC MORE INFORMATION !!! |
---|
| 1712 | |
---|
| 1713 | } // closing else otherwise |
---|
| 1714 | |
---|
| 1715 | } // closing else case I!=1 |
---|
| 1716 | |
---|
| 1717 | } // closing else for EordI=0 |
---|
| 1718 | |
---|
| 1719 | if (EordI<0) {print("ERROR in chart"); print(idchart); ~; break;} |
---|
| 1720 | |
---|
| 1721 | //----------------------- guardar de momento-------- |
---|
| 1722 | // if (EordI==0) {auxfchart[1]=chart[idchart]; |
---|
| 1723 | // finalchart=finalchart+auxfchart; |
---|
| 1724 | // L2=Gamma(expJ,c,n); // HAY QUE APLICARLO AL M NO AL J |
---|
| 1725 | // invariant[idchart]=L2[2]; |
---|
| 1726 | // auxinv[idchart]=list(0); |
---|
| 1727 | // nsons[idchart]=0; |
---|
| 1728 | // idchart=idchart+1;} |
---|
| 1729 | //------------------------------------------------ |
---|
| 1730 | |
---|
| 1731 | |
---|
| 1732 | } // END ELSE |
---|
| 1733 | //--------------------------------------------------- |
---|
| 1734 | |
---|
| 1735 | } // END LOOP WHILE |
---|
| 1736 | |
---|
| 1737 | tnum=step[nchart]; |
---|
| 1738 | |
---|
| 1739 | totalinv=resfunction(invariant,auxinv,nchart,n); |
---|
| 1740 | |
---|
| 1741 | return(chart,finalchart,invariant,nchart,step,nsons,auxinv,mobile,totalinv); |
---|
| 1742 | } |
---|
| 1743 | example |
---|
| 1744 | {"EXAMPLE:"; echo = 2; |
---|
| 1745 | ring r = 0,(x(1..2)),dp; |
---|
| 1746 | ideal J=x(1)^2-x(2)^3; |
---|
| 1747 | list L=Eresol(J); |
---|
| 1748 | "Please press return after each break point to see the next element of the output list"; |
---|
| 1749 | L[1][1]; // information of the first chart, L[1] list of charts |
---|
| 1750 | ~; |
---|
| 1751 | L[2]; // list of charts with information about sons |
---|
| 1752 | ~; |
---|
| 1753 | L[3]; // invariant, "#" means solved chart |
---|
| 1754 | ~; |
---|
| 1755 | L[4]; // number of charts, 7 in this example |
---|
| 1756 | ~; |
---|
| 1757 | L[5]; // height corresponding to each chart |
---|
| 1758 | ~; |
---|
| 1759 | L[6]; // number of sons |
---|
| 1760 | ~; |
---|
| 1761 | L[7]; // auxiliary invariant |
---|
| 1762 | ~; |
---|
| 1763 | L[8]; // H exceptional divisors and more information |
---|
| 1764 | ~; |
---|
| 1765 | L[9]; // complete resolution function |
---|
| 1766 | |
---|
| 1767 | "Second example, write L[i] to see the i-th component of the list"; |
---|
| 1768 | ring r = 0,(x(1..3)),dp; |
---|
| 1769 | ideal J=x(1)^2*x(2),x(3)^3; // SOLVED! |
---|
| 1770 | list L=Eresol(J); |
---|
| 1771 | L[4]; // 16 charts |
---|
| 1772 | L[9]; // complete resolution function |
---|
| 1773 | ~; |
---|
| 1774 | |
---|
| 1775 | "Third example, write L[i] to see the i-th component of the list"; |
---|
| 1776 | ring r = 0,(x(1..2)),dp; |
---|
| 1777 | ideal J=x(1)^3-x(1)*x(2)^3; |
---|
| 1778 | list L=Eresol(J); |
---|
| 1779 | L[4]; // 8 charts, rational exponents |
---|
| 1780 | L[9]; // complete resolution function |
---|
| 1781 | ~; |
---|
| 1782 | } |
---|
| 1783 | |
---|
| 1784 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1785 | |
---|
| 1786 | proc resfunction(list invariant, list auxinv, int nchart,int n) |
---|
| 1787 | "USAGE: resfunction(invariant,auxinv,nchart,n); |
---|
| 1788 | invariant, auxinv lists, nchart, n integers |
---|
| 1789 | COMPUTE: Patch the resolution function |
---|
| 1790 | RETURN: The complete resolution function |
---|
| 1791 | EXAMPLE: example resfunction; shows an example |
---|
| 1792 | " |
---|
| 1793 | { |
---|
| 1794 | int i,j,l,k; |
---|
| 1795 | list patchfun,aux; |
---|
| 1796 | |
---|
| 1797 | for (i=1;i<=nchart;i++){patchfun[i]=invariant[i];} |
---|
| 1798 | |
---|
| 1799 | for (i=1;i<=nchart;i++){if (auxinv[i][1]!=0 and size(auxinv[i])==3){l=size(invariant[i]); |
---|
| 1800 | for (j=1;j<=l;j++){ |
---|
| 1801 | if (invariant[i][j]==0){aux=auxinv[i]; |
---|
| 1802 | patchfun[i][j]=aux; |
---|
| 1803 | if (l<n){for (k=j+1;k<=n;k++){patchfun[i][k]="*";}}}} |
---|
| 1804 | |
---|
| 1805 | } |
---|
| 1806 | else{ |
---|
| 1807 | if (auxinv[i][1]==1 and size(auxinv[i])==1){l=size(invariant[i]); |
---|
| 1808 | if (l<n){for (k=l+1;k<=n;k++){patchfun[i][k]="*";}} |
---|
| 1809 | } |
---|
| 1810 | } |
---|
| 1811 | } |
---|
| 1812 | |
---|
| 1813 | return(patchfun); |
---|
| 1814 | } |
---|
| 1815 | example |
---|
| 1816 | {"EXAMPLE:"; echo = 2; |
---|
| 1817 | ring r = 0,(x(1..2)),dp; |
---|
| 1818 | ideal J=x(1)^2-x(2)^3; |
---|
| 1819 | list L=Eresol(J); |
---|
| 1820 | L[3]; // incomplete resolution function |
---|
| 1821 | resfunction(L[3],L[7],7,2); // complete resolution function |
---|
| 1822 | } |
---|
| 1823 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1824 | // // |
---|
| 1825 | // MAIN PROCEDURE // |
---|
| 1826 | // // |
---|
| 1827 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
| 1828 | |
---|
| 1829 | proc BINresol(ideal J) |
---|
| 1830 | "USAGE: BINresol(J); J ideal |
---|
| 1831 | RETURN: E-resolution of singularities of a binomial ideal J in terms of the affine charts, see example |
---|
| 1832 | EXAMPLE: example BINresol; shows an example |
---|
| 1833 | " |
---|
| 1834 | { |
---|
| 1835 | |
---|
| 1836 | int p,n; |
---|
| 1837 | |
---|
| 1838 | p=char(basering); |
---|
| 1839 | n=nvars(basering); // already computed in Eresol, it can be improved |
---|
| 1840 | |
---|
| 1841 | def rr=basering; |
---|
| 1842 | |
---|
| 1843 | // INTERNAL CHANGE: changing the name of the variables, only if it is necessary |
---|
| 1844 | |
---|
| 1845 | list Mout=changeoriginalvar(); |
---|
| 1846 | |
---|
| 1847 | if (Mout[2]==1){ |
---|
| 1848 | def r=Mout[1]; |
---|
| 1849 | setring r; |
---|
| 1850 | ideal chy=maxideal(1); |
---|
| 1851 | map frr=rr,chy; |
---|
| 1852 | ideal J=frr(J); |
---|
| 1853 | } |
---|
| 1854 | // else{def r=basering;} // CHECK THAT IS NECESSARY !!! |
---|
| 1855 | |
---|
| 1856 | // IF WE ARE IN POSTIVE CHAR |
---|
| 1857 | |
---|
| 1858 | if (p>0){list Lring=ringlist(basering); |
---|
| 1859 | Lring[1]=0; |
---|
| 1860 | // def r=basering; |
---|
| 1861 | def Rnew=ring(Lring); |
---|
| 1862 | setring Rnew; |
---|
| 1863 | ideal chy=maxideal(1); |
---|
| 1864 | map fRnew=r,chy; |
---|
| 1865 | ideal J=fRnew(J); |
---|
| 1866 | |
---|
| 1867 | // E-RESOLUTION, Computations in char 0 |
---|
| 1868 | |
---|
| 1869 | list L=Eresol(J); |
---|
| 1870 | |
---|
| 1871 | // STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL |
---|
| 1872 | |
---|
| 1873 | // not implemented yet, CHAR p !!!! |
---|
| 1874 | |
---|
| 1875 | // STEP 3: DO THE E-RESOLUTION AGAIN (char 0 again) |
---|
| 1876 | |
---|
| 1877 | |
---|
| 1878 | // generating output in char p |
---|
| 1879 | |
---|
| 1880 | int q=lcmofall(L[4],L[8]); // lcm of the denominators |
---|
| 1881 | |
---|
| 1882 | list B=genoutput(L[1],L[8],L[4],L[6],n,q,p); // generate output needed for visualization |
---|
| 1883 | |
---|
| 1884 | |
---|
| 1885 | // setring r; // Back to the basering |
---|
| 1886 | // ideal chy=maxideal(1); |
---|
| 1887 | // map fr=Rnew,chy; |
---|
| 1888 | // list L=fr(L); |
---|
| 1889 | // list B=fr(B); |
---|
| 1890 | |
---|
| 1891 | } |
---|
| 1892 | |
---|
| 1893 | else{ |
---|
| 1894 | |
---|
| 1895 | // E-RESOLUTION |
---|
| 1896 | |
---|
| 1897 | list L=Eresol(J); |
---|
| 1898 | |
---|
| 1899 | // STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL |
---|
| 1900 | |
---|
| 1901 | // not implemented yet |
---|
| 1902 | |
---|
| 1903 | // STEP 3: DO THE E-RESOLUTION AGAIN |
---|
| 1904 | |
---|
| 1905 | |
---|
| 1906 | // generating output |
---|
| 1907 | |
---|
| 1908 | int q=lcmofall(L[4],L[8]); |
---|
| 1909 | |
---|
| 1910 | list B=genoutput(L[1],L[8],L[4],L[6],n,q,p); |
---|
| 1911 | |
---|
| 1912 | } |
---|
| 1913 | |
---|
| 1914 | return(B); |
---|
| 1915 | } |
---|
| 1916 | example |
---|
| 1917 | {"EXAMPLE:"; echo = 2; |
---|
| 1918 | ring r = 0,(x(1..2)),dp; |
---|
| 1919 | ideal J=x(1)^2-x(2)^3; |
---|
| 1920 | list B=BINresol(J); |
---|
| 1921 | B[1]; // list of final charts |
---|
| 1922 | B[2]; // list of all charts |
---|
| 1923 | |
---|
| 1924 | ring r = 2,(x(1..3)),dp; |
---|
| 1925 | ideal J=x(1)^2-x(2)^2*x(3)^2; |
---|
| 1926 | list B=BINresol(J); |
---|
| 1927 | B[2]; // list of all charts |
---|
| 1928 | } |
---|
| 1929 | /////////////////////////////////////////////////////// |
---|
| 1930 | |
---|
| 1931 | proc Maxord(list L,int n) |
---|
| 1932 | "USAGE: Maxord(L,n); L list, n integer |
---|
| 1933 | COMPUTE: Find the maximal entry of a list, input is a list defining a monomial |
---|
| 1934 | RETURN: maximum entry of a list and its position |
---|
| 1935 | EXAMPLE: example Maxord; shows an example |
---|
| 1936 | " |
---|
| 1937 | {int i,can; |
---|
| 1938 | number canmax; |
---|
| 1939 | list aux; |
---|
| 1940 | |
---|
| 1941 | canmax=1; |
---|
| 1942 | can=1; |
---|
| 1943 | for (i=1;i<=n;i++) |
---|
| 1944 | { if (L[i]>=canmax and i>=can) |
---|
| 1945 | {canmax=L[i]; can=i;}} |
---|
| 1946 | |
---|
| 1947 | return(canmax,can); |
---|
| 1948 | } |
---|
| 1949 | example |
---|
| 1950 | {"EXAMPLE:"; echo = 2; |
---|
| 1951 | ring r = 0,(x(1..3)),dp; |
---|
| 1952 | ideal J=x(1)^2*x(2)*x(3)^5; |
---|
| 1953 | list L=data(J,1,3); |
---|
| 1954 | L[2]; // list of exponents |
---|
| 1955 | Maxord(L[2][1][1],3); |
---|
| 1956 | } |
---|
| 1957 | /////////////////////////////////////////////////////// |
---|
| 1958 | |
---|
| 1959 | proc Gamma(list expM,number c,int n) |
---|
| 1960 | "USAGE: Gamma(L,c,n); L list, c number, n integer |
---|
| 1961 | COMPUTE: The Gamma function, resolution function corresponding to the monomial case |
---|
| 1962 | RETURN: lists of maximum exponents in L, value of Gamma function, center of blow up |
---|
| 1963 | EXAMPLE: example Gamma; shows an example |
---|
| 1964 | " |
---|
| 1965 | {int i,j,k,l,cont,can; |
---|
| 1966 | intvec upla; |
---|
| 1967 | number canmax; |
---|
| 1968 | list expM2,gamma,L,aux,maxlist,center,aux2; |
---|
| 1969 | |
---|
| 1970 | i=1; |
---|
| 1971 | cont=0; |
---|
| 1972 | expM2=expM; |
---|
| 1973 | |
---|
| 1974 | while (cont==0 and i<=n) |
---|
| 1975 | { |
---|
| 1976 | |
---|
| 1977 | L=Maxord(expM2,n); |
---|
| 1978 | aux=L[1]; |
---|
| 1979 | maxlist=maxlist + aux; |
---|
| 1980 | can=L[2]; |
---|
| 1981 | |
---|
| 1982 | if (i==1) {upla=can; center=can;} |
---|
| 1983 | else {upla=upla,can; aux2=can; center=center+aux2;} |
---|
| 1984 | |
---|
| 1985 | canmax=sum(maxlist); |
---|
| 1986 | if (canmax>=c) |
---|
| 1987 | {gamma[1]=-i; gamma[2]=canmax/c; gamma[3]=upla; cont=1;} |
---|
| 1988 | else {expM2[can]=0;} |
---|
| 1989 | i=i+1; |
---|
| 1990 | } |
---|
| 1991 | return(maxlist,gamma,center); |
---|
| 1992 | } |
---|
| 1993 | example |
---|
| 1994 | {"EXAMPLE:"; echo = 2; |
---|
| 1995 | ring r = 0,(x(1..5)),dp; |
---|
| 1996 | ideal J=x(1)^2*x(2)*x(3)^5*x(4)^2*x(5)^3; |
---|
| 1997 | list L=data(J,1,5); |
---|
| 1998 | list G=Gamma(L[2][1][1],9,5); // critical value c=9 |
---|
| 1999 | G[1]; // maximum exponents in the ideal |
---|
| 2000 | G[2]; // maximal value of Gamma function |
---|
| 2001 | G[3]; // center given by Gamma |
---|
| 2002 | } |
---|
| 2003 | /////////////////////////////////////////////////////// |
---|
| 2004 | |
---|
| 2005 | proc convertdata(list C,list L, int n, list flaglist) |
---|
| 2006 | "USAGE: convertdata(C,L,n,flaglist); |
---|
| 2007 | C, L, flaglist lists, n integer |
---|
| 2008 | COMPUTE: Compute the ideal corresponding to the given lists C,L |
---|
| 2009 | RETURN: an ideal whose coefficients are given by C, exponents given by L |
---|
| 2010 | EXAMPLE: example convertdata; shows an example |
---|
| 2011 | " |
---|
| 2012 | {int i,j,k,a,b,lon; |
---|
| 2013 | poly aux,aux1,aux2,aux3,f; |
---|
| 2014 | ideal J; |
---|
| 2015 | |
---|
| 2016 | aux=poly(0); |
---|
| 2017 | aux1=poly(1); |
---|
| 2018 | aux2=poly(0); |
---|
| 2019 | aux3=poly(1); |
---|
| 2020 | |
---|
| 2021 | |
---|
| 2022 | k=size(L); |
---|
| 2023 | for (i=1;i<=k;i++){lon=size(C[i]); |
---|
| 2024 | if (lon==1){ // variables in the monomial |
---|
| 2025 | for (j=1;j<=n;j++){a=int(poly(L[i][1][j])); |
---|
| 2026 | if (a!=0){ |
---|
| 2027 | if (flaglist[j]==0){aux=poly(x(j)^a); |
---|
| 2028 | aux1=aux1*aux;} |
---|
| 2029 | else {aux=poly(y(j)^a); |
---|
| 2030 | aux1=aux1*aux;} |
---|
| 2031 | } |
---|
| 2032 | } |
---|
| 2033 | if (C[i][1]!=0){aux1=C[i][1]*aux1;} // we add the coefficient |
---|
| 2034 | else {aux1=0;} |
---|
| 2035 | |
---|
| 2036 | J[i]=aux1; |
---|
| 2037 | aux1=poly(1); |
---|
| 2038 | } |
---|
| 2039 | |
---|
| 2040 | else{ // variables in the binomial |
---|
| 2041 | |
---|
| 2042 | for (j=1;j<=n;j++){a=int(poly(L[i][1][j])); b=int(poly(L[i][2][j])); |
---|
| 2043 | |
---|
| 2044 | if (a!=0){ |
---|
| 2045 | if (flaglist[j]==0){aux=poly(x(j)^a); |
---|
| 2046 | aux1=aux1*aux;} |
---|
| 2047 | else {aux=poly(y(j)^a); |
---|
| 2048 | aux1=aux1*aux;} |
---|
| 2049 | } |
---|
| 2050 | |
---|
| 2051 | if (b!=0){ |
---|
| 2052 | if (flaglist[j]==0){aux2=poly(x(j)^b); |
---|
| 2053 | aux3=aux3*aux2;} |
---|
| 2054 | else {aux2=poly(y(j)^b); |
---|
| 2055 | aux3=aux3*aux2;} |
---|
| 2056 | } |
---|
| 2057 | } |
---|
| 2058 | // we add the coefficients |
---|
| 2059 | if (C[i][1]!=0){aux1=C[i][1]*aux1;} |
---|
| 2060 | else {aux1=0;} |
---|
| 2061 | if (C[i][2]!=0){aux3=C[i][2]*aux3;} |
---|
| 2062 | else {aux3=0;} |
---|
| 2063 | |
---|
| 2064 | f=aux1+aux3; |
---|
| 2065 | J[i]=f; |
---|
| 2066 | aux1=poly(1); |
---|
| 2067 | aux3=poly(1); |
---|
| 2068 | |
---|
| 2069 | } |
---|
| 2070 | } |
---|
| 2071 | return(J); |
---|
| 2072 | } |
---|
| 2073 | example |
---|
| 2074 | {"EXAMPLE:"; echo = 2; |
---|
| 2075 | ring r = 0,(x(1..4),y(5)),dp; |
---|
| 2076 | list M=identifyvar(); |
---|
| 2077 | ideal J=x(1)^2*y(5)^2-x(2)^2*x(3)^2,6*x(4)^2; |
---|
| 2078 | list L=data(J,2,5); |
---|
| 2079 | L[1]; // Coefficients |
---|
| 2080 | L[2]; // Exponents |
---|
| 2081 | ideal J2=convertdata(L[1],L[2],5,M); |
---|
| 2082 | J2; |
---|
| 2083 | } |
---|
| 2084 | |
---|
| 2085 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 2086 | |
---|
| 2087 | proc lcmofall(int nchart,list mobile) |
---|
| 2088 | "USAGE: lcmofall(nchart,mobile); |
---|
| 2089 | nchart integer, mobile list of lists |
---|
| 2090 | COMPUTE: Compute the lcm of the denominators of the E-orders of all the charts |
---|
| 2091 | RETURN: an integer given the lcm |
---|
| 2092 | NOTE: CALL BEFORE salida |
---|
| 2093 | EXAMPLE: example lcmofall; shows an example |
---|
| 2094 | " |
---|
| 2095 | |
---|
| 2096 | { |
---|
| 2097 | int i,m,tip,mcmall; |
---|
| 2098 | intvec numall; |
---|
| 2099 | |
---|
| 2100 | for (i=2;i<=nchart+1;i++){ |
---|
| 2101 | tip=mobile[i][1]; |
---|
| 2102 | if (tip!=1){numall=numall,tip;} |
---|
| 2103 | } |
---|
| 2104 | m=size(numall); |
---|
| 2105 | |
---|
| 2106 | if (m==1){mcmall=1;} |
---|
| 2107 | else{ |
---|
| 2108 | if (numall[1]==0){numall=numall[2..m];} |
---|
| 2109 | mcmall=lcm(numall);} |
---|
| 2110 | |
---|
| 2111 | return(mcmall); |
---|
| 2112 | } |
---|
| 2113 | example |
---|
| 2114 | {"EXAMPLE:"; echo = 2; |
---|
| 2115 | ring r = 0,(x(1..2)),dp; |
---|
| 2116 | ideal J=x(1)^3-x(1)*x(2)^3; |
---|
| 2117 | list L=Eresol(J); |
---|
| 2118 | L[4]; // 8 charts, rational exponents |
---|
| 2119 | L[8][2][2]; // E-orders at the first chart |
---|
| 2120 | lcmofall(8,L[8]); |
---|
| 2121 | } |
---|
| 2122 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 2123 | |
---|
| 2124 | proc salida(int idchart,list chart,list mobile,int numson,intvec previousa,int n,int q,int p) |
---|
| 2125 | "USAGE: salida(idchart,chart,mobile,numson,previousa,n,q,p); |
---|
| 2126 | idchart, numson, n, q, p integers, chart, mobile, lists, previousa intvec |
---|
| 2127 | COMPUTE: CONVERT THE OUTPUT OF A CHART IN A RING, WHERE DEFINE A BASIC OBJECT (BO) |
---|
| 2128 | RETURN: the ring corresponding to the chart |
---|
| 2129 | EXAMPLE: example salida; shows an example |
---|
| 2130 | " |
---|
| 2131 | { |
---|
| 2132 | int l,i,m,aux,parent,m4,j; |
---|
| 2133 | intvec Hhist,EOhist,aux7,aux9; |
---|
| 2134 | list expJ,Coef,BO,blwhist,Eolist,hipercoef,hiperexp; |
---|
| 2135 | list flag; |
---|
| 2136 | |
---|
| 2137 | // chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp |
---|
| 2138 | // mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2]; |
---|
| 2139 | |
---|
| 2140 | // we need to define the suitable ring at this chart |
---|
| 2141 | |
---|
| 2142 | list Lring=ringlist(basering); |
---|
| 2143 | def RR2=basering; |
---|
| 2144 | |
---|
| 2145 | flag=chart[6]; |
---|
| 2146 | string newl; |
---|
| 2147 | |
---|
| 2148 | for (l=1;l<=n; l++){if (flag[l]==1){newl=string(l); |
---|
| 2149 | Lring[2][l]="y("+newl+")";} } |
---|
| 2150 | |
---|
| 2151 | |
---|
| 2152 | def RRnew=ring(Lring); |
---|
| 2153 | setring RRnew; |
---|
| 2154 | ideal chy=maxideal(1); |
---|
| 2155 | map fRnew=RR2,chy; |
---|
| 2156 | |
---|
| 2157 | list chart=fRnew(chart); |
---|
| 2158 | |
---|
| 2159 | list mobile2=fRnew(mobile); |
---|
| 2160 | |
---|
| 2161 | |
---|
| 2162 | flag=chart[6]; |
---|
| 2163 | |
---|
| 2164 | // we need to convert expJ and Coef to an ideal |
---|
| 2165 | |
---|
| 2166 | expJ=chart[4]; |
---|
| 2167 | Coef=chart[5]; |
---|
| 2168 | Hhist=chart[7]; |
---|
| 2169 | blwhist=chart[8]; |
---|
| 2170 | |
---|
| 2171 | // now the ideal will be correctly defined in the ring Rnew |
---|
| 2172 | |
---|
| 2173 | ideal J2=convertdata(Coef,expJ,n,flag); // Computations in RRnew |
---|
| 2174 | |
---|
| 2175 | //------------------------------------------------------------------------------ |
---|
| 2176 | // START TO CREATE THE BO corresponding to this chart |
---|
| 2177 | |
---|
| 2178 | BO=createBO(J2); |
---|
| 2179 | |
---|
| 2180 | // MODIFY BO WITH THE INFORMATION OF THE CHART |
---|
| 2181 | |
---|
| 2182 | // BO[1] an ideal, say W_i, defining the ambient space of the i-th chart of the blowing up |
---|
| 2183 | // If there are hyperbolic equations, we put them here |
---|
| 2184 | |
---|
| 2185 | hipercoef=chart[10]; |
---|
| 2186 | hiperexp=chart[11]; |
---|
| 2187 | |
---|
| 2188 | if (size(hipercoef)!=0){ |
---|
| 2189 | ideal ambJ=convertdata(hipercoef,hiperexp,n,flag); |
---|
| 2190 | BO[1]=ambJ; |
---|
| 2191 | } |
---|
| 2192 | |
---|
| 2193 | // BO[2] an ideal defining the controlled transform |
---|
| 2194 | |
---|
| 2195 | BO[2]=J2; |
---|
| 2196 | |
---|
| 2197 | // BO[3] intvec, tupla containing the maximal E-order of BO[2] |
---|
| 2198 | |
---|
| 2199 | if (numson==0){BO[3]=1;} // we write 1 if the chart is a final chart |
---|
| 2200 | else{ |
---|
| 2201 | Eolist=mobile2[2]; // otherwise, convert the list of E-orders in an intvec |
---|
| 2202 | m=size(Eolist); |
---|
| 2203 | aux=int(Eolist[1]*q); |
---|
| 2204 | EOhist=aux; |
---|
| 2205 | |
---|
| 2206 | if (m>1){for (i=2;i<=m;i++){aux=int(Eolist[i]*q); EOhist=EOhist,aux;}} |
---|
| 2207 | |
---|
| 2208 | BO[3]=EOhist; |
---|
| 2209 | } |
---|
| 2210 | |
---|
| 2211 | // BO[4] the list of exceptional divisors given by Hhist |
---|
| 2212 | |
---|
| 2213 | BO[4]=constructH(Hhist,n,flag); |
---|
| 2214 | |
---|
| 2215 | // BO[5] an ideal defining the map K[W] ----> K[Wi] given by blwhist |
---|
| 2216 | |
---|
| 2217 | BO[5]=constructblwup(blwhist,n,chy,flag); |
---|
| 2218 | |
---|
| 2219 | // BO[6] an intvec, BO[6][j]=1 indicates that <BO[4][j],BO[2]>=1, i.e. the |
---|
| 2220 | // strict transform does not meet the j-th exceptional divisor |
---|
| 2221 | |
---|
| 2222 | m4=size(BO[4]); |
---|
| 2223 | ideal auxydeal; |
---|
| 2224 | ideal Jint; |
---|
| 2225 | |
---|
| 2226 | for (j=1;j<=m4;j++){ |
---|
| 2227 | |
---|
| 2228 | auxydeal=BO[4][j]+J2; |
---|
| 2229 | Jint=std(auxydeal); |
---|
| 2230 | |
---|
| 2231 | if (size(Jint)==1 and Jint[1]==1){BO[6][j]=1;} |
---|
| 2232 | else{BO[6][j]=0;} |
---|
| 2233 | } |
---|
| 2234 | |
---|
| 2235 | // BO[7] intvec, the index of the first blown-up object in the resolution process |
---|
| 2236 | // leading to this object for which the value of b was BO[3] |
---|
| 2237 | // the subsequent ones are the indices for the Coeff-Objects |
---|
| 2238 | // of BO[2] used when determining the center |
---|
| 2239 | // index of last element of H^- in H |
---|
| 2240 | |
---|
| 2241 | |
---|
| 2242 | if (numson!=0){BO[7]=mobile2[8];} // it is always -1 at the final charts |
---|
| 2243 | |
---|
| 2244 | // BO[8] a matrix indicating that BO[4][i] meets BO[4][j] by BO[8][i,j]=1 for i < j |
---|
| 2245 | |
---|
| 2246 | if (m4>0){ |
---|
| 2247 | matrix aux8[m4][m4]; |
---|
| 2248 | |
---|
| 2249 | BO[8]=aux8; |
---|
| 2250 | |
---|
| 2251 | ideal auxydeal2; |
---|
| 2252 | ideal Jint2; |
---|
| 2253 | |
---|
| 2254 | for (i=1;i<=m4;i++){ |
---|
| 2255 | for (j=i+1;j<=m4;j++){ |
---|
| 2256 | auxydeal2=BO[4][i]+BO[4][j]; |
---|
| 2257 | Jint2=std(auxydeal2); |
---|
| 2258 | |
---|
| 2259 | if (size(Jint2)==1 and Jint2[1]==1){BO[8][i,j]=0;} |
---|
| 2260 | else{ for (l=1;l<j;l++){BO[8][l,j]=1;} } |
---|
| 2261 | } |
---|
| 2262 | |
---|
| 2263 | } |
---|
| 2264 | } |
---|
| 2265 | else{ matrix aux8[1][1]; |
---|
| 2266 | BO[8]=aux8;} |
---|
| 2267 | |
---|
| 2268 | |
---|
| 2269 | // BO[9] INTERNAL DATA, second component of Villamayor resolution function, |
---|
| 2270 | // only needed to use the visualization procedures |
---|
| 2271 | |
---|
| 2272 | int m3=size(BO[3]); |
---|
| 2273 | |
---|
| 2274 | if (m3==1){aux9=intvec(0);} |
---|
| 2275 | else{ aux9[1]=0; |
---|
| 2276 | for (i=2;i<=m3;i++){aux9=aux9,0;} |
---|
| 2277 | } |
---|
| 2278 | |
---|
| 2279 | BO[9]=aux9; |
---|
| 2280 | |
---|
| 2281 | //------------------------------------------------------------------------------ |
---|
| 2282 | |
---|
| 2283 | // START TO CREATE THE extra information corresponding to this chart |
---|
| 2284 | |
---|
| 2285 | /////////////// Short description of data in a chart /////////////////// |
---|
| 2286 | // All chart data is stored in an object of type ring, the following |
---|
| 2287 | // variables are always present in such a ring: |
---|
| 2288 | |
---|
| 2289 | // BO: already created |
---|
| 2290 | |
---|
| 2291 | // cent: ideal, describing the upcoming center determined by the algorithm |
---|
| 2292 | |
---|
| 2293 | ideal cent=tradtoideal(previousa,J2,flag); |
---|
| 2294 | |
---|
| 2295 | |
---|
| 2296 | // path= module (autoconverted to matrix) |
---|
| 2297 | // path[1][idchart]=parent[idchart] index of the parent-chart in resolution history of this chart |
---|
| 2298 | // path[2][idchart]=index of this chart in relation with its brother-charts |
---|
| 2299 | |
---|
| 2300 | module path=chart[9]; |
---|
| 2301 | |
---|
| 2302 | |
---|
| 2303 | // lastMap: ideal, describing the preceding blow up leading to this chart |
---|
| 2304 | |
---|
| 2305 | ideal lastMap=constructlastblwup(blwhist,n,chy,flag); |
---|
| 2306 | |
---|
| 2307 | |
---|
| 2308 | //------------------------------------------------------------------------------ |
---|
| 2309 | |
---|
| 2310 | // EXTRA INFORMATION NEEDED |
---|
| 2311 | |
---|
| 2312 | list invSat=ideal(0),aux9; |
---|
| 2313 | |
---|
| 2314 | |
---|
| 2315 | // BACK TO THE CHAR OF THE ORIGINAL RING, IF IT HAD p>0 |
---|
| 2316 | |
---|
| 2317 | if (p>0){ |
---|
| 2318 | |
---|
| 2319 | list Lring; |
---|
| 2320 | Lring=ringlist(RRnew); |
---|
| 2321 | Lring[1]=p; |
---|
| 2322 | def auxRnew=ring(Lring); |
---|
| 2323 | |
---|
| 2324 | kill Lring; |
---|
| 2325 | setring auxRnew; |
---|
| 2326 | ideal chy=maxideal(1); |
---|
| 2327 | map frnew=RRnew,chy; |
---|
| 2328 | def BO=frnew(BO); |
---|
| 2329 | |
---|
| 2330 | // def chart=frr(chart); |
---|
| 2331 | def invSat=frnew(invSat); |
---|
| 2332 | def lastMap=frnew(lastMap); |
---|
| 2333 | def cent=frnew(cent); |
---|
| 2334 | def path=frnew(path); |
---|
| 2335 | |
---|
| 2336 | } |
---|
| 2337 | |
---|
| 2338 | // export everything needed |
---|
| 2339 | |
---|
| 2340 | export BO; |
---|
| 2341 | export(invSat); |
---|
| 2342 | export lastMap; |
---|
| 2343 | export path; |
---|
| 2344 | export cent; |
---|
| 2345 | |
---|
| 2346 | if (p==0){return(RRnew);} |
---|
| 2347 | else{ |
---|
| 2348 | return(auxRnew);} |
---|
| 2349 | } |
---|
| 2350 | example |
---|
| 2351 | {"EXAMPLE:"; echo = 2; |
---|
| 2352 | ring r = 0,(x(1..2)),dp; |
---|
| 2353 | ideal J=x(1)^2-x(2)^3; |
---|
| 2354 | list L=Eresol(J); |
---|
| 2355 | list B=salida(5,L[1][5],L[8][6],2,L[1][3][3],2,1,0); // chart 5 |
---|
| 2356 | def RR=B[1]; |
---|
| 2357 | setring RR; |
---|
| 2358 | BO; |
---|
| 2359 | "press return to see next example"; ~; |
---|
| 2360 | |
---|
| 2361 | ring r = 0,(x(1..2)),dp; |
---|
| 2362 | ideal J=x(1)^2-x(2)^3; |
---|
| 2363 | list L=Eresol(J); |
---|
| 2364 | list B=salida(7,L[1][7],L[8][8],0,L[1][5][3],2,1,0); // chart 7 |
---|
| 2365 | def RR=B[1]; |
---|
| 2366 | setring RR; |
---|
| 2367 | BO; |
---|
| 2368 | showBO(BO); |
---|
| 2369 | "press return to see next example"; ~; |
---|
| 2370 | |
---|
| 2371 | ring r = 0,(x(1..2)),dp; |
---|
| 2372 | ideal J=x(1)^3-x(1)*x(2)^3; |
---|
| 2373 | list L=Eresol(J); // 8 charts, rational exponents |
---|
| 2374 | list B=salida(1,L[1][1],L[8][2],2,0,2,2,0); // CHART 1 |
---|
| 2375 | def RR=B[1]; |
---|
| 2376 | setring RR; |
---|
| 2377 | BO; |
---|
| 2378 | |
---|
| 2379 | } |
---|
| 2380 | |
---|
| 2381 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 2382 | // CONVERT THE OUTPUT OF Eresol IN A LIST OF RINGS, WHERE A BASIC OBJECT (BO) IS DEFINED |
---|
| 2383 | // IN ORDER TO INTEGRATE THIS LIBRARY INSIDE THE LIBRARY resolve.lib |
---|
| 2384 | |
---|
| 2385 | proc genoutput(list chart,list mobile,int nchart,list nsons,int n,int q, int p) |
---|
| 2386 | "USAGE: genoutput(chart,mobile,nchart,nsons,n,q,p); |
---|
| 2387 | chart, mobile, nsons lists, nchart, n,q, p integers |
---|
| 2388 | RETURN: two lists, the first one gives the rings corresponding to the final charts, |
---|
| 2389 | the second one is the list of all rings corresponding to the affine charts of the resolution process |
---|
| 2390 | EXAMPLE: example genoutput; shows an example |
---|
| 2391 | " |
---|
| 2392 | { |
---|
| 2393 | int idchart,parent; |
---|
| 2394 | list auxlist,solvedrings,totalringlist,previousa; |
---|
| 2395 | list auxlistenp,solvedringsenp,totalringenp; |
---|
| 2396 | |
---|
| 2397 | // chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp |
---|
| 2398 | // mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2]; |
---|
| 2399 | |
---|
| 2400 | idchart=1; |
---|
| 2401 | |
---|
| 2402 | // first loop, construct list previousa |
---|
| 2403 | |
---|
| 2404 | while (idchart<=nchart) |
---|
| 2405 | { |
---|
| 2406 | if (idchart==1){previousa[1]=chart[2][3];} |
---|
| 2407 | |
---|
| 2408 | else{ |
---|
| 2409 | |
---|
| 2410 | // if there are no sons, the next center is nothing |
---|
| 2411 | |
---|
| 2412 | if (nsons[idchart]==0){previousa[idchart]=0;} |
---|
| 2413 | |
---|
| 2414 | // always fill the parent |
---|
| 2415 | |
---|
| 2416 | parent=chart[idchart][1]; |
---|
| 2417 | previousa[parent]=chart[idchart][3]; |
---|
| 2418 | } |
---|
| 2419 | idchart=idchart+1; |
---|
| 2420 | } |
---|
| 2421 | |
---|
| 2422 | // HERE BEGIN THE LOOP |
---|
| 2423 | |
---|
| 2424 | idchart=1; |
---|
| 2425 | |
---|
| 2426 | while (idchart<=nchart) |
---|
| 2427 | { |
---|
| 2428 | |
---|
| 2429 | def auxexit=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q,p); |
---|
| 2430 | |
---|
| 2431 | if (p>0){ // we need the computations in char 0 too |
---|
| 2432 | def auxexitenp=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q,0); |
---|
| 2433 | } |
---|
| 2434 | else{def auxexitenp=auxexit;} |
---|
| 2435 | |
---|
| 2436 | // we add the ring to the list of all rings |
---|
| 2437 | |
---|
| 2438 | auxlist[1]=auxexit; |
---|
| 2439 | totalringlist=totalringlist+auxlist; |
---|
| 2440 | |
---|
| 2441 | auxlistenp[1]=auxexitenp; |
---|
| 2442 | totalringenp=totalringenp+auxlistenp; |
---|
| 2443 | |
---|
| 2444 | // if the chart has no sons, add it to the list of final charts |
---|
| 2445 | |
---|
| 2446 | if (nsons[idchart]==0){solvedrings=solvedrings+auxlist; |
---|
| 2447 | solvedringsenp=solvedringsenp+auxlistenp; |
---|
| 2448 | } |
---|
| 2449 | |
---|
| 2450 | auxlist=list(); |
---|
| 2451 | auxlistenp=list(); |
---|
| 2452 | kill auxexit; |
---|
| 2453 | kill auxexitenp; |
---|
| 2454 | |
---|
| 2455 | idchart=idchart+1; |
---|
| 2456 | |
---|
| 2457 | } // EXIT WHILE |
---|
| 2458 | |
---|
| 2459 | return(solvedrings,totalringlist,solvedringsenp,totalringenp); |
---|
| 2460 | } |
---|
| 2461 | example |
---|
| 2462 | {"EXAMPLE:"; echo = 2; |
---|
| 2463 | ring r = 0,(x(1..2)),dp; |
---|
| 2464 | ideal J=x(1)^3-x(1)*x(2)^3; |
---|
| 2465 | list L=Eresol(J); // 8 charts, rational exponents |
---|
| 2466 | list B=genoutput(L[1],L[8],L[4],L[6],2,2,0); // generates the output |
---|
| 2467 | presentTree(B); |
---|
| 2468 | list iden0=collectDiv(B); |
---|
| 2469 | ResTree(B,iden0[1]); // generates the resolution tree |
---|
| 2470 | |
---|
| 2471 | // Use presentTree(B); to see the final charts |
---|
| 2472 | // To see the tree type in another shell |
---|
| 2473 | // dot -Tjpg ResTree.dot -o ResTree.jpg |
---|
| 2474 | // /usr/bin/X11/xv ResTree.jpg |
---|
| 2475 | |
---|
| 2476 | } |
---|
| 2477 | ///////////////////////////////////////////////////////////////////// |
---|
| 2478 | |
---|
| 2479 | proc computemcm(list Eolist) |
---|
| 2480 | "USAGE: computemcm(Eolist); Eolist list |
---|
| 2481 | RETURN: an integer, the least common multiple of the denominators of the E-orders |
---|
| 2482 | NOTE: Make the same as lcmofall but for one chart. NECESSARY BECAUSE THE E-ORDERS ARE OF TYPE NUMBER!! |
---|
| 2483 | EXAMPLE: example computemcm; shows an example |
---|
| 2484 | " |
---|
| 2485 | { |
---|
| 2486 | int m,i,aux,mcmchart; |
---|
| 2487 | intvec num; |
---|
| 2488 | |
---|
| 2489 | m=size(Eolist); |
---|
| 2490 | |
---|
| 2491 | if (m==1){mcmchart=int(denominator(Eolist[1])); return(mcmchart);} |
---|
| 2492 | |
---|
| 2493 | if (m>1){num=int(denominator(Eolist[1])); |
---|
| 2494 | for (i=2;i<=m;i++){aux=int(denominator(Eolist[i])); |
---|
| 2495 | num=num,aux; }} |
---|
| 2496 | |
---|
| 2497 | mcmchart=lcm(num); |
---|
| 2498 | |
---|
| 2499 | return(mcmchart); |
---|
| 2500 | } |
---|
| 2501 | example |
---|
| 2502 | {"EXAMPLE:"; echo = 2; |
---|
| 2503 | ring r = 0,(x(1..2)),dp; |
---|
| 2504 | ideal J=x(1)^3-x(1)*x(2)^3; |
---|
| 2505 | list L=Eresol(J); // 8 charts, rational exponents |
---|
| 2506 | L[8][2][2]; // maximal E-order at the first chart |
---|
| 2507 | computemcm(L[8][2][2]); |
---|
| 2508 | |
---|
| 2509 | } |
---|
| 2510 | ///////////////////////////////////////////////////////////////////// |
---|
| 2511 | |
---|
| 2512 | proc constructH(intvec Hhist,int n,list flag) |
---|
| 2513 | "USAGE: constructH(Hhist,n,flag); |
---|
| 2514 | Hhist intvec, n integer, flag list |
---|
| 2515 | RETURN: the list of exceptional divisors accumulated at this chart |
---|
| 2516 | EXAMPLE: example constructH; shows an example |
---|
| 2517 | " |
---|
| 2518 | { |
---|
| 2519 | int i,j,m,l; |
---|
| 2520 | list exceplist; |
---|
| 2521 | ideal aux; |
---|
| 2522 | |
---|
| 2523 | m=size(Hhist); |
---|
| 2524 | if (Hhist[1]==0 and m>1){Hhist=Hhist[2..m]; m=m-1; |
---|
| 2525 | |
---|
| 2526 | for (i=1;i<=m;i++){ |
---|
| 2527 | l=Hhist[i]; |
---|
| 2528 | if (flag[l]==0){aux=ideal(poly(x(l))); } |
---|
| 2529 | else {aux=ideal(poly(y(l))); } |
---|
| 2530 | |
---|
| 2531 | exceplist[i]=aux; |
---|
| 2532 | } |
---|
| 2533 | // eliminate repeated variables |
---|
| 2534 | for (i=1;i<=m;i++){for (j=1;j<=m;j++){ |
---|
| 2535 | if (Hhist[i]==Hhist[j] and i!=j){ |
---|
| 2536 | if (i<j){exceplist[i]=ideal(1);} |
---|
| 2537 | if (i>j){exceplist[j]=ideal(1);} |
---|
| 2538 | } |
---|
| 2539 | } |
---|
| 2540 | } |
---|
| 2541 | |
---|
| 2542 | } |
---|
| 2543 | else {exceplist=list();} |
---|
| 2544 | |
---|
| 2545 | // else {exceplist=list(ideal(0));} // IF IT FAILS USE THIS |
---|
| 2546 | |
---|
| 2547 | return(exceplist); |
---|
| 2548 | } |
---|
| 2549 | example |
---|
| 2550 | {"EXAMPLE:"; echo = 2; |
---|
| 2551 | ring r = 0,(x(1..3)),dp; |
---|
| 2552 | list flag=identifyvar(); |
---|
| 2553 | ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; |
---|
| 2554 | list L=Eresol(J); // 7 charts |
---|
| 2555 | // history of the exceptional divisors at the 7-th chart |
---|
| 2556 | L[1][7][7]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts |
---|
| 2557 | constructH(L[1][7][7],3,flag); |
---|
| 2558 | |
---|
| 2559 | } |
---|
| 2560 | ///////////////////////////////////////////////////////////////////// |
---|
| 2561 | |
---|
| 2562 | proc constructblwup(list blwhist,int n,ideal chy,list flag) |
---|
| 2563 | "USAGE: constructblwup(blwhist,n,chy,flag); |
---|
| 2564 | blwhist, flag lists, n integer, chy ideal |
---|
| 2565 | RETURN: the ideal defining the map K[W] --> K[Wi], |
---|
| 2566 | which gives the composition map of all the blowing up leading to this chart |
---|
| 2567 | NOTE: NECESSARY START WITH COLUMNS |
---|
| 2568 | EXAMPLE: example constructblwup; shows an example |
---|
| 2569 | " |
---|
| 2570 | { |
---|
| 2571 | int i,j,m,m2; |
---|
| 2572 | poly aux2; |
---|
| 2573 | |
---|
| 2574 | m=size(blwhist[1]); |
---|
| 2575 | |
---|
| 2576 | for (j=1;j<=m;j++){ |
---|
| 2577 | for (i=1;i<=n;i++){ m2=blwhist[i][j]; |
---|
| 2578 | |
---|
| 2579 | // If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not |
---|
| 2580 | |
---|
| 2581 | if (m2!=0){ |
---|
| 2582 | if (flag[m2]==0){aux2=poly(x(m2));} |
---|
| 2583 | else {aux2=poly(y(m2));} |
---|
| 2584 | |
---|
| 2585 | // And then substitute this variable for the corresponding product in the whole ideal |
---|
| 2586 | |
---|
| 2587 | if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);} |
---|
| 2588 | else {chy=subst(chy,y(i),y(i)*aux2);} |
---|
| 2589 | |
---|
| 2590 | } |
---|
| 2591 | } |
---|
| 2592 | |
---|
| 2593 | } |
---|
| 2594 | |
---|
| 2595 | return(chy); |
---|
| 2596 | } |
---|
| 2597 | example |
---|
| 2598 | {"EXAMPLE:"; echo = 2; |
---|
| 2599 | ring r = 0,(x(1..3)),dp; |
---|
| 2600 | list flag=identifyvar(); |
---|
| 2601 | ideal chy=maxideal(1); |
---|
| 2602 | ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; |
---|
| 2603 | list L=Eresol(J); // 7 charts |
---|
| 2604 | // history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time |
---|
| 2605 | L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts |
---|
| 2606 | constructblwup(L[1][7][8],3,chy,flag); |
---|
| 2607 | } |
---|
| 2608 | ///////////////////////////////////////////////////////////////////// |
---|
| 2609 | |
---|
| 2610 | proc constructlastblwup(list blwhist,int n,ideal chy,list flag) |
---|
| 2611 | "USAGE: constructlastblwup(blwhist,n,chy,flag); |
---|
| 2612 | blwhist, flag lists, n integer, chy ideal |
---|
| 2613 | RETURN: the ideal defining the last blow up |
---|
| 2614 | NOTE: NECESSARY START WITH COLUMNS |
---|
| 2615 | EXAMPLE: example constructlastblwup; shows an example |
---|
| 2616 | " |
---|
| 2617 | { |
---|
| 2618 | int i,j,m,m2; |
---|
| 2619 | poly aux2; |
---|
| 2620 | |
---|
| 2621 | m=size(blwhist[1]); |
---|
| 2622 | |
---|
| 2623 | if (m>0){ |
---|
| 2624 | for (i=1;i<=n;i++){ m2=blwhist[i][m]; |
---|
| 2625 | |
---|
| 2626 | // If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not |
---|
| 2627 | |
---|
| 2628 | if (m2!=0){ |
---|
| 2629 | if (flag[m2]==0){aux2=poly(x(m2));} |
---|
| 2630 | else {aux2=poly(y(m2));} |
---|
| 2631 | |
---|
| 2632 | // And then substitute this variable for the corresponding product in the whole ideal |
---|
| 2633 | |
---|
| 2634 | if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);} |
---|
| 2635 | else {chy=subst(chy,y(i),y(i)*aux2);} |
---|
| 2636 | |
---|
| 2637 | } |
---|
| 2638 | } |
---|
| 2639 | } |
---|
| 2640 | |
---|
| 2641 | return(chy); |
---|
| 2642 | } |
---|
| 2643 | example |
---|
| 2644 | {"EXAMPLE:"; echo = 2; |
---|
| 2645 | ring r = 0,(x(1..3)),dp; |
---|
| 2646 | list flag=identifyvar(); |
---|
| 2647 | ideal chy=maxideal(1); |
---|
| 2648 | ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; |
---|
| 2649 | list L=Eresol(J); // 7 charts |
---|
| 2650 | // history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time |
---|
| 2651 | L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts |
---|
| 2652 | constructlastblwup(L[1][7][8],3,chy,flag); |
---|
| 2653 | } |
---|
| 2654 | ///////////////////////////////////////////////////////////////////// |
---|
| 2655 | |
---|
| 2656 | proc tradtoideal(intvec a,ideal J2,list flag) |
---|
| 2657 | "USAGE: tradtoideal(a,J2,flag); |
---|
| 2658 | a intvec, J2 ideal, flag list |
---|
| 2659 | COMPUTE: traslate to an ideal the intvec defining the center |
---|
| 2660 | RETURN: the ideal of the center, given by the intvec a, or J2 if a=0 |
---|
| 2661 | EXAMPLE: example tradtoideal; shows an example |
---|
| 2662 | " |
---|
| 2663 | { |
---|
| 2664 | int i,m; |
---|
| 2665 | ideal acenter,aux2; |
---|
| 2666 | |
---|
| 2667 | if (a==0){acenter=J2;} |
---|
| 2668 | else{ |
---|
| 2669 | m=size(a); |
---|
| 2670 | for (i=1;i<=m;i++){ |
---|
| 2671 | if (flag[a[i]]==0){aux2=poly(x(a[i]));} |
---|
| 2672 | else {aux2=poly(y(a[i]));} |
---|
| 2673 | |
---|
| 2674 | acenter=acenter+aux2; |
---|
| 2675 | } |
---|
| 2676 | } |
---|
| 2677 | return(acenter); |
---|
| 2678 | } |
---|
| 2679 | example |
---|
| 2680 | {"EXAMPLE:"; echo = 2; |
---|
| 2681 | ring r = 0,(x(1..3)),dp; |
---|
| 2682 | list flag=identifyvar(); |
---|
| 2683 | ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3; |
---|
| 2684 | intvec a=1,3; // first center of blowing up |
---|
| 2685 | tradtoideal(a,J,flag); |
---|
| 2686 | } |
---|
| 2687 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
| 2688 | // OPERATIONS WITH LISTS |
---|
| 2689 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
| 2690 | |
---|
| 2691 | proc iniD(int n) |
---|
| 2692 | "USAGE: iniD(n); n integer |
---|
| 2693 | RETURN: list of lists of zeros of size n |
---|
| 2694 | EXAMPLE: example iniD; shows an example |
---|
| 2695 | " |
---|
| 2696 | {int i,j; |
---|
| 2697 | list D,auxD; |
---|
| 2698 | for (j=1;j<=n; j++) {auxD[j]=0;} |
---|
| 2699 | for (i=1;i<=n; i++) {D[i]=auxD;} |
---|
| 2700 | return(D); |
---|
| 2701 | } |
---|
| 2702 | example |
---|
| 2703 | {"EXAMPLE:"; echo = 2; |
---|
| 2704 | iniD(3); |
---|
| 2705 | } |
---|
| 2706 | ///////////////////////////////////////////////////////// |
---|
| 2707 | |
---|
| 2708 | proc sumlist(list L1,list L2) |
---|
| 2709 | "USAGE: sumlist(L1,L2); L1,L2 lists, (size(L1)==size(L2)) |
---|
| 2710 | RETURN: a list, sum of L1 and L2 |
---|
| 2711 | EXAMPLE: example sumlist; shows an example |
---|
| 2712 | " |
---|
| 2713 | { |
---|
| 2714 | int i,k; |
---|
| 2715 | list sumL; |
---|
| 2716 | k=size(L1); |
---|
| 2717 | if (size(L2)!=k) {return("ERROR en sumlist, lists must have the same size");} |
---|
| 2718 | for (i=1;i<=k;i++) {sumL[i]=L1[i]+L2[i];} |
---|
| 2719 | return(sumL); |
---|
| 2720 | } |
---|
| 2721 | example |
---|
| 2722 | {"EXAMPLE:"; echo = 2; |
---|
| 2723 | list L1=1,2,3; |
---|
| 2724 | list L2=5,9,7; |
---|
| 2725 | sumlist(L1,L2); |
---|
| 2726 | } |
---|
| 2727 | /////////////////////////////////////////////////////// |
---|
| 2728 | |
---|
| 2729 | proc reslist(list L1,list L2) |
---|
| 2730 | "USAGE: reslist(L1,L2); L1,L2 lists, (size(L1)==size(L2)) |
---|
| 2731 | RETURN: a list, subtraction of L1 and L2 |
---|
| 2732 | EXAMPLE: example reslist; shows an example |
---|
| 2733 | " |
---|
| 2734 | { |
---|
| 2735 | int i,k; |
---|
| 2736 | list resL; |
---|
| 2737 | k=size(L1); |
---|
| 2738 | if (size(L2)!=k) {return("ERROR en reslist, lists must have the same size");} |
---|
| 2739 | for (i=1;i<=k;i++) {resL[i]=L1[i]-L2[i];} |
---|
| 2740 | return(resL); |
---|
| 2741 | } |
---|
| 2742 | example |
---|
| 2743 | {"EXAMPLE:"; echo = 2; |
---|
| 2744 | list L1=1,2,3; |
---|
| 2745 | list L2=5,9,7; |
---|
| 2746 | reslist(L1,L2); |
---|
| 2747 | } |
---|
| 2748 | ////////////////////////////////////////////////////// |
---|
| 2749 | |
---|
| 2750 | proc multiplylist(list L,number a) |
---|
| 2751 | "USAGE: multiplylist(L,a); L list, a number |
---|
| 2752 | RETURN: list of elements of type number, multiplication of L times a |
---|
| 2753 | EXAMPLE: example multiplylist; shows an example |
---|
| 2754 | " |
---|
| 2755 | {int i,k; |
---|
| 2756 | list newL,bb; |
---|
| 2757 | number b; |
---|
| 2758 | k=size(L); |
---|
| 2759 | for (i=1;i<=k;i++) {b=L[i]*a; bb=b; newL=newL+bb;} |
---|
| 2760 | return(newL); |
---|
| 2761 | } |
---|
| 2762 | example |
---|
| 2763 | {"EXAMPLE:"; echo = 2; |
---|
| 2764 | ring r = 0,(x(1..3)),dp; |
---|
| 2765 | list L=1,2,3; |
---|
| 2766 | multiplylist(L,1/5); |
---|
| 2767 | } |
---|
| 2768 | /////////////////////////////////////////////////////// |
---|
| 2769 | |
---|
| 2770 | proc dividelist(list L1,list L2) |
---|
| 2771 | "USAGE: dividelist(L1,L2); L1,L2 lists |
---|
| 2772 | RETURN: list of elements of type number, division of L1 by L2 |
---|
| 2773 | EXAMPLE: example dividelist; shows an example |
---|
| 2774 | " |
---|
| 2775 | {int i,k,k1,k2; |
---|
| 2776 | list LL,bb; |
---|
| 2777 | number a1,a2,b; |
---|
| 2778 | k1=size(L1); |
---|
| 2779 | k2=size(L2); |
---|
| 2780 | if (k2!=k1) {print("ERROR en dividelist, lists must have the same size");} |
---|
| 2781 | if (k1<=k2) {k=k1;} |
---|
| 2782 | else {k=k2;} |
---|
| 2783 | for (i=1;i<=k;i++) |
---|
| 2784 | {a1=L1[i]; a2=L2[i]; b=a1/a2; bb=b; LL=LL+bb;} |
---|
| 2785 | return(LL); |
---|
| 2786 | } |
---|
| 2787 | example |
---|
| 2788 | {"EXAMPLE:"; echo = 2; |
---|
| 2789 | ring r = 0,(x(1..3)),dp; |
---|
| 2790 | list L1=1,2,3; |
---|
| 2791 | list L2=5,9,7; |
---|
| 2792 | dividelist(L1,L2); |
---|
| 2793 | } |
---|
| 2794 | /////////////////////////////////////////////////////// |
---|
| 2795 | |
---|
| 2796 | proc createlist(list L1,list L2) |
---|
| 2797 | "USAGE: createlist(L1,L2); L1,L2 lists, (size(L1)==size(L2)) |
---|
| 2798 | RETURN: list of lists of two elements, the first one of L1 and the second of L2 |
---|
| 2799 | EXAMPLE: example createlist; shows an example |
---|
| 2800 | " |
---|
| 2801 | {int i,k; |
---|
| 2802 | list L,aux; |
---|
| 2803 | k=size(L1); |
---|
| 2804 | if (size(L2)!=k) {return("ERROR en createlist, lists must have the same size");} |
---|
| 2805 | L=list0(k); |
---|
| 2806 | for (i=1;i<=k;i++) {if (L1[i]!=0) {aux=L1[i],L2[i]; L[i]=aux;} |
---|
| 2807 | else {L=delete(L,i);}} |
---|
| 2808 | return(L); |
---|
| 2809 | } |
---|
| 2810 | example |
---|
| 2811 | {"EXAMPLE:"; echo = 2; |
---|
| 2812 | list L1=1,2,3; |
---|
| 2813 | list L2=5,9,7; |
---|
| 2814 | createlist(L1,L2); |
---|
| 2815 | } |
---|
| 2816 | /////////////////////////////////////////////////////// |
---|
| 2817 | proc list0(int n) |
---|
| 2818 | "USAGE: list0(n); n integer |
---|
| 2819 | RETURN: list of n zeros |
---|
| 2820 | EXAMPLE: example list0; shows an example |
---|
| 2821 | " |
---|
| 2822 | {int i; |
---|
| 2823 | list L0; |
---|
| 2824 | for (i=1;i<=n;i++) {L0[i]=0;} |
---|
| 2825 | return(L0); |
---|
| 2826 | } |
---|
| 2827 | example |
---|
| 2828 | {"EXAMPLE:"; echo = 2; |
---|
| 2829 | list0(4); |
---|
| 2830 | } |
---|
| 2831 | //////////////////////////////////////////////////////////// |
---|