1 | ////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id: reszeta.lib,v 1.2 2005-05-03 14:22:08 Singular Exp $"; |
---|
3 | category="Commutative Algebra"; |
---|
4 | info=" |
---|
5 | LIBRARY: zeta.lib Zeta-function of Denef and Loeser |
---|
6 | |
---|
7 | AUTHORS: A. Fruehbis-Krueger, anne@mathematik.uni-kl.de, |
---|
8 | @* G. Pfister, pfister@mathematik.uni-kl.de |
---|
9 | |
---|
10 | MAIN PROCEDURES: |
---|
11 | intersectionDiv(L) computes intersection form and genera of exceptional |
---|
12 | divisors (isolated singularities of surfaces) |
---|
13 | spectralNeg(L) computes negative spectral numbers |
---|
14 | (isolated hypersurface singularity) |
---|
15 | zetaDL computes Denef-Loeser zeta function |
---|
16 | (hypersurface singularity of dimension 2) |
---|
17 | |
---|
18 | AUXILLIARY PROCEDURES: |
---|
19 | collectDiv(L[,iv]) identify exceptional divisors in different charts |
---|
20 | (embedded and non-embedded case) |
---|
21 | prepEmbDiv(L[,b]) prepare list of divisors (including components |
---|
22 | of strict transform, embedded case) |
---|
23 | abstractR(L) pass from embedded to non-embedded resolution |
---|
24 | "; |
---|
25 | LIB "resol.lib"; |
---|
26 | LIB "solve.lib"; |
---|
27 | LIB "normal.lib"; |
---|
28 | /////////////////////////////////////////////////////////////////////////////// |
---|
29 | // the repunit: ring R=0,x,dp;number n=10;(n^1031-1)/9; is prim // |
---|
30 | // record: number m=2;number k=m^20996011-1; // |
---|
31 | /////////////////////////////////////////////////////////////////////////////// |
---|
32 | static |
---|
33 | proc spectral1(poly h,list re, list DL,intvec v, intvec n) |
---|
34 | "Internal procedure - no help and no example available |
---|
35 | " |
---|
36 | { |
---|
37 | //--- compute one spectral number |
---|
38 | //--- DL is output of prepEmbDiv |
---|
39 | int i; |
---|
40 | intvec w=computeH(h,re,DL); |
---|
41 | number gw=number(w[1]+v[1])/number(n[1]); |
---|
42 | for(i=2;i<=size(v);i++) |
---|
43 | { |
---|
44 | if(gw>number(w[i]+v[i])/number(n[i])) |
---|
45 | { |
---|
46 | gw=number(w[i]+v[i])/number(n[i]); |
---|
47 | } |
---|
48 | } |
---|
49 | return(gw-1); |
---|
50 | } |
---|
51 | |
---|
52 | /////////////////////////////////////////////////////////////////////////////// |
---|
53 | |
---|
54 | proc spectralNeg(list re,list #) |
---|
55 | "USAGE: spectralNeg(L); |
---|
56 | @* L = list of rings |
---|
57 | ASSUME: L is output of resolution of singularities |
---|
58 | RETURN: list of numbers, each a spectral number in (-1,0] |
---|
59 | EXAMPLE: example spectralNeg; shows an example |
---|
60 | " |
---|
61 | { |
---|
62 | //----------------------------------------------------------------------------- |
---|
63 | // Initialization and Sanity Checks |
---|
64 | //----------------------------------------------------------------------------- |
---|
65 | int i,j,l; |
---|
66 | number bound; |
---|
67 | list resu; |
---|
68 | if(size(#)>0) |
---|
69 | { |
---|
70 | //--- undocumented feature: |
---|
71 | //--- if # is not empty it computes numbers up to this bound, |
---|
72 | //--- not necessarily spectral numbers |
---|
73 | bound=number(#[1]); |
---|
74 | } |
---|
75 | //--- get list of embedded divisors |
---|
76 | list DL=prepEmbDiv(re,1); |
---|
77 | int k=1; |
---|
78 | ideal I,delI; |
---|
79 | number g; |
---|
80 | int m=nvars(basering); |
---|
81 | //--- prepare the multiplicities of exceptional divisors N and nu |
---|
82 | intvec v=computeV(re,DL); // nu |
---|
83 | intvec n=computeN(re,DL); // N |
---|
84 | //--------------------------------------------------------------------------- |
---|
85 | // start computation, first case separately, then loop |
---|
86 | //--------------------------------------------------------------------------- |
---|
87 | resu[1]=spectral1(1,re,DL,v,n); // first number, corresponding to |
---|
88 | // volume form itself |
---|
89 | if(resu[1]>=bound) |
---|
90 | { |
---|
91 | //--- exceeds bound ==> not a spectral number |
---|
92 | resu=delete(resu,1); |
---|
93 | return(resu); |
---|
94 | } |
---|
95 | delI=std(ideal(0)); |
---|
96 | while(k) |
---|
97 | { |
---|
98 | //--- now run through all monomial x volume form, degree by degree |
---|
99 | j++; |
---|
100 | k=0; |
---|
101 | I=maxideal(j); |
---|
102 | I=reduce(I,delI); |
---|
103 | for(i=1;i<=size(I);i++) |
---|
104 | { |
---|
105 | //--- all monomials in degree j |
---|
106 | g=spectral1(I[i],re,DL,v,n); |
---|
107 | if(g<bound) |
---|
108 | { |
---|
109 | //--- otherwise g exceeds bound ==> not a spectral number |
---|
110 | k=1; |
---|
111 | l=1; |
---|
112 | while(resu[l]<g) |
---|
113 | { |
---|
114 | l++; |
---|
115 | if(l==size(resu)+1){break;} |
---|
116 | } |
---|
117 | if(l==size(resu)+1){resu[size(resu)+1]=g;} |
---|
118 | if(resu[l]!=g){resu=insert(resu,g,l-1);} |
---|
119 | } |
---|
120 | else |
---|
121 | { |
---|
122 | delI[size(delI)+1]=I[i]; |
---|
123 | } |
---|
124 | } |
---|
125 | attrib(delI,"isSB",1); |
---|
126 | } |
---|
127 | return(resu); |
---|
128 | } |
---|
129 | example |
---|
130 | {"EXAMPLE:"; |
---|
131 | echo = 2; |
---|
132 | ring R=0,(x,y,z),dp; |
---|
133 | ideal I=x3+y4+z5; |
---|
134 | list L=resolve(I,"K"); |
---|
135 | spectralNeg(L); |
---|
136 | LIB"gaussman.lib"; |
---|
137 | ring r=0,(x,y,z),ds; |
---|
138 | poly f=x3+y4+z5; |
---|
139 | spectrum(f); |
---|
140 | } |
---|
141 | |
---|
142 | /////////////////////////////////////////////////////////////////////////////// |
---|
143 | static |
---|
144 | proc ordE(ideal J,ideal E,ideal W) |
---|
145 | "Internal procedure - no help and no example available |
---|
146 | " |
---|
147 | { |
---|
148 | //--- compute multiplicity of E in J -- embedded in W |
---|
149 | int s; |
---|
150 | if(size(J)==0){~;ERROR("ordE: J=0");} |
---|
151 | ideal Estd=std(E+W); |
---|
152 | ideal Epow=1; |
---|
153 | ideal Jquot=1; |
---|
154 | while(size(reduce(Jquot,Estd))!=0) |
---|
155 | { |
---|
156 | s++; |
---|
157 | Epow=Epow*E; |
---|
158 | Jquot=quotient(Epow+W,J); |
---|
159 | } |
---|
160 | return(s-1); |
---|
161 | } |
---|
162 | |
---|
163 | /////////////////////////////////////////////////////////////////////////////// |
---|
164 | static |
---|
165 | proc computeV(list re, list DL) |
---|
166 | "Internal procedure - no help and no example available |
---|
167 | " |
---|
168 | { |
---|
169 | //--- computes for every divisor E_i its multiplicity + 1 in pi^*(w) |
---|
170 | //--- w a non-vanishing 1-form |
---|
171 | //--- note: DL is output of prepEmbDiv |
---|
172 | |
---|
173 | //----------------------------------------------------------------------------- |
---|
174 | // Initialization |
---|
175 | //----------------------------------------------------------------------------- |
---|
176 | def R=basering; |
---|
177 | int i,j,k,n; |
---|
178 | intvec v,w; |
---|
179 | list iden=DL; |
---|
180 | v[size(iden)]=0; |
---|
181 | |
---|
182 | //---------------------------------------------------------------------------- |
---|
183 | // Run through all exceptional divisors |
---|
184 | //---------------------------------------------------------------------------- |
---|
185 | for(k=1;k<=size(iden);k++) |
---|
186 | { |
---|
187 | for(i=1;i<=size(iden[k]);i++) |
---|
188 | { |
---|
189 | if(defined(S)){kill S;} |
---|
190 | def S=re[2][iden[k][i][1]]; |
---|
191 | setring S; |
---|
192 | if((!v[k])&&(defined(EList))) |
---|
193 | { |
---|
194 | if(defined(II)){kill II;} |
---|
195 | //--- we might be embedded in a non-trivial BO[1] |
---|
196 | //--- take this into account when forming the jacobi-determinant |
---|
197 | ideal II=jacobDet(BO[5],BO[1]); |
---|
198 | if(size(II)!=0) |
---|
199 | { |
---|
200 | v[k]=ordE(II,EList[iden[k][i][2]],BO[1])+1; |
---|
201 | } |
---|
202 | } |
---|
203 | setring R; |
---|
204 | } |
---|
205 | } |
---|
206 | return(v); |
---|
207 | } |
---|
208 | example |
---|
209 | {"EXAMPLE:"; |
---|
210 | echo = 2; |
---|
211 | ring R=0,(x,y,z),dp; |
---|
212 | ideal I=(x-y)*(x-z)*(y-z)-z4; |
---|
213 | list re=resolve(I,1); |
---|
214 | intvec v=computeV(re); |
---|
215 | v; |
---|
216 | } |
---|
217 | |
---|
218 | /////////////////////////////////////////////////////////////////////////////// |
---|
219 | static |
---|
220 | proc jacobDet(ideal I, ideal J) |
---|
221 | "Internal procedure - no help and no example available |
---|
222 | " |
---|
223 | { |
---|
224 | //--- Returns the Jacobian determinant of the morphism |
---|
225 | //--- K[x_1,...,x_m]--->K[y_1,...,y_n]/J defined by x_i ---> I_i. |
---|
226 | //--- Let basering=K[y_1,...,y_n], l=n-dim(basering/J), |
---|
227 | //--- I=<I_1,...,I_m>, J=<J_1,...,J_r> |
---|
228 | //--- For each subset v in {1,...,n} of l elements and |
---|
229 | //--- w in {1,...,r} of l elements |
---|
230 | //--- let K_v,w be the ideal generated by the n-l-minors of the matrix |
---|
231 | //--- (diff(I_i,y_j)+ |
---|
232 | //--- \sum_k diff(I_i,y_v[k])*diff(J_w[k],y_j))_{j not in v multiplied with |
---|
233 | //--- the determinant of (diff(J_w[i],y_v[j])) |
---|
234 | //--- the sum of all such ideals K_v,w plus J is returned. |
---|
235 | |
---|
236 | //---------------------------------------------------------------------------- |
---|
237 | // Initialization |
---|
238 | //---------------------------------------------------------------------------- |
---|
239 | int n=nvars(basering); |
---|
240 | int i,j,k; |
---|
241 | intvec u,v,w,x; |
---|
242 | matrix MI[ncols(I)][n]=jacob(I); |
---|
243 | matrix N=unitmat(n); |
---|
244 | matrix L; |
---|
245 | ideal K=J; |
---|
246 | if(size(J)==0) |
---|
247 | { |
---|
248 | K=minor(MI,n); |
---|
249 | } |
---|
250 | |
---|
251 | //--------------------------------------------------------------------------- |
---|
252 | // Do calculation as described above. |
---|
253 | // separately for case size(J)=1 |
---|
254 | //--------------------------------------------------------------------------- |
---|
255 | if(size(J)==1) |
---|
256 | { |
---|
257 | matrix MJ[ncols(J)][n]=jacob(J); |
---|
258 | N=concat(N,transpose(MJ)); |
---|
259 | v=1..n; |
---|
260 | for(i=1;i<=n;i++) |
---|
261 | { |
---|
262 | L=transpose(permcol(N,i,n+1)); |
---|
263 | if(i==1){w=2..n;} |
---|
264 | if(i==n){w=1..n-1;} |
---|
265 | if((i!=1)&&(i!=n)){w=1..i-1,i+1..n;} |
---|
266 | L=submat(L,v,w); |
---|
267 | L=MI*L; |
---|
268 | K=K+minor(L,n-1)*MJ[1,i]; |
---|
269 | } |
---|
270 | } |
---|
271 | if(size(J)>1) |
---|
272 | { |
---|
273 | matrix MJ[ncols(J)][n]=jacob(J); |
---|
274 | matrix SMJ; |
---|
275 | N=concat(N,transpose(MJ)); |
---|
276 | ideal Jstd=std(J); |
---|
277 | int l=n-dim(Jstd); |
---|
278 | int r=ncols(J); |
---|
279 | list L1=indexSet(n,l); |
---|
280 | list L2=indexSet(r,l); |
---|
281 | for(i=1;i<=size(L1);i++) |
---|
282 | { |
---|
283 | for(j=1;j<=size(L2);j++) |
---|
284 | { |
---|
285 | for(k=1;k<=size(L1[i]);k++) |
---|
286 | { |
---|
287 | if(L1[i][k]){v[size(v)+1]=k;} |
---|
288 | } |
---|
289 | v=v[2..size(v)]; |
---|
290 | for(k=1;k<=size(L2[j]);k++) |
---|
291 | { |
---|
292 | if(L2[j][k]){w[size(w)+1]=k;} |
---|
293 | } |
---|
294 | w=w[2..size(w)]; |
---|
295 | SMJ=submat(MJ,w,v); |
---|
296 | L=N; |
---|
297 | for(k=1;k<=l;k++) |
---|
298 | { |
---|
299 | L=permcol(L,v[k],n+w[k]); |
---|
300 | } |
---|
301 | u=1..n; |
---|
302 | x=1..n; |
---|
303 | v=sort(v)[1]; |
---|
304 | for(k=l;k>=1;k--) |
---|
305 | { |
---|
306 | if(v[k]) |
---|
307 | { |
---|
308 | u=deleteInt(u,v[k],1); |
---|
309 | } |
---|
310 | } |
---|
311 | L=transpose(submat(L,u,x)); |
---|
312 | L=MI*L; |
---|
313 | K=K+minor(L,n-l)*det(SMJ); |
---|
314 | } |
---|
315 | } |
---|
316 | } |
---|
317 | return(K); |
---|
318 | } |
---|
319 | |
---|
320 | /////////////////////////////////////////////////////////////////////////////// |
---|
321 | static |
---|
322 | proc computeH(ideal h,list re,list DL) |
---|
323 | "Internal procedure - no help and no example available |
---|
324 | " |
---|
325 | { |
---|
326 | //--- additional procedure to computeV, allows |
---|
327 | //--- computation for polynomial x volume form |
---|
328 | //--- by computing the contribution of the polynomial h |
---|
329 | //--- Note: DL is output of prepEmbDiv |
---|
330 | |
---|
331 | //---------------------------------------------------------------------------- |
---|
332 | // Initialization |
---|
333 | //---------------------------------------------------------------------------- |
---|
334 | def R=basering; |
---|
335 | ideal II=h; |
---|
336 | list iden=DL; |
---|
337 | def T=re[2][1]; |
---|
338 | setring T; |
---|
339 | int i,k; |
---|
340 | intvec v; |
---|
341 | v[size(iden)]=0; |
---|
342 | if(deg(II[1])==0){return(v);} |
---|
343 | //---------------------------------------------------------------------------- |
---|
344 | // Run through all exceptional divisors |
---|
345 | //---------------------------------------------------------------------------- |
---|
346 | for(k=1;k<=size(iden);k++) |
---|
347 | { |
---|
348 | for(i=1;i<=size(iden[k]);i++) |
---|
349 | { |
---|
350 | if(defined(S)){kill S;} |
---|
351 | def S=re[2][iden[k][i][1]]; |
---|
352 | setring S; |
---|
353 | if((!v[k])&&(defined(EList))) |
---|
354 | { |
---|
355 | if(defined(JJ)){kill JJ;} |
---|
356 | if(defined(phi)){kill phi;} |
---|
357 | map phi=T,BO[5]; |
---|
358 | ideal JJ=phi(II); |
---|
359 | if(size(JJ)!=0) |
---|
360 | { |
---|
361 | v[k]=ordE(JJ,EList[iden[k][i][2]],BO[1]); |
---|
362 | } |
---|
363 | } |
---|
364 | setring R; |
---|
365 | } |
---|
366 | } |
---|
367 | return(v); |
---|
368 | } |
---|
369 | |
---|
370 | ////////////////////////////////////////////////////////////////////////////// |
---|
371 | static |
---|
372 | proc computeN(list re,list DL) |
---|
373 | "Internal procedure - no help and no example available |
---|
374 | " |
---|
375 | { |
---|
376 | //--- computes for every (Q-irred.) divisor E_i its multiplicity in f \circ pi |
---|
377 | //--- DL is output of prepEmbDiv |
---|
378 | |
---|
379 | //---------------------------------------------------------------------------- |
---|
380 | // Initialization |
---|
381 | //---------------------------------------------------------------------------- |
---|
382 | def R=basering; |
---|
383 | list iden=DL; |
---|
384 | def T=re[2][1]; |
---|
385 | setring T; |
---|
386 | ideal J=BO[2]; |
---|
387 | int i,k; |
---|
388 | intvec v; |
---|
389 | v[size(iden)]=0; |
---|
390 | //---------------------------------------------------------------------------- |
---|
391 | // Run through all exceptional divisors |
---|
392 | //---------------------------------------------------------------------------- |
---|
393 | for(k=1;k<=size(iden);k++) |
---|
394 | { |
---|
395 | for(i=1;i<=size(iden[k]);i++) |
---|
396 | { |
---|
397 | if(defined(S)){kill S;} |
---|
398 | def S=re[2][iden[k][i][1]]; |
---|
399 | setring S; |
---|
400 | if((!v[k])&&(defined(EList))) |
---|
401 | { |
---|
402 | if(defined(II)){kill II;} |
---|
403 | if(defined(phi)){kill phi;} |
---|
404 | map phi=T,BO[5]; |
---|
405 | ideal II=phi(J); |
---|
406 | if(size(II)!=0) |
---|
407 | { |
---|
408 | v[k]=ordE(II,EList[iden[k][i][2]],BO[1]); |
---|
409 | } |
---|
410 | } |
---|
411 | setring R; |
---|
412 | } |
---|
413 | } |
---|
414 | return(v); |
---|
415 | } |
---|
416 | example |
---|
417 | {"EXAMPLE:"; |
---|
418 | echo = 2; |
---|
419 | ring R=0,(x,y,z),dp; |
---|
420 | ideal I=(x-y)*(x-z)*(y-z)-z4; |
---|
421 | list re=resolve(I,1); |
---|
422 | intvec v=computeN(re); |
---|
423 | v; |
---|
424 | } |
---|
425 | ////////////////////////////////////////////////////////////////////////////// |
---|
426 | static |
---|
427 | proc countEijk(list re,list iden,intvec iv,list #) |
---|
428 | "Internal procedure - no help and no example available |
---|
429 | " |
---|
430 | { |
---|
431 | //--- count the number of points in the intersection of 3 exceptional |
---|
432 | //--- hyperplanes (of dimension 2) - one of them is allowed to be a component |
---|
433 | //--- of the strict transform |
---|
434 | |
---|
435 | //---------------------------------------------------------------------------- |
---|
436 | // Initialization |
---|
437 | //---------------------------------------------------------------------------- |
---|
438 | int i,j,k,comPa,numPts,localCase; |
---|
439 | intvec ituple,jtuple,ktuple; |
---|
440 | list chList,tmpList; |
---|
441 | def R=basering; |
---|
442 | if(size(#)>0) |
---|
443 | { |
---|
444 | if(string(#[1])=="local") |
---|
445 | { |
---|
446 | localCase=1; |
---|
447 | } |
---|
448 | } |
---|
449 | //---------------------------------------------------------------------------- |
---|
450 | // Find common charts |
---|
451 | //---------------------------------------------------------------------------- |
---|
452 | for(i=1;i<=size(iden[iv[1]]);i++) |
---|
453 | { |
---|
454 | //--- find common charts - only for final charts |
---|
455 | if(defined(S)) {kill S;} |
---|
456 | def S=re[2][iden[iv[1]][i][1]]; |
---|
457 | setring S; |
---|
458 | if(!defined(EList)) |
---|
459 | { |
---|
460 | i++; |
---|
461 | setring R; |
---|
462 | continue; |
---|
463 | } |
---|
464 | setring R; |
---|
465 | kill ituple,jtuple,ktuple; |
---|
466 | intvec ituple=iden[iv[1]][i]; |
---|
467 | intvec jtuple=findInIVList(1,ituple[1],iden[iv[2]]); |
---|
468 | intvec ktuple=findInIVList(1,ituple[1],iden[iv[3]]); |
---|
469 | if((size(jtuple)!=1)&&(size(ktuple)!=1)) |
---|
470 | { |
---|
471 | //--- chList contains all information about the common charts, |
---|
472 | //--- each entry represents a chart and contains three intvecs from iden |
---|
473 | //--- one for each E_l |
---|
474 | kill tmpList; |
---|
475 | list tmpList=ituple,jtuple,ktuple; |
---|
476 | chList[size(chList)+1]=tmpList; |
---|
477 | i++; |
---|
478 | if(i<=size(iden[iv[1]])) |
---|
479 | { |
---|
480 | continue; |
---|
481 | } |
---|
482 | else |
---|
483 | { |
---|
484 | break; |
---|
485 | } |
---|
486 | } |
---|
487 | } |
---|
488 | if(size(chList)==0) |
---|
489 | { |
---|
490 | //--- no common chart !!! |
---|
491 | return(int(0)); |
---|
492 | } |
---|
493 | //---------------------------------------------------------------------------- |
---|
494 | // Count points in common charts |
---|
495 | //---------------------------------------------------------------------------- |
---|
496 | for(i=1;i<=size(chList);i++) |
---|
497 | { |
---|
498 | //--- run through all common charts |
---|
499 | if(defined(S)) { kill S;} |
---|
500 | def S=re[2][chList[i][1][1]]; |
---|
501 | setring S; |
---|
502 | //--- intersection in this chart |
---|
503 | if(defined(interId)){kill interId;} |
---|
504 | if(localCase==1) |
---|
505 | { |
---|
506 | //--- in this case we need to intersect with \pi^-1(0) |
---|
507 | ideal interId=EList[chList[i][1][2]]+EList[chList[i][2][2]] |
---|
508 | +EList[chList[i][3][2]]+BO[5]; |
---|
509 | } |
---|
510 | else |
---|
511 | { |
---|
512 | ideal interId=EList[chList[i][1][2]]+EList[chList[i][2][2]] |
---|
513 | +EList[chList[i][3][2]]; |
---|
514 | } |
---|
515 | interId=std(interId); |
---|
516 | if(defined(otherId)) {kill otherId;} |
---|
517 | ideal otherId=1; |
---|
518 | for(j=1;j<i;j++) |
---|
519 | { |
---|
520 | //--- run through the previously computed ones |
---|
521 | if(defined(opath)){kill opath;} |
---|
522 | def opath=imap(re[2][chList[j][1][1]],path); |
---|
523 | comPa=1; |
---|
524 | while(opath[1,comPa]==path[1,comPa]) |
---|
525 | { |
---|
526 | comPa++; |
---|
527 | if((comPa>ncols(path))||(comPa>ncols(opath))) break; |
---|
528 | } |
---|
529 | comPa=int(leadcoef(path[1,comPa-1])); |
---|
530 | otherId=otherId+interId; |
---|
531 | otherId=intersect(otherId, |
---|
532 | fetchInTree(re,chList[j][1][1], |
---|
533 | comPa,chList[i][1][1],"interId",iden)); |
---|
534 | } |
---|
535 | otherId=std(otherId); |
---|
536 | //--- do not count each point more than once |
---|
537 | interId=sat(interId,otherId)[1]; |
---|
538 | export(interId); |
---|
539 | if(dim(interId)>0) |
---|
540 | { |
---|
541 | ERROR("CountEijk: intersection not zerodimensional"); |
---|
542 | } |
---|
543 | //--- add the remaining number of points to the total point count numPts |
---|
544 | numPts=numPts+vdim(interId); |
---|
545 | } |
---|
546 | return(numPts); |
---|
547 | } |
---|
548 | ////////////////////////////////////////////////////////////////////////////// |
---|
549 | static |
---|
550 | proc chiEij(list re, list iden, intvec iv) |
---|
551 | "Internal procedure - no help and no example available |
---|
552 | " |
---|
553 | { |
---|
554 | //!!! Copy of chiEij_local adjusted for non-local case |
---|
555 | //!!! changes must be made in both copies |
---|
556 | |
---|
557 | //--- compute the Euler characteristic of the intersection |
---|
558 | //--- curve of two exceptional hypersurfaces (of dimension 2) |
---|
559 | //--- one of which is allowed to be a component of the strict transform |
---|
560 | //--- using the formula chi(Eij)=2-2g(Eij) |
---|
561 | |
---|
562 | //---------------------------------------------------------------------------- |
---|
563 | // Initialization |
---|
564 | //---------------------------------------------------------------------------- |
---|
565 | int i,j,k,chi,g; |
---|
566 | intvec ituple,jtuple,inters; |
---|
567 | def R=basering; |
---|
568 | //---------------------------------------------------------------------------- |
---|
569 | // Find a common chart in which they intersect |
---|
570 | //---------------------------------------------------------------------------- |
---|
571 | for(i=1;i<=size(iden[iv[1]]);i++) |
---|
572 | { |
---|
573 | //--- find a common chart in which they intersect: only for final charts |
---|
574 | if(defined(S)) {kill S;} |
---|
575 | def S=re[2][iden[iv[1]][i][1]]; |
---|
576 | setring S; |
---|
577 | if(!defined(EList)) |
---|
578 | { |
---|
579 | i++; |
---|
580 | setring R; |
---|
581 | continue; |
---|
582 | } |
---|
583 | setring R; |
---|
584 | kill ituple,jtuple; |
---|
585 | intvec ituple=iden[iv[1]][i]; |
---|
586 | intvec jtuple=findInIVList(1,ituple[1],iden[iv[2]]); |
---|
587 | if(size(jtuple)==1) |
---|
588 | { |
---|
589 | if(i<size(iden[iv[1]])) |
---|
590 | { |
---|
591 | //--- not in this chart |
---|
592 | i++; |
---|
593 | continue; |
---|
594 | } |
---|
595 | else |
---|
596 | { |
---|
597 | if(size(inters)==1) |
---|
598 | { |
---|
599 | //--- E_i and E_j do not meet at all |
---|
600 | return("leer"); |
---|
601 | } |
---|
602 | else |
---|
603 | { |
---|
604 | return(chi); |
---|
605 | } |
---|
606 | } |
---|
607 | } |
---|
608 | //---------------------------------------------------------------------------- |
---|
609 | // Run through common charts and compute the Euler characteristic of |
---|
610 | // each component of Eij. |
---|
611 | // As soon as a component has been treated in a chart, it will not be used in |
---|
612 | // any subsequent charts. |
---|
613 | //---------------------------------------------------------------------------- |
---|
614 | if(defined(S)) {kill S;} |
---|
615 | def S=re[2][ituple[1]]; |
---|
616 | setring S; |
---|
617 | //--- interId: now all components in this chart, |
---|
618 | //--- but we want only new components |
---|
619 | if(defined(interId)){kill interId;} |
---|
620 | ideal interId=EList[ituple[2]]+EList[jtuple[2]]; |
---|
621 | interId=std(interId); |
---|
622 | //--- doneId: already considered components |
---|
623 | if(defined(doneId)){kill doneId;} |
---|
624 | ideal doneId=1; |
---|
625 | for(j=2;j<=size(inters);j++) |
---|
626 | { |
---|
627 | //--- fetch the components which have already been dealt with via fetchInTree |
---|
628 | if(defined(opath)) {kill opath;} |
---|
629 | def opath=imap(re[2][inters[j]],path); |
---|
630 | k=1; |
---|
631 | while((k<ncols(opath))&&(k<ncols(path))) |
---|
632 | { |
---|
633 | if(path[1,k+1]!=opath[1,k+1]) break; |
---|
634 | k++; |
---|
635 | } |
---|
636 | if(defined(comPa)) {kill comPa;} |
---|
637 | int comPa=int(leadcoef(path[1,k])); |
---|
638 | if(defined(tempId)){kill tempId;} |
---|
639 | ideal tempId=fetchInTree(re,inters[j],comPa, |
---|
640 | iden[iv[1]][i][1],"interId",iden); |
---|
641 | doneId=intersect(doneId,tempId); |
---|
642 | kill tempId; |
---|
643 | } |
---|
644 | //--- only consider new components in interId |
---|
645 | interId=sat(interId,doneId)[1]; |
---|
646 | if(dim(interId)>1) |
---|
647 | { |
---|
648 | ERROR("genus_Eij: higher dimensional intersection"); |
---|
649 | } |
---|
650 | if(dim(interId)>=0) |
---|
651 | { |
---|
652 | //--- save the index of the current chart for future use |
---|
653 | export(interId); |
---|
654 | inters[size(inters)+1]=iden[iv[1]][i][1]; |
---|
655 | } |
---|
656 | BO[1]=std(BO[1]); |
---|
657 | if(((dim(interId)<=0)&&(dim(BO[1])>2))|| |
---|
658 | ((dim(interId)<0)&&(dim(BO[1])==2))) |
---|
659 | { |
---|
660 | if(i<size(iden[iv[1]])) |
---|
661 | { |
---|
662 | //--- not in this chart |
---|
663 | setring R; |
---|
664 | i++; |
---|
665 | continue; |
---|
666 | } |
---|
667 | else |
---|
668 | { |
---|
669 | if(size(inters)==1) |
---|
670 | { |
---|
671 | //--- E_i and E_j do not meet at all |
---|
672 | return("leer"); |
---|
673 | } |
---|
674 | else |
---|
675 | { |
---|
676 | return(chi); |
---|
677 | } |
---|
678 | } |
---|
679 | } |
---|
680 | g=genus(interId); |
---|
681 | |
---|
682 | //--- chi is the Euler characteristic of the (disjoint !!!) union of the |
---|
683 | //--- considered components |
---|
684 | //--- remark: components are disjoint, because the E_i are normal crossing!!! |
---|
685 | |
---|
686 | chi=chi+(2-2*g); |
---|
687 | } |
---|
688 | return(chi); |
---|
689 | } |
---|
690 | ////////////////////////////////////////////////////////////////////////////// |
---|
691 | ////////////////////////////////////////////////////////////////////////////// |
---|
692 | static |
---|
693 | proc chiEij_local(list re, list iden, intvec iv) |
---|
694 | "Internal procedure - no help and no example available |
---|
695 | " |
---|
696 | { |
---|
697 | //!!! Copy of chiEij adjusted for local case |
---|
698 | //!!! changes must be made in both copies |
---|
699 | |
---|
700 | //--- we have to consider two different cases: |
---|
701 | //--- case1: E_i \cap E_j \cap \pi^-1(0) is a curve |
---|
702 | //--- compute the Euler characteristic of the intersection |
---|
703 | //--- curve of two exceptional hypersurfaces (of dimension 2) |
---|
704 | //--- one of which is allowed to be a component of the strict transform |
---|
705 | //--- using the formula chi(Eij)=2-2g(Eij) |
---|
706 | //--- case2: E_i \cap E_j \cap \pi^-1(0) is a set of points |
---|
707 | //--- count the points |
---|
708 | |
---|
709 | //---------------------------------------------------------------------------- |
---|
710 | // Initialization |
---|
711 | //---------------------------------------------------------------------------- |
---|
712 | int i,j,k,chi,g,points; |
---|
713 | intvec ituple,jtuple,inters; |
---|
714 | def R=basering; |
---|
715 | //---------------------------------------------------------------------------- |
---|
716 | // Find a common chart in which they intersect |
---|
717 | //---------------------------------------------------------------------------- |
---|
718 | for(i=1;i<=size(iden[iv[1]]);i++) |
---|
719 | { |
---|
720 | //--- find a common chart in which they intersect: only for final charts |
---|
721 | if(defined(S)) {kill S;} |
---|
722 | def S=re[2][iden[iv[1]][i][1]]; |
---|
723 | setring S; |
---|
724 | if(!defined(EList)) |
---|
725 | { |
---|
726 | i++; |
---|
727 | setring R; |
---|
728 | continue; |
---|
729 | } |
---|
730 | setring R; |
---|
731 | kill ituple,jtuple; |
---|
732 | intvec ituple=iden[iv[1]][i]; |
---|
733 | intvec jtuple=findInIVList(1,ituple[1],iden[iv[2]]); |
---|
734 | if(size(jtuple)==1) |
---|
735 | { |
---|
736 | if(i<size(iden[iv[1]])) |
---|
737 | { |
---|
738 | //--- not in this chart |
---|
739 | i++; |
---|
740 | continue; |
---|
741 | } |
---|
742 | else |
---|
743 | { |
---|
744 | if(size(inters)==1) |
---|
745 | { |
---|
746 | //--- E_i and E_j do not meet at all |
---|
747 | return("leer"); |
---|
748 | } |
---|
749 | else |
---|
750 | { |
---|
751 | return(chi); |
---|
752 | } |
---|
753 | } |
---|
754 | } |
---|
755 | //---------------------------------------------------------------------------- |
---|
756 | // Run through common charts and compute the Euler characteristic of |
---|
757 | // each component of Eij. |
---|
758 | // As soon as a component has been treated in a chart, it will not be used in |
---|
759 | // any subsequent charts. |
---|
760 | //---------------------------------------------------------------------------- |
---|
761 | if(defined(S)) {kill S;} |
---|
762 | def S=re[2][ituple[1]]; |
---|
763 | setring S; |
---|
764 | //--- interId: now all components in this chart, |
---|
765 | //--- but we want only new components |
---|
766 | if(defined(interId)){kill interId;} |
---|
767 | ideal interId=EList[ituple[2]]+EList[jtuple[2]]+BO[5]; |
---|
768 | interId=std(interId); |
---|
769 | //--- doneId: already considered components |
---|
770 | if(defined(doneId)){kill doneId;} |
---|
771 | ideal doneId=1; |
---|
772 | for(j=2;j<=size(inters);j++) |
---|
773 | { |
---|
774 | //--- fetch the components which have already been dealt with via fetchInTree |
---|
775 | if(defined(opath)) {kill opath;} |
---|
776 | def opath=imap(re[2][inters[j]],path); |
---|
777 | k=1; |
---|
778 | while((k<ncols(opath))&&(k<ncols(path))) |
---|
779 | { |
---|
780 | if(path[1,k+1]!=opath[1,k+1]) break; |
---|
781 | k++; |
---|
782 | } |
---|
783 | if(defined(comPa)) {kill comPa;} |
---|
784 | int comPa=int(leadcoef(path[1,k])); |
---|
785 | if(defined(tempId)){kill tempId;} |
---|
786 | ideal tempId=fetchInTree(re,inters[j],comPa, |
---|
787 | iden[iv[1]][i][1],"interId",iden); |
---|
788 | doneId=intersect(doneId,tempId); |
---|
789 | kill tempId; |
---|
790 | } |
---|
791 | //--- only consider new components in interId |
---|
792 | interId=sat(interId,doneId)[1]; |
---|
793 | if(dim(interId)>1) |
---|
794 | { |
---|
795 | ERROR("genus_Eij: higher dimensional intersection"); |
---|
796 | } |
---|
797 | if(dim(interId)>=0) |
---|
798 | { |
---|
799 | //--- save the index of the current chart for future use |
---|
800 | export(interId); |
---|
801 | inters[size(inters)+1]=iden[iv[1]][i][1]; |
---|
802 | } |
---|
803 | BO[1]=std(BO[1]); |
---|
804 | if(dim(interId)<0) |
---|
805 | { |
---|
806 | if(i<size(iden[iv[1]])) |
---|
807 | { |
---|
808 | //--- not in this chart |
---|
809 | setring R; |
---|
810 | i++; |
---|
811 | continue; |
---|
812 | } |
---|
813 | else |
---|
814 | { |
---|
815 | if(size(inters)==1) |
---|
816 | { |
---|
817 | //--- E_i and E_j do not meet at all |
---|
818 | return("leer"); |
---|
819 | } |
---|
820 | else |
---|
821 | { |
---|
822 | return(chi); |
---|
823 | } |
---|
824 | } |
---|
825 | } |
---|
826 | if((dim(interId)==0)&&(dim(std(BO[1]))>2)) |
---|
827 | { |
---|
828 | //--- for sets of points the Euler characteristic is just |
---|
829 | //--- the number of points |
---|
830 | //--- fat points are impossible, since everything is smooth and n.c. |
---|
831 | chi=chi+vdim(interId); |
---|
832 | points=1; |
---|
833 | } |
---|
834 | else |
---|
835 | { |
---|
836 | if(points==1) |
---|
837 | { |
---|
838 | ERROR("components of intersection do not have same dimension"); |
---|
839 | } |
---|
840 | g=genus(interId); |
---|
841 | //--- chi is the Euler characteristic of the (disjoint !!!) union of the |
---|
842 | //--- considered components |
---|
843 | //--- remark: components are disjoint, because the E_i are normal crossing!!! |
---|
844 | chi=chi+(2-2*g); |
---|
845 | } |
---|
846 | } |
---|
847 | return(chi); |
---|
848 | } |
---|
849 | ////////////////////////////////////////////////////////////////////////////// |
---|
850 | static |
---|
851 | proc computeChiE(list re, list iden) |
---|
852 | "Internal procedure - no help and no example available |
---|
853 | " |
---|
854 | { |
---|
855 | //--- compute the Euler characteristic of the exceptional hypersurfaces |
---|
856 | //--- (of dimension 2), not considering the components of the strict |
---|
857 | //--- transform |
---|
858 | |
---|
859 | //---------------------------------------------------------------------------- |
---|
860 | // Initialization |
---|
861 | //---------------------------------------------------------------------------- |
---|
862 | int i,j,k,m,thisE,otherE; |
---|
863 | def R=basering; |
---|
864 | intvec nulliv,chi_temp,kvec; |
---|
865 | nulliv[size(iden)]=0; |
---|
866 | list chi_E; |
---|
867 | for(i=1;i<=size(iden);i++) |
---|
868 | { |
---|
869 | chi_E[i]=list(); |
---|
870 | } |
---|
871 | //--------------------------------------------------------------------------- |
---|
872 | // Run through the list of charts and compute the Euler characteristic of |
---|
873 | // the new exceptional hypersurface and change the values for the old ones |
---|
874 | // according to the blow-up which has just been performed |
---|
875 | // For initialization reasons, treat the case of the first blow-up separately |
---|
876 | //--------------------------------------------------------------------------- |
---|
877 | for(i=2;i<=size(re[2]);i++) |
---|
878 | { |
---|
879 | //--- run through all charts |
---|
880 | if(defined(S)){kill S;} |
---|
881 | def S=re[2][i]; |
---|
882 | setring S; |
---|
883 | m=int(leadcoef(path[1,ncols(path)])); |
---|
884 | if(defined(Spa)){kill Spa;} |
---|
885 | def Spa=re[2][m]; |
---|
886 | if(size(BO[4])==1) |
---|
887 | { |
---|
888 | //--- just one exceptional divisor |
---|
889 | thisE=1; |
---|
890 | setring Spa; |
---|
891 | if(i==2) |
---|
892 | { |
---|
893 | //--- have not set the initial value of chi(E_1) yet |
---|
894 | if(dim(std(cent))==0) |
---|
895 | { |
---|
896 | //--- center was point ==> new except. div. is a P^2 |
---|
897 | list templist=3*vdim(std(BO[1]+cent)),nulliv; |
---|
898 | } |
---|
899 | else |
---|
900 | { |
---|
901 | //--- center was curve ==> new except. div. is curve x P^1 |
---|
902 | list templist=4-4*genus(BO[1]+cent),nulliv; |
---|
903 | } |
---|
904 | chi_E[1]=templist; |
---|
905 | kill templist; |
---|
906 | } |
---|
907 | setring S; |
---|
908 | i++; |
---|
909 | if(i<size(re[2])) |
---|
910 | { |
---|
911 | continue; |
---|
912 | } |
---|
913 | else |
---|
914 | { |
---|
915 | break; |
---|
916 | } |
---|
917 | } |
---|
918 | for(j=1;j<=size(iden);j++) |
---|
919 | { |
---|
920 | //--- find out which exceptional divisor has just been born |
---|
921 | if(inIVList(intvec(i,size(BO[4])),iden[j])) |
---|
922 | { |
---|
923 | //--- found it |
---|
924 | thisE=j; |
---|
925 | break; |
---|
926 | } |
---|
927 | } |
---|
928 | //--- now setup new chi and change the previous ones appropriately |
---|
929 | setring Spa; |
---|
930 | if(size(chi_E[thisE])==0) |
---|
931 | { |
---|
932 | //--- have not set the initial value of chi(E_thisE) yet |
---|
933 | if(dim(std(cent))==0) |
---|
934 | { |
---|
935 | //--- center was point ==> new except. div. is a P^2 |
---|
936 | list templist=3*vdim(std(BO[1]+cent)),nulliv; |
---|
937 | } |
---|
938 | else |
---|
939 | { |
---|
940 | //--- center was curve ==> new except. div. is a C x P^1 |
---|
941 | list templist=4-4*genus(BO[1]+cent),nulliv; |
---|
942 | } |
---|
943 | chi_E[thisE]=templist; |
---|
944 | kill templist; |
---|
945 | } |
---|
946 | for(j=1;j<=size(BO[4]);j++) |
---|
947 | { |
---|
948 | //--- we are in the parent ring ==> thisE is not yet born |
---|
949 | //--- all the other E_i have already been initialized, but the chi |
---|
950 | //--- might change with the current blow-up at cent |
---|
951 | if(BO[6][j]==1) |
---|
952 | { |
---|
953 | //--- ignore empty sets |
---|
954 | j++; |
---|
955 | if(j<=size(BO[4])) |
---|
956 | { |
---|
957 | continue; |
---|
958 | } |
---|
959 | else |
---|
960 | { |
---|
961 | break; |
---|
962 | } |
---|
963 | } |
---|
964 | for(k=1;k<=size(iden);k++) |
---|
965 | { |
---|
966 | //--- find global index of BO[4][j] |
---|
967 | if(inIVList(intvec(m,j),iden[k])) |
---|
968 | { |
---|
969 | otherE=k; |
---|
970 | break; |
---|
971 | } |
---|
972 | } |
---|
973 | if(chi_E[otherE][2][thisE]==1) |
---|
974 | { |
---|
975 | //--- already considered this one |
---|
976 | j++; |
---|
977 | if(j<=size(BO[4])) |
---|
978 | { |
---|
979 | continue; |
---|
980 | } |
---|
981 | else |
---|
982 | { |
---|
983 | break; |
---|
984 | } |
---|
985 | } |
---|
986 | //--------------------------------------------------------------------------- |
---|
987 | // update chi according to the formula |
---|
988 | // chi(E_k^transf)=chi(E_k) - chi(C \cap E_k) + chi(E_k \cap E_new) |
---|
989 | // for convenience of implementation, we first compute |
---|
990 | // chi(E_k) - chi(C \cap E_k) |
---|
991 | // and afterwards add the last term chi(E_k \cap E_new) |
---|
992 | //--------------------------------------------------------------------------- |
---|
993 | ideal CinE=std(cent+BO[4][j]+BO[1]); // this is C \cap E_k |
---|
994 | if(dim(CinE)==1) |
---|
995 | { |
---|
996 | //--- center meets E_k in a curve |
---|
997 | chi_temp[otherE]=chi_E[otherE][1]-(2-2*genus(CinE)); |
---|
998 | } |
---|
999 | if(dim(CinE)==0) |
---|
1000 | { |
---|
1001 | //--- center meets E_k in points |
---|
1002 | chi_temp[otherE]=chi_E[otherE][1]-vdim(std(CinE)); |
---|
1003 | } |
---|
1004 | kill CinE; |
---|
1005 | setring S; |
---|
1006 | //--- now we are back in the i-th ring in the list |
---|
1007 | ideal CinE=std(BO[4][j]+BO[4][size(BO[4])]+BO[1]); |
---|
1008 | // this is E_k \cap E_new |
---|
1009 | if(dim(CinE)==1) |
---|
1010 | { |
---|
1011 | //--- if the two divisors meet, they meet in a curve |
---|
1012 | chi_E[otherE][1]=chi_temp[otherE]+(2-2*genus(CinE)); |
---|
1013 | chi_E[otherE][2][thisE]=1; // this blow-up of E_k is done |
---|
1014 | } |
---|
1015 | kill CinE; |
---|
1016 | setring Spa; |
---|
1017 | } |
---|
1018 | } |
---|
1019 | setring R; |
---|
1020 | return(chi_E); |
---|
1021 | } |
---|
1022 | ////////////////////////////////////////////////////////////////////////////// |
---|
1023 | static |
---|
1024 | proc computeChiE_local(list re, list iden) |
---|
1025 | "Internal procedure - no help and no example available |
---|
1026 | " |
---|
1027 | { |
---|
1028 | //--- compute the Euler characteristic of the intersection of the |
---|
1029 | //--- exceptional hypersurfaces with \pi^-1(0) which can be of |
---|
1030 | //--- dimension 1 or 2 - not considering the components of the strict |
---|
1031 | //--- transform |
---|
1032 | |
---|
1033 | //---------------------------------------------------------------------------- |
---|
1034 | // Initialization |
---|
1035 | //---------------------------------------------------------------------------- |
---|
1036 | int i,j,k,aa,m,n,thisE,otherE; |
---|
1037 | def R=basering; |
---|
1038 | intvec nulliv,chi_temp,kvec,dimEi,endiv; |
---|
1039 | nulliv[size(iden)]=0; |
---|
1040 | dimEi[size(iden)]=0; |
---|
1041 | endiv[size(re[2])]=0; |
---|
1042 | list chi_E; |
---|
1043 | for(i=1;i<=size(iden);i++) |
---|
1044 | { |
---|
1045 | chi_E[i]=list(); |
---|
1046 | } |
---|
1047 | //--------------------------------------------------------------------------- |
---|
1048 | // Run through the list of charts and compute the Euler characteristic of |
---|
1049 | // the new exceptional hypersurface and change the values for the old ones |
---|
1050 | // according to the blow-up which has just been performed |
---|
1051 | // For initialization reasons, treat the case of the first blow-up separately |
---|
1052 | //--------------------------------------------------------------------------- |
---|
1053 | for(i=2;i<=size(re[2]);i++) |
---|
1054 | { |
---|
1055 | //--- run through all charts |
---|
1056 | if(defined(S)){kill S;} |
---|
1057 | def S=re[2][i]; |
---|
1058 | setring S; |
---|
1059 | if(defined(EList)) |
---|
1060 | { |
---|
1061 | endiv[i]=1; |
---|
1062 | } |
---|
1063 | m=int(leadcoef(path[1,ncols(path)])); |
---|
1064 | if(defined(Spa)){kill Spa;} |
---|
1065 | def Spa=re[2][m]; |
---|
1066 | if(size(BO[4])==1) |
---|
1067 | { |
---|
1068 | //--- just one exceptional divisor |
---|
1069 | thisE=1; |
---|
1070 | setring Spa; |
---|
1071 | if(i==2) |
---|
1072 | { |
---|
1073 | //--- have not set the initial value of chi(E_1) yet |
---|
1074 | //--- in the local case, we need to know whether the center contains 0 |
---|
1075 | if(size(reduce(cent,std(maxideal(1))))!=0) |
---|
1076 | { |
---|
1077 | //--- first center does not meet 0 |
---|
1078 | list templist=0,nulliv; |
---|
1079 | dimEi[1]=-1; |
---|
1080 | } |
---|
1081 | else |
---|
1082 | { |
---|
1083 | if(dim(std(cent))==0) |
---|
1084 | { |
---|
1085 | //--- center was point ==> new except. div. is a P^2 |
---|
1086 | list templist=3*vdim(std(BO[1]+cent)),nulliv; |
---|
1087 | dimEi[1]=2; |
---|
1088 | } |
---|
1089 | else |
---|
1090 | { |
---|
1091 | //--- center was curve ==> intersection of new exceptional divisor |
---|
1092 | //--- with \pi^-1(0) is a curve, namely P^1 |
---|
1093 | setring S; |
---|
1094 | list templist=2,nulliv; |
---|
1095 | dimEi[1]=1; |
---|
1096 | } |
---|
1097 | } |
---|
1098 | chi_E[1]=templist; |
---|
1099 | kill templist; |
---|
1100 | } |
---|
1101 | setring S; |
---|
1102 | i++; |
---|
1103 | if(i<size(re[2])) |
---|
1104 | { |
---|
1105 | continue; |
---|
1106 | } |
---|
1107 | else |
---|
1108 | { |
---|
1109 | break; |
---|
1110 | } |
---|
1111 | } |
---|
1112 | for(j=1;j<=size(iden);j++) |
---|
1113 | { |
---|
1114 | //--- find out which exceptional divisor has just been born |
---|
1115 | if(inIVList(intvec(i,size(BO[4])),iden[j])) |
---|
1116 | { |
---|
1117 | //--- found it |
---|
1118 | thisE=j; |
---|
1119 | break; |
---|
1120 | } |
---|
1121 | } |
---|
1122 | //--- now setup new chi and change the previous ones appropriately |
---|
1123 | setring Spa; |
---|
1124 | if(size(chi_E[thisE])==0) |
---|
1125 | { |
---|
1126 | //--- have not set the initial value of chi(E_thisE) yet |
---|
1127 | if(deg(std(cent+BO[5])[1])==0) |
---|
1128 | { |
---|
1129 | if(dim(std(cent))==0) |
---|
1130 | { |
---|
1131 | //--- \pi^-1(0) does not meet the Q-point cent |
---|
1132 | list templist=0,nulliv; |
---|
1133 | dimEi[thisE]=-1; |
---|
1134 | } |
---|
1135 | //--- if cent is a curve, the intersection point might simply be outside |
---|
1136 | //--- of this chart!!! |
---|
1137 | } |
---|
1138 | else |
---|
1139 | { |
---|
1140 | if(dim(std(cent))==0) |
---|
1141 | { |
---|
1142 | //--- center was point ==> new except. div. is a P^2 |
---|
1143 | list templist=3*vdim(std(BO[1]+cent)),nulliv; |
---|
1144 | dimEi[thisE]=2; |
---|
1145 | } |
---|
1146 | else |
---|
1147 | { |
---|
1148 | //--- center was curve ==> new except. div. is a C x P^1 |
---|
1149 | if(dim(std(cent+BO[5]))==1) |
---|
1150 | { |
---|
1151 | //--- whole curve is in \pi^-1(0) |
---|
1152 | list templist=4-4*genus(BO[1]+cent),nulliv; |
---|
1153 | dimEi[thisE]=2; |
---|
1154 | } |
---|
1155 | else |
---|
1156 | { |
---|
1157 | //--- curve meets \pi^-1(0) in points |
---|
1158 | //--- in S, the intersection will be a curve!!! |
---|
1159 | setring S; |
---|
1160 | list templist=2-2*genus(BO[1]+BO[4][size(BO[4])]+BO[5]),nulliv; |
---|
1161 | dimEi[thisE]=1; |
---|
1162 | setring Spa; |
---|
1163 | } |
---|
1164 | } |
---|
1165 | } |
---|
1166 | if(defined(templist)) |
---|
1167 | { |
---|
1168 | chi_E[thisE]=templist; |
---|
1169 | } |
---|
1170 | kill templist; |
---|
1171 | } |
---|
1172 | for(j=1;j<=size(BO[4]);j++) |
---|
1173 | { |
---|
1174 | //--- we are in the parent ring ==> thisE is not yet born |
---|
1175 | //--- all the other E_i have already been initialized, but the chi |
---|
1176 | //--- might change with the current blow-up at cent |
---|
1177 | if(BO[6][j]==1) |
---|
1178 | { |
---|
1179 | //--- ignore empty sets |
---|
1180 | j++; |
---|
1181 | if(j<=size(BO[4])) |
---|
1182 | { |
---|
1183 | continue; |
---|
1184 | } |
---|
1185 | else |
---|
1186 | { |
---|
1187 | break; |
---|
1188 | } |
---|
1189 | } |
---|
1190 | for(k=1;k<=size(iden);k++) |
---|
1191 | { |
---|
1192 | //--- find global index of BO[4][j] |
---|
1193 | if(inIVList(intvec(m,j),iden[k])) |
---|
1194 | { |
---|
1195 | otherE=k; |
---|
1196 | break; |
---|
1197 | } |
---|
1198 | } |
---|
1199 | if(dimEi[otherE]<=1) |
---|
1200 | { |
---|
1201 | //--- dimEi[otherE]==-1: center leading to this E does not meet \pi^-1(0) |
---|
1202 | //--- dimEi[otherE]== 0: center leading to this E does not meet \pi^-1(0) |
---|
1203 | //--- in any previously visited charts |
---|
1204 | //--- maybe in some other branch later, but has nothing |
---|
1205 | //--- to do with this center |
---|
1206 | //--- dimEi[otherE]== 1: E \cap \pi^-1(0) is curve |
---|
1207 | //--- ==> chi is birational invariant |
---|
1208 | j++; |
---|
1209 | if(j<=size(BO[4])) |
---|
1210 | { |
---|
1211 | continue; |
---|
1212 | } |
---|
1213 | break; |
---|
1214 | } |
---|
1215 | if(chi_E[otherE][2][thisE]==1) |
---|
1216 | { |
---|
1217 | //--- already considered this one |
---|
1218 | j++; |
---|
1219 | if(j<=size(BO[4])) |
---|
1220 | { |
---|
1221 | continue; |
---|
1222 | } |
---|
1223 | else |
---|
1224 | { |
---|
1225 | break; |
---|
1226 | } |
---|
1227 | } |
---|
1228 | //--------------------------------------------------------------------------- |
---|
1229 | // update chi according to the formula |
---|
1230 | // chi(E_k^transf)=chi(E_k) - chi(C \cap E_k) + chi(E_k \cap E_new) |
---|
1231 | // for convenience of implementation, we first compute |
---|
1232 | // chi(E_k) - chi(C \cap E_k) |
---|
1233 | // and afterwards add the last term chi(E_k \cap E_new) |
---|
1234 | //--------------------------------------------------------------------------- |
---|
1235 | ideal CinE=std(cent+BO[4][j]+BO[1]); // this is C \cap E_k |
---|
1236 | if(dim(CinE)==1) |
---|
1237 | { |
---|
1238 | //--- center meets E_k in a curve |
---|
1239 | chi_temp[otherE]=chi_E[otherE][1]-(2-2*genus(CinE)); |
---|
1240 | } |
---|
1241 | if(dim(CinE)==0) |
---|
1242 | { |
---|
1243 | //--- center meets E_k in points |
---|
1244 | chi_temp[otherE]=chi_E[otherE][1]-vdim(std(CinE)); |
---|
1245 | } |
---|
1246 | kill CinE; |
---|
1247 | setring S; |
---|
1248 | //--- now we are back in the i-th ring in the list |
---|
1249 | ideal CinE=std(BO[4][j]+BO[4][size(BO[4])]+BO[1]); |
---|
1250 | // this is E_k \cap E_new |
---|
1251 | if(dim(CinE)==1) |
---|
1252 | { |
---|
1253 | //--- if the two divisors meet, they meet in a curve |
---|
1254 | chi_E[otherE][1]=chi_temp[otherE]+(2-2*genus(CinE)); |
---|
1255 | chi_E[otherE][2][thisE]=1; // this blow-up of E_k is done |
---|
1256 | } |
---|
1257 | kill CinE; |
---|
1258 | setring Spa; |
---|
1259 | } |
---|
1260 | } |
---|
1261 | //--- we still need to clean-up the 1-dimensional E_i \cap \pi^-1(0) |
---|
1262 | for(i=1;i<=size(dimEi);i++) |
---|
1263 | { |
---|
1264 | if(dimEi[i]!=1) |
---|
1265 | { |
---|
1266 | //--- not 1-dimensional ==> skip |
---|
1267 | i++; |
---|
1268 | if(i>size(dimEi)) break; |
---|
1269 | continue; |
---|
1270 | } |
---|
1271 | if(defined(myCharts)) {kill myCharts;} |
---|
1272 | intvec myCharts; |
---|
1273 | chi_E[i]=0; |
---|
1274 | for(j=1;j<=size(re[2]);j++) |
---|
1275 | { |
---|
1276 | if(endiv[j]==0) |
---|
1277 | { |
---|
1278 | //--- not an endChart ==> skip |
---|
1279 | j++; |
---|
1280 | if(j>size(re[2])) break; |
---|
1281 | continue; |
---|
1282 | } |
---|
1283 | if(defined(mtuple)) {kill mtuple;} |
---|
1284 | intvec mtuple=findInIVList(1,j,iden[i]); |
---|
1285 | if(size(mtuple)==1) |
---|
1286 | { |
---|
1287 | //-- nothing to do with this Ei ==> skip |
---|
1288 | j++; |
---|
1289 | if(j>size(re[2])) break; |
---|
1290 | continue; |
---|
1291 | } |
---|
1292 | myCharts[size(myCharts)+1]=j; |
---|
1293 | if(defined(S)){kill S;} |
---|
1294 | def S=re[2][j]; |
---|
1295 | setring S; |
---|
1296 | if(defined(interId)){kill interId;} |
---|
1297 | //--- all components |
---|
1298 | ideal interId=std(BO[4][mtuple[2]]+BO[5]); |
---|
1299 | if(defined(myPts)){kill myPts;} |
---|
1300 | ideal myPts=1; |
---|
1301 | export(myPts); |
---|
1302 | export(interId); |
---|
1303 | if(defined(doneId)){kill doneId;} |
---|
1304 | if(defined(donePts)){kill donePts;} |
---|
1305 | ideal donePts=1; |
---|
1306 | ideal doneId=1; |
---|
1307 | for(k=2;k<size(myCharts);k++) |
---|
1308 | { |
---|
1309 | //--- fetch the components which have already been dealt with via fetchInTree |
---|
1310 | if(defined(opath)) {kill opath;} |
---|
1311 | def opath=imap(re[2][myCharts[k][1]],path); |
---|
1312 | aa=1; |
---|
1313 | while((aa<ncols(opath))&&(aa<ncols(path))) |
---|
1314 | { |
---|
1315 | if(path[1,aa+1]!=opath[1,aa+1]) break; |
---|
1316 | aa++; |
---|
1317 | } |
---|
1318 | if(defined(comPa)) {kill comPa;} |
---|
1319 | int comPa=int(leadcoef(path[1,aa])); |
---|
1320 | if(defined(tempId)){kill tempId;} |
---|
1321 | ideal tempId=fetchInTree(re,myCharts[k][1],comPa,j,"interId",iden); |
---|
1322 | doneId=intersect(doneId,tempId); |
---|
1323 | kill tempId; |
---|
1324 | ideal tempId=fetchInTree(re,myCharts[k][1],comPa,j,"myPts",iden); |
---|
1325 | donePts=intersect(donePts,tempId); |
---|
1326 | kill tempId; |
---|
1327 | } |
---|
1328 | //--- drop components which have already been dealt with |
---|
1329 | interId=sat(interId,doneId)[1]; |
---|
1330 | list pr=minAssGTZ(interId); |
---|
1331 | myPts=std(interId+doneId); |
---|
1332 | for(k=1;k<=size(pr);k++) |
---|
1333 | { |
---|
1334 | for(n=k+1;n<=size(pr);n++) |
---|
1335 | { |
---|
1336 | myPts=intersect(myPts,std(pr[k]+pr[n])); |
---|
1337 | } |
---|
1338 | if(deg(std(pr[k])[1])>0) |
---|
1339 | { |
---|
1340 | chi_E[i]=chi_E[i]+(2-2*genus(pr[k])); |
---|
1341 | } |
---|
1342 | } |
---|
1343 | myPts=sat(myPts,donePts)[1]; |
---|
1344 | chi_E[i]=chi_E[i]-vdim(myPts); |
---|
1345 | } |
---|
1346 | } |
---|
1347 | setring R; |
---|
1348 | return(chi_E); |
---|
1349 | } |
---|
1350 | ////////////////////////////////////////////////////////////////////////////// |
---|
1351 | static |
---|
1352 | proc chi_ast(list re,list iden,list #) |
---|
1353 | "Internal procedure - no help and no example available |
---|
1354 | " |
---|
1355 | { |
---|
1356 | //--- compute the Euler characteristic of the Ei,Eij,Eijk and the |
---|
1357 | //--- corresponding Ei^*,Eij^*,Eijk^* by preparing the input to the |
---|
1358 | //--- specialized auxilliary procedures and then recombining the results |
---|
1359 | |
---|
1360 | //---------------------------------------------------------------------------- |
---|
1361 | // Initialization |
---|
1362 | //---------------------------------------------------------------------------- |
---|
1363 | int i,j,k,g; |
---|
1364 | intvec tiv; |
---|
1365 | list chi_ijk,chi_ij,chi_i,ast_ijk,ast_ij,ast_i,tmplist,g_ij,emptylist; |
---|
1366 | list leererSchnitt; |
---|
1367 | def R=basering; |
---|
1368 | ring Rhelp=0,@t,dp; |
---|
1369 | setring R; |
---|
1370 | //---------------------------------------------------------------------------- |
---|
1371 | // first compute the chi(Eij) and at the same time |
---|
1372 | // check whether E_i \cap E_j is empty |
---|
1373 | // the formula is |
---|
1374 | // chi_ij=2-2*genus(E_i \cap E_j) |
---|
1375 | //---------------------------------------------------------------------------- |
---|
1376 | if(size(#)>0) |
---|
1377 | { |
---|
1378 | "Entering chi_ast"; |
---|
1379 | } |
---|
1380 | for(i=1;i<=size(iden)-1;i++) |
---|
1381 | { |
---|
1382 | for(j=i+1;j<=size(iden);j++) |
---|
1383 | { |
---|
1384 | if(defined(blub)){kill blub;} |
---|
1385 | def blub=chiEij(re,iden,intvec(i,j)); |
---|
1386 | if(typeof(blub)=="int") |
---|
1387 | { |
---|
1388 | tmplist=intvec(i,j),blub; |
---|
1389 | } |
---|
1390 | else |
---|
1391 | { |
---|
1392 | leererSchnitt[size(leererSchnitt)+1]=intvec(i,j); |
---|
1393 | tmplist=intvec(i,j),0; |
---|
1394 | } |
---|
1395 | chi_ij[size(chi_ij)+1]=tmplist; |
---|
1396 | } |
---|
1397 | } |
---|
1398 | if(size(#)>0) |
---|
1399 | { |
---|
1400 | "chi_ij computed"; |
---|
1401 | } |
---|
1402 | //----------------------------------------------------------------------------- |
---|
1403 | // compute chi(Eijk)=chi^*(Eijk) by counting the points in the intersection |
---|
1404 | // chi_ijk=#(E_i \cap E_j \cap E_k) |
---|
1405 | // ast_ijk=chi_ijk |
---|
1406 | //----------------------------------------------------------------------------- |
---|
1407 | for(i=1;i<=size(iden)-2;i++) |
---|
1408 | { |
---|
1409 | for(j=i+1;j<=size(iden)-1;j++) |
---|
1410 | { |
---|
1411 | for(k=j+1;k<=size(iden);k++) |
---|
1412 | { |
---|
1413 | if(inIVList(intvec(i,j),leererSchnitt)) |
---|
1414 | { |
---|
1415 | tmplist=intvec(i,j,k),0; |
---|
1416 | } |
---|
1417 | else |
---|
1418 | { |
---|
1419 | tmplist=intvec(i,j,k),countEijk(re,iden,intvec(i,j,k)); |
---|
1420 | } |
---|
1421 | chi_ijk[size(chi_ijk)+1]=tmplist; |
---|
1422 | } |
---|
1423 | } |
---|
1424 | } |
---|
1425 | ast_ijk=chi_ijk; |
---|
1426 | if(size(#)>0) |
---|
1427 | { |
---|
1428 | "chi_ijk computed"; |
---|
1429 | } |
---|
1430 | //---------------------------------------------------------------------------- |
---|
1431 | // construct chi(Eij^*) by the formula |
---|
1432 | // ast_ij=chi_ij - sum_ijk chi_ijk, |
---|
1433 | // where k runs over all indices != i,j |
---|
1434 | //---------------------------------------------------------------------------- |
---|
1435 | for(i=1;i<=size(chi_ij);i++) |
---|
1436 | { |
---|
1437 | ast_ij[i]=chi_ij[i]; |
---|
1438 | for(k=1;k<=size(chi_ijk);k++) |
---|
1439 | { |
---|
1440 | if(((chi_ijk[k][1][1]==chi_ij[i][1][1])|| |
---|
1441 | (chi_ijk[k][1][2]==chi_ij[i][1][1]))&& |
---|
1442 | ((chi_ijk[k][1][2]==chi_ij[i][1][2])|| |
---|
1443 | (chi_ijk[k][1][3]==chi_ij[i][1][2]))) |
---|
1444 | { |
---|
1445 | ast_ij[i][2]=ast_ij[i][2]-chi_ijk[k][2]; |
---|
1446 | } |
---|
1447 | } |
---|
1448 | } |
---|
1449 | if(size(#)>0) |
---|
1450 | { |
---|
1451 | "ast_ij computed"; |
---|
1452 | } |
---|
1453 | //---------------------------------------------------------------------------- |
---|
1454 | // construct ast_i according to the following formulae |
---|
1455 | // ast_i=0 if E_i is (Q- resp. C-)component of the strict transform |
---|
1456 | // chi_i=3*n if E_i originates from blowing up a Q-point, |
---|
1457 | // which consists of n (different) C-points |
---|
1458 | // chi_i=2-2g(C) if E_i originates from blowing up a (Q-)curve C |
---|
1459 | // (chi_i=n*(2-2g(C_i))=2-2g(C), |
---|
1460 | // where C=\cup C_i, C_i \cap C_j = \emptyset) |
---|
1461 | // if E_i is not a component of the strict transform, then |
---|
1462 | // ast_i=chi_i - sum_{j!=i} ast_ij |
---|
1463 | //---------------------------------------------------------------------------- |
---|
1464 | for(i=1;i<=size(iden);i++) |
---|
1465 | { |
---|
1466 | if(defined(S)) {kill S;} |
---|
1467 | def S=re[2][iden[i][1][1]]; |
---|
1468 | setring S; |
---|
1469 | if(iden[i][1][2]>size(BO[4])) |
---|
1470 | { |
---|
1471 | i--; |
---|
1472 | break; |
---|
1473 | } |
---|
1474 | } |
---|
1475 | list idenE=iden; |
---|
1476 | while(size(idenE)>i) |
---|
1477 | { |
---|
1478 | idenE=delete(idenE,size(idenE)); |
---|
1479 | } |
---|
1480 | list cl=computeChiE(re,idenE); |
---|
1481 | for(i=1;i<=size(idenE);i++) |
---|
1482 | { |
---|
1483 | chi_i[i]=list(intvec(i),cl[i][1]); |
---|
1484 | } |
---|
1485 | if(size(#)>0) |
---|
1486 | { |
---|
1487 | "chi_i computed"; |
---|
1488 | } |
---|
1489 | for(i=1;i<=size(idenE);i++) |
---|
1490 | { |
---|
1491 | ast_i[i]=chi_i[i]; |
---|
1492 | for(j=1;j<=size(ast_ij);j++) |
---|
1493 | { |
---|
1494 | if((ast_ij[j][1][1]==i)||(ast_ij[j][1][2]==i)) |
---|
1495 | { |
---|
1496 | ast_i[i][2]=ast_i[i][2]-chi_ij[j][2]; |
---|
1497 | } |
---|
1498 | } |
---|
1499 | for(j=1;j<=size(ast_ijk);j++) |
---|
1500 | { |
---|
1501 | if((ast_ijk[j][1][1]==i)||(ast_ijk[j][1][2]==i) |
---|
1502 | ||(ast_ijk[j][1][3]==i)) |
---|
1503 | { |
---|
1504 | ast_i[i][2]=ast_i[i][2]+chi_ijk[j][2]; |
---|
1505 | } |
---|
1506 | } |
---|
1507 | } |
---|
1508 | for(i=size(idenE)+1;i<=size(iden);i++) |
---|
1509 | { |
---|
1510 | ast_i[i]=list(intvec(i),0); |
---|
1511 | } |
---|
1512 | //--- results are in ast_i, ast_ij and ast_ijk |
---|
1513 | //--- all are of the form intvec(indices),int(value) |
---|
1514 | list result=ast_i,ast_ij,ast_ijk; |
---|
1515 | return(result); |
---|
1516 | } |
---|
1517 | ////////////////////////////////////////////////////////////////////////////// |
---|
1518 | static |
---|
1519 | proc chi_ast_local(list re,list iden,list #) |
---|
1520 | "Internal procedure - no help and no example available |
---|
1521 | " |
---|
1522 | { |
---|
1523 | //--- compute the Euler characteristic of the Ei,Eij,Eijk and the |
---|
1524 | //--- corresponding Ei^*,Eij^*,Eijk^* by preparing the input to the |
---|
1525 | //--- specialized auxilliary procedures and then recombining the results |
---|
1526 | |
---|
1527 | //---------------------------------------------------------------------------- |
---|
1528 | // Initialization |
---|
1529 | //---------------------------------------------------------------------------- |
---|
1530 | int i,j,k,g; |
---|
1531 | intvec tiv; |
---|
1532 | list chi_ijk,chi_ij,chi_i,ast_ijk,ast_ij,ast_i,tmplist,g_ij,emptylist; |
---|
1533 | list leererSchnitt; |
---|
1534 | def R=basering; |
---|
1535 | ring Rhelp=0,@t,dp; |
---|
1536 | setring R; |
---|
1537 | //---------------------------------------------------------------------------- |
---|
1538 | // first compute |
---|
1539 | // if E_i \cap E_j \cap \pi^-1(0) is a curve: |
---|
1540 | // chi(Eij) and at the same time |
---|
1541 | // check whether E_i \cap E_j is empty |
---|
1542 | // the formula is |
---|
1543 | // chi_ij=2-2*genus(E_i \cap E_j) |
---|
1544 | // otherwise (points): |
---|
1545 | // chi(E_ij) by counting the points |
---|
1546 | //---------------------------------------------------------------------------- |
---|
1547 | if(size(#)>0) |
---|
1548 | { |
---|
1549 | "Entering chi_ast_local"; |
---|
1550 | } |
---|
1551 | for(i=1;i<=size(iden)-1;i++) |
---|
1552 | { |
---|
1553 | for(j=i+1;j<=size(iden);j++) |
---|
1554 | { |
---|
1555 | if(defined(blub)){kill blub;} |
---|
1556 | def blub=chiEij_local(re,iden,intvec(i,j)); |
---|
1557 | if(typeof(blub)=="int") |
---|
1558 | { |
---|
1559 | tmplist=intvec(i,j),blub; |
---|
1560 | } |
---|
1561 | else |
---|
1562 | { |
---|
1563 | leererSchnitt[size(leererSchnitt)+1]=intvec(i,j); |
---|
1564 | tmplist=intvec(i,j),0; |
---|
1565 | } |
---|
1566 | chi_ij[size(chi_ij)+1]=tmplist; |
---|
1567 | } |
---|
1568 | } |
---|
1569 | if(size(#)>0) |
---|
1570 | { |
---|
1571 | "chi_ij computed"; |
---|
1572 | } |
---|
1573 | //----------------------------------------------------------------------------- |
---|
1574 | // compute chi(Eijk)=chi^*(Eijk) by counting the points in the intersection |
---|
1575 | // chi_ijk=#(E_i \cap E_j \cap E_k \cap \pi^-1(0)) |
---|
1576 | // ast_ijk=chi_ijk |
---|
1577 | //----------------------------------------------------------------------------- |
---|
1578 | for(i=1;i<=size(iden)-2;i++) |
---|
1579 | { |
---|
1580 | for(j=i+1;j<=size(iden)-1;j++) |
---|
1581 | { |
---|
1582 | for(k=j+1;k<=size(iden);k++) |
---|
1583 | { |
---|
1584 | if(inIVList(intvec(i,j),leererSchnitt)) |
---|
1585 | { |
---|
1586 | tmplist=intvec(i,j,k),0; |
---|
1587 | } |
---|
1588 | else |
---|
1589 | { |
---|
1590 | tmplist=intvec(i,j,k),countEijk(re,iden,intvec(i,j,k),"local"); |
---|
1591 | } |
---|
1592 | chi_ijk[size(chi_ijk)+1]=tmplist; |
---|
1593 | } |
---|
1594 | } |
---|
1595 | } |
---|
1596 | ast_ijk=chi_ijk; |
---|
1597 | if(size(#)>0) |
---|
1598 | { |
---|
1599 | "chi_ijk computed"; |
---|
1600 | } |
---|
1601 | //---------------------------------------------------------------------------- |
---|
1602 | // construct chi(Eij^*) by the formula |
---|
1603 | // ast_ij=chi_ij - sum_ijk chi_ijk, |
---|
1604 | // where k runs over all indices != i,j |
---|
1605 | //---------------------------------------------------------------------------- |
---|
1606 | for(i=1;i<=size(chi_ij);i++) |
---|
1607 | { |
---|
1608 | ast_ij[i]=chi_ij[i]; |
---|
1609 | for(k=1;k<=size(chi_ijk);k++) |
---|
1610 | { |
---|
1611 | if(((chi_ijk[k][1][1]==chi_ij[i][1][1])|| |
---|
1612 | (chi_ijk[k][1][2]==chi_ij[i][1][1]))&& |
---|
1613 | ((chi_ijk[k][1][2]==chi_ij[i][1][2])|| |
---|
1614 | (chi_ijk[k][1][3]==chi_ij[i][1][2]))) |
---|
1615 | { |
---|
1616 | ast_ij[i][2]=ast_ij[i][2]-chi_ijk[k][2]; |
---|
1617 | } |
---|
1618 | } |
---|
1619 | } |
---|
1620 | if(size(#)>0) |
---|
1621 | { |
---|
1622 | "ast_ij computed"; |
---|
1623 | } |
---|
1624 | //---------------------------------------------------------------------------- |
---|
1625 | // construct ast_i according to the following formulae |
---|
1626 | // ast_i=0 if E_i is (Q- resp. C-)component of the strict transform |
---|
1627 | // if E_i \cap \pi^-1(0) is of dimension 2: |
---|
1628 | // chi_i=3*n if E_i originates from blowing up a Q-point, |
---|
1629 | // which consists of n (different) C-points |
---|
1630 | // chi_i=2-2g(C) if E_i originates from blowing up a (Q-)curve C |
---|
1631 | // (chi_i=n*(2-2g(C_i))=2-2g(C), |
---|
1632 | // where C=\cup C_i, C_i \cap C_j = \emptyset) |
---|
1633 | // if E_i \cap \pi^-1(0) is a curve: |
---|
1634 | // use the formula chi_i=2-2*genus(E_i \cap \pi^-1(0)) |
---|
1635 | // |
---|
1636 | // for E_i not a component of the strict transform we have |
---|
1637 | // ast_i=chi_i - sum_{j!=i} ast_ij |
---|
1638 | //---------------------------------------------------------------------------- |
---|
1639 | for(i=1;i<=size(iden);i++) |
---|
1640 | { |
---|
1641 | if(defined(S)) {kill S;} |
---|
1642 | def S=re[2][iden[i][1][1]]; |
---|
1643 | setring S; |
---|
1644 | if(iden[i][1][2]>size(BO[4])) |
---|
1645 | { |
---|
1646 | i--; |
---|
1647 | break; |
---|
1648 | } |
---|
1649 | } |
---|
1650 | list idenE=iden; |
---|
1651 | while(size(idenE)>i) |
---|
1652 | { |
---|
1653 | idenE=delete(idenE,size(idenE)); |
---|
1654 | } |
---|
1655 | list cl=computeChiE_local(re,idenE); |
---|
1656 | for(i=1;i<=size(idenE);i++) |
---|
1657 | { |
---|
1658 | chi_i[i]=list(intvec(i),cl[i][1]); |
---|
1659 | } |
---|
1660 | if(size(#)>0) |
---|
1661 | { |
---|
1662 | "chi_i computed"; |
---|
1663 | } |
---|
1664 | for(i=1;i<=size(idenE);i++) |
---|
1665 | { |
---|
1666 | ast_i[i]=chi_i[i]; |
---|
1667 | for(j=1;j<=size(ast_ij);j++) |
---|
1668 | { |
---|
1669 | if((ast_ij[j][1][1]==i)||(ast_ij[j][1][2]==i)) |
---|
1670 | { |
---|
1671 | ast_i[i][2]=ast_i[i][2]-chi_ij[j][2]; |
---|
1672 | } |
---|
1673 | } |
---|
1674 | for(j=1;j<=size(ast_ijk);j++) |
---|
1675 | { |
---|
1676 | if((ast_ijk[j][1][1]==i)||(ast_ijk[j][1][2]==i) |
---|
1677 | ||(ast_ijk[j][1][3]==i)) |
---|
1678 | { |
---|
1679 | ast_i[i][2]=ast_i[i][2]+chi_ijk[j][2]; |
---|
1680 | } |
---|
1681 | } |
---|
1682 | } |
---|
1683 | for(i=size(idenE)+1;i<=size(iden);i++) |
---|
1684 | { |
---|
1685 | ast_i[i]=list(intvec(i),0); |
---|
1686 | } |
---|
1687 | //--- results are in ast_i, ast_ij and ast_ijk |
---|
1688 | //--- all are of the form intvec(indices),int(value) |
---|
1689 | //"End of chi_ast_local"; |
---|
1690 | //~; |
---|
1691 | list result=ast_i,ast_ij,ast_ijk; |
---|
1692 | return(result); |
---|
1693 | } |
---|
1694 | ////////////////////////////////////////////////////////////////////////////// |
---|
1695 | |
---|
1696 | proc zetaDL(list re,int d,list #) |
---|
1697 | "USAGE: zetaDL(L,d[,s][,a]); |
---|
1698 | @* L = list of rings |
---|
1699 | @* d = integer |
---|
1700 | @* s = string |
---|
1701 | @* a = integer |
---|
1702 | ASSUME: L is the output of resolution of singularities |
---|
1703 | COMPUTE: local Denef-Loeser zeta function, if string s is present and |
---|
1704 | has the value 'lokal'; global Denef-Loeser zeta function otherwise |
---|
1705 | RETURN: string specifying the Denef-Loeser zeta function if a not present |
---|
1706 | list l if a is present, where |
---|
1707 | l[1]: string specifying the Denef-Loeser zeta function |
---|
1708 | l[2]: list ast, |
---|
1709 | ast[1]=chi(Ei^*) |
---|
1710 | ast[2]=chi(Eij^*) |
---|
1711 | ast[3]=chi(Eijk^*) |
---|
1712 | l[3]: intvec nu of multiplicites as needed in computation of zeta |
---|
1713 | function |
---|
1714 | l[4]: intvec N of multiplicities as needed in compuation of zeta |
---|
1715 | function |
---|
1716 | EXAMPLE: example zetaDL; shows an example |
---|
1717 | " |
---|
1718 | { |
---|
1719 | //---------------------------------------------------------------------------- |
---|
1720 | // Initialization |
---|
1721 | //---------------------------------------------------------------------------- |
---|
1722 | def R=basering; |
---|
1723 | int show_all,i; |
---|
1724 | if(size(#)>0) |
---|
1725 | { |
---|
1726 | if((typeof(#[1])=="int")||(size(#)>1)) |
---|
1727 | { |
---|
1728 | show_all=1; |
---|
1729 | } |
---|
1730 | if(typeof(#[1])=="string") |
---|
1731 | { |
---|
1732 | if((#[1]=="local")||(#[1]=="lokal")) |
---|
1733 | { |
---|
1734 | // ERROR("Local case not implemented yet"); |
---|
1735 | "Local Case: Assuming that no (!) charts were dropped"; |
---|
1736 | "during calculation of the resolution (option \"A\")"; |
---|
1737 | int localComp=1; |
---|
1738 | } |
---|
1739 | else |
---|
1740 | { |
---|
1741 | "Computing global zeta function"; |
---|
1742 | } |
---|
1743 | } |
---|
1744 | } |
---|
1745 | //---------------------------------------------------------------------------- |
---|
1746 | // Identify the embedded divisors and chi(Ei^*), chi(Eij^*) and chi(Eijk^*) |
---|
1747 | // as well as the integer vector V(=nu) and N |
---|
1748 | //---------------------------------------------------------------------------- |
---|
1749 | list iden=prepEmbDiv(re); //--- identify the E_i |
---|
1750 | //!!! TIMING: E8 takes 520 sec ==> needs speed up |
---|
1751 | if(!defined(localComp)) |
---|
1752 | { |
---|
1753 | list ast_list=chi_ast(re,iden); //--- compute chi(E^*) |
---|
1754 | } |
---|
1755 | else |
---|
1756 | { |
---|
1757 | list ast_list=chi_ast_local(re,iden); |
---|
1758 | } |
---|
1759 | intvec Vvec=computeV(re,iden); //--- nu |
---|
1760 | intvec Nvec=computeN(re,iden); //--- N |
---|
1761 | //---------------------------------------------------------------------------- |
---|
1762 | // Build a new ring with one parameter s |
---|
1763 | // and compute Zeta_top^(d) in its ground field |
---|
1764 | //---------------------------------------------------------------------------- |
---|
1765 | ring Qs=(0,s),x,dp; |
---|
1766 | number zetaTop=0; |
---|
1767 | number enum,denom; |
---|
1768 | denom=1; |
---|
1769 | for(i=1;i<=size(Nvec);i++) |
---|
1770 | { |
---|
1771 | denom=denom*(Vvec[i]+s*Nvec[i]); |
---|
1772 | } |
---|
1773 | //--- factors for which index set J consists of one element |
---|
1774 | //--- (do something only if d divides N_j) |
---|
1775 | for(i=1;i<=size(ast_list[1]);i++) |
---|
1776 | { |
---|
1777 | if((((Nvec[ast_list[1][i][1][1]]/d)*d)-Nvec[ast_list[1][i][1][1]]==0)&& |
---|
1778 | (ast_list[1][i][2]!=0)) |
---|
1779 | { |
---|
1780 | enum=enum+ast_list[1][i][2]*(denom/(Vvec[ast_list[1][i][1][1]]+s*Nvec[ast_list[1][i][1][1]])); |
---|
1781 | } |
---|
1782 | } |
---|
1783 | //--- factors for which index set J consists of two elements |
---|
1784 | //--- (do something only if d divides both N_i and N_j) |
---|
1785 | //!!! TIMING: E8 takes 690 sec and has 703 elements |
---|
1786 | //!!! ==> need to implement a smarter way to do this |
---|
1787 | //!!! e.g. build up enumerator and denominator separately, thus not |
---|
1788 | //!!! searching for common factors in each step |
---|
1789 | for(i=1;i<=size(ast_list[2]);i++) |
---|
1790 | { |
---|
1791 | if((((Nvec[ast_list[2][i][1][1]]/d)*d)-Nvec[ast_list[2][i][1][1]]==0)&& |
---|
1792 | (((Nvec[ast_list[2][i][1][2]]/d)*d)-Nvec[ast_list[2][i][1][2]]==0)&& |
---|
1793 | (ast_list[2][i][2]!=0)) |
---|
1794 | { |
---|
1795 | enum=enum+ast_list[2][i][2]*(denom/((Vvec[ast_list[2][i][1][1]]+s*Nvec[ast_list[2][i][1][1]])*(Vvec[ast_list[2][i][1][2]]+s*Nvec[ast_list[2][i][1][2]]))); |
---|
1796 | } |
---|
1797 | } |
---|
1798 | //--- factors for which index set J consists of three elements |
---|
1799 | //--- (do something only if d divides N_i, N_j and N_k) |
---|
1800 | //!!! TIMING: E8 takes 490 sec and has 8436 elements |
---|
1801 | //!!! ==> same kind of improvements as in the previous case needed |
---|
1802 | for(i=1;i<=size(ast_list[3]);i++) |
---|
1803 | { |
---|
1804 | if((((Nvec[ast_list[3][i][1][1]]/d)*d)-Nvec[ast_list[3][i][1][1]]==0)&& |
---|
1805 | (((Nvec[ast_list[3][i][1][2]]/d)*d)-Nvec[ast_list[3][i][1][2]]==0)&& |
---|
1806 | (((Nvec[ast_list[3][i][1][3]]/d)*d)-Nvec[ast_list[3][i][1][3]]==0)&& |
---|
1807 | (ast_list[3][i][2]!=0)) |
---|
1808 | { |
---|
1809 | enum=enum+ast_list[3][i][2]*(denom/((Vvec[ast_list[3][i][1][1]]+s*Nvec[ast_list[3][i][1][1]])*(Vvec[ast_list[3][i][1][2]]+s*Nvec[ast_list[3][i][1][2]])*(Vvec[ast_list[3][i][1][3]]+s*Nvec[ast_list[3][i][1][3]]))); |
---|
1810 | } |
---|
1811 | } |
---|
1812 | zetaTop=enum/denom; |
---|
1813 | zetaTop=numerator(zetaTop)/denominator(zetaTop); |
---|
1814 | string zetaStr=string(zetaTop); |
---|
1815 | |
---|
1816 | if(show_all) |
---|
1817 | { |
---|
1818 | list result=zetaStr,ast_list,Vvec,Nvec; |
---|
1819 | } |
---|
1820 | else |
---|
1821 | { |
---|
1822 | list result=zetaStr; |
---|
1823 | } |
---|
1824 | setring R; |
---|
1825 | return(result); |
---|
1826 | } |
---|
1827 | example |
---|
1828 | {"EXAMPLE:"; |
---|
1829 | echo = 2; |
---|
1830 | ring R=0,(x,y,z),dp; |
---|
1831 | ideal I=x2+y2+z3; |
---|
1832 | list re=resolve(I,"K"); |
---|
1833 | zetaDL(re,1); |
---|
1834 | I=(xz+y2)*(xz+y2+x2)+z5; |
---|
1835 | list L=resolve(I,"K"); |
---|
1836 | zetaDL(L,1); |
---|
1837 | |
---|
1838 | //===== expected zeta function ========= |
---|
1839 | // (130s+20s^2+87)/(1+s)/(3+4s)/(29+40s) |
---|
1840 | //====================================== |
---|
1841 | } |
---|
1842 | ////////////////////////////////////////////////////////////////////////////// |
---|
1843 | |
---|
1844 | proc abstractR(list re) |
---|
1845 | "USAGE: abstractR(L); |
---|
1846 | @* L = list of rings |
---|
1847 | ASSUME: L is output of resolution of singularities |
---|
1848 | NOTE: currently only implemented for isolated surface singularities |
---|
1849 | RETURN: list l |
---|
1850 | l[1]: intvec, where |
---|
1851 | l[1][i]=1 if the corresponding ring is a final chart |
---|
1852 | of non-embedded resolution |
---|
1853 | l[1][i]=0 otherwise |
---|
1854 | l[2]: intvec, where |
---|
1855 | l[2][i]=1 if the corresponding ring does not occur |
---|
1856 | in the non-embedded resolution |
---|
1857 | l[2][i]=0 otherwise |
---|
1858 | l[3]: list L |
---|
1859 | EXAMPLE: example abstractR; shows an example |
---|
1860 | " |
---|
1861 | { |
---|
1862 | //--------------------------------------------------------------------------- |
---|
1863 | // Initialization and sanity checks |
---|
1864 | //--------------------------------------------------------------------------- |
---|
1865 | def R=basering; |
---|
1866 | //---Test whether we are in the irreducible surface case |
---|
1867 | def S=re[2][1]; |
---|
1868 | setring S; |
---|
1869 | BO[2]=BO[2]+BO[1]; |
---|
1870 | if(dim(std(BO[2]))!=2) |
---|
1871 | { |
---|
1872 | ERROR("NOT A SURFACE"); |
---|
1873 | } |
---|
1874 | if(dim(std(slocus(BO[2])))>0) |
---|
1875 | { |
---|
1876 | ERROR("NOT AN ISOLATED SINGULARITY"); |
---|
1877 | } |
---|
1878 | setring R; |
---|
1879 | int i,j,k,l,i0; |
---|
1880 | intvec deleted; |
---|
1881 | intvec endiv; |
---|
1882 | endiv[size(re[2])]=0; |
---|
1883 | deleted[size(re[2])]=0; |
---|
1884 | //----------------------------------------------------------------------------- |
---|
1885 | // run through all rings, only consider final charts |
---|
1886 | // for each final chart follow the list of charts leading up to it until |
---|
1887 | // we encounter a chart which is not finished in the non-embedded case |
---|
1888 | //----------------------------------------------------------------------------- |
---|
1889 | for(i=1;i<=size(re[2]);i++) |
---|
1890 | { |
---|
1891 | if(defined(S)){kill S;} |
---|
1892 | def S=re[2][i]; |
---|
1893 | setring S; |
---|
1894 | if(size(reduce(cent,std(BO[2]+BO[1])))!=0) |
---|
1895 | { |
---|
1896 | //--- only consider endrings |
---|
1897 | i++; |
---|
1898 | continue; |
---|
1899 | } |
---|
1900 | i0=i; |
---|
1901 | for(j=ncols(path);j>=2;j--) |
---|
1902 | { |
---|
1903 | //--- walk backwards through history |
---|
1904 | if(j==2) |
---|
1905 | { |
---|
1906 | endiv[i0]=1; |
---|
1907 | break; |
---|
1908 | } |
---|
1909 | k=int(leadcoef(path[1,j])); |
---|
1910 | if((deleted[k]==1)||(endiv[k]==1)) |
---|
1911 | { |
---|
1912 | deleted[i0]=1; |
---|
1913 | break; |
---|
1914 | } |
---|
1915 | if(defined(SPa)){kill SPa;} |
---|
1916 | def SPa=re[2][k]; |
---|
1917 | setring SPa; |
---|
1918 | l=int(leadcoef(path[1,ncols(path)])); |
---|
1919 | if(defined(SPa2)){kill SPa2;} |
---|
1920 | def SPa2=re[2][l]; |
---|
1921 | setring SPa2; |
---|
1922 | if((deleted[l]==1)||(endiv[l]==1)) |
---|
1923 | { |
---|
1924 | //--- parent was already treated via different final chart |
---|
1925 | //--- we may safely inherit the data |
---|
1926 | deleted[i0]=1; |
---|
1927 | setring S; |
---|
1928 | i0=k; |
---|
1929 | j--; |
---|
1930 | continue; |
---|
1931 | } |
---|
1932 | setring SPa; |
---|
1933 | //!!! Idea of Improvement: |
---|
1934 | //!!! BESSER: rueckwaerts gehend nur testen ob glatt |
---|
1935 | //!!! danach vorwaerts bis zum ersten Mal abstractNC |
---|
1936 | //!!! ACHTUNG: rueckweg unterwegs notieren - wir haben nur vergangenheit! |
---|
1937 | if((deg(std(slocus(BO[2]))[1])!=0)||(!abstractNC(BO))) |
---|
1938 | { |
---|
1939 | //--- not finished in the non-embedded case |
---|
1940 | endiv[i0]=1; |
---|
1941 | break; |
---|
1942 | } |
---|
1943 | //--- unnecessary chart in non-embedded case |
---|
1944 | setring S; |
---|
1945 | deleted[i0]=1; |
---|
1946 | i0=k; |
---|
1947 | } |
---|
1948 | } |
---|
1949 | //----------------------------------------------------------------------------- |
---|
1950 | // Clean up the intvec deleted and return the result |
---|
1951 | //----------------------------------------------------------------------------- |
---|
1952 | setring R; |
---|
1953 | for(i=1;i<=size(endiv);i++) |
---|
1954 | { |
---|
1955 | if(endiv[i]==1) |
---|
1956 | { |
---|
1957 | if(defined(S)) {kill S;} |
---|
1958 | def S=re[2][i]; |
---|
1959 | setring S; |
---|
1960 | for(j=3;j<ncols(path);j++) |
---|
1961 | { |
---|
1962 | if((endiv[int(leadcoef(path[1,j]))]==1)|| |
---|
1963 | (deleted[int(leadcoef(path[1,j]))]==1)) |
---|
1964 | { |
---|
1965 | deleted[int(leadcoef(path[1,j+1]))]=1; |
---|
1966 | endiv[int(leadcoef(path[1,j+1]))]=0; |
---|
1967 | } |
---|
1968 | } |
---|
1969 | if((endiv[int(leadcoef(path[1,ncols(path)]))]==1)|| |
---|
1970 | (deleted[int(leadcoef(path[1,ncols(path)]))]==1)) |
---|
1971 | { |
---|
1972 | deleted[i]=1; |
---|
1973 | endiv[i]=0; |
---|
1974 | } |
---|
1975 | } |
---|
1976 | } |
---|
1977 | list resu=endiv,deleted,re; |
---|
1978 | return(resu); |
---|
1979 | } |
---|
1980 | example |
---|
1981 | {"EXAMPLE:"; |
---|
1982 | echo = 2; |
---|
1983 | ring r = 0,(x,y,z),dp; |
---|
1984 | ideal I=x2+y2+z11; |
---|
1985 | list L=resolve(I); |
---|
1986 | abstractR(L); |
---|
1987 | } |
---|
1988 | ////////////////////////////////////////////////////////////////////////////// |
---|
1989 | static |
---|
1990 | proc decompE(list BO) |
---|
1991 | "Internal procedure - no help and no example available |
---|
1992 | " |
---|
1993 | { |
---|
1994 | //--- compute the list of exceptional divisors, including the components |
---|
1995 | //--- of the strict transform in the non-embedded case |
---|
1996 | //--- (computation over Q !!!) |
---|
1997 | def R=basering; |
---|
1998 | list Elist,prList; |
---|
1999 | int i; |
---|
2000 | for(i=1;i<=size(BO[4]);i++) |
---|
2001 | { |
---|
2002 | Elist[i]=BO[4][i]; |
---|
2003 | } |
---|
2004 | /* practical speed up (part 1 of 3) -- no theoretical relevance |
---|
2005 | ideal M=maxideal(1); |
---|
2006 | M[1]=var(nvars(basering)); |
---|
2007 | M[nvars(basering)]=var(1); |
---|
2008 | map phi=R,M; |
---|
2009 | */ |
---|
2010 | ideal KK=BO[2]; |
---|
2011 | |
---|
2012 | /* practical speed up (part 2 of 3) |
---|
2013 | KK=phi(KK); |
---|
2014 | */ |
---|
2015 | prList=minAssGTZ(KK); |
---|
2016 | |
---|
2017 | /* practical speed up (part 3 of 3) |
---|
2018 | prList=phi(prList); |
---|
2019 | */ |
---|
2020 | |
---|
2021 | for(i=1;i<=size(prList);i++) |
---|
2022 | { |
---|
2023 | Elist[size(BO[4])+i]=prList[i]; |
---|
2024 | } |
---|
2025 | return(Elist); |
---|
2026 | } |
---|
2027 | ////////////////////////////////////////////////////////////////////////////// |
---|
2028 | |
---|
2029 | proc prepEmbDiv(list re, list #) |
---|
2030 | "USAGE: prepEmbDiv(L[,a]); |
---|
2031 | @* L = list of rings |
---|
2032 | @* a = integer |
---|
2033 | ASSUME: L is output of resolution of singularities |
---|
2034 | COMPUTE: if a is not present: exceptional divisors including components |
---|
2035 | of the strict transform |
---|
2036 | otherwise: only exceptional divisors |
---|
2037 | RETURN: list of Q-irreducible exceptional divisors (embedded case) |
---|
2038 | EXAMPLE: example prepEmbDiv; shows an example |
---|
2039 | " |
---|
2040 | { |
---|
2041 | //--- 1) in each final chart, a list of (decomposed) exceptional divisors |
---|
2042 | //--- is created (and exported) |
---|
2043 | //--- 2) the strict transform is decomposed |
---|
2044 | //--- 3) the exceptional divisors (including the strict transform) |
---|
2045 | //--- in the different charts are compared, identified and this |
---|
2046 | //--- information collected into a list which is then returned |
---|
2047 | //--------------------------------------------------------------------------- |
---|
2048 | // Initialization |
---|
2049 | //--------------------------------------------------------------------------- |
---|
2050 | int i,j,k,ncomps,offset,found,a,b,c,d; |
---|
2051 | list tmpList; |
---|
2052 | def R=basering; |
---|
2053 | //--- identify identical exceptional divisors |
---|
2054 | //--- (note: we are in the embedded case) |
---|
2055 | list iden=collectDiv(re)[2]; |
---|
2056 | //--------------------------------------------------------------------------- |
---|
2057 | // Go to each final chart and create the EList |
---|
2058 | //--------------------------------------------------------------------------- |
---|
2059 | for(i=1;i<=size(iden[size(iden)]);i++) |
---|
2060 | { |
---|
2061 | if(defined(S)){kill S;} |
---|
2062 | def S=re[2][iden[size(iden)][i][1]]; |
---|
2063 | setring S; |
---|
2064 | if(defined(EList)){kill EList;} |
---|
2065 | list EList=decompE(BO); |
---|
2066 | export(EList); |
---|
2067 | setring R; |
---|
2068 | kill S; |
---|
2069 | } |
---|
2070 | //--- save original iden for further use and then drop |
---|
2071 | //--- strict transform from it |
---|
2072 | list iden0=iden; |
---|
2073 | iden=delete(iden,size(iden)); |
---|
2074 | if(size(#)>0) |
---|
2075 | { |
---|
2076 | //--- we are not interested in the strict transform of X |
---|
2077 | return(iden); |
---|
2078 | } |
---|
2079 | //---------------------------------------------------------------------------- |
---|
2080 | // Run through all final charts and collect and identify all components of |
---|
2081 | // the strict transform |
---|
2082 | //---------------------------------------------------------------------------- |
---|
2083 | //--- first final chart - to be used for initialization |
---|
2084 | def S=re[2][iden0[size(iden0)][1][1]]; |
---|
2085 | setring S; |
---|
2086 | ncomps=size(EList)-size(BO[4]); |
---|
2087 | if((ncomps==1)&&(deg(std(EList[size(EList)])[1])==0)) |
---|
2088 | { |
---|
2089 | ncomps=0; |
---|
2090 | } |
---|
2091 | offset=size(BO[4]); |
---|
2092 | for(i=1;i<=ncomps;i++) |
---|
2093 | { |
---|
2094 | //--- add components of strict transform |
---|
2095 | tmpList[1]=intvec(iden0[size(iden0)][1][1],size(BO[4])+i); |
---|
2096 | iden[size(iden)+1]=tmpList; |
---|
2097 | } |
---|
2098 | //--- now run through the other final charts |
---|
2099 | for(i=2;i<=size(iden0[size(iden0)]);i++) |
---|
2100 | { |
---|
2101 | if(defined(S2)){kill S2;} |
---|
2102 | def S2=re[2][iden0[size(iden0)][i][1]]; |
---|
2103 | setring S2; |
---|
2104 | //--- determine common parent of this ring and re[2][iden0[size(iden0)][1][1]] |
---|
2105 | if(defined(opath)){kill opath;} |
---|
2106 | def opath=imap(S,path); |
---|
2107 | j=1; |
---|
2108 | while(opath[1,j]==path[1,j]) |
---|
2109 | { |
---|
2110 | j++; |
---|
2111 | if((j>ncols(path))||(j>ncols(opath))) break; |
---|
2112 | } |
---|
2113 | if(defined(li1)){kill li1;} |
---|
2114 | list li1; |
---|
2115 | //--- fetch the components we have considered in |
---|
2116 | //--- re[2][iden0[size(iden0)][1][1]] |
---|
2117 | //--- via the resolution tree |
---|
2118 | for(k=1;k<=ncomps;k++) |
---|
2119 | { |
---|
2120 | if(defined(id1)){kill id1;} |
---|
2121 | string tempstr="EList["+string(eval(k+offset))+"]"; |
---|
2122 | ideal id1=fetchInTree(re,iden0[size(iden0)][1][1], |
---|
2123 | int(leadcoef(path[1,j-1])), |
---|
2124 | iden0[size(iden0)][i][1],tempstr,iden0,1); |
---|
2125 | kill tempstr; |
---|
2126 | li1[k]=id1; |
---|
2127 | kill id1; |
---|
2128 | } |
---|
2129 | //--- do the comparison |
---|
2130 | for(k=size(BO[4])+1;k<=size(EList);k++) |
---|
2131 | { |
---|
2132 | //--- only components of the strict transform are interesting |
---|
2133 | if((size(BO[4])+1==size(EList))&&(deg(std(EList[size(EList)])[1])==0)) |
---|
2134 | { |
---|
2135 | break; |
---|
2136 | } |
---|
2137 | found=0; |
---|
2138 | for(j=1;j<=size(li1);j++) |
---|
2139 | { |
---|
2140 | if((size(reduce(li1[j],std(EList[k])))==0)&& |
---|
2141 | (size(reduce(EList[k],std(li1[j])))==0)) |
---|
2142 | { |
---|
2143 | //--- found a match |
---|
2144 | li1[j]=ideal(1); |
---|
2145 | iden[size(iden0)-1+j][size(iden[size(iden0)-1+j])+1]= |
---|
2146 | intvec(iden0[size(iden0)][i][1],k); |
---|
2147 | found=1; |
---|
2148 | break; |
---|
2149 | } |
---|
2150 | } |
---|
2151 | if(!found) |
---|
2152 | { |
---|
2153 | //--- no match yet, maybe there are entries not corresponding to the |
---|
2154 | //--- initialization of the list -- collected in list repair |
---|
2155 | if(!defined(repair)) |
---|
2156 | { |
---|
2157 | //--- no entries in repair, we add the very first one |
---|
2158 | list repair; |
---|
2159 | repair[1]=list(intvec(iden0[size(iden0)][i][1],k)); |
---|
2160 | } |
---|
2161 | else |
---|
2162 | { |
---|
2163 | //--- compare against repair, and add the item appropriately |
---|
2164 | //--- steps of comparison as before |
---|
2165 | for(c=1;c<=size(repair);c++) |
---|
2166 | { |
---|
2167 | for(d=1;d<=size(repair[c]);d++) |
---|
2168 | { |
---|
2169 | if(defined(opath)) {kill opath;} |
---|
2170 | def opath=imap(re[2][repair[c][d][1]],path); |
---|
2171 | b=0; |
---|
2172 | while(path[1,b+1]==opath[1,b+1]) |
---|
2173 | { |
---|
2174 | b++; |
---|
2175 | if((b>ncols(path)-1)||(b>ncols(opath)-1)) break; |
---|
2176 | } |
---|
2177 | b=int(leadcoef(path[1,b])); |
---|
2178 | string tempstr="EList["+string(eval(repair[c][d][2])) |
---|
2179 | +"]"; |
---|
2180 | if(defined(id1)){kill id1;} |
---|
2181 | ideal id1=fetchInTree(re,repair[c][d][1],b, |
---|
2182 | iden0[size(iden0)][i][1],tempstr,iden0,1); |
---|
2183 | kill tempstr; |
---|
2184 | if((size(reduce(EList[k],std(id1)))==0)&& |
---|
2185 | (size(reduce(id1,std(EList[k])))==0)) |
---|
2186 | { |
---|
2187 | repair[c][size(repair[c])+1]=intvec(iden0[size(iden0)][i][1],k); |
---|
2188 | break; |
---|
2189 | } |
---|
2190 | } |
---|
2191 | if(d<=size(repair[c])) |
---|
2192 | { |
---|
2193 | break; |
---|
2194 | } |
---|
2195 | } |
---|
2196 | if(c>size(repair)) |
---|
2197 | { |
---|
2198 | repair[size(repair)+1]=list(intvec(iden0[size(iden0)][i][1],k)); |
---|
2199 | } |
---|
2200 | } |
---|
2201 | } |
---|
2202 | } |
---|
2203 | } |
---|
2204 | if(defined(repair)) |
---|
2205 | { |
---|
2206 | //--- there were further components, add them |
---|
2207 | for(c=1;c<=size(repair);c++) |
---|
2208 | { |
---|
2209 | iden[size(iden)+1]=repair[c]; |
---|
2210 | } |
---|
2211 | kill repair; |
---|
2212 | } |
---|
2213 | //--- up to now only Q-irred components - not C-irred components !!! |
---|
2214 | return(iden); |
---|
2215 | } |
---|
2216 | example |
---|
2217 | {"EXAMPLE:"; |
---|
2218 | echo = 2; |
---|
2219 | ring R=0,(x,y,z),dp; |
---|
2220 | ideal I=x2+y2+z11; |
---|
2221 | list L=resolve(I); |
---|
2222 | prepEmbDiv(L); |
---|
2223 | } |
---|
2224 | /////////////////////////////////////////////////////////////////////////////// |
---|
2225 | static |
---|
2226 | proc decompEinX(list BO) |
---|
2227 | "Internal procedure - no help and no example available |
---|
2228 | " |
---|
2229 | { |
---|
2230 | //--- decomposition of exceptional divisor, non-embedded resolution. |
---|
2231 | //--- even a single exceptional divisor may be Q-reducible when considered |
---|
2232 | //--- as divisor on the strict transform |
---|
2233 | |
---|
2234 | //---------------------------------------------------------------------------- |
---|
2235 | // Initialization |
---|
2236 | //---------------------------------------------------------------------------- |
---|
2237 | int i,j,k,de,contact; |
---|
2238 | intmat interMat; |
---|
2239 | list dcE,tmpList,prList,sa,nullList; |
---|
2240 | string mpol,compList; |
---|
2241 | def R=basering; |
---|
2242 | ideal I; |
---|
2243 | //---------------------------------------------------------------------------- |
---|
2244 | // pass to divisors on V(J) and throw away components already present as |
---|
2245 | // previous exceptional divisors |
---|
2246 | //---------------------------------------------------------------------------- |
---|
2247 | for(i=1;i<=size(BO[4]);i++) |
---|
2248 | { |
---|
2249 | I=BO[4][i]+BO[2]; |
---|
2250 | for(j=i+1;j<=size(BO[4]);j++) |
---|
2251 | { |
---|
2252 | sa=sat(I,BO[4][j]+BO[2]); |
---|
2253 | if(sa[2]) |
---|
2254 | { |
---|
2255 | I=sa[1]; |
---|
2256 | } |
---|
2257 | } |
---|
2258 | //!!! Practical improvement - not yet implemented: |
---|
2259 | //!!!hier den Input besser aufbereiten (cf. J. Wahl's example) |
---|
2260 | //!!!I[1]=x(2)^15*y(2)^9+3*x(2)^10*y(2)^6+3*x(2)^5*y(2)^3+x(2)+1; |
---|
2261 | //!!!I[2]=x(2)^8*y(2)^6+y(0); |
---|
2262 | //!!!heuristisch die Ordnung so waehlen, dass y(0) im Prinzip eliminiert |
---|
2263 | //!!!wird. |
---|
2264 | //----------------------------------------------------------------------------- |
---|
2265 | // 1) decompose exceptional divisor (over Q) |
---|
2266 | // 2) check whether there are C-reducible Q-components |
---|
2267 | // 3) if necessary, find appropriate field extension of Q to decompose |
---|
2268 | // 4) in each chart collect information in list dcE and export it |
---|
2269 | //----------------------------------------------------------------------------- |
---|
2270 | prList=primdecGTZ(I); |
---|
2271 | for(j=1;j<=size(prList);j++) |
---|
2272 | { |
---|
2273 | tmpList=grad(prList[j][2]); |
---|
2274 | de=tmpList[1]; |
---|
2275 | interMat=tmpList[2]; |
---|
2276 | mpol=tmpList[3]; |
---|
2277 | compList=tmpList[4]; |
---|
2278 | nullList=tmpList[5]; |
---|
2279 | contact=Kontakt(prList[j][1],BO[2]); |
---|
2280 | tmpList=prList[j][2],de,contact,interMat,mpol,compList,nullList; |
---|
2281 | prList[j]=tmpList; |
---|
2282 | } |
---|
2283 | dcE[i]=prList; |
---|
2284 | } |
---|
2285 | return(dcE); |
---|
2286 | } |
---|
2287 | ////////////////////////////////////////////////////////////////////////////// |
---|
2288 | static |
---|
2289 | proc getMinpoly(poly p) |
---|
2290 | "Internal procedure - no help and no example available |
---|
2291 | " |
---|
2292 | { |
---|
2293 | //---assume that p is a polynomial in 2 variables and irreducible |
---|
2294 | //---over Q. Computes an irreducible polynomial mp in one variable |
---|
2295 | //---over Q such that p splits completely over the splitting field of mp |
---|
2296 | //---returns mp as a string |
---|
2297 | //---use a variant of the algorithm of S. Gao |
---|
2298 | def R=basering; |
---|
2299 | int i,j,k,a,b,m,n; |
---|
2300 | intvec v; |
---|
2301 | string mp="poly p=t-1;"; |
---|
2302 | list Li=string(1); |
---|
2303 | list re=mp,Li,1; |
---|
2304 | |
---|
2305 | //---check which variables occur in p |
---|
2306 | for(i=1;i<=nvars(basering);i++) |
---|
2307 | { |
---|
2308 | if(p!=subst(p,var(i),0)){v[size(v)+1]=i;} |
---|
2309 | } |
---|
2310 | |
---|
2311 | //---the poly is constant |
---|
2312 | if(size(v)==1){return(re);} |
---|
2313 | |
---|
2314 | //---the poly depends only on one variable or is homogeneous |
---|
2315 | //---in 2 variables |
---|
2316 | if((size(v)==2)||((size(v)==3)&&(homog(p)))) |
---|
2317 | { |
---|
2318 | if((size(v)==3)&&(homog(p))) |
---|
2319 | { |
---|
2320 | p=subst(p,var(v[3]),1); |
---|
2321 | } |
---|
2322 | ring Rhelp=0,var(v[2]),dp; |
---|
2323 | poly p=imap(R,p); |
---|
2324 | ring Shelp=0,t,dp; |
---|
2325 | poly p=fetch(Rhelp,p); |
---|
2326 | int de=deg(p); |
---|
2327 | p=simplifyMinpoly(p); |
---|
2328 | Li=getNumZeros(p); |
---|
2329 | short=0; |
---|
2330 | mp="poly p="+string(p)+";"; |
---|
2331 | re=mp,Li,de; |
---|
2332 | setring R; |
---|
2333 | return(re); |
---|
2334 | } |
---|
2335 | v=v[2..size(v)]; |
---|
2336 | if(size(v)>2){ERROR("getMinpoly:input depends on more then 2 variables");} |
---|
2337 | |
---|
2338 | //---the general case, the poly is considered as poly in x an y now |
---|
2339 | ring T=0,(x,y),lp; |
---|
2340 | ideal M,N; |
---|
2341 | M[nvars(R)]=0; |
---|
2342 | N[nvars(R)]=0; |
---|
2343 | M[v[1]]=x; |
---|
2344 | N[v[1]]=y; |
---|
2345 | M[v[2]]=y; |
---|
2346 | N[v[2]]=x; |
---|
2347 | map phi=R,M; |
---|
2348 | map psi=R,N; |
---|
2349 | poly p=phi(p); |
---|
2350 | poly q=psi(p); |
---|
2351 | ring Thelp=(0,x),y,dp; |
---|
2352 | poly p=imap(T,p); |
---|
2353 | poly q=imap(T,q); |
---|
2354 | n=deg(p); //---the degree with respect to y |
---|
2355 | m=deg(q); //---the degree with respect to x |
---|
2356 | setring T; |
---|
2357 | ring A=0,(u(1..m*(n+1)),v(1..(m+1)*n),x,y,t),dp; |
---|
2358 | poly f=imap(T,p); |
---|
2359 | poly g; |
---|
2360 | poly h; |
---|
2361 | for(i=0;i<=m-1;i++) |
---|
2362 | { |
---|
2363 | for(j=0;j<=n;j++) |
---|
2364 | { |
---|
2365 | g=g+u(i*(n+1)+j+1)*x^i*y^j; |
---|
2366 | } |
---|
2367 | } |
---|
2368 | for(i=0;i<=m;i++) |
---|
2369 | { |
---|
2370 | for(j=0;j<=n-1;j++) |
---|
2371 | { |
---|
2372 | h=h+v(i*n+j+1)*x^i*y^j; |
---|
2373 | } |
---|
2374 | } |
---|
2375 | poly L=f*(diff(g,y)-diff(h,x))+h*diff(f,x)-g*diff(f,y); |
---|
2376 | //---according to the theory f is absolutely irreducible iff |
---|
2377 | //---L(g,h)=0 has no non-trivial solution g,h |
---|
2378 | //---(g=diff(f,x),h=diff(f,y) is always a solution) |
---|
2379 | //---therefore we compute a vector space basis of G |
---|
2380 | //---G={g in Q[x,y],deg_x(g)<m,|exist h, such that L(g,h)=0} |
---|
2381 | //---dim(G)=a is the number of factors of f in C[x,y] |
---|
2382 | matrix M=coef(L,xy); |
---|
2383 | ideal J=M[2,1..ncols(M)]; |
---|
2384 | option(redSB); |
---|
2385 | J=std(J); |
---|
2386 | option(noredSB); |
---|
2387 | poly gred=reduce(g,J); |
---|
2388 | ideal G; |
---|
2389 | for(i=1;i<=m*(n+1);i++) |
---|
2390 | { |
---|
2391 | if(gred!=subst(gred,u(i),0)) |
---|
2392 | { |
---|
2393 | G[size(G)+1]=subst(gred,u(i),1); |
---|
2394 | } |
---|
2395 | } |
---|
2396 | for(i=1;i<=n*(m+1);i++) |
---|
2397 | { |
---|
2398 | if(gred!=subst(gred,v(i),0)) |
---|
2399 | { |
---|
2400 | G[size(G)+1]=subst(gred,v(i),1); |
---|
2401 | } |
---|
2402 | } |
---|
2403 | for(i=1;i<=m*(n+1);i++) |
---|
2404 | { |
---|
2405 | G=subst(G,u(i),0); |
---|
2406 | } |
---|
2407 | for(i=1;i<=n*(m+1);i++) |
---|
2408 | { |
---|
2409 | G=subst(G,v(i),0); |
---|
2410 | } |
---|
2411 | //---the number of factors in C[x,y] |
---|
2412 | a=size(G); |
---|
2413 | for(i=1;i<=a;i++) |
---|
2414 | { |
---|
2415 | G[i]=simplify(G[i],1); |
---|
2416 | } |
---|
2417 | if(a==1) |
---|
2418 | { |
---|
2419 | //---f is absolutely irreducible |
---|
2420 | setring R; |
---|
2421 | return(re); |
---|
2422 | } |
---|
2423 | //---let g in G be any non-trivial element (g not in <diff(f,x)>) |
---|
2424 | //---according to the theory f=product over all c in C of the |
---|
2425 | //---gcd(f,g-c*diff(f,x)) |
---|
2426 | //---let g_1,...,g_a be a basis of G and write |
---|
2427 | //---g*g_i=sum a_ij*g_j*diff(f,x) mod f |
---|
2428 | //---let B=(a_ij) and ch=det(t*unitmat(a)-B) the characteristic |
---|
2429 | //---polynomial then the number of distinct irreducible factors |
---|
2430 | //---of gcd(f,g-c*diff(f,x)) in C[x,y] is equal to the multiplicity |
---|
2431 | //---of c as a root of ch. |
---|
2432 | //---in our special situation (f is irreducible over Q) ch should |
---|
2433 | //---be irreducible and the different roots of ch lead to the |
---|
2434 | //---factors of f, i.e. ch is the minpoly we are looking for |
---|
2435 | |
---|
2436 | poly fh=homog(f,t); |
---|
2437 | //---homogenization is used to obtain a constant matrix using lift |
---|
2438 | ideal Gh=homog(G,t); |
---|
2439 | int dh,df; |
---|
2440 | df=deg(fh); |
---|
2441 | for(i=1;i<=a;i++) |
---|
2442 | { |
---|
2443 | if(deg(Gh[i])>dh){dh=deg(Gh[i]);} |
---|
2444 | } |
---|
2445 | for(i=1;i<=a;i++) |
---|
2446 | { |
---|
2447 | Gh[i]=t^(dh-deg(Gh[i]))*Gh[i]; |
---|
2448 | } |
---|
2449 | ideal GF=simplify(diff(fh,x),1)*Gh,fh; |
---|
2450 | poly ch; |
---|
2451 | matrix LI; |
---|
2452 | matrix B[a][a]; |
---|
2453 | matrix E=unitmat(a); |
---|
2454 | poly gran; |
---|
2455 | ideal fac; |
---|
2456 | for(i=1;i<=a;i++) |
---|
2457 | { |
---|
2458 | LI=lift(GF,t^(df-1-dh)*Gh[i]*Gh); |
---|
2459 | B=LI[1..a,1..a]; |
---|
2460 | ch=det(t*E-B); |
---|
2461 | //---irreducibility test |
---|
2462 | fac=factorize(ch,1); |
---|
2463 | if(deg(fac[1])==a) |
---|
2464 | { |
---|
2465 | ch=simplifyMinpoly(ch); |
---|
2466 | Li=getNumZeros(ch); |
---|
2467 | int de=deg(ch); |
---|
2468 | short=0; |
---|
2469 | mp="poly p="+string(ch)+";"; |
---|
2470 | re=mp,Li,de; |
---|
2471 | setring R; |
---|
2472 | return(re); |
---|
2473 | } |
---|
2474 | } |
---|
2475 | ERROR("getMinpoly:not found:please send the example to the authors"); |
---|
2476 | } |
---|
2477 | ////////////////////////////////////////////////////////////////////////////// |
---|
2478 | static |
---|
2479 | proc getNumZeros(poly p) |
---|
2480 | "Internal procedure - no help and no example available |
---|
2481 | " |
---|
2482 | { |
---|
2483 | //--- compute numerically (!!!) the zeros of the minimal polynomial |
---|
2484 | def R=basering; |
---|
2485 | ring S=0,t,dp; |
---|
2486 | poly p=imap(R,p); |
---|
2487 | def L=laguerre_solve(p,30); |
---|
2488 | //!!! practical improvement: |
---|
2489 | //!!! testen ob die Nullstellen signifikant verschieden sind |
---|
2490 | //!!! und im Notfall Genauigkeit erhoehen |
---|
2491 | list re; |
---|
2492 | int i; |
---|
2493 | for(i=1;i<=size(L);i++) |
---|
2494 | { |
---|
2495 | re[i]=string(L[i]); |
---|
2496 | } |
---|
2497 | setring R; |
---|
2498 | return(re); |
---|
2499 | } |
---|
2500 | ////////////////////////////////////////////////////////////////////////////// |
---|
2501 | static |
---|
2502 | proc simplifyMinpoly(poly p) |
---|
2503 | "Internal procedure - no help and no example available |
---|
2504 | " |
---|
2505 | { |
---|
2506 | //--- describe field extension in a simple way |
---|
2507 | p=cleardenom(p); |
---|
2508 | int n=int(leadcoef(p)); |
---|
2509 | int d=deg(p); |
---|
2510 | int i,k; |
---|
2511 | int re=1; |
---|
2512 | number s=1; |
---|
2513 | |
---|
2514 | list L=primefactors(n); |
---|
2515 | |
---|
2516 | for(i=1;i<=size(L[1]);i++) |
---|
2517 | { |
---|
2518 | k=L[2][i] mod d; |
---|
2519 | s=1/number((L[1][i])^(L[2][i]/d)); |
---|
2520 | if(!k){p=subst(p,t,s*t);} |
---|
2521 | } |
---|
2522 | p=cleardenom(p); |
---|
2523 | n=int(leadcoef(subst(p,t,0))); |
---|
2524 | L=primefactors(n); |
---|
2525 | for(i=1;i<=size(L[1]);i++) |
---|
2526 | { |
---|
2527 | k=L[2][i] mod d; |
---|
2528 | s=(L[1][i])^(L[2][i]/d); |
---|
2529 | if(!k){p=subst(p,t,s*t);} |
---|
2530 | } |
---|
2531 | p=cleardenom(p); |
---|
2532 | return(p); |
---|
2533 | } |
---|
2534 | /////////////////////////////////////////////////////////////////////////////// |
---|
2535 | static |
---|
2536 | proc grad(ideal I) |
---|
2537 | "Internal procedure - no help and no example available |
---|
2538 | " |
---|
2539 | { |
---|
2540 | //--- computes the number of components over C |
---|
2541 | //--- for a prime ideal of height 1 over Q |
---|
2542 | def R=basering; |
---|
2543 | int n=nvars(basering); |
---|
2544 | string mp="poly p=t-1;"; |
---|
2545 | string str=string(1); |
---|
2546 | list zeroList=string(1); |
---|
2547 | int i,j,k,l,d,e,c,mi; |
---|
2548 | ideal Istd=std(I); |
---|
2549 | intmat interMat; |
---|
2550 | d=dim(Istd); |
---|
2551 | if(d==-1){return(list(0,0,mp,str,zeroList));} |
---|
2552 | if(d!=1){ERROR("ideal is not one-dimensional");} |
---|
2553 | ideal Sloc=std(slocus(I)); |
---|
2554 | if(deg(Sloc[1])>0) |
---|
2555 | { |
---|
2556 | //---This is only to test that in case of singularities we have only |
---|
2557 | //---one singular point which is a normal crossing |
---|
2558 | //---consider the different singular points |
---|
2559 | ideal M; |
---|
2560 | list pr=minAssGTZ(Sloc); |
---|
2561 | if(size(pr)>1){ERROR("grad:more then one singular point");} |
---|
2562 | for(l=1;l<=size(pr);l++) |
---|
2563 | { |
---|
2564 | M=std(pr[l]); |
---|
2565 | d=vdim(M); |
---|
2566 | if(d!=1) |
---|
2567 | { |
---|
2568 | //---now we have to extend the field |
---|
2569 | if(defined(S)){kill S;} |
---|
2570 | ring S=0,x(1..n),lp; |
---|
2571 | ideal M=fetch(R,M); |
---|
2572 | ideal I=fetch(R,I); |
---|
2573 | ideal jmap; |
---|
2574 | map phi=S,maxideal(1);; |
---|
2575 | ideal Mstd=std(M); |
---|
2576 | //---M has to be in general position with respect to lp, i.e. |
---|
2577 | //---vdim(M)=deg(M[1]) |
---|
2578 | poly p=Mstd[1]; |
---|
2579 | e=vdim(Mstd); |
---|
2580 | while(e!=deg(p)) |
---|
2581 | { |
---|
2582 | jmap=randomLast(100); |
---|
2583 | phi=S,jmap; |
---|
2584 | Mstd=std(phi(M)); |
---|
2585 | p=Mstd[1]; |
---|
2586 | } |
---|
2587 | I=phi(I); |
---|
2588 | kill phi; |
---|
2589 | //---now it is in general position an M[1] defines the field extension |
---|
2590 | //---Q[x]/M over Q |
---|
2591 | ring Shelp=0,t,dp; |
---|
2592 | ideal helpmap; |
---|
2593 | helpmap[n]=t; |
---|
2594 | map psi=S,helpmap; |
---|
2595 | poly p=psi(p); |
---|
2596 | ring T=(0,t),x(1..n),lp; |
---|
2597 | poly p=imap(Shelp,p); |
---|
2598 | //---we are now in the polynomial ring over the field Q[x]/M |
---|
2599 | minpoly=leadcoef(p); |
---|
2600 | ideal M=imap(S,Mstd); |
---|
2601 | M=M,var(n)-t; |
---|
2602 | ideal I=fetch(S,I); |
---|
2603 | } |
---|
2604 | //---we construct a map phi which maps M to maxideal(1) |
---|
2605 | option(redSB); |
---|
2606 | ideal Mstd=-simplify(std(M),1); |
---|
2607 | option(noredSB); |
---|
2608 | for(i=1;i<=n;i++) |
---|
2609 | { |
---|
2610 | Mstd=subst(Mstd,var(i),-var(i)); |
---|
2611 | M[n-i+1]=Mstd[i]; |
---|
2612 | } |
---|
2613 | M=M[1..n]; |
---|
2614 | //---go to the localization with respect to <x> |
---|
2615 | if(d!=1) |
---|
2616 | { |
---|
2617 | ring Tloc=(0,t),x(1..n),ds; |
---|
2618 | poly p=imap(Shelp,p); |
---|
2619 | minpoly=leadcoef(p); |
---|
2620 | ideal M=fetch(T,M); |
---|
2621 | map phi=T,M; |
---|
2622 | } |
---|
2623 | else |
---|
2624 | { |
---|
2625 | ring Tloc=0,x(1..n),ds; |
---|
2626 | ideal M=fetch(R,M); |
---|
2627 | map phi=R,M; |
---|
2628 | } |
---|
2629 | ideal I=phi(I); |
---|
2630 | ideal Istd=std(I); |
---|
2631 | mi=mi+milnor(Istd); |
---|
2632 | if(mi>l) |
---|
2633 | { |
---|
2634 | ERROR("grad:divisor is really singular"); |
---|
2635 | } |
---|
2636 | setring R; |
---|
2637 | } |
---|
2638 | } |
---|
2639 | intvec ind=indepSet(Istd,1)[1]; |
---|
2640 | for(i=1;i<=n;i++){if(ind[i]) break;} |
---|
2641 | //---the i-th variable is the independent one |
---|
2642 | ring Shelp=0,x(1..n),dp; |
---|
2643 | ideal I=fetch(R,I); |
---|
2644 | if(defined(S)){kill S;} |
---|
2645 | if(i==1){ring S=(0,x(1)),x(2..n),lp;} |
---|
2646 | if(i==n){ring S=(0,x(n)),x(1..n-1),lp;} |
---|
2647 | if((i!=1)&&(i!=n)){ring S=(0,x(i)),(x(1..i-1),x(i+1..n)),lp;} |
---|
2648 | //---I is zero-dimensional now |
---|
2649 | ideal I=imap(Shelp,I); |
---|
2650 | ideal Istd=std(I); |
---|
2651 | ideal jmap; |
---|
2652 | map phi; |
---|
2653 | poly p=Istd[1]; |
---|
2654 | e=vdim(Istd); |
---|
2655 | if(e==1) |
---|
2656 | { |
---|
2657 | setring R; |
---|
2658 | str=string(I); |
---|
2659 | list resi=1,interMat,mp,str,zeroList; |
---|
2660 | return(resi); |
---|
2661 | } |
---|
2662 | //---move I to general position with respect to lp |
---|
2663 | if(e!=deg(p)) |
---|
2664 | { |
---|
2665 | jmap=randomLast(5); |
---|
2666 | phi=S,jmap; |
---|
2667 | Istd=std(phi(I)); |
---|
2668 | p=Istd[1]; |
---|
2669 | } |
---|
2670 | while(e!=deg(p)) |
---|
2671 | { |
---|
2672 | jmap=randomLast(100); |
---|
2673 | phi=S,jmap; |
---|
2674 | Istd=std(phi(I)); |
---|
2675 | p=Istd[1]; |
---|
2676 | } |
---|
2677 | setring Shelp; |
---|
2678 | poly p=imap(S,p); |
---|
2679 | list Q=getMinpoly(p); |
---|
2680 | int de=Q[3]; |
---|
2681 | mp=Q[1]; |
---|
2682 | //!!!diese Stelle effizienter machen |
---|
2683 | //!!!minAssGTZ vermeiden durch direkte Betrachtung von |
---|
2684 | //!!!p und mp und evtl. Quotientenbildung |
---|
2685 | //!!!bisher nicht zeitkritisch |
---|
2686 | string Tesr="ring Tes=(0,t),("+varstr(R)+"),dp;"; |
---|
2687 | execute(Tesr); |
---|
2688 | execute(mp); |
---|
2689 | minpoly=leadcoef(p); |
---|
2690 | ideal I=fetch(R,I); |
---|
2691 | list pr=minAssGTZ(I); |
---|
2692 | ideal allgEbene=randomLast(100)[nvars(basering)]; |
---|
2693 | int minpts=vdim(std(I+allgEbene)); |
---|
2694 | ideal tempi; |
---|
2695 | j=1; |
---|
2696 | for(i=1;i<=size(pr);i++) |
---|
2697 | { |
---|
2698 | tempi=std(pr[i]+allgEbene); |
---|
2699 | if(vdim(tempi)<minpts) |
---|
2700 | { |
---|
2701 | minpts=vdim(tempi); |
---|
2702 | j=i; |
---|
2703 | } |
---|
2704 | } |
---|
2705 | tempi=pr[j]; |
---|
2706 | str=string(tempi); |
---|
2707 | kill interMat; |
---|
2708 | setring R; |
---|
2709 | intmat interMat[de][de]=intersComp(str,mp,Q[2],str,mp,Q[2]); |
---|
2710 | list resi=de,interMat,mp,str,Q[2]; |
---|
2711 | return(resi); |
---|
2712 | } |
---|
2713 | //////////////////////////////////////////////////////////////////////////// |
---|
2714 | static |
---|
2715 | proc Kontakt(ideal I, ideal K) |
---|
2716 | "Internal procedure - no help and no example available |
---|
2717 | " |
---|
2718 | { |
---|
2719 | //---Let K be a prime ideal and I an ideal not contained in K |
---|
2720 | //---computes a maximalideal M=<x(1)-a1,...,x(n)-an>, ai in a field |
---|
2721 | //---extension of Q, containing I+K and an integer a |
---|
2722 | //---such that in the localization of the polynomial ring with |
---|
2723 | //---respect to M the ideal I is not contained in K+M^a+1 but in M^a in |
---|
2724 | def R=basering; |
---|
2725 | int n=nvars(basering); |
---|
2726 | int i,j,k,d,e; |
---|
2727 | ideal J=std(I+K); |
---|
2728 | if(dim(J)==-1){return(0);} |
---|
2729 | ideal W; |
---|
2730 | //---choice of the maximal ideal M |
---|
2731 | for(i=1;i<=n;i++) |
---|
2732 | { |
---|
2733 | W=std(J,var(i)); |
---|
2734 | d=dim(W); |
---|
2735 | if(d==0) break; |
---|
2736 | } |
---|
2737 | i=1;k=2; |
---|
2738 | while((d)&&(i<n)) |
---|
2739 | { |
---|
2740 | W=std(J,var(i)+var(k)); |
---|
2741 | d=dim(W); |
---|
2742 | if(k==n){i++;k=i;} |
---|
2743 | if(k<n){k++;} |
---|
2744 | } |
---|
2745 | while(d) |
---|
2746 | { |
---|
2747 | W=std(J,randomid(maxideal(1))[1]); |
---|
2748 | d=dim(W); |
---|
2749 | } |
---|
2750 | //---now we have a collection om maximalideals and choose one with dim Q[x]/M |
---|
2751 | //---minimal |
---|
2752 | list pr=minAssGTZ(W); |
---|
2753 | d=vdim(std(pr[1])); |
---|
2754 | k=1; |
---|
2755 | for(i=2;i<=size(pr);i++) |
---|
2756 | { |
---|
2757 | if(d==1) break; |
---|
2758 | e=vdim(std(pr[i])); |
---|
2759 | if(e<d){k=i;d=e;} |
---|
2760 | } |
---|
2761 | //---M is fixed now |
---|
2762 | //---if dim Q[x]/M =1 we localize at M |
---|
2763 | ideal M=pr[k]; |
---|
2764 | if(d!=1) |
---|
2765 | { |
---|
2766 | //---now we have to extend the field |
---|
2767 | if(defined(S)){kill S;} |
---|
2768 | ring S=0,x(1..n),lp; |
---|
2769 | ideal M=fetch(R,M); |
---|
2770 | ideal I=fetch(R,I); |
---|
2771 | ideal K=fetch(R,K); |
---|
2772 | ideal jmap; |
---|
2773 | map phi=S,maxideal(1);; |
---|
2774 | ideal Mstd=std(M); |
---|
2775 | //---M has to be in general position with respect to lp, i.e. |
---|
2776 | //---vdim(M)=deg(M[1]) |
---|
2777 | poly p=Mstd[1]; |
---|
2778 | e=vdim(Mstd); |
---|
2779 | while(e!=deg(p)) |
---|
2780 | { |
---|
2781 | jmap=randomLast(100); |
---|
2782 | phi=S,jmap; |
---|
2783 | Mstd=std(phi(M)); |
---|
2784 | p=Mstd[1]; |
---|
2785 | } |
---|
2786 | I=phi(I); |
---|
2787 | K=phi(K); |
---|
2788 | kill phi; |
---|
2789 | //---now it is in general position an M[1] defines the field extension |
---|
2790 | //---Q[x]/M over Q |
---|
2791 | ring Shelp=0,t,dp; |
---|
2792 | ideal helpmap; |
---|
2793 | helpmap[n]=t; |
---|
2794 | map psi=S,helpmap; |
---|
2795 | poly p=psi(p); |
---|
2796 | ring T=(0,t),x(1..n),lp; |
---|
2797 | poly p=imap(Shelp,p); |
---|
2798 | //---we are now in the polynomial ring over the field Q[x]/M |
---|
2799 | minpoly=leadcoef(p); |
---|
2800 | ideal M=imap(S,Mstd); |
---|
2801 | M=M,var(n)-t; |
---|
2802 | ideal I=fetch(S,I); |
---|
2803 | ideal K=fetch(S,K); |
---|
2804 | } |
---|
2805 | //---we construct a map phi which maps M to maxideal(1) |
---|
2806 | option(redSB); |
---|
2807 | ideal Mstd=-simplify(std(M),1); |
---|
2808 | option(noredSB); |
---|
2809 | for(i=1;i<=n;i++) |
---|
2810 | { |
---|
2811 | Mstd=subst(Mstd,var(i),-var(i)); |
---|
2812 | M[n-i+1]=Mstd[i]; |
---|
2813 | } |
---|
2814 | M=M[1..n]; |
---|
2815 | //---go to the localization with respect to <x> |
---|
2816 | if(d!=1) |
---|
2817 | { |
---|
2818 | ring Tloc=(0,t),x(1..n),ds; |
---|
2819 | poly p=imap(Shelp,p); |
---|
2820 | minpoly=leadcoef(p); |
---|
2821 | ideal M=fetch(T,M); |
---|
2822 | map phi=T,M; |
---|
2823 | } |
---|
2824 | else |
---|
2825 | { |
---|
2826 | ring Tloc=0,x(1..n),ds; |
---|
2827 | ideal M=fetch(R,M); |
---|
2828 | map phi=R,M; |
---|
2829 | } |
---|
2830 | ideal K=phi(K); |
---|
2831 | ideal I=phi(I); |
---|
2832 | //---compute the order of I in (Q[x]/M)[[x]]/K |
---|
2833 | k=1;d=0; |
---|
2834 | while(!d) |
---|
2835 | { |
---|
2836 | k++; |
---|
2837 | d=size(reduce(I,std(maxideal(k)+K))); |
---|
2838 | } |
---|
2839 | setring R; |
---|
2840 | return(k-1); |
---|
2841 | } |
---|
2842 | example |
---|
2843 | {"EXAMPLE:"; |
---|
2844 | echo = 2; |
---|
2845 | ring r = 0,(x,y,z),dp; |
---|
2846 | ideal I=x4+z4+1; |
---|
2847 | ideal K=x+y2+z2; |
---|
2848 | Kontakt(I,K); |
---|
2849 | } |
---|
2850 | ////////////////////////////////////////////////////////////////////////////// |
---|
2851 | static |
---|
2852 | proc abstractNC(list BO) |
---|
2853 | "Internal procedure - no help and no example available |
---|
2854 | " |
---|
2855 | { |
---|
2856 | //--- check normal crossing property |
---|
2857 | //--- used for passing from embedded to non-embedded resolution |
---|
2858 | //---------------------------------------------------------------------------- |
---|
2859 | // Initialization |
---|
2860 | //---------------------------------------------------------------------------- |
---|
2861 | int i,k,j,flag; |
---|
2862 | list L; |
---|
2863 | ideal J; |
---|
2864 | if(dim(std(cent))>0){return(1);} |
---|
2865 | //---------------------------------------------------------------------------- |
---|
2866 | // check each exceptional divisor on V(J) |
---|
2867 | //---------------------------------------------------------------------------- |
---|
2868 | for(i=1;i<=size(BO[4]);i++) |
---|
2869 | { |
---|
2870 | if(dim(std(BO[2]+BO[4][i]))>0) |
---|
2871 | { |
---|
2872 | //--- really something to do |
---|
2873 | J=radical(BO[4][i]+BO[2]); |
---|
2874 | if(deg(std(slocus(J))[1])!=0) |
---|
2875 | { |
---|
2876 | if(!nodes(J)) |
---|
2877 | { |
---|
2878 | //--- really singular, not only nodes ==> not normal crossing |
---|
2879 | return(0); |
---|
2880 | } |
---|
2881 | } |
---|
2882 | for(k=1;k<=size(L);k++) |
---|
2883 | { |
---|
2884 | //--- run through previously considered divisors |
---|
2885 | //--- we do not want to bother with the same one twice |
---|
2886 | if((size(reduce(J,std(L[k])))==0)&&(size(reduce(L[k],std(J)))==0)) |
---|
2887 | { |
---|
2888 | //--- already considered this one |
---|
2889 | flag=1;break; |
---|
2890 | } |
---|
2891 | //--- drop previously considered exceptional divisors from the current one |
---|
2892 | J=sat(J,L[k])[1]; |
---|
2893 | if(deg(std(J)[1])==0) |
---|
2894 | { |
---|
2895 | //--- nothing remaining |
---|
2896 | flag=1;break; |
---|
2897 | } |
---|
2898 | } |
---|
2899 | if(flag==0) |
---|
2900 | { |
---|
2901 | //--- add exceptional divisor to the list |
---|
2902 | L[size(L)+1]=J; |
---|
2903 | } |
---|
2904 | flag=0; |
---|
2905 | } |
---|
2906 | } |
---|
2907 | //--------------------------------------------------------------------------- |
---|
2908 | // check intersection properties between different exceptional divisors |
---|
2909 | //--------------------------------------------------------------------------- |
---|
2910 | for(k=1;k<size(L);k++) |
---|
2911 | { |
---|
2912 | for(i=k+1;i<=size(L);i++) |
---|
2913 | { |
---|
2914 | if(!nodes(intersect(L[k],L[i]))) |
---|
2915 | { |
---|
2916 | //--- divisors Ek and Ei do not meet in a node but in a singularity |
---|
2917 | //--- which is not allowed to occur ==> not normal crossing |
---|
2918 | return(0); |
---|
2919 | } |
---|
2920 | for(j=i+1;j<=size(L);j++) |
---|
2921 | { |
---|
2922 | if(deg(std(L[i]+L[j]+L[k])[1])>0) |
---|
2923 | { |
---|
2924 | //--- three divisors meet simultaneously ==> not normal crossing |
---|
2925 | return(0); |
---|
2926 | } |
---|
2927 | } |
---|
2928 | } |
---|
2929 | } |
---|
2930 | //--- we reached this point ==> normal crossing |
---|
2931 | return(1); |
---|
2932 | } |
---|
2933 | ////////////////////////////////////////////////////////////////////////////// |
---|
2934 | static |
---|
2935 | proc nodes(ideal J) |
---|
2936 | "Internal procedure - no help and no example available |
---|
2937 | " |
---|
2938 | { |
---|
2939 | //--- check whether at most nodes occur as singularities |
---|
2940 | ideal K=std(slocus(J)); |
---|
2941 | if(deg(K[1])==0){return(1);} |
---|
2942 | if(dim(K)>0){return(0);} |
---|
2943 | if(vdim(K)!=vdim(std(radical(K)))){return(0);} |
---|
2944 | return(1); |
---|
2945 | } |
---|
2946 | ////////////////////////////////////////////////////////////////////////////// |
---|
2947 | |
---|
2948 | proc intersectionDiv(list re) |
---|
2949 | "USAGE: intersectionDiv(L); |
---|
2950 | @* L = list of rings |
---|
2951 | ASSUME: L is output of resolution of singularities |
---|
2952 | (only case of isolated surface singularities) |
---|
2953 | COMPUTE: intersection matrix and genera of the exceptional divisors |
---|
2954 | (considered as curves on the strict transform) |
---|
2955 | RETURN: list l, where |
---|
2956 | l[1]: intersection matrix of exceptional divisors |
---|
2957 | l[2]: intvec, genera of exceptional divisors |
---|
2958 | l[3]: divisorList, encoding the identification of the divisors |
---|
2959 | EXAMPLE: example intersectionDiv; shows an example |
---|
2960 | " |
---|
2961 | { |
---|
2962 | //---------------------------------------------------------------------------- |
---|
2963 | //--- Computes in case of surface singularities (non-embedded resolution): |
---|
2964 | //--- the intersection of the divisors (on the surface) |
---|
2965 | //--- assuming that re=resolve(J) |
---|
2966 | //---------------------------------------------------------------------------- |
---|
2967 | def R=basering; |
---|
2968 | //---Test whether we are in the irreducible surface case |
---|
2969 | def S=re[2][1]; |
---|
2970 | setring S; |
---|
2971 | BO[2]=BO[2]+BO[1]; // make sure we are living in the smooth W |
---|
2972 | if(dim(std(BO[2]))!=2) |
---|
2973 | { |
---|
2974 | ERROR("The given original object is not a surface"); |
---|
2975 | } |
---|
2976 | if(dim(std(slocus(BO[2])))>0) |
---|
2977 | { |
---|
2978 | ERROR("The given original object has non-isolated singularities."); |
---|
2979 | } |
---|
2980 | setring R; |
---|
2981 | //---------------------------------------------------------------------------- |
---|
2982 | // Compute a non-embedded resolution from the given embedded one by |
---|
2983 | // dropping redundant trailing blow-ups |
---|
2984 | //---------------------------------------------------------------------------- |
---|
2985 | list resu,tmpiden,templist; |
---|
2986 | intvec divcomp; |
---|
2987 | int i,j,k,offset1,offset2,a,b,c,d,q,found; |
---|
2988 | //--- compute non-embedded resolution |
---|
2989 | list abst=abstractR(re); |
---|
2990 | intvec endiv=abst[1]; |
---|
2991 | intvec deleted=abst[2]; |
---|
2992 | //--- identify the divisors in the various final charts |
---|
2993 | list iden=collectDiv(re,deleted)[2]; |
---|
2994 | // list of final divisors |
---|
2995 | list iden0=iden; // backup copy of iden for later use |
---|
2996 | |
---|
2997 | iden=delete(iden,size(iden)); // drop list of endRings from iden |
---|
2998 | //--------------------------------------------------------------------------- |
---|
2999 | // In iden, only the final charts should be listed, whereas iden0 contains |
---|
3000 | // everything. |
---|
3001 | //--------------------------------------------------------------------------- |
---|
3002 | for(i=1;i<=size(iden);i++) |
---|
3003 | { |
---|
3004 | k=size(iden[i]); |
---|
3005 | tmpiden=iden[i]; |
---|
3006 | for(j=k;j>0;j--) |
---|
3007 | { |
---|
3008 | if(!endiv[iden[i][j][1]]) |
---|
3009 | { |
---|
3010 | //---not a final chart |
---|
3011 | tmpiden=delete(tmpiden,j); |
---|
3012 | } |
---|
3013 | } |
---|
3014 | if(size(tmpiden)==0) |
---|
3015 | { |
---|
3016 | //--- oops, this divisor does not appear in final charts |
---|
3017 | iden=delete(iden,i); |
---|
3018 | continue; |
---|
3019 | } |
---|
3020 | else |
---|
3021 | { |
---|
3022 | iden[i]=tmpiden; |
---|
3023 | } |
---|
3024 | } |
---|
3025 | //--------------------------------------------------------------------------- |
---|
3026 | // Even though the exceptional divisors were irreducible in the embedded |
---|
3027 | // case, they may very well have become reducible after intersection with |
---|
3028 | // the strict transform of the original object. |
---|
3029 | // ===> compute a decomposition for each divisor in each of the final charts |
---|
3030 | // and change the entries of iden accordingly |
---|
3031 | // In particular, it is important to keep track of the identification of the |
---|
3032 | // components of the divisors in each of the charts |
---|
3033 | //--------------------------------------------------------------------------- |
---|
3034 | int n=size(iden); |
---|
3035 | for(i=1;i<=size(re[2]);i++) |
---|
3036 | { |
---|
3037 | if(endiv[i]) |
---|
3038 | { |
---|
3039 | def SN=re[2][i]; |
---|
3040 | setring SN; |
---|
3041 | if(defined(dcE)){kill dcE;} |
---|
3042 | list dcE=decompEinX(BO); // decomposition of exceptional divisors |
---|
3043 | export(dcE); |
---|
3044 | setring R; |
---|
3045 | kill SN; |
---|
3046 | } |
---|
3047 | } |
---|
3048 | if(defined(tmpiden)){kill tmpiden;} |
---|
3049 | list tmpiden=iden; |
---|
3050 | for(i=1;i<=size(iden);i++) |
---|
3051 | { |
---|
3052 | for(j=size(iden[i]);j>0;j--) |
---|
3053 | { |
---|
3054 | def SN=re[2][iden[i][j][1]]; |
---|
3055 | setring SN; |
---|
3056 | if(size(dcE[iden[i][j][2]])==1) |
---|
3057 | { |
---|
3058 | if(dcE[iden[i][j][2]][1][2]==0) |
---|
3059 | { |
---|
3060 | tmpiden[i]=delete(tmpiden[i],j); |
---|
3061 | } |
---|
3062 | } |
---|
3063 | setring R; |
---|
3064 | kill SN; |
---|
3065 | } |
---|
3066 | } |
---|
3067 | for(i=size(tmpiden);i>0;i--) |
---|
3068 | { |
---|
3069 | if(size(tmpiden[i])==0) |
---|
3070 | { |
---|
3071 | tmpiden=delete(tmpiden,i); |
---|
3072 | } |
---|
3073 | } |
---|
3074 | iden=tmpiden; |
---|
3075 | kill tmpiden; |
---|
3076 | list tmpiden; |
---|
3077 | //--- change entries of iden accordingly |
---|
3078 | for(i=1;i<=size(iden);i++) |
---|
3079 | { |
---|
3080 | //--- first set up new entries in iden if necessary - using the first chart |
---|
3081 | //--- in which we see the respective exceptional divisor |
---|
3082 | if(defined(S)){kill S;} |
---|
3083 | def S=re[2][iden[i][1][1]]; |
---|
3084 | //--- considering first entry for i-th divisor |
---|
3085 | setring S; |
---|
3086 | a=size(dcE[iden[i][1][2]]); |
---|
3087 | for(j=1;j<=a;j++) |
---|
3088 | { |
---|
3089 | //--- reducible - add to the list considering each component as an exceptional |
---|
3090 | //--- divisor in its own right |
---|
3091 | list tl; |
---|
3092 | tl[1]=intvec(iden[i][1][1],iden[i][1][2],j); |
---|
3093 | tmpiden[size(tmpiden)+1]=tl; |
---|
3094 | kill tl; |
---|
3095 | } |
---|
3096 | //--- now identify the components in the other charts w.r.t. the ones in the |
---|
3097 | //--- first chart which have already been added to the list |
---|
3098 | for(j=2;j<=size(iden[i]);j++) |
---|
3099 | { |
---|
3100 | //--- considering remaining entries for the same original divisor |
---|
3101 | if(defined(S2)){kill S2;} |
---|
3102 | def S2=re[2][iden[i][j][1]]; |
---|
3103 | setring S2; |
---|
3104 | //--- determine common parent of this ring and re[2][iden[i][1][1]] |
---|
3105 | if(defined(opath)){kill opath;} |
---|
3106 | def opath=imap(S,path); |
---|
3107 | b=1; |
---|
3108 | while(opath[1,b]==path[1,b]) |
---|
3109 | { |
---|
3110 | b++; |
---|
3111 | if((b>ncols(path))||(b>ncols(opath))) break; |
---|
3112 | } |
---|
3113 | if(defined(li1)){kill li1;} |
---|
3114 | list li1; |
---|
3115 | //--- fetch the components we have considered in re[2][iden[i][1][1]] |
---|
3116 | //--- via the resolution tree |
---|
3117 | for(k=1;k<=a;k++) |
---|
3118 | { |
---|
3119 | string tempstr="dcE["+string(eval(iden[i][1][2]))+"]["+string(k)+"][1]"; |
---|
3120 | if(defined(id1)){kill id1;} |
---|
3121 | ideal id1=fetchInTree(re,iden[i][1][1],int(leadcoef(path[1,b-1])), |
---|
3122 | iden[i][j][1],tempstr,iden0,1); |
---|
3123 | kill tempstr; |
---|
3124 | li1[k]=radical(id1); // for comparison only the geometric |
---|
3125 | // object matters |
---|
3126 | kill id1; |
---|
3127 | } |
---|
3128 | //--- compare the components we have fetched with the components in the |
---|
3129 | //--- current ring |
---|
3130 | for(k=1;k<=size(dcE[iden[i][j][2]]);k++) |
---|
3131 | { |
---|
3132 | found=0; |
---|
3133 | for(b=1;b<=size(li1);b++) |
---|
3134 | { |
---|
3135 | if((size(reduce(li1[b],std(dcE[iden[i][j][2]][k][1])))==0)&& |
---|
3136 | (size(reduce(dcE[iden[i][j][2]][k][1],std(li1[b]+BO[2])))==0)) |
---|
3137 | { |
---|
3138 | li1[b]=ideal(1); |
---|
3139 | tmpiden[size(tmpiden)-a+b][size(tmpiden[size(tmpiden)-a+b])+1]= |
---|
3140 | intvec(iden[i][j][1],iden[i][j][2],k); |
---|
3141 | found=1; |
---|
3142 | break; |
---|
3143 | } |
---|
3144 | } |
---|
3145 | if(!found) |
---|
3146 | { |
---|
3147 | if(!defined(repair)) |
---|
3148 | { |
---|
3149 | list repair; |
---|
3150 | repair[1]=list(intvec(iden[i][j][1],iden[i][j][2],k)); |
---|
3151 | } |
---|
3152 | else |
---|
3153 | { |
---|
3154 | for(c=1;c<=size(repair);c++) |
---|
3155 | { |
---|
3156 | for(d=1;d<=size(repair[c]);d++) |
---|
3157 | { |
---|
3158 | if(defined(opath)) {kill opath;} |
---|
3159 | def opath=imap(re[2][repair[c][d][1]],path); |
---|
3160 | q=0; |
---|
3161 | while(path[1,q+1]==opath[1,q+1]) |
---|
3162 | { |
---|
3163 | q++; |
---|
3164 | if((q>ncols(path)-1)||(q>ncols(opath)-1)) break; |
---|
3165 | } |
---|
3166 | q=int(leadcoef(path[1,q])); |
---|
3167 | string tempstr="dcE["+string(eval(repair[c][d][2]))+"]["+string(eval(repair[c][d][3]))+"][1]"; |
---|
3168 | if(defined(id1)){kill id1;} |
---|
3169 | ideal id1=fetchInTree(re,repair[c][d][1],q, |
---|
3170 | iden[i][j][1],tempstr,iden0,1); |
---|
3171 | kill tempstr; |
---|
3172 | //!!! sind die nicht schon radical? |
---|
3173 | id1=radical(id1); // for comparison |
---|
3174 | // only the geometric |
---|
3175 | // object matters |
---|
3176 | if((size(reduce(dcE[iden[i][j][2]][k][1],std(id1+BO[2])))==0)&& |
---|
3177 | (size(reduce(id1+BO[2],std(dcE[iden[i][j][2]][k][1])))==0)) |
---|
3178 | { |
---|
3179 | repair[c][size(repair[c])+1]=intvec(iden[i][j][1],iden[i][j][2],k); |
---|
3180 | break; |
---|
3181 | } |
---|
3182 | } |
---|
3183 | if(d<=size(repair[c])) |
---|
3184 | { |
---|
3185 | break; |
---|
3186 | } |
---|
3187 | } |
---|
3188 | if(c>size(repair)) |
---|
3189 | { |
---|
3190 | repair[size(repair)+1]=list(intvec(iden[i][j][1],iden[i][j][2],k)); |
---|
3191 | } |
---|
3192 | } |
---|
3193 | } |
---|
3194 | } |
---|
3195 | } |
---|
3196 | if(defined(repair)) |
---|
3197 | { |
---|
3198 | for(c=1;c<=size(repair);c++) |
---|
3199 | { |
---|
3200 | tmpiden[size(tmpiden)+1]=repair[c]; |
---|
3201 | } |
---|
3202 | kill repair; |
---|
3203 | } |
---|
3204 | } |
---|
3205 | setring R; |
---|
3206 | for(i=size(tmpiden);i>0;i--) |
---|
3207 | { |
---|
3208 | if(size(tmpiden[i])==0) |
---|
3209 | { |
---|
3210 | tmpiden=delete(tmpiden,i); |
---|
3211 | continue; |
---|
3212 | } |
---|
3213 | } |
---|
3214 | iden=tmpiden; // store the modified divisor list |
---|
3215 | kill tmpiden; // and clean up temporary objects |
---|
3216 | //--------------------------------------------------------------------------- |
---|
3217 | // Now we have decomposed everything into irreducible components over Q, |
---|
3218 | // but over C there might still be some reducible ones left: |
---|
3219 | // Determine the number of components over C. |
---|
3220 | //--------------------------------------------------------------------------- |
---|
3221 | n=0; |
---|
3222 | for(i=1;i<=size(iden);i++) |
---|
3223 | { |
---|
3224 | if(defined(S)) {kill S;} |
---|
3225 | def S=re[2][iden[i][1][1]]; |
---|
3226 | setring S; |
---|
3227 | divcomp[i]=ncols(dcE[iden[i][1][2]][iden[i][1][3]][4]); |
---|
3228 | // number of components of the Q-irreducible curve dcE[iden[i][1][2]] |
---|
3229 | n=n+divcomp[i]; |
---|
3230 | setring R; |
---|
3231 | } |
---|
3232 | //--------------------------------------------------------------------------- |
---|
3233 | // set up the entries Inters[i,j] , i!=j, in the intersection matrix: |
---|
3234 | // we have to compute the intersection of the exceptional divisors (over C) |
---|
3235 | // i.e. we have to work in over appropriate algebraic extension of Q. |
---|
3236 | // (1) plug the intersection matrices of the components of the same Q-irred. |
---|
3237 | // divisor into the correct position in the intersection matrix |
---|
3238 | // (2) for comparison of Ei,k and Ej,l move to a chart where both divisors |
---|
3239 | // are present, fetch the components from the very first chart containing |
---|
3240 | // the respective divisor and then compare by using intersComp |
---|
3241 | // (4) put the result into the correct position in the integer matrix Inters |
---|
3242 | //--------------------------------------------------------------------------- |
---|
3243 | //--- some initialization |
---|
3244 | int comPai,comPaj; |
---|
3245 | intvec v,w; |
---|
3246 | intmat Inters[n][n]; |
---|
3247 | //--- run through all Q-irreducible exceptional divisors |
---|
3248 | for(i=1;i<=size(iden);i++) |
---|
3249 | { |
---|
3250 | if(divcomp[i]>1) |
---|
3251 | { |
---|
3252 | //--- (1) put the intersection matrix for Ei,k with Ei,l into the correct place |
---|
3253 | for(k=1;k<=size(iden[i]);k++) |
---|
3254 | { |
---|
3255 | if(defined(tempmat)){kill tempmat;} |
---|
3256 | intmat tempmat=imap(re[2][iden[i][k][1]],dcE)[iden[i][k][2]][iden[i][k][3]][4]; |
---|
3257 | if(size(ideal(tempmat))!=0) |
---|
3258 | { |
---|
3259 | Inters[i+offset1..(i+offset1+divcomp[i]-1), |
---|
3260 | i+offset1..(i+offset1+divcomp[i]-1)]= |
---|
3261 | tempmat[1..nrows(tempmat),1..ncols(tempmat)]; |
---|
3262 | break; |
---|
3263 | } |
---|
3264 | kill tempmat; |
---|
3265 | } |
---|
3266 | } |
---|
3267 | offset2=offset1+divcomp[i]-1; |
---|
3268 | //--- set up the components over C of the i-th exceptional divisor |
---|
3269 | if(defined(S)){kill S;} |
---|
3270 | def S=re[2][iden[i][1][1]]; |
---|
3271 | setring S; |
---|
3272 | if(defined(idlisti)) {kill idlisti;} |
---|
3273 | list idlisti; |
---|
3274 | idlisti[1]=dcE[iden[i][1][2]][iden[i][1][3]][6]; |
---|
3275 | export(idlisti); |
---|
3276 | setring R; |
---|
3277 | //--- run through the remaining exceptional divisors and check whether they |
---|
3278 | //--- have a chart in common with the i-th divisor |
---|
3279 | for(j=i+1;j<=size(iden);j++) |
---|
3280 | { |
---|
3281 | kill templist; |
---|
3282 | list templist; |
---|
3283 | for(k=1;k<=size(iden[i]);k++) |
---|
3284 | { |
---|
3285 | intvec tiv2=findInIVList(1,iden[i][k][1],iden[j]); |
---|
3286 | if(size(tiv2)!=1) |
---|
3287 | { |
---|
3288 | //--- tiv2[1] is a common chart for the divisors i and j |
---|
3289 | tiv2[4..6]=iden[i][k]; |
---|
3290 | templist[size(templist)+1]=tiv2; |
---|
3291 | } |
---|
3292 | kill tiv2; |
---|
3293 | } |
---|
3294 | if(size(templist)==0) |
---|
3295 | { |
---|
3296 | //--- the two (Q-irred) divisors do not appear in any chart simultaneously |
---|
3297 | offset2=offset2+divcomp[j]-1; |
---|
3298 | j++; |
---|
3299 | continue; |
---|
3300 | } |
---|
3301 | for(k=1;k<=size(templist);k++) |
---|
3302 | { |
---|
3303 | if(defined(S)) {kill S;} |
---|
3304 | //--- set up the components over C of the j-th exceptional divisor |
---|
3305 | def S=re[2][iden[j][1][1]]; |
---|
3306 | setring S; |
---|
3307 | if(defined(idlistj)) {kill idlistj;} |
---|
3308 | list idlistj; |
---|
3309 | idlistj[1]=dcE[iden[j][1][2]][iden[j][1][3]][6]; |
---|
3310 | export(idlistj); |
---|
3311 | if(defined(opath)){kill opath;} |
---|
3312 | def opath=imap(re[2][templist[k][1]],path); |
---|
3313 | comPaj=1; |
---|
3314 | while(opath[1,comPaj]==path[1,comPaj]) |
---|
3315 | { |
---|
3316 | comPaj++; |
---|
3317 | if((comPaj>ncols(opath))||(comPaj>ncols(path))) break; |
---|
3318 | } |
---|
3319 | comPaj=int(leadcoef(path[1,comPaj-1])); |
---|
3320 | setring R; |
---|
3321 | kill S; |
---|
3322 | def S=re[2][iden[i][1][1]]; |
---|
3323 | setring S; |
---|
3324 | if(defined(opath)){kill opath;} |
---|
3325 | def opath=imap(re[2][templist[k][1]],path); |
---|
3326 | comPai=1; |
---|
3327 | while(opath[1,comPai]==path[1,comPai]) |
---|
3328 | { |
---|
3329 | comPai++; |
---|
3330 | if((comPai>ncols(opath))||(comPai>ncols(path))) break; |
---|
3331 | } |
---|
3332 | comPai=int(leadcoef(opath[1,comPai-1])); |
---|
3333 | setring R; |
---|
3334 | kill S; |
---|
3335 | def S=re[2][templist[k][1]]; |
---|
3336 | setring S; |
---|
3337 | if(defined(il)) {kill il;} |
---|
3338 | if(defined(jl)) {kill jl;} |
---|
3339 | if(defined(str1)) {kill str1;} |
---|
3340 | if(defined(str2)) {kill str2;} |
---|
3341 | string str1="idlisti"; |
---|
3342 | string str2="idlistj"; |
---|
3343 | attrib(str1,"algext",imap(re[2][iden[i][1][1]],dcE)[iden[i][1][2]][iden[i][1][3]][5]); |
---|
3344 | attrib(str2,"algext",imap(re[2][iden[j][1][1]],dcE)[iden[j][1][2]][iden[j][1][3]][5]); |
---|
3345 | list il=fetchInTree(re,iden[i][1][1],comPai, |
---|
3346 | templist[k][1],str1,iden0,1); |
---|
3347 | list jl=fetchInTree(re,iden[j][1][1],comPaj, |
---|
3348 | templist[k][1],str2,iden0,1); |
---|
3349 | list nulli=imap(re[2][iden[i][1][1]],dcE)[iden[i][1][2]][iden[i][1][3]][7]; |
---|
3350 | list nullj=imap(re[2][iden[j][1][1]],dcE)[iden[j][1][2]][iden[j][1][3]][7]; |
---|
3351 | string mpi=imap(re[2][iden[i][1][1]],dcE)[iden[i][1][2]][iden[i][1][3]][5]; |
---|
3352 | string mpj=imap(re[2][iden[j][1][1]],dcE)[iden[j][1][2]][iden[j][1][3]][5]; |
---|
3353 | if(defined(tintMat)){kill tintMat;} |
---|
3354 | intmat tintMat=intersComp(il[1],mpi,nulli,jl[1],mpj,nullj); |
---|
3355 | kill mpi; |
---|
3356 | kill mpj; |
---|
3357 | kill nulli; |
---|
3358 | kill nullj; |
---|
3359 | for(a=1;a<=divcomp[i];a++) |
---|
3360 | { |
---|
3361 | for(b=1;b<=divcomp[j];b++) |
---|
3362 | { |
---|
3363 | if(tintMat[a,b]!=0) |
---|
3364 | { |
---|
3365 | Inters[i+offset1+a-1,j+offset2+b-1]=tintMat[a,b]; |
---|
3366 | Inters[j+offset2+b-1,i+offset1+a-1]=tintMat[a,b]; |
---|
3367 | } |
---|
3368 | } |
---|
3369 | } |
---|
3370 | } |
---|
3371 | offset2=offset2+divcomp[j]-1; |
---|
3372 | } |
---|
3373 | offset1=offset1+divcomp[i]-1; |
---|
3374 | } |
---|
3375 | Inters=addSelfInter(re,Inters,iden,iden0,endiv); |
---|
3376 | intvec GenusIden; |
---|
3377 | |
---|
3378 | list tl_genus; |
---|
3379 | a=1; |
---|
3380 | for(i=1;i<=size(iden);i++) |
---|
3381 | { |
---|
3382 | tl_genus=genus_E(re,iden0,iden[i][1]); |
---|
3383 | for(j=1;j<=tl_genus[2];j++) |
---|
3384 | { |
---|
3385 | GenusIden[a]=tl_genus[1]; |
---|
3386 | a++; |
---|
3387 | } |
---|
3388 | } |
---|
3389 | |
---|
3390 | list retlist=Inters,GenusIden,iden,divcomp; |
---|
3391 | return(retlist); |
---|
3392 | } |
---|
3393 | example |
---|
3394 | {"EXAMPLE:"; |
---|
3395 | echo = 2; |
---|
3396 | ring r = 0,(x(1..3)),dp(3); |
---|
3397 | ideal J=x(3)^5+x(2)^4+x(1)^3+x(1)*x(2)*x(3); |
---|
3398 | list re=resolve(J); |
---|
3399 | list di=intersectionDiv(re); |
---|
3400 | di; |
---|
3401 | |
---|
3402 | } |
---|
3403 | ////////////////////////////////////////////////////////////////////////////// |
---|
3404 | static |
---|
3405 | proc intersComp(string str1, |
---|
3406 | string mp1, |
---|
3407 | list null1, |
---|
3408 | string str2, |
---|
3409 | string mp2, |
---|
3410 | list null2) |
---|
3411 | "Internal procedure - no help and no example available |
---|
3412 | " |
---|
3413 | { |
---|
3414 | //--- format of input |
---|
3415 | //--- str1 : ideal (over field extension 1) |
---|
3416 | //--- mp1 : minpoly of field extension 1 |
---|
3417 | //--- null1: numerical zeros of minpoly |
---|
3418 | //--- str2 : ideal (over field extension 2) |
---|
3419 | //--- mp2 : minpoly of field extension 2 |
---|
3420 | //--- null2: numerical zeros of minpoly |
---|
3421 | |
---|
3422 | //--- determine intersection matrix of the C-components defined by the input |
---|
3423 | |
---|
3424 | //--------------------------------------------------------------------------- |
---|
3425 | // Initialization |
---|
3426 | //--------------------------------------------------------------------------- |
---|
3427 | int ii,jj,same; |
---|
3428 | def R=basering; |
---|
3429 | intmat InterMat[size(null1)][size(null2)]; |
---|
3430 | ring ringst=0,(t,s),dp; |
---|
3431 | //--------------------------------------------------------------------------- |
---|
3432 | // Add new variables s and t and compare the minpolys and ideals |
---|
3433 | // to find out whether they are identical |
---|
3434 | //--------------------------------------------------------------------------- |
---|
3435 | def S=R+ringst; |
---|
3436 | setring S; |
---|
3437 | if((mp1==mp2)&&(str1==str2)) |
---|
3438 | { |
---|
3439 | same=1; |
---|
3440 | } |
---|
3441 | //--- define first Q-component/C-components, substitute t by s |
---|
3442 | string tempstr="ideal id1="+str1+";"; |
---|
3443 | execute(tempstr); |
---|
3444 | execute(mp1); |
---|
3445 | id1=subst(id1,t,s); |
---|
3446 | poly q=subst(p,t,s); |
---|
3447 | kill p; |
---|
3448 | //--- define second Q-component/C-components |
---|
3449 | tempstr="ideal id2="+str2+";"; |
---|
3450 | execute(tempstr); |
---|
3451 | execute(mp2); |
---|
3452 | //--- do the intersection |
---|
3453 | ideal interId=id1+id2+ideal(p)+ideal(q); |
---|
3454 | if(same) |
---|
3455 | { |
---|
3456 | interId=quotient(interId,t-s); |
---|
3457 | } |
---|
3458 | interId=std(interId); |
---|
3459 | //--- refine the comparison by passing to each of the numerical zeros |
---|
3460 | //--- of the two minpolys |
---|
3461 | ideal stid=nselect(interId,1,nvars(R)); |
---|
3462 | ring compl_st=complex,(s,t),dp; |
---|
3463 | def stid=imap(S,stid); |
---|
3464 | ideal tempid,tempid2; |
---|
3465 | for(ii=1;ii<=size(null1);ii++) |
---|
3466 | { |
---|
3467 | tempstr="number numi="+null1[ii]+";"; |
---|
3468 | execute(tempstr); |
---|
3469 | tempid=subst(stid,s,numi); |
---|
3470 | kill numi; |
---|
3471 | for(jj=1;jj<=size(null2);jj++) |
---|
3472 | { |
---|
3473 | tempstr="number numj="+null2[jj]+";"; |
---|
3474 | execute(tempstr); |
---|
3475 | tempid2=subst(tempid,t,numj); |
---|
3476 | kill numj; |
---|
3477 | if(size(tempid2)==0) |
---|
3478 | { |
---|
3479 | InterMat[ii,jj]=1; |
---|
3480 | } |
---|
3481 | } |
---|
3482 | } |
---|
3483 | //--- sanity check; as both Q-components were Q-irreducible, |
---|
3484 | //--- summation over all entries of a single row must lead to the same |
---|
3485 | //--- result, no matter which row is chosen |
---|
3486 | //--- dito for the columns |
---|
3487 | int cou,cou1; |
---|
3488 | for(ii=1;ii<=ncols(InterMat);ii++) |
---|
3489 | { |
---|
3490 | cou=0; |
---|
3491 | for(jj=1;jj<=nrows(InterMat);jj++) |
---|
3492 | { |
---|
3493 | cou=cou+InterMat[jj,ii]; |
---|
3494 | } |
---|
3495 | if(ii==1){cou1=cou;} |
---|
3496 | if(cou1!=cou){ERROR("intersComp:matrix has wrong entries");} |
---|
3497 | } |
---|
3498 | for(ii=1;ii<=nrows(InterMat);ii++) |
---|
3499 | { |
---|
3500 | cou=0; |
---|
3501 | for(jj=1;jj<=ncols(InterMat);jj++) |
---|
3502 | { |
---|
3503 | cou=cou+InterMat[ii,jj]; |
---|
3504 | } |
---|
3505 | if(ii==1){cou1=cou;} |
---|
3506 | if(cou1!=cou){ERROR("intersComp:matrix has wrong entries");} |
---|
3507 | } |
---|
3508 | return(InterMat); |
---|
3509 | } |
---|
3510 | ///////////////////////////////////////////////////////////////////////////// |
---|
3511 | static |
---|
3512 | proc addSelfInter(list re,intmat Inters,list iden,list iden0,intvec endiv) |
---|
3513 | "Internal procedure - no help and no example available |
---|
3514 | " |
---|
3515 | { |
---|
3516 | //--------------------------------------------------------------------------- |
---|
3517 | // Initialization |
---|
3518 | //--------------------------------------------------------------------------- |
---|
3519 | def R=basering; |
---|
3520 | int i,j,k,l,a,b; |
---|
3521 | int n=size(iden); |
---|
3522 | intvec v,w; |
---|
3523 | list satlist; |
---|
3524 | def T=re[2][1]; |
---|
3525 | setring T; |
---|
3526 | poly p; |
---|
3527 | p=var(1); //any linear form will do, |
---|
3528 | //but this one is most convenient |
---|
3529 | ideal F=ideal(p); |
---|
3530 | //---------------------------------------------------------------------------- |
---|
3531 | // lift linear form to every end ring, determine the multiplicity of |
---|
3532 | // the exceptional divisors and store it in Flist |
---|
3533 | //---------------------------------------------------------------------------- |
---|
3534 | list templist; |
---|
3535 | intvec tiv; |
---|
3536 | for(i=1;i<=size(endiv);i++) |
---|
3537 | { |
---|
3538 | if(endiv[i]==1) |
---|
3539 | { |
---|
3540 | kill v; |
---|
3541 | intvec v; |
---|
3542 | a=0; |
---|
3543 | if(defined(S)) {kill S;} |
---|
3544 | def S=re[2][i]; |
---|
3545 | setring S; |
---|
3546 | map resi=T,BO[5]; |
---|
3547 | ideal F=resi(F)+BO[2]; |
---|
3548 | ideal Ftemp=F; |
---|
3549 | list Flist; |
---|
3550 | if(defined(satlist)){kill satlist;} |
---|
3551 | list satlist; |
---|
3552 | for(a=1;a<=size(dcE);a++) |
---|
3553 | { |
---|
3554 | for(b=1;b<=size(dcE[a]);b++) |
---|
3555 | { |
---|
3556 | Ftemp=sat(Ftemp,dcE[a][b][1])[1]; |
---|
3557 | } |
---|
3558 | } |
---|
3559 | F=sat(F,Ftemp)[1]; |
---|
3560 | Flist[1]=Ftemp; |
---|
3561 | Ftemp=1; |
---|
3562 | list pr=primdecGTZ(F); |
---|
3563 | v[size(pr)]=0; |
---|
3564 | for(j=1;j<=size(pr);j++) |
---|
3565 | { |
---|
3566 | for(a=1;a<=size(dcE);a++) |
---|
3567 | { |
---|
3568 | if(j==1) |
---|
3569 | { |
---|
3570 | kill tiv; |
---|
3571 | intvec tiv; |
---|
3572 | tiv[size(dcE[a])]=0; |
---|
3573 | templist[a]=tiv; |
---|
3574 | if(v[j]==1) |
---|
3575 | { |
---|
3576 | a++; |
---|
3577 | continue; |
---|
3578 | } |
---|
3579 | } |
---|
3580 | if(dcE[a][1][2]==0) |
---|
3581 | { |
---|
3582 | a++; |
---|
3583 | continue; |
---|
3584 | } |
---|
3585 | for(b=1;b<=size(dcE[a]);b++) |
---|
3586 | { |
---|
3587 | if((size(reduce(dcE[a][b][1],std(pr[j][2])))==0)&& |
---|
3588 | (size(reduce(pr[j][2],std(dcE[a][b][1])))==0)) |
---|
3589 | { |
---|
3590 | templist[a][b]=Vielfachheit(pr[j][1],pr[j][2]); |
---|
3591 | v[j]=1; |
---|
3592 | break; |
---|
3593 | } |
---|
3594 | } |
---|
3595 | if((v[j]==1)&&(j>1)) break; |
---|
3596 | } |
---|
3597 | } |
---|
3598 | kill v; |
---|
3599 | intvec v; |
---|
3600 | Flist[2]=templist; |
---|
3601 | } |
---|
3602 | } |
---|
3603 | //----------------------------------------------------------------------------- |
---|
3604 | // Now set up all the data: |
---|
3605 | // 1. run through all exceptional divisors in iden and determine the |
---|
3606 | // coefficients c_i of the divisor of F. ===> civ |
---|
3607 | // 2. determine the intersection locus of F^bar and the Ei and from this data |
---|
3608 | // the F^bar.Ei . ===> intF |
---|
3609 | //----------------------------------------------------------------------------- |
---|
3610 | intvec civ; |
---|
3611 | intvec intF; |
---|
3612 | intF[ncols(Inters)]=0; |
---|
3613 | int offset,comPa,ncomp,vd; |
---|
3614 | for(i=1;i<=size(iden);i++) |
---|
3615 | { |
---|
3616 | ncomp=0; |
---|
3617 | for(j=1;j<=size(iden[i]);j++) |
---|
3618 | { |
---|
3619 | if(defined(S)) {kill S;} |
---|
3620 | def S=re[2][iden[i][j][1]]; |
---|
3621 | setring S; |
---|
3622 | if((size(civ)<i+offset+1)&& |
---|
3623 | (((Flist[2][iden[i][j][2]][iden[i][j][3]])!=0)||(j==size(iden[i])))) |
---|
3624 | { |
---|
3625 | ncomp=ncols(dcE[iden[i][j][2]][iden[i][j][3]][4]); |
---|
3626 | for(k=1;k<=ncomp;k++) |
---|
3627 | { |
---|
3628 | civ[i+offset+k]=Flist[2][iden[i][j][2]][iden[i][j][3]]; |
---|
3629 | if(deg(std(slocus(dcE[iden[i][j][2]][iden[i][j][3]][1]))[1])>0) |
---|
3630 | { |
---|
3631 | civ[i+offset+k]=civ[i+k]; |
---|
3632 | } |
---|
3633 | } |
---|
3634 | } |
---|
3635 | if(defined(interId)) {kill interId;} |
---|
3636 | ideal interId=dcE[iden[i][j][2]][iden[i][j][3]][1]+Flist[1]; |
---|
3637 | if(defined(interList)) {kill interList;} |
---|
3638 | list interList; |
---|
3639 | interList[1]=string(interId); |
---|
3640 | interList[2]=ideal(0); |
---|
3641 | export(interList); |
---|
3642 | if(defined(doneId)) {kill doneId;} |
---|
3643 | if(defined(tempId)) {kill tempId;} |
---|
3644 | ideal doneId=ideal(1); |
---|
3645 | if(defined(dl)) {kill dl;} |
---|
3646 | list dl; |
---|
3647 | for(k=1;k<j;k++) |
---|
3648 | { |
---|
3649 | if(defined(St)) {kill St;} |
---|
3650 | def St=re[2][iden[i][k][1]]; |
---|
3651 | setring St; |
---|
3652 | if(defined(str)){kill str;} |
---|
3653 | string str="interId="+interList[1]+";"; |
---|
3654 | execute(str); |
---|
3655 | if(deg(std(interId)[1])==0) |
---|
3656 | { |
---|
3657 | setring S; |
---|
3658 | k++; |
---|
3659 | continue; |
---|
3660 | } |
---|
3661 | setring S; |
---|
3662 | if(defined(opath)) {kill opath;} |
---|
3663 | def opath=imap(re[2][iden[i][k][1]],path); |
---|
3664 | comPa=1; |
---|
3665 | while(opath[1,comPa]==path[1,comPa]) |
---|
3666 | { |
---|
3667 | comPa++; |
---|
3668 | if((comPa>ncols(path))||(comPa>ncols(opath))) break; |
---|
3669 | } |
---|
3670 | comPa=int(leadcoef(path[1,comPa-1])); |
---|
3671 | if(defined(str)) {kill str;} |
---|
3672 | string str="interList"; |
---|
3673 | attrib(str,"algext","poly p=t-1;"); |
---|
3674 | dl=fetchInTree(re,iden[i][k][1],comPa,iden[i][j][1],str,iden0,1); |
---|
3675 | if(defined(tempId)){kill tempId;} |
---|
3676 | str="ideal tempId="+dl[1]+";"; |
---|
3677 | execute(str); |
---|
3678 | doneId=intersect(doneId,tempId); |
---|
3679 | str="interId="+interList[1]+";"; |
---|
3680 | execute(str); |
---|
3681 | interId=sat(interId,doneId)[1]; |
---|
3682 | interList[1]=string(interId); |
---|
3683 | } |
---|
3684 | interId=std(interId); |
---|
3685 | if(dim(interId)>0) |
---|
3686 | { |
---|
3687 | "oops, intersection not a set of points"; |
---|
3688 | ~; |
---|
3689 | } |
---|
3690 | vd=vdim(interId); |
---|
3691 | if(vd>0) |
---|
3692 | { |
---|
3693 | for(k=i+offset;k<=i+offset+ncomp-1;k++) |
---|
3694 | { |
---|
3695 | intF[k]=intF[k]+(vd/ncomp); |
---|
3696 | } |
---|
3697 | } |
---|
3698 | } |
---|
3699 | offset=size(civ)-i-1; |
---|
3700 | } |
---|
3701 | if(defined(tiv)){kill tiv;} |
---|
3702 | intvec tiv=civ[2..size(civ)]; |
---|
3703 | civ=tiv; |
---|
3704 | kill tiv; |
---|
3705 | //----------------------------------------------------------------------------- |
---|
3706 | // Using the F_total= sum c_i Ei + F^bar, the intersection matrix Inters and |
---|
3707 | // the f^bar.Ei, determine the selfintersection numbers of the Ei from the |
---|
3708 | // equation F_total.Ei=0 and store it in the diagonal of Inters. |
---|
3709 | //----------------------------------------------------------------------------- |
---|
3710 | intvec diag=Inters*civ+intF; |
---|
3711 | for(i=1;i<=size(diag);i++) |
---|
3712 | { |
---|
3713 | Inters[i,i]=-diag[i]/civ[i]; |
---|
3714 | } |
---|
3715 | return(Inters); |
---|
3716 | } |
---|
3717 | ////////////////////////////////////////////////////////////////////////////// |
---|
3718 | static |
---|
3719 | proc invSort(list re, list #) |
---|
3720 | "Internal procedure - no help and no example available |
---|
3721 | " |
---|
3722 | { |
---|
3723 | int i,j,k,markier,EZeiger,offset; |
---|
3724 | intvec v,e; |
---|
3725 | intvec deleted; |
---|
3726 | if(size(#)>0) |
---|
3727 | { |
---|
3728 | deleted=#[1]; |
---|
3729 | } |
---|
3730 | else |
---|
3731 | { |
---|
3732 | deleted[size(re[2])]=0; |
---|
3733 | } |
---|
3734 | list LE,HI; |
---|
3735 | def R=basering; |
---|
3736 | //---------------------------------------------------------------------------- |
---|
3737 | // Go through all rings |
---|
3738 | //---------------------------------------------------------------------------- |
---|
3739 | for(i=1;i<=size(re[2]);i++) |
---|
3740 | { |
---|
3741 | if(deleted[i]){i++;continue} |
---|
3742 | def S=re[2][i]; |
---|
3743 | setring S; |
---|
3744 | //---------------------------------------------------------------------------- |
---|
3745 | // Determine Invariant |
---|
3746 | //---------------------------------------------------------------------------- |
---|
3747 | if((size(BO[3])==size(BO[9]))||(size(BO[3])==size(BO[9])+1)) |
---|
3748 | { |
---|
3749 | if(defined(merk2)){kill merk2;} |
---|
3750 | intvec merk2; |
---|
3751 | EZeiger=0; |
---|
3752 | for(j=1;j<=size(BO[9]);j++) |
---|
3753 | { |
---|
3754 | offset=0; |
---|
3755 | if(BO[7][j]==-1) |
---|
3756 | { |
---|
3757 | BO[7][j]=size(BO[4])-EZeiger; |
---|
3758 | } |
---|
3759 | for(k=EZeiger+1;(k<=EZeiger+BO[7][j])&&(k<=size(BO[4]));k++) |
---|
3760 | { |
---|
3761 | if(BO[6][k]==2) |
---|
3762 | { |
---|
3763 | offset++; |
---|
3764 | } |
---|
3765 | } |
---|
3766 | EZeiger=EZeiger+BO[7][1]; |
---|
3767 | merk2[3*j-2]=BO[3][j]; |
---|
3768 | merk2[3*j-1]=BO[9][j]-offset; |
---|
3769 | if(size(invSat[2])>j) |
---|
3770 | { |
---|
3771 | merk2[3*j]=-invSat[2][j]; |
---|
3772 | } |
---|
3773 | else |
---|
3774 | { |
---|
3775 | if(j<size(BO[9])) |
---|
3776 | { |
---|
3777 | "!!!!!problem with invSat";~; |
---|
3778 | } |
---|
3779 | } |
---|
3780 | } |
---|
3781 | if((size(BO[3])>size(BO[9]))) |
---|
3782 | { |
---|
3783 | merk2[size(merk2)+1]=BO[3][size(BO[3])]; |
---|
3784 | } |
---|
3785 | if((size(merk2)%3)==0) |
---|
3786 | { |
---|
3787 | intvec tintvec=merk2[1..size(merk2)-1]; |
---|
3788 | merk2=tintvec; |
---|
3789 | kill tintvec; |
---|
3790 | } |
---|
3791 | } |
---|
3792 | else |
---|
3793 | { |
---|
3794 | ERROR("This situation should not occur, please send the example |
---|
3795 | to the authors."); |
---|
3796 | } |
---|
3797 | //---------------------------------------------------------------------------- |
---|
3798 | // Save invariant describing current center as an object in this ring |
---|
3799 | // We also store information on the intersection with the center and the |
---|
3800 | // exceptional divisors |
---|
3801 | //---------------------------------------------------------------------------- |
---|
3802 | cent=std(cent); |
---|
3803 | kill e; |
---|
3804 | intvec e; |
---|
3805 | for(j=1;j<=size(BO[4]);j++) |
---|
3806 | { |
---|
3807 | if(size(reduce(BO[4][j],std(cent+BO[1])))==0) |
---|
3808 | { |
---|
3809 | e[j]=1; |
---|
3810 | } |
---|
3811 | else |
---|
3812 | { |
---|
3813 | e[j]=0; |
---|
3814 | } |
---|
3815 | } |
---|
3816 | if(size(ideal(merk2))==0) |
---|
3817 | { |
---|
3818 | markier=1; |
---|
3819 | } |
---|
3820 | if((size(merk2)%3==0)&&(merk2[size(merk2)]==0)) |
---|
3821 | { |
---|
3822 | intvec blabla=merk2[1..size(merk2)-1]; |
---|
3823 | merk2=blabla; |
---|
3824 | kill blabla; |
---|
3825 | } |
---|
3826 | if(defined(invCenter)){kill invCenter;} |
---|
3827 | list invCenter=cent,merk2,e; |
---|
3828 | export invCenter; |
---|
3829 | //---------------------------------------------------------------------------- |
---|
3830 | // Insert it into correct place in the list |
---|
3831 | //---------------------------------------------------------------------------- |
---|
3832 | if(i==1) |
---|
3833 | { |
---|
3834 | if(!markier) |
---|
3835 | { |
---|
3836 | HI=intvec(merk2[1]+1),intvec(1); |
---|
3837 | } |
---|
3838 | else |
---|
3839 | { |
---|
3840 | HI=intvec(778),intvec(1); // some really large integer |
---|
3841 | // will be changed at the end!!! |
---|
3842 | } |
---|
3843 | LE[1]=HI; |
---|
3844 | i++; |
---|
3845 | setring R; |
---|
3846 | kill S; |
---|
3847 | continue; |
---|
3848 | } |
---|
3849 | if(markier==1) |
---|
3850 | { |
---|
3851 | if(i==2) |
---|
3852 | { |
---|
3853 | HI=intvec(777),intvec(2); // same really large integer-1 |
---|
3854 | LE[2]=HI; |
---|
3855 | i++; |
---|
3856 | setring R; |
---|
3857 | kill S; |
---|
3858 | continue; |
---|
3859 | } |
---|
3860 | else |
---|
3861 | { |
---|
3862 | if(ncols(path)==2) |
---|
3863 | { |
---|
3864 | LE[2][2][size(LE[2][2])+1]=i; |
---|
3865 | i++; |
---|
3866 | setring R; |
---|
3867 | kill S; |
---|
3868 | continue; |
---|
3869 | } |
---|
3870 | else |
---|
3871 | { |
---|
3872 | markier=0; |
---|
3873 | } |
---|
3874 | } |
---|
3875 | } |
---|
3876 | j=1; |
---|
3877 | def SOld=re[2][int(leadcoef(path[1,ncols(path)]))]; |
---|
3878 | setring SOld; |
---|
3879 | merk2=invCenter[2]; |
---|
3880 | setring S; |
---|
3881 | kill SOld; |
---|
3882 | while(merk2<LE[j][1]) |
---|
3883 | { |
---|
3884 | j++; |
---|
3885 | if(j>size(LE)) break; |
---|
3886 | } |
---|
3887 | HI=merk2,intvec(i); |
---|
3888 | if(j<=size(LE)) |
---|
3889 | { |
---|
3890 | if(merk2>LE[j][1]) |
---|
3891 | { |
---|
3892 | LE=insert(LE,HI,j-1); |
---|
3893 | } |
---|
3894 | else |
---|
3895 | { |
---|
3896 | while((merk2==LE[j][1])&&(size(merk2)<size(LE[j][1]))) |
---|
3897 | { |
---|
3898 | j++; |
---|
3899 | if(j>size(LE)) break; |
---|
3900 | } |
---|
3901 | if(j<=size(LE)) |
---|
3902 | { |
---|
3903 | if((merk2!=LE[j][1])||(size(merk2)!=size(LE[j][1]))) |
---|
3904 | { |
---|
3905 | LE=insert(LE,HI,j-1); |
---|
3906 | } |
---|
3907 | else |
---|
3908 | { |
---|
3909 | LE[j][2][size(LE[j][2])+1]=i; |
---|
3910 | } |
---|
3911 | } |
---|
3912 | else |
---|
3913 | { |
---|
3914 | LE[size(LE)+1]=HI; |
---|
3915 | } |
---|
3916 | } |
---|
3917 | } |
---|
3918 | else |
---|
3919 | { |
---|
3920 | LE[size(LE)+1]=HI; |
---|
3921 | } |
---|
3922 | setring R; |
---|
3923 | kill S; |
---|
3924 | } |
---|
3925 | if((LE[1][1]==intvec(778)) && (size(LE)>2)) |
---|
3926 | { |
---|
3927 | LE[1][1]=intvec(LE[3][1][1]+2); // by now we know what 'sufficiently |
---|
3928 | LE[2][1]=intvec(LE[3][1][1]+1); // large' is |
---|
3929 | } |
---|
3930 | return(LE); |
---|
3931 | } |
---|
3932 | example |
---|
3933 | {"EXAMPLE:"; |
---|
3934 | echo = 2; |
---|
3935 | ring r = 0,(x(1..3)),dp(3); |
---|
3936 | ideal J=x(1)^3-x(1)*x(2)^3+x(3)^2; |
---|
3937 | list re=resolve(J,1); |
---|
3938 | list di=invSort(re); |
---|
3939 | di; |
---|
3940 | } |
---|
3941 | ///////////////////////////////////////////////////////////////////////////// |
---|
3942 | static |
---|
3943 | proc addToRE(intvec v,int x,list RE) |
---|
3944 | "Internal procedure - no help and no example available |
---|
3945 | " |
---|
3946 | { |
---|
3947 | //--- auxilliary procedure for collectDiv, |
---|
3948 | //--- inserting an entry at the correct place |
---|
3949 | int i=1; |
---|
3950 | while(i<=size(RE)) |
---|
3951 | { |
---|
3952 | if(v==RE[i][1]) |
---|
3953 | { |
---|
3954 | RE[i][2][size(RE[i][2])+1]=x; |
---|
3955 | return(RE); |
---|
3956 | } |
---|
3957 | if(v>RE[i][1]) |
---|
3958 | { |
---|
3959 | list templist=v,intvec(x); |
---|
3960 | RE=insert(RE,templist,i-1); |
---|
3961 | return(RE); |
---|
3962 | } |
---|
3963 | i++; |
---|
3964 | } |
---|
3965 | list templist=v,intvec(x); |
---|
3966 | RE=insert(RE,templist,size(RE)); |
---|
3967 | return(RE); |
---|
3968 | } |
---|
3969 | //////////////////////////////////////////////////////////////////////////// |
---|
3970 | |
---|
3971 | proc collectDiv(list re,list #) |
---|
3972 | "USAGE: collectDiv(L); |
---|
3973 | @* L = list of rings |
---|
3974 | ASSUME: L is output of resolution of singularities |
---|
3975 | COMPUTE: list representing the identification of the exceptional divisors |
---|
3976 | in the various charts |
---|
3977 | RETURN: list l, where |
---|
3978 | l[1]: intmat, entry k in position i,j implies BO[4][j] of chart i |
---|
3979 | is divisor k (if k!=0) |
---|
3980 | if k==0, no divisor corresponding to i,j |
---|
3981 | l[2]: list ll, where each entry of ll is a list of intvecs |
---|
3982 | entry i,j in list ll[k] implies BO[4][j] of chart i |
---|
3983 | is divisor k |
---|
3984 | l[3]: list L |
---|
3985 | EXAMPLE: example collectDiv; shows an example |
---|
3986 | " |
---|
3987 | { |
---|
3988 | //------------------------------------------------------------------------ |
---|
3989 | // Initialization |
---|
3990 | //------------------------------------------------------------------------ |
---|
3991 | int i,j,k,l,m,maxk,maxj,mPa,oPa,interC,pa,ignoreL,iTotal; |
---|
3992 | int mLast,oLast=1,1; |
---|
3993 | intvec deleted; |
---|
3994 | //--- sort the rings by the invariant which controlled the last of the |
---|
3995 | //--- exceptional divisors |
---|
3996 | if(size(#)>0) |
---|
3997 | { |
---|
3998 | deleted=#[1]; |
---|
3999 | } |
---|
4000 | else |
---|
4001 | { |
---|
4002 | deleted[size(re[2])]=0; |
---|
4003 | } |
---|
4004 | list LE=invSort(re,deleted); |
---|
4005 | list LEtotal=LE; |
---|
4006 | intmat M[size(re[2])][size(re[2])]; |
---|
4007 | intvec invar,tempiv; |
---|
4008 | def R=basering; |
---|
4009 | list divList; |
---|
4010 | list RE,SE; |
---|
4011 | intvec myEi,otherEi,tempe; |
---|
4012 | int co=2; |
---|
4013 | |
---|
4014 | while(size(LE)>0) |
---|
4015 | { |
---|
4016 | //------------------------------------------------------------------------ |
---|
4017 | // Run through the sorted list LE whose entries are lists containing |
---|
4018 | // the invariant and the numbers of all rings corresponding to it |
---|
4019 | //------------------------------------------------------------------------ |
---|
4020 | for(i=co;i<=size(LE);i++) |
---|
4021 | { |
---|
4022 | //--- i==1 in first iteration: |
---|
4023 | //--- the original ring which did not arise from a blow-up |
---|
4024 | //--- hence there are no exceptional divisors to be identified there ; |
---|
4025 | |
---|
4026 | //------------------------------------------------------------------------ |
---|
4027 | // For each fixed value of the invariant, run through all corresponding |
---|
4028 | // rings |
---|
4029 | //------------------------------------------------------------------------ |
---|
4030 | for(l=1;l<=size(LE[i][2]);l++) |
---|
4031 | { |
---|
4032 | if(defined(S)){kill S;} |
---|
4033 | def S=re[2][LE[i][2][l]]; |
---|
4034 | setring S; |
---|
4035 | if(size(BO[4])>maxj){maxj=size(BO[4]);} |
---|
4036 | //--- all exceptional divisors, except the last one, were previously |
---|
4037 | //--- identified - hence we can simply inherit the data from the parent ring |
---|
4038 | for(j=1;j<size(BO[4]);j++) |
---|
4039 | { |
---|
4040 | if(deg(std(BO[4][j])[1])>0) |
---|
4041 | { |
---|
4042 | k=int(leadcoef(path[1,ncols(path)])); |
---|
4043 | k=M[k,j]; |
---|
4044 | if(k==0) |
---|
4045 | { |
---|
4046 | RE=addToRE(LE[i][1],LE[i][2][l],RE); |
---|
4047 | ignoreL=1; |
---|
4048 | break; |
---|
4049 | } |
---|
4050 | M[LE[i][2][l],j]=k; |
---|
4051 | tempiv=LE[i][2][l],j; |
---|
4052 | divList[k][size(divList[k])+1]=tempiv; |
---|
4053 | } |
---|
4054 | } |
---|
4055 | if(ignoreL){ignoreL=0;l++;continue;} |
---|
4056 | //---------------------------------------------------------------------------- |
---|
4057 | // In the remaining part of the procedure, the identification of the last |
---|
4058 | // exceptional divisor takes place. |
---|
4059 | // Step 1: check whether there is a previously considered ring with the |
---|
4060 | // same parent; if this is the case, we can again inherit the data |
---|
4061 | // Step 1':check whether the parent had a stored center which it then used |
---|
4062 | // in this case, we are dealing with an additional component of this |
---|
4063 | // divisor: store it in the integer otherComp |
---|
4064 | // Step 2: if no appropriate ring was found in step 1, we check whether |
---|
4065 | // there is a previously considered ring, in the parent of which |
---|
4066 | // the center intersects the same exceptional divisors as the center |
---|
4067 | // in our parent. |
---|
4068 | // if such a ring does not exist: new exceptional divisor |
---|
4069 | // if it exists: see below |
---|
4070 | //---------------------------------------------------------------------------- |
---|
4071 | if(path[1,ncols(path)-1]==0) |
---|
4072 | { |
---|
4073 | //--- current ring originated from very first blow-up |
---|
4074 | //--- hence exceptional divisor is the first one |
---|
4075 | M[LE[i][2][l],1]=1; |
---|
4076 | if(size(divList)>0) |
---|
4077 | { |
---|
4078 | divList[1][size(divList[1])+1]=intvec(LE[i][2][l],j); |
---|
4079 | } |
---|
4080 | else |
---|
4081 | { |
---|
4082 | divList[1]=list(intvec(LE[i][2][l],j)); |
---|
4083 | } |
---|
4084 | l++; |
---|
4085 | continue; |
---|
4086 | } |
---|
4087 | if(l==1) |
---|
4088 | { |
---|
4089 | list TE=addToRE(LE[i][1],1,SE); |
---|
4090 | if(size(TE)!=size(SE)) |
---|
4091 | { |
---|
4092 | //--- new value of invariant hence new exceptional divisor |
---|
4093 | SE=TE; |
---|
4094 | divList[size(divList)+1]=list(intvec(LE[i][2][l],j)); |
---|
4095 | M[LE[i][2][l],j]=size(divList); |
---|
4096 | } |
---|
4097 | kill TE; |
---|
4098 | } |
---|
4099 | for(k=1;k<=size(LEtotal);k++) |
---|
4100 | { |
---|
4101 | if(LE[i][1]==LEtotal[k][1]) |
---|
4102 | { |
---|
4103 | iTotal=k; |
---|
4104 | break; |
---|
4105 | } |
---|
4106 | } |
---|
4107 | //--- Step 1 |
---|
4108 | k=1; |
---|
4109 | while(LEtotal[iTotal][2][k]<LE[i][2][l]) |
---|
4110 | { |
---|
4111 | if(defined(tempPath)){kill tempPath;} |
---|
4112 | def tempPath=imap(re[2][LEtotal[iTotal][2][k]],path); |
---|
4113 | if(tempPath[1,ncols(tempPath)]==path[1,ncols(path)]) |
---|
4114 | { |
---|
4115 | //--- Same parent, hence we inherit our data |
---|
4116 | m=size(imap(re[2][LEtotal[iTotal][2][k]],BO)[4]); |
---|
4117 | m=M[LEtotal[iTotal][2][k],m]; |
---|
4118 | if(m==0) |
---|
4119 | { |
---|
4120 | RE=addToRE(LE[i][1],LE[i][2][l],RE); |
---|
4121 | ignoreL=1; |
---|
4122 | break; |
---|
4123 | } |
---|
4124 | M[LE[i][2][l],j]=m; |
---|
4125 | tempiv=LE[i][2][l],j; |
---|
4126 | divList[m][size(divList[m])+1]=tempiv; |
---|
4127 | break; |
---|
4128 | } |
---|
4129 | k++; |
---|
4130 | if(k>size(LEtotal[iTotal][2])) {break;} |
---|
4131 | } |
---|
4132 | if(ignoreL){ignoreL=0;l++;continue;} |
---|
4133 | //--- Step 1', if necessary |
---|
4134 | if(M[LE[i][2][l],j]==0) |
---|
4135 | { |
---|
4136 | int savedCent; |
---|
4137 | def SPa1=re[2][int(leadcoef(path[1,ncols(path)]))]; |
---|
4138 | // parent ring |
---|
4139 | setring SPa1; |
---|
4140 | if(size(BO)>9) |
---|
4141 | { |
---|
4142 | if(size(BO[10])>0) |
---|
4143 | { |
---|
4144 | savedCent=1; |
---|
4145 | } |
---|
4146 | } |
---|
4147 | if(!savedCent) |
---|
4148 | { |
---|
4149 | def SPa2=re[2][int(leadcoef(path[1,ncols(path)]))]; |
---|
4150 | map lMa=SPa2,lastMap; |
---|
4151 | // map leading from grandparent to parent |
---|
4152 | list transBO=lMa(BO); |
---|
4153 | // actually we only need BO[10], but this is an |
---|
4154 | // object not a name |
---|
4155 | list tempsat; |
---|
4156 | if(size(transBO)>9) |
---|
4157 | { |
---|
4158 | //--- there were saved centers |
---|
4159 | while((k<=size(transBO[10])) & (savedCent==0)) |
---|
4160 | { |
---|
4161 | tempsat=sat(transBO[10][k][1],BO[4][size(BO[4])]); |
---|
4162 | if(deg(tempsat[1][1])!=0) |
---|
4163 | { |
---|
4164 | //--- saved center can be seen in this affine chart |
---|
4165 | if((size(reduce(tempsat[1],std(cent)))==0) && |
---|
4166 | (size(reduce(cent,tempsat[1]))==0)) |
---|
4167 | { |
---|
4168 | //--- this was the saved center which was used |
---|
4169 | savedCent=1; |
---|
4170 | } |
---|
4171 | } |
---|
4172 | k++; |
---|
4173 | } |
---|
4174 | } |
---|
4175 | kill lMa; // clean up temporary objects |
---|
4176 | kill tempsat; |
---|
4177 | kill transBO; |
---|
4178 | } |
---|
4179 | setring S; // back to the ring which we want to consider |
---|
4180 | if(savedCent==1) |
---|
4181 | { |
---|
4182 | vector otherComp; |
---|
4183 | otherComp[M[int(leadcoef(path[1,ncols(path)])),size(BO[4])-1]] |
---|
4184 | =1; |
---|
4185 | } |
---|
4186 | kill savedCent; |
---|
4187 | if (defined(SPa2)){kill SPa2;} |
---|
4188 | kill SPa1; |
---|
4189 | } |
---|
4190 | //--- Step 2, if necessary |
---|
4191 | if(M[LE[i][2][l],j]==0) |
---|
4192 | { |
---|
4193 | //--- we are not done after step 1 and 2 |
---|
4194 | pa=int(leadcoef(path[1,ncols(path)])); // parent ring |
---|
4195 | tempe=imap(re[2][pa],invCenter)[3]; // intersection there |
---|
4196 | kill myEi; |
---|
4197 | intvec myEi; |
---|
4198 | for(k=1;k<=size(tempe);k++) |
---|
4199 | { |
---|
4200 | if(tempe[k]==1) |
---|
4201 | { |
---|
4202 | //--- center meets this exceptional divisor |
---|
4203 | myEi[size(myEi)+1]=M[pa,k]; |
---|
4204 | mLast=k; |
---|
4205 | } |
---|
4206 | } |
---|
4207 | //--- ring in which the last divisor we meet is new-born |
---|
4208 | mPa=int(leadcoef(path[1,mLast+2])); |
---|
4209 | k=1; |
---|
4210 | while(LEtotal[iTotal][2][k]<LE[i][2][l]) |
---|
4211 | { |
---|
4212 | //--- perform the same preparations for the ring we want to compare with |
---|
4213 | if(defined(tempPath)){kill tempPath;} |
---|
4214 | def tempPath=imap(re[2][LEtotal[iTotal][2][k]],path); |
---|
4215 | // its ancestors |
---|
4216 | tempe=imap(re[2][int(leadcoef(tempPath[1,ncols(tempPath)]))], |
---|
4217 | invCenter)[3]; // its intersections |
---|
4218 | kill otherEi; |
---|
4219 | intvec otherEi; |
---|
4220 | for(m=1;m<=size(tempe);m++) |
---|
4221 | { |
---|
4222 | if(tempe[m]==1) |
---|
4223 | { |
---|
4224 | //--- its center meets this exceptional divisor |
---|
4225 | otherEi[size(otherEi)+1] |
---|
4226 | =M[int(leadcoef(tempPath[1,ncols(tempPath)])),m]; |
---|
4227 | oLast=m; |
---|
4228 | } |
---|
4229 | } |
---|
4230 | if(myEi!=otherEi) |
---|
4231 | { |
---|
4232 | //--- not the same center because of intersection properties with the |
---|
4233 | //--- exceptional divisor |
---|
4234 | k++; |
---|
4235 | if(k>size(LEtotal[iTotal][2])) |
---|
4236 | { |
---|
4237 | break; |
---|
4238 | } |
---|
4239 | else |
---|
4240 | { |
---|
4241 | continue; |
---|
4242 | } |
---|
4243 | } |
---|
4244 | //---------------------------------------------------------------------------- |
---|
4245 | // Current situation: |
---|
4246 | // 1. the last exceptional divisor could not be identified by simply |
---|
4247 | // considering its parent |
---|
4248 | // 2. it could not be proved to be a new one by considering its intersections |
---|
4249 | // with previous exceptional divisors |
---|
4250 | //---------------------------------------------------------------------------- |
---|
4251 | if(defined(bool1)) { kill bool1;} |
---|
4252 | int bool1= |
---|
4253 | compareE(re,LE[i][2][l],LEtotal[iTotal][2][k],divList); |
---|
4254 | if(bool1) |
---|
4255 | { |
---|
4256 | //--- found some non-empty intersection |
---|
4257 | if(bool1==1) |
---|
4258 | { |
---|
4259 | //--- it is really the same exceptional divisor |
---|
4260 | m=size(imap(re[2][LEtotal[iTotal][2][k]],BO)[4]); |
---|
4261 | m=M[LEtotal[iTotal][2][k],m]; |
---|
4262 | if(m==0) |
---|
4263 | { |
---|
4264 | RE=addToRE(LE[i][1],LE[i][2][l],RE); |
---|
4265 | ignoreL=1; |
---|
4266 | break; |
---|
4267 | } |
---|
4268 | M[LE[i][2][l],j]=m; |
---|
4269 | tempiv=LE[i][2][l],j; |
---|
4270 | divList[m][size(divList[m])+1]=tempiv; |
---|
4271 | break; |
---|
4272 | } |
---|
4273 | else |
---|
4274 | { |
---|
4275 | m=size(imap(re[2][LEtotal[iTotal][2][k]],BO)[4]); |
---|
4276 | m=M[LEtotal[iTotal][2][k],m]; |
---|
4277 | if(m!=0) |
---|
4278 | { |
---|
4279 | otherComp[m]=1; |
---|
4280 | } |
---|
4281 | } |
---|
4282 | } |
---|
4283 | k++; |
---|
4284 | if(k>size(LEtotal[iTotal][2])) |
---|
4285 | { |
---|
4286 | break; |
---|
4287 | } |
---|
4288 | } |
---|
4289 | if(ignoreL){ignoreL=0;l++;continue;} |
---|
4290 | if( M[LE[i][2][l],j]==0) |
---|
4291 | { |
---|
4292 | divList[size(divList)+1]=list(intvec(LE[i][2][l],j)); |
---|
4293 | M[LE[i][2][l],j]=size(divList); |
---|
4294 | } |
---|
4295 | } |
---|
4296 | setring R; |
---|
4297 | kill S; |
---|
4298 | } |
---|
4299 | } |
---|
4300 | LE=RE; |
---|
4301 | co=1; |
---|
4302 | kill RE; |
---|
4303 | list RE; |
---|
4304 | } |
---|
4305 | //---------------------------------------------------------------------------- |
---|
4306 | // Add the strict transform to the list of divisors at the last place |
---|
4307 | // and clean up M |
---|
4308 | //---------------------------------------------------------------------------- |
---|
4309 | //--- add strict transform |
---|
4310 | for(i=1;i<=size(re[2]);i++) |
---|
4311 | { |
---|
4312 | if(defined(S)){kill S;} |
---|
4313 | def S=re[2][i]; |
---|
4314 | setring S; |
---|
4315 | if(size(reduce(cent,std(BO[2])))==0) |
---|
4316 | { |
---|
4317 | tempiv=i,0; |
---|
4318 | RE[size(RE)+1]=tempiv; |
---|
4319 | } |
---|
4320 | setring R; |
---|
4321 | } |
---|
4322 | divList[size(divList)+1]=RE; |
---|
4323 | //--- drop trailing zero-columns of M |
---|
4324 | intvec iv0; |
---|
4325 | iv0[nrows(M)]=0; |
---|
4326 | for(i=ncols(M);i>0;i--) |
---|
4327 | { |
---|
4328 | if(intvec(M[1..nrows(M),i])!=iv0) break; |
---|
4329 | } |
---|
4330 | intmat N[nrows(M)][i]; |
---|
4331 | for(i=1;i<=ncols(N);i++) |
---|
4332 | { |
---|
4333 | N[1..nrows(M),i]=M[1..nrows(M),i]; |
---|
4334 | } |
---|
4335 | kill M; |
---|
4336 | intmat M=N; |
---|
4337 | list retlist=cleanUpDiv(re,M,divList); |
---|
4338 | return(retlist); |
---|
4339 | } |
---|
4340 | example |
---|
4341 | {"EXAMPLE:"; |
---|
4342 | echo = 2; |
---|
4343 | ring R=0,(x,y,z),dp; |
---|
4344 | ideal I=xyz+x4+y4+z4; |
---|
4345 | //we really need to blow up curves even if the generic point of |
---|
4346 | //the curve the total transform is n.c. |
---|
4347 | //this occurs here in r[2][5] |
---|
4348 | list re=resolve(I,1); |
---|
4349 | list di=collectDiv(re); |
---|
4350 | di; |
---|
4351 | } |
---|
4352 | ////////////////////////////////////////////////////////////////////////////// |
---|
4353 | static |
---|
4354 | proc cleanUpDiv(list re,intmat M,list divList) |
---|
4355 | "Internal procedure - no help and no example available |
---|
4356 | " |
---|
4357 | { |
---|
4358 | //--- It may occur that two different entries of invSort coincide on the |
---|
4359 | //--- first part up to the last entry of the shorter one. In this case |
---|
4360 | //--- exceptional divisors may appear in both entries of the invSort-list. |
---|
4361 | //--- To correct this, we now compare the final collection of Divisors |
---|
4362 | //--- for coinciding ones. |
---|
4363 | int i,j,k,a,oPa,mPa,comPa,mdim,odim; |
---|
4364 | def R=basering; |
---|
4365 | for(i=1;i<=size(divList)-2;i++) |
---|
4366 | { |
---|
4367 | if(defined(Sm)){kill Sm;} |
---|
4368 | def Sm=re[2][divList[i][1][1]]; |
---|
4369 | setring Sm; |
---|
4370 | mPa=int(leadcoef(path[1,ncols(path)])); |
---|
4371 | if(defined(SmPa)){kill SmPa;} |
---|
4372 | def SmPa=re[2][mPa]; |
---|
4373 | setring SmPa; |
---|
4374 | mdim=dim(std(BO[1]+cent)); |
---|
4375 | setring Sm; |
---|
4376 | if(mPa==1) |
---|
4377 | { |
---|
4378 | //--- very first divisor originates exactly from the first blow-up |
---|
4379 | //--- there cannot be any mistake here |
---|
4380 | i++; |
---|
4381 | continue; |
---|
4382 | } |
---|
4383 | for(j=i+1;j<=size(divList)-1;j++) |
---|
4384 | { |
---|
4385 | setring Sm; |
---|
4386 | for(k=1;k<=size(divList[j]);k++) |
---|
4387 | { |
---|
4388 | if(size(findInIVList(1,divList[j][k][1],divList[i]))>1) |
---|
4389 | { |
---|
4390 | //--- same divisor cannot appear twice in the same chart |
---|
4391 | k=-1; |
---|
4392 | break; |
---|
4393 | } |
---|
4394 | } |
---|
4395 | if(k==-1) |
---|
4396 | { |
---|
4397 | j++; |
---|
4398 | if(j>size(divList)-1) break; |
---|
4399 | continue; |
---|
4400 | } |
---|
4401 | if(defined(opath)){kill opath;} |
---|
4402 | def opath=imap(re[2][divList[j][1][1]],path); |
---|
4403 | oPa=int(leadcoef(opath[1,ncols(opath)])); |
---|
4404 | if(defined(SoPa)){kill SoPa;} |
---|
4405 | def SoPa=re[2][oPa]; |
---|
4406 | setring SoPa; |
---|
4407 | odim=dim(std(BO[1]+cent)); |
---|
4408 | setring Sm; |
---|
4409 | if(mdim!=odim) |
---|
4410 | { |
---|
4411 | //--- different dimension ==> cannot be same center |
---|
4412 | j++; |
---|
4413 | if(j>size(divList)-1) break; |
---|
4414 | continue; |
---|
4415 | } |
---|
4416 | comPa=1; |
---|
4417 | while(path[1,comPa]==opath[1,comPa]) |
---|
4418 | { |
---|
4419 | comPa++; |
---|
4420 | if((comPa>ncols(path))||(comPa>ncols(opath))) break; |
---|
4421 | } |
---|
4422 | comPa=int(leadcoef(path[1,comPa-1])); |
---|
4423 | if(defined(SPa)){kill SPa;} |
---|
4424 | def SPa=re[2][mPa]; |
---|
4425 | setring SPa; |
---|
4426 | if(defined(tempIdE)){kill tempIdE;} |
---|
4427 | ideal tempIdE=fetchInTree(re,oPa,comPa,mPa,"cent",divList); |
---|
4428 | if((size(reduce(cent,std(tempIdE)))!=0)|| |
---|
4429 | (size(reduce(tempIdE,std(cent)))!=0)) |
---|
4430 | { |
---|
4431 | //--- it is not the same divisor! |
---|
4432 | j++; |
---|
4433 | if(j>size(divList)) |
---|
4434 | { |
---|
4435 | break; |
---|
4436 | } |
---|
4437 | else |
---|
4438 | { |
---|
4439 | continue; |
---|
4440 | } |
---|
4441 | } |
---|
4442 | for(k=1;k<=size(divList[j]);k++) |
---|
4443 | { |
---|
4444 | //--- append the entries of the j-th divisor (which is actually also the i-th) |
---|
4445 | //--- to the i-th divisor |
---|
4446 | divList[i][size(divList[i])+1]=divList[j][k]; |
---|
4447 | } |
---|
4448 | divList=delete(divList,j); //kill obsolete entry from the list |
---|
4449 | for(k=1;k<=nrows(M);k++) |
---|
4450 | { |
---|
4451 | for(a=1;a<=ncols(M);a++) |
---|
4452 | { |
---|
4453 | if(M[k,a]==j) |
---|
4454 | { |
---|
4455 | //--- j-th divisor is actually the i-th one |
---|
4456 | M[k,a]=i; |
---|
4457 | } |
---|
4458 | if(M[k,a]>j) |
---|
4459 | { |
---|
4460 | //--- index j was deleted from the list ==> all subsequent indices dropped by |
---|
4461 | //--- one |
---|
4462 | M[k,a]=M[k,a]-1; |
---|
4463 | } |
---|
4464 | } |
---|
4465 | } |
---|
4466 | j--; //do not forget to consider new j-th entry |
---|
4467 | } |
---|
4468 | } |
---|
4469 | setring R; |
---|
4470 | list retlist=M,divList; |
---|
4471 | return(retlist); |
---|
4472 | } |
---|
4473 | ///////////////////////////////////////////////////////////////////////////// |
---|
4474 | static |
---|
4475 | proc findTrans(ideal Z, ideal E, list notE, list #) |
---|
4476 | "Internal procedure - no help and no example available |
---|
4477 | " |
---|
4478 | { |
---|
4479 | //---Auxilliary procedure for fetchInTree! |
---|
4480 | //---Assume E prime ideal, Z+E eqidimensional, |
---|
4481 | //---ht(E)+r=ht(Z+E). Compute P=<p[1],...,p[r]> in Z+E, and poly f, |
---|
4482 | //---such that radical(Z+E)=radical((E+P):f) |
---|
4483 | int i,j,d,e; |
---|
4484 | ideal Estd=std(E); |
---|
4485 | //!!! alternative to subsequent line: |
---|
4486 | //!!! ideal Zstd=std(radical(Z+E)); |
---|
4487 | ideal Zstd=std(Z+E); |
---|
4488 | ideal J=1; |
---|
4489 | if(size(#)>0) |
---|
4490 | { |
---|
4491 | J=#[1]; |
---|
4492 | } |
---|
4493 | if(deg(Zstd[1])==0){return(list(ideal(1),poly(1)));} |
---|
4494 | for(i=1;i<=size(notE);i++) |
---|
4495 | { |
---|
4496 | notE[i]=std(notE[i]); |
---|
4497 | } |
---|
4498 | ideal Zred=simplify(reduce(Z,Estd),2); |
---|
4499 | if(size(Zred)==0){Z,Estd;~;ERROR("Z is contained in E");} |
---|
4500 | ideal P,Q,Qstd; |
---|
4501 | Q=Estd; |
---|
4502 | attrib(Q,"isSB",1); |
---|
4503 | d=dim(Estd); |
---|
4504 | e=dim(Zstd); |
---|
4505 | for(i=1;i<=size(Zred);i++) |
---|
4506 | { |
---|
4507 | Qstd=std(Q,Zred[i]); |
---|
4508 | if(dim(Qstd)<d) |
---|
4509 | { |
---|
4510 | d=dim(Qstd); |
---|
4511 | P[size(P)+1]=Zred[i]; |
---|
4512 | Q=Qstd; |
---|
4513 | attrib(Q,"isSB",1); |
---|
4514 | if(d==e) break; |
---|
4515 | } |
---|
4516 | } |
---|
4517 | list pr=minAssGTZ(E+P); |
---|
4518 | list sr=minAssGTZ(J+P); |
---|
4519 | i=0; |
---|
4520 | Q=1; |
---|
4521 | list qr; |
---|
4522 | |
---|
4523 | while(i<size(pr)) |
---|
4524 | { |
---|
4525 | i++; |
---|
4526 | Qstd=std(pr[i]); |
---|
4527 | Zred=simplify(reduce(Zstd,Qstd),2); |
---|
4528 | if(size(Zred)==0) |
---|
4529 | { |
---|
4530 | qr[size(qr)+1]=pr[i]; |
---|
4531 | pr=delete(pr,i); |
---|
4532 | i--; |
---|
4533 | } |
---|
4534 | else |
---|
4535 | { |
---|
4536 | Q=intersect(Q,pr[i]); |
---|
4537 | } |
---|
4538 | } |
---|
4539 | i=0; |
---|
4540 | while(i<size(sr)) |
---|
4541 | { |
---|
4542 | i++; |
---|
4543 | Qstd=std(sr[i]+E); |
---|
4544 | Zred=simplify(reduce(Zstd,Qstd),2); |
---|
4545 | if((size(Zred)!=0)||(dim(Qstd)!=dim(Zstd))) |
---|
4546 | { |
---|
4547 | Q=intersect(Q,sr[i]); |
---|
4548 | } |
---|
4549 | } |
---|
4550 | poly f; |
---|
4551 | for(i=1;i<=size(Q);i++) |
---|
4552 | { |
---|
4553 | f=Q[i]; |
---|
4554 | for(e=1;e<=size(qr);e++) |
---|
4555 | { |
---|
4556 | if(reduce(f,std(qr[e]))==0){f=0;break;} |
---|
4557 | } |
---|
4558 | for(j=1;j<=size(notE);j++) |
---|
4559 | { |
---|
4560 | if(reduce(f,notE[j])==0){f=0; break;} |
---|
4561 | } |
---|
4562 | if(f!=0) break; |
---|
4563 | } |
---|
4564 | i=0; |
---|
4565 | while(f==0) |
---|
4566 | { |
---|
4567 | i++; |
---|
4568 | f=randomid(Q)[1]; |
---|
4569 | for(e=1;e<=size(qr);e++) |
---|
4570 | { |
---|
4571 | if(reduce(f,std(qr[e]))==0){f=0;break;} |
---|
4572 | } |
---|
4573 | for(j=1;j<=size(notE);j++) |
---|
4574 | { |
---|
4575 | if(reduce(f,notE[j])==0){f=0; break;} |
---|
4576 | } |
---|
4577 | if(f!=0) break; |
---|
4578 | if(i>20) |
---|
4579 | { |
---|
4580 | ~; |
---|
4581 | ERROR("findTrans:Hier ist was faul"); |
---|
4582 | } |
---|
4583 | } |
---|
4584 | |
---|
4585 | list resu=P,f; |
---|
4586 | return(resu); |
---|
4587 | } |
---|
4588 | ///////////////////////////////////////////////////////////////////////////// |
---|
4589 | static |
---|
4590 | proc compareE(list L, int m, int o, list DivL) |
---|
4591 | "Internal procedure - no help and no example available |
---|
4592 | " |
---|
4593 | { |
---|
4594 | //---------------------------------------------------------------------------- |
---|
4595 | // We want to compare the divisors BO[4][size(BO[4])] of the rings |
---|
4596 | // L[2][m] and L[2][o]. |
---|
4597 | // In the initialization step, we collect all necessary data from those |
---|
4598 | // those rings. In particular, we determine at what point (in the resolution |
---|
4599 | // history) the branches for L[2][m] and L[2][o] were separated, denoting |
---|
4600 | // the corresponding ring indices by mPa, oPa and comPa. |
---|
4601 | //---------------------------------------------------------------------------- |
---|
4602 | def R=basering; |
---|
4603 | int i,j,k,len; |
---|
4604 | |
---|
4605 | //-- find direct parents and branching point in resolution history |
---|
4606 | matrix tpm=imap(L[2][m],path); |
---|
4607 | matrix tpo=imap(L[2][o],path); |
---|
4608 | int m1,o1=int(leadcoef(tpm[1,ncols(tpm)])), |
---|
4609 | int(leadcoef(tpo[1,ncols(tpo)])); |
---|
4610 | while((i<ncols(tpo)) && (i<ncols(tpm))) |
---|
4611 | { |
---|
4612 | if(tpm[1,i+1]!=tpo[1,i+1]) break; |
---|
4613 | i++; |
---|
4614 | } |
---|
4615 | int branchpos=i; |
---|
4616 | int comPa=int(leadcoef(tpm[1,branchpos])); // last common ancestor |
---|
4617 | //---------------------------------------------------------------------------- |
---|
4618 | // simple checks to save us some work in obvious cases |
---|
4619 | //---------------------------------------------------------------------------- |
---|
4620 | if((comPa==m1)||(comPa==o1)) |
---|
4621 | { |
---|
4622 | //--- one is in the history of the other ==> they cannot give rise |
---|
4623 | //--- to the same divisor |
---|
4624 | return(0); |
---|
4625 | } |
---|
4626 | def T=L[2][o1]; |
---|
4627 | setring T; |
---|
4628 | int dimCo1=dim(std(cent+BO[1])); |
---|
4629 | def S=L[2][m1]; |
---|
4630 | setring S; |
---|
4631 | int dimCm1=dim(std(cent+BO[1])); |
---|
4632 | if(dimCm1!=dimCo1) |
---|
4633 | { |
---|
4634 | //--- centers do not have same dimension ==> they cannot give rise |
---|
4635 | //--- to the same divisor |
---|
4636 | return(0); |
---|
4637 | } |
---|
4638 | //---------------------------------------------------------------------------- |
---|
4639 | // fetch the center via the tree for comparison |
---|
4640 | //---------------------------------------------------------------------------- |
---|
4641 | if(defined(invLocus0)) {kill invLocus0;} |
---|
4642 | ideal invLocus0=fetchInTree(L,o1,comPa,m1,"cent",DivL); |
---|
4643 | // blow down from L[2][o1] to L[2][comPa] and then up to L[2][m1] |
---|
4644 | if(deg(std(invLocus0+invCenter[1]+BO[1])[1])!=0) |
---|
4645 | { |
---|
4646 | setring R; |
---|
4647 | return(int(1)); |
---|
4648 | } |
---|
4649 | if(size(BO)>9) |
---|
4650 | { |
---|
4651 | for(i=1;i<=size(BO[10]);i++) |
---|
4652 | { |
---|
4653 | if(deg(std(invLocus0+BO[10][i][1]+BO[1])[1])!=0) |
---|
4654 | { |
---|
4655 | if(dim(std(BO[10][i][1]+BO[1])) > |
---|
4656 | dim(std(invLocus0+BO[10][i][1]+BO[1]))) |
---|
4657 | { |
---|
4658 | ERROR("Internal Error: Please send this example to the authors."); |
---|
4659 | } |
---|
4660 | setring R; |
---|
4661 | return(int(2)); |
---|
4662 | } |
---|
4663 | } |
---|
4664 | } |
---|
4665 | setring R; |
---|
4666 | return(int(0)); |
---|
4667 | //---------------------------------------------------------------------------- |
---|
4668 | // Return-Values: |
---|
4669 | // TRUE (=1) if the exceptional divisors coincide, |
---|
4670 | // TRUE (=2) if the exceptional divisors originate from different |
---|
4671 | // components of the same center |
---|
4672 | // FALSE (=0) otherwise |
---|
4673 | //---------------------------------------------------------------------------- |
---|
4674 | } |
---|
4675 | ////////////////////////////////////////////////////////////////////////////// |
---|
4676 | |
---|
4677 | proc fetchInTree(list L, |
---|
4678 | int o1, |
---|
4679 | int comPa, |
---|
4680 | int m1, |
---|
4681 | string idname, |
---|
4682 | list DivL, |
---|
4683 | list #); |
---|
4684 | "Internal procedure - no help and no example available |
---|
4685 | " |
---|
4686 | { |
---|
4687 | //---------------------------------------------------------------------------- |
---|
4688 | // Initialization and Sanity Checks |
---|
4689 | //---------------------------------------------------------------------------- |
---|
4690 | int i,j,k,m,branchPos,inJ,exception; |
---|
4691 | string algext; |
---|
4692 | //--- we need to be in L[2][m1] |
---|
4693 | def R=basering; |
---|
4694 | ideal test_for_the_same_ring=-77; |
---|
4695 | def Sm1=L[2][m1]; |
---|
4696 | setring Sm1; |
---|
4697 | if(!defined(test_for_the_same_ring)) |
---|
4698 | { |
---|
4699 | //--- we are not in L[2][m1] |
---|
4700 | ERROR("basering has to coincide with L[2][m1]"); |
---|
4701 | } |
---|
4702 | else |
---|
4703 | { |
---|
4704 | //--- we are in L[2][m1] |
---|
4705 | kill test_for_the_same_ring; |
---|
4706 | } |
---|
4707 | //--- non-embedded case? |
---|
4708 | if(size(#)>0) |
---|
4709 | { |
---|
4710 | inJ=1; |
---|
4711 | } |
---|
4712 | //--- do parameter values make sense? |
---|
4713 | if(comPa<1) |
---|
4714 | { |
---|
4715 | ERROR("Common Parent should at least be the first ring!"); |
---|
4716 | } |
---|
4717 | //--- do we need to pass to an algebraic field extension of Q? |
---|
4718 | if(typeof(attrib(idname,"algext"))=="string") |
---|
4719 | { |
---|
4720 | algext=attrib(idname,"algext"); |
---|
4721 | } |
---|
4722 | //--- check wheter comPa is in the history of m1 |
---|
4723 | //--- same test for o1 can be done later on (on the fly) |
---|
4724 | if(m1==comPa) |
---|
4725 | { |
---|
4726 | j=1; |
---|
4727 | i=ncols(path)+1; |
---|
4728 | } |
---|
4729 | else |
---|
4730 | { |
---|
4731 | for(i=1;i<=ncols(path);i++) |
---|
4732 | { |
---|
4733 | if(int(leadcoef(path[1,i]))==comPa) |
---|
4734 | { |
---|
4735 | //--- comPa occurs in the history |
---|
4736 | j=1; |
---|
4737 | break; |
---|
4738 | } |
---|
4739 | } |
---|
4740 | } |
---|
4741 | branchPos=i; |
---|
4742 | if(j==0) |
---|
4743 | { |
---|
4744 | ERROR("L[2][comPa] not in history of L[2][m1]!"); |
---|
4745 | } |
---|
4746 | //---------------------------------------------------------------------------- |
---|
4747 | // Blow down ideal "idname" from L[2][o1] to L[2][comPa], where the latter |
---|
4748 | // is assumed to be the common parent of L[2][o1] and L[2][m1] |
---|
4749 | //---------------------------------------------------------------------------- |
---|
4750 | if(size(algext)>0) |
---|
4751 | { |
---|
4752 | if(defined(tstr)){kill tstr;} |
---|
4753 | string tstr="ring So1=(0,t),("+varstr(L[2][o1])+"),("+ordstr(L[2][o1])+");"; |
---|
4754 | execute(tstr); |
---|
4755 | setring So1; |
---|
4756 | execute(algext); |
---|
4757 | minpoly=leadcoef(p); |
---|
4758 | if(defined(id1)) { kill id1; } |
---|
4759 | if(defined(id2)) { kill id2; } |
---|
4760 | if(defined(idlist)) { kill idlist; } |
---|
4761 | execute("int bool2=defined("+idname+");"); |
---|
4762 | if(bool2==0) |
---|
4763 | { |
---|
4764 | execute("list ttlist=imap(L[2][o1],"+idname+");"); |
---|
4765 | } |
---|
4766 | else |
---|
4767 | { |
---|
4768 | execute("list ttlist="+idname+";"); |
---|
4769 | } |
---|
4770 | kill bool2; |
---|
4771 | def BO=imap(L[2][o1],BO); |
---|
4772 | def path=imap(L[2][o1],path); |
---|
4773 | def lastMap=imap(L[2][o1],lastMap); |
---|
4774 | ideal id2=1; |
---|
4775 | if(defined(notE)){kill notE;} |
---|
4776 | list notE; |
---|
4777 | intvec nE; |
---|
4778 | list idlist; |
---|
4779 | for(i=1;i<=size(ttlist);i++) |
---|
4780 | { |
---|
4781 | if((i==size(ttlist))&&(typeof(ttlist[i])!="string")) break; |
---|
4782 | execute("ideal tid="+ttlist[i]+";"); |
---|
4783 | idlist[i]=list(tid,ideal(1),nE); |
---|
4784 | kill tid; |
---|
4785 | } |
---|
4786 | } |
---|
4787 | else |
---|
4788 | { |
---|
4789 | def So1=L[2][o1]; |
---|
4790 | setring So1; |
---|
4791 | if(defined(id1)) { kill id1; } |
---|
4792 | if(defined(id2)) { kill id2; } |
---|
4793 | if(defined(idlist)) { kill idlist; } |
---|
4794 | execute("ideal id1="+idname+";"); |
---|
4795 | if(deg(std(id1)[1])==0) |
---|
4796 | { |
---|
4797 | //--- problems with findTrans if id1 is empty set |
---|
4798 | //!!! todo: also correct in if branch!!! |
---|
4799 | setring R; |
---|
4800 | return(ideal(1)); |
---|
4801 | } |
---|
4802 | // id1=radical(id1); |
---|
4803 | |
---|
4804 | ideal id2=1; |
---|
4805 | list idlist; |
---|
4806 | if(defined(notE)){kill notE;} |
---|
4807 | list notE; |
---|
4808 | intvec nE; |
---|
4809 | idlist[1]=list(id1,id2,nE); |
---|
4810 | } |
---|
4811 | if(defined(tli)){kill tli;} |
---|
4812 | list tli; |
---|
4813 | if(defined(id1)) { kill id1; } |
---|
4814 | if(defined(id2)) { kill id2; } |
---|
4815 | ideal id1; |
---|
4816 | ideal id2; |
---|
4817 | if(defined(Etemp)){kill Etemp;} |
---|
4818 | ideal Etemp; |
---|
4819 | for(m=1;m<=size(idlist);m++) |
---|
4820 | { |
---|
4821 | id1=idlist[m][1]; |
---|
4822 | id2=idlist[m][2]; |
---|
4823 | nE=idlist[m][3]; |
---|
4824 | for(i=branchPos-1;i<=size(BO[4]);i++) |
---|
4825 | { |
---|
4826 | if(size(reduce(BO[4][i],std(id1+BO[1])))==0) |
---|
4827 | { |
---|
4828 | if(size(reduce(id1,std(BO[4][i])))!=0) |
---|
4829 | { |
---|
4830 | Etemp=BO[4][i]; |
---|
4831 | if(npars(basering)>0) |
---|
4832 | { |
---|
4833 | if(defined(prtemp)){kill prtemp;} |
---|
4834 | list prtemp=minAssGTZ(BO[4][i]); |
---|
4835 | if(size(prtemp)>1) |
---|
4836 | { |
---|
4837 | Etemp=ideal(1); |
---|
4838 | for(j=1;j<=size(prtemp);j++) |
---|
4839 | { |
---|
4840 | if(size(reduce(prtemp[j],std(id1)))==0) |
---|
4841 | { |
---|
4842 | Etemp=prtemp[j]; |
---|
4843 | break; |
---|
4844 | } |
---|
4845 | } |
---|
4846 | if(deg(std(Etemp)[1])==0) |
---|
4847 | { |
---|
4848 | ERROR("fetchInTree:something wrong in field extension"); |
---|
4849 | } |
---|
4850 | } |
---|
4851 | } |
---|
4852 | if(inJ) |
---|
4853 | { |
---|
4854 | tli=findTrans(id1+BO[2],Etemp,notE,BO[2]); |
---|
4855 | } |
---|
4856 | else |
---|
4857 | { |
---|
4858 | tli=findTrans(id1,Etemp,notE); |
---|
4859 | } |
---|
4860 | } |
---|
4861 | else |
---|
4862 | { |
---|
4863 | tli[1]=ideal(0); |
---|
4864 | tli[2]=ideal(1); |
---|
4865 | } |
---|
4866 | id1=tli[1]; |
---|
4867 | id2=id2*tli[2]; |
---|
4868 | notE[size(notE)+1]=BO[4][i]; |
---|
4869 | for(j=1;j<=size(DivL);j++) |
---|
4870 | { |
---|
4871 | if(inIVList(intvec(o1,i),DivL[j])) |
---|
4872 | { |
---|
4873 | nE[size(nE)+1]=j; |
---|
4874 | break; |
---|
4875 | } |
---|
4876 | } |
---|
4877 | if(size(nE)<size(notE)) |
---|
4878 | { |
---|
4879 | ERROR("fetchInTree: divisor not found in divL"); |
---|
4880 | } |
---|
4881 | } |
---|
4882 | idlist[m][1]=id1; |
---|
4883 | idlist[m][2]=id2; |
---|
4884 | idlist[m][3]=nE; |
---|
4885 | } |
---|
4886 | } |
---|
4887 | if(o1>1) |
---|
4888 | { |
---|
4889 | while(int(leadcoef(path[1,ncols(path)]))>=comPa) |
---|
4890 | { |
---|
4891 | if((int(leadcoef(path[1,ncols(path)]))>comPa)&& |
---|
4892 | (int(leadcoef(path[1,ncols(path)-1]))<comPa)) |
---|
4893 | { |
---|
4894 | ERROR("L[2][comPa] not in history of L[2][o1]!"); |
---|
4895 | } |
---|
4896 | def S=basering; |
---|
4897 | if(int(leadcoef(path[1,ncols(path)]))==1) |
---|
4898 | { |
---|
4899 | //--- that's the very first ring!!! |
---|
4900 | int und_jetzt_raus; |
---|
4901 | } |
---|
4902 | if(defined(T)){kill T;} |
---|
4903 | if(size(algext)>0) |
---|
4904 | { |
---|
4905 | if(defined(T0)){kill T0;} |
---|
4906 | def T0=L[2][int(leadcoef(path[1,ncols(path)]))]; |
---|
4907 | if(defined(tstr)){kill tstr;} |
---|
4908 | string tstr="ring T=(0,t),(" |
---|
4909 | +varstr(L[2][int(leadcoef(path[1,ncols(path)]))])+"),(" |
---|
4910 | +ordstr(L[2][int(leadcoef(path[1,ncols(path)]))])+");"; |
---|
4911 | execute(tstr); |
---|
4912 | setring T; |
---|
4913 | execute(algext); |
---|
4914 | minpoly=leadcoef(p); |
---|
4915 | kill tstr; |
---|
4916 | def BO=imap(T0,BO); |
---|
4917 | if(!defined(und_jetzt_raus)) |
---|
4918 | { |
---|
4919 | def path=imap(T0,path); |
---|
4920 | def lastMap=imap(T0,lastMap); |
---|
4921 | } |
---|
4922 | if(defined(idlist)){kill idlist;} |
---|
4923 | list idlist=list(list(ideal(1),ideal(1))); |
---|
4924 | } |
---|
4925 | else |
---|
4926 | { |
---|
4927 | def T=L[2][int(leadcoef(path[1,ncols(path)]))]; |
---|
4928 | setring T; |
---|
4929 | if(defined(id1)) { kill id1; } |
---|
4930 | if(defined(id2)) { kill id2; } |
---|
4931 | if(defined(idlist)){kill idlist;} |
---|
4932 | list idlist=list(list(ideal(1),ideal(1))); |
---|
4933 | } |
---|
4934 | setring S; |
---|
4935 | if(defined(phi)) { kill phi; } |
---|
4936 | map phi=T,lastMap; |
---|
4937 | for(m=1;m<=size(idlist);m++) |
---|
4938 | { |
---|
4939 | if(defined(id1)){kill id1;} |
---|
4940 | if(defined(id2)){kill id2;} |
---|
4941 | ideal id1=idlist[m][1]+BO[1]; |
---|
4942 | ideal id2=idlist[m][2]+BO[1]; |
---|
4943 | nE=idlist[m][3]; |
---|
4944 | setring T; |
---|
4945 | ideal id1=preimage(S,phi,id1); |
---|
4946 | ideal id2=preimage(S,phi,id2); |
---|
4947 | idlist[m]=list(id1,id2,nE); |
---|
4948 | kill id1,id2; |
---|
4949 | setring S; |
---|
4950 | } |
---|
4951 | setring T; |
---|
4952 | //--- after blowing down we might again be sitting inside a relevant |
---|
4953 | //--- exceptional divisor |
---|
4954 | for(m=1;m<=size(idlist);m++) |
---|
4955 | { |
---|
4956 | if(defined(id1)) {kill id1;} |
---|
4957 | if(defined(id2)) {kill id2;} |
---|
4958 | if(defined(notE)) {kill notE;} |
---|
4959 | if(defined(notE)) {kill notE;} |
---|
4960 | list notE; |
---|
4961 | ideal id1=idlist[m][1]; |
---|
4962 | ideal id2=idlist[m][2]; |
---|
4963 | nE=idlist[m][3]; |
---|
4964 | for(i=branchPos-1;i<=size(BO[4]);i++) |
---|
4965 | { |
---|
4966 | if(size(reduce(BO[4][i],std(id1)))==0) |
---|
4967 | { |
---|
4968 | if(size(reduce(id1,std(BO[4][i])))!=0) |
---|
4969 | { |
---|
4970 | if(defined(Etemp)) {kill Etemp;} |
---|
4971 | ideal Etemp=BO[4][i]; |
---|
4972 | if(npars(basering)>0) |
---|
4973 | { |
---|
4974 | if(defined(prtemp)){kill prtemp;} |
---|
4975 | list prtemp=minAssGTZ(BO[4][i]); |
---|
4976 | if(size(prtemp)>1) |
---|
4977 | { |
---|
4978 | Etemp=ideal(1); |
---|
4979 | for(j=1;j<=size(prtemp);j++) |
---|
4980 | { |
---|
4981 | if(size(reduce(prtemp[j],std(id1)))==0) |
---|
4982 | { |
---|
4983 | Etemp=prtemp[j]; |
---|
4984 | break; |
---|
4985 | } |
---|
4986 | } |
---|
4987 | if(deg(std(Etemp)[1])==0) |
---|
4988 | { |
---|
4989 | ERROR("fetchInTree:something wrong in field extension"); |
---|
4990 | } |
---|
4991 | } |
---|
4992 | } |
---|
4993 | if(defined(tli)) {kill tli;} |
---|
4994 | if(inJ) |
---|
4995 | { |
---|
4996 | def tli=findTrans(id1+BO[2],Etemp,notE,BO[2]); |
---|
4997 | } |
---|
4998 | else |
---|
4999 | { |
---|
5000 | def tli=findTrans(id1,Etemp,notE); |
---|
5001 | } |
---|
5002 | } |
---|
5003 | else |
---|
5004 | { |
---|
5005 | tli[1]=ideal(0); |
---|
5006 | tli[2]=ideal(1); |
---|
5007 | } |
---|
5008 | id1=tli[1]; |
---|
5009 | id2=id2*tli[2]; |
---|
5010 | notE[size(notE)+1]=BO[4][i]; |
---|
5011 | for(j=1;j<=size(DivL);j++) |
---|
5012 | { |
---|
5013 | if(inIVList(intvec(o1,i),DivL[j])) |
---|
5014 | { |
---|
5015 | nE[size(nE)+1]=j; |
---|
5016 | break; |
---|
5017 | } |
---|
5018 | } |
---|
5019 | if(size(nE)<size(notE)) |
---|
5020 | { |
---|
5021 | ERROR("fetchInTree: divisor not found in divL"); |
---|
5022 | } |
---|
5023 | } |
---|
5024 | idlist[m][1]=id1; |
---|
5025 | idlist[m][2]=id2; |
---|
5026 | idlist[m][3]=nE; |
---|
5027 | } |
---|
5028 | } |
---|
5029 | kill S; |
---|
5030 | if(defined(und_jetzt_raus)) |
---|
5031 | { |
---|
5032 | kill und_jetzt_raus; |
---|
5033 | break; |
---|
5034 | } |
---|
5035 | } |
---|
5036 | } |
---|
5037 | //---------------------------------------------------------------------------- |
---|
5038 | // Blow up ideal id1 from L[2][comPa] to L[2][m1]. To this end, first |
---|
5039 | // determine the path to follow and save it in path_togo. |
---|
5040 | //---------------------------------------------------------------------------- |
---|
5041 | if(m1==comPa) |
---|
5042 | { |
---|
5043 | if(size(algext)==0) |
---|
5044 | { |
---|
5045 | return(idlist[1][1]); |
---|
5046 | } |
---|
5047 | else |
---|
5048 | { |
---|
5049 | list retlist; |
---|
5050 | for(m=1;m<=size(idlist);m++) |
---|
5051 | { |
---|
5052 | retlist[m]=string(idlist[m][1]); |
---|
5053 | } |
---|
5054 | return(retlist); |
---|
5055 | } |
---|
5056 | } |
---|
5057 | if(defined(path_m1)) { kill path_m1; } |
---|
5058 | matrix path_m1=imap(Sm1,path); |
---|
5059 | intvec path_togo; |
---|
5060 | for(i=1;i<=ncols(path_m1);i++) |
---|
5061 | { |
---|
5062 | if(path_m1[1,i]>=comPa) |
---|
5063 | { |
---|
5064 | path_togo=path_togo,int(leadcoef(path_m1[1,i])); |
---|
5065 | } |
---|
5066 | } |
---|
5067 | path_togo=path_togo[2..size(path_togo)],m1; |
---|
5068 | i=1; |
---|
5069 | while(i<size(path_togo)) |
---|
5070 | { |
---|
5071 | def S=basering; |
---|
5072 | if(defined(T)){kill T;} |
---|
5073 | if(size(algext)>0) |
---|
5074 | { |
---|
5075 | if(defined(T0)){kill T0;} |
---|
5076 | def T0=L[2][path_togo[i+1]]; |
---|
5077 | if(defined(tstr)){kill tstr;} |
---|
5078 | string tstr="ring T=(0,t),(" +varstr(T0)+"),(" +ordstr(T0)+");"; |
---|
5079 | execute(tstr); |
---|
5080 | setring T; |
---|
5081 | execute(algext); |
---|
5082 | minpoly=leadcoef(p); |
---|
5083 | kill tstr; |
---|
5084 | def path=imap(T0,path); |
---|
5085 | def BO=imap(T0,BO); |
---|
5086 | def lastMap=imap(T0,lastMap); |
---|
5087 | if(defined(phi)){kill phi;} |
---|
5088 | map phi=S,lastMap; |
---|
5089 | list idlist=phi(idlist); |
---|
5090 | } |
---|
5091 | else |
---|
5092 | { |
---|
5093 | def T=L[2][path_togo[i+1]]; |
---|
5094 | setring T; |
---|
5095 | if(defined(phi)) { kill phi; } |
---|
5096 | map phi=S,lastMap; |
---|
5097 | if(defined(idlist)) {kill idlist;} |
---|
5098 | list idlist=phi(idlist); |
---|
5099 | idlist[1][1]=radical(idlist[1][1]); |
---|
5100 | idlist[1][2]=radical(idlist[1][2]); |
---|
5101 | } |
---|
5102 | for(m=1;m<=size(idlist);m++) |
---|
5103 | { |
---|
5104 | idlist[m][1]=sat(idlist[m][1]+BO[1],BO[4][size(BO[4])])[1]; |
---|
5105 | idlist[m][2]=sat(idlist[m][2],BO[4][size(BO[4])])[1]; |
---|
5106 | } |
---|
5107 | if((size(algext)==0)&&(deg(std(idlist[1][1])[1])==0)) |
---|
5108 | { |
---|
5109 | //--- strict transform empty in this chart, it will stay empty till the end |
---|
5110 | setring Sm1; |
---|
5111 | return(ideal(1)); |
---|
5112 | } |
---|
5113 | kill S; |
---|
5114 | i++; |
---|
5115 | } |
---|
5116 | ideal E,bla; |
---|
5117 | intvec kv; |
---|
5118 | list retlist; |
---|
5119 | for(m=1;m<=size(idlist);m++) |
---|
5120 | { |
---|
5121 | for(j=2;j<=size(idlist[m][3]);j++) |
---|
5122 | { |
---|
5123 | kv=findInIVList(1,path_togo[size(path_togo)],DivL[idlist[m][3][j]]); |
---|
5124 | if(kv!=intvec(0)) |
---|
5125 | { |
---|
5126 | E=E+BO[4][kv[2]]; |
---|
5127 | } |
---|
5128 | } |
---|
5129 | bla=quotient(idlist[m][1]+E,idlist[m][2]); |
---|
5130 | retlist[m]=string(bla); |
---|
5131 | } |
---|
5132 | if(size(algext)==0) |
---|
5133 | { |
---|
5134 | return(bla); |
---|
5135 | } |
---|
5136 | return(retlist); |
---|
5137 | } |
---|
5138 | ///////////////////////////////////////////////////////////////////////////// |
---|
5139 | static |
---|
5140 | proc findInIVList(int pos, int val, list ivl) |
---|
5141 | "Internal procedure - no help and no example available |
---|
5142 | " |
---|
5143 | { |
---|
5144 | //--- find entry with value val at position pos in list of intvecs |
---|
5145 | //--- and return the corresponding entry |
---|
5146 | int i; |
---|
5147 | for(i=1;i<=size(ivl);i++) |
---|
5148 | { |
---|
5149 | if(ivl[i][pos]==val) |
---|
5150 | { |
---|
5151 | return(ivl[i]); |
---|
5152 | } |
---|
5153 | } |
---|
5154 | return(intvec(0)); |
---|
5155 | } |
---|
5156 | ///////////////////////////////////////////////////////////////////////////// |
---|
5157 | static |
---|
5158 | proc inIVList(intvec iv, list li) |
---|
5159 | "Internal procedure - no help and no example available |
---|
5160 | " |
---|
5161 | { |
---|
5162 | //--- if intvec iv is contained in list li return 1, 0 otherwise |
---|
5163 | int i; |
---|
5164 | int s=size(iv); |
---|
5165 | for(i=1;i<=size(li);i++) |
---|
5166 | { |
---|
5167 | if(typeof(li[i])!="intvec"){ERROR("Not integer vector in the list");} |
---|
5168 | if(s==size(li[i])) |
---|
5169 | { |
---|
5170 | if(iv==li[i]){return(1);} |
---|
5171 | } |
---|
5172 | } |
---|
5173 | return(0); |
---|
5174 | } |
---|
5175 | ////////////////////////////////////////////////////////////////////////////// |
---|
5176 | static |
---|
5177 | proc Vielfachheit(ideal J,ideal I) |
---|
5178 | "Internal procedure - no help and no example available |
---|
5179 | " |
---|
5180 | { |
---|
5181 | //--- auxilliary procedure for addSelfInter |
---|
5182 | //--- compute multiplicity, suitable for the special situation there |
---|
5183 | int d=1; |
---|
5184 | int vd; |
---|
5185 | int c; |
---|
5186 | poly p; |
---|
5187 | ideal Ip,Jp; |
---|
5188 | while((d>0)||(!vd)) |
---|
5189 | { |
---|
5190 | p=randomLast(100)[nvars(basering)]; |
---|
5191 | Ip=std(I+ideal(p)); |
---|
5192 | c++; |
---|
5193 | if(c>20){ERROR("Vielfachheit: Dimension is wrong");} |
---|
5194 | d=dim(Ip); |
---|
5195 | vd=vdim(Ip); |
---|
5196 | } |
---|
5197 | Jp=std(J+ideal(p)); |
---|
5198 | return(vdim(Jp)/vdim(Ip)); |
---|
5199 | } |
---|
5200 | ///////////////////////////////////////////////////////////////////////////// |
---|
5201 | static |
---|
5202 | proc genus_E(list re, list iden0, intvec Eindex) |
---|
5203 | "Internal procedure - no help and no example available |
---|
5204 | " |
---|
5205 | { |
---|
5206 | int i,ge,gel,num; |
---|
5207 | def R=basering; |
---|
5208 | ring Rhelp=0,@t,dp; |
---|
5209 | def S=re[2][Eindex[1]]; |
---|
5210 | setring S; |
---|
5211 | def Sh=S+Rhelp; |
---|
5212 | //---------------------------------------------------------------------------- |
---|
5213 | //--- The Q-component X is reducible over C, decomposes into s=num components |
---|
5214 | //--- X_i, we assume they have n.c. |
---|
5215 | //--- s*g(X_i)=g(X)+s-1. |
---|
5216 | //---------------------------------------------------------------------------- |
---|
5217 | if(defined(I2)){kill I2;} |
---|
5218 | ideal I2=dcE[Eindex[2]][Eindex[3]][1]; |
---|
5219 | num=ncols(dcE[Eindex[2]][Eindex[3]][4]); |
---|
5220 | setring Sh; |
---|
5221 | if(defined(I2)){kill I2;} |
---|
5222 | ideal I2=imap(S,I2); |
---|
5223 | I2=homog(I2,@t); |
---|
5224 | ge=genus(I2); |
---|
5225 | gel=(ge+(num-1))/num; |
---|
5226 | if(gel*num-ge-num+1!=0){ERROR("genus_E: not divisible by num");} |
---|
5227 | setring R; |
---|
5228 | return(gel,num); |
---|
5229 | } |
---|
5230 | |
---|