source: git/Singular/LIB/reszeta.lib @ 4bde6b

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