source: git/Singular/LIB/resolve.lib @ cb8103a

spielwiese
Last change on this file since cb8103a was 66d68c, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13499 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 135.2 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Algebraic Geometry";
4info="
5LIBRARY:  resolve.lib      Resolution of singularities (Desingularization)
6                           Algorithm of Villamayor
7AUTHORS:  A. Fruehbis-Krueger,  anne@mathematik.uni-kl.de,
8@*        G. Pfister,           pfister@mathematik.uni-kl.de
9
10REFERENCES:
11[1] J.Kollar: Lectures on Resolution of Singularities, Princeton University Press (2007)
12@*  (contains large overview over various known methods for curves and surfaces as well as
13@*   a detailed description of the approach in the general case)
14[2] A.Bravo, S.Encinas, O.Villamayor: A Simplified Proof of Desingularisation and
15@*  Applications, Rev. Math. Iberoamericana 21 (2005), 349-458
16@*  (description of the algorithmic proof of desingularization in characteristic zero
17@*  which underlies this implementation)
18[3] A.Fruehbis-Krueger: Computational Aspects of Singularities, in J.-P. Brasselet,
19@*  J.Damon et al.: Singularities in Geometry and Topology, World Scientific
20@*  Publishing, 253--327 (2007)
21@*  (chapter 4 contains a detailed discussion on algorithmic desingularization and
22@*  efficiency aspects thereof)
23
24PROCEDURES:
25 blowUp(J,C[,W,E]) computes the blowing up of the variety V(J) (considered
26                   as embedded in V(W)) in the (smooth) center V(C),
27 blowUp2(J,C)      computes the blowing up of the variety V(J) in the
28                   (possibly singular) center V(C)
29 Center(J[,W,E])   computes 'Villamayor'-center for blow up
30 resolve(J)        computes the desingularization of the variety V(J)
31
32 showBO(BO)      prints the content of a BO in more human readable form
33 presentTree(L)  prints the final charts in more human readable form
34 showDataTypes() prints help text for output data types
35
36 blowUpBO(BO,C)  computes the blowing up of the variety V(BO[1]) in the
37                 center V(C). BO is a list (basic object), C is an ideal
38 createBO(J,W,E) creates basic object from input data
39 CenterBO(BO)    computes the center for the next blow-up of the
40                 given basic object
41 Delta(BO)       apply the Delta-operator of [Bravo,Encinas,Villamayor]
42 DeltaList(BO)   list of results of Delta^0 to Delta^bmax
43";
44LIB "elim.lib";
45LIB "primdec.lib";
46LIB "presolve.lib";
47LIB "linalg.lib";
48LIB "sing.lib";
49///////////////////////////////////////////////////////////////////////////////
50//                      Tasks:
51//                     1) optimization of the local case
52//                     2) optimization in Coeff
53//                     3) change invariant to represent coeff=1 case
54///////////////////////////////////////////////////////////////////////////////
55
56proc showDataTypes()
57"USAGE:  showDataTypes();
58RETURN:  nothing, only pretty printing of extended version of help text
59EXAMPLE: none
60"
61{
62"                               ";
63"//////////////// Short description of data type BO ///////////////////";
64"BO[1] an ideal, say Wi, defining the ambient space of the i-th chart";
65"      of the blowing up";
66"BO[2] an ideal defining the strict transform";
67"BO[3] intvec, the first integer b such that in the original object";
68"      (Delta^b(BO[2]))==1";
69"      the subsequent integers have the same property for Coeff-Objects";
70"      of BO[2] used when determining the center";
71"BO[4] the list of exceptional divisors";
72"BO[5] an ideal defining the map K[W] ----> K[Wi]";
73"BO[6] an intvec BO[6][j]=1 indicates that <BO[4][j],BO[2]>=1, i.e. the";
74"      strict transform does not meet the j-th exceptional divisor";
75"BO[7] intvec,";
76"      the index of the first blown-up object in the resolution process";
77"      leading to this object for which the value of b was BO[3]";
78"      the subsequent ones are the indices for the Coeff-Objects";
79"      of BO[2] used when determining the center";
80"BO[i], i>7: internal data";
81"                      ";
82pause();
83"              ";
84"///////////// Short description of data in a chart ///////////////////";
85"All chart data is stored in an object of type ring, the following ";
86"variables are always present in such a ring:";
87"BO:      list of type basic object containing the variety, ambient space,";
88"         exceptional divisors and further data (see previous page)";
89"cent:    ideal, describing the upcoming center determined by the algorithm";
90"path:    module (autoconverted to matrix)";
91"         path[1][i]= (i-1)st chart in resolution history of this chart";
92"         path[2][i]= index of chart of the blow up leading to i-th chart";
93"                     in resolution history";
94"lastMap: ideal, describing the preceding blow up leading to this chart";
95"                      ";
96pause();
97"                  ";
98"///////////// Short description of type resolution data //////////////";
99"list L containing two lists:";
100"L[1]: list of rings, containing the final charts";
101"L[2]: list of rings, containing all charts created in the resolution";
102"      process";
103"The most convenient way to view these data is the procedure presentTree()";
104"from this library. Alternatively, it can be digested using tools from";
105"the libraries reszeta.lib and resgraph.lib";
106"//////////////////////////////////////////////////////////////////////";
107}
108///////////////////////////////////////////////////////////////////////////////
109proc createBO(ideal J,list #)
110"USAGE:  createBO(J[,W][,E]);
111@*       J,W = ideals
112@*       E     = list
113ASSUME:  J  = ideal containing W ( W = 0 if not specified)
114@*       E  = list of smooth hypersurfaces (e.g. exceptional divisors)
115RETURN:  list BO representing a basic object :
116         BO[1] ideal W, if W has been specified; ideal(0) otherwise
117         BO[2] ideal J
118         BO[3] intvec
119         BO[4] the list E of exceptional divisors if specified;
120               empty list otherwise
121         BO[5] an ideal defining the identity map
122         BO[6] an intvec
123         BO[7] intvec
124         BO[8] a matrix
125         entries 3,5,6,7,8 are initialized appropriately for use of CenterBO
126         and blowUpBO
127EXAMPLE: example createBO;  shows an example
128"
129{
130  ideal W;
131  list E;
132  ideal abb=maxideal(1);
133  intvec v;
134  intvec bvec;
135  intvec w=-1;
136  matrix intE;
137  if(size(#)>0)
138  {
139    if(typeof(#[1])=="ideal")
140    {
141      W=#[1];
142    }
143    if(typeof(#[1])=="list")
144    {
145      E=#[1];
146    }
147    if(size(#)>1)
148    {
149      if((typeof(#[2])=="list") && (size(E)==0))
150      {
151        E=#[2];
152      }
153      if((typeof(#[2])=="ideal") && (size(W)==0))
154      {
155        W=#[2];
156      }
157    }
158  }
159  list BO=W,J,bvec,E,abb,v,w,intE;
160  return(BO);
161}
162example
163{"EXAMPLE:";
164   echo = 2;
165   ring R=0,(x,y,z),dp;
166   ideal J=x2-y3;
167   createBO(J,ideal(z));
168}
169///////////////////////////////////////////////////////////////////////////////
170
171proc blowUp(ideal J,ideal C,list #)
172"USAGE:  blowUp(J,C[,W][,E]);
173         W,J,C = ideals,
174         E     = list
175ASSUME:  J  = ideal containing W ( W = 0 if not specified)
176@*       C  = ideal containing J
177@*       E  = list of smooth hypersurfaces (e.g. exceptional divisors)
178NOTE:    W the ideal of the ambient space, C the ideal of the center of
179         the blowup and J the ideal of the variety
180         Important difference to blowUp2:
181         - the ambient space V(W) is blown up and V(J) transformed in it
182@*         - V(C) is assumed to be non-singular
183COMPUTE: the blowing up of W in C, the exceptional locus, the strict
184         transform of J and the blowup map
185RETURN:  list, say l, of size at most size(C),
186
187         l[i] is the affine ring corresponding to the i-th chart
188         each l[i] contains the ideals
189         - aS, ideal of the blownup ambient space
190         - sT, ideal of the strict transform
191         - eD, ideal of the exceptional divisor
192         - bM, ideal corresponding to the blowup map
193
194         l[i] also contains a list BO, which can best be viewed with showBO(BO)
195         detailed information on the data type BO can be viewed via the
196         command showDataTypes();
197EXAMPLE: example blowUp;  shows an example
198"
199{
200  def S=basering;
201  ideal W;
202  list E;
203  ideal abb=maxideal(1);
204  intvec v;
205  intvec bvec;
206  intvec w=-1;
207  matrix intE;
208  if(size(#)>0)
209  {
210    if(typeof(#[1])=="ideal")
211    {
212      W=#[1];
213    }
214    if(typeof(#[1])=="list")
215    {
216      E=#[1];
217    }
218    if(size(#)>1)
219    {
220      if((typeof(#[2])=="list") && (size(E)==0))
221      {
222        E=#[2];
223      }
224      if((typeof(#[2])=="ideal") && (size(W)==0))
225      {
226        W=#[2];
227      }
228    }
229  }
230  list BO=W,J,bvec,E,abb,v,w,intE;
231  int locaT;
232  export locaT;
233  list blow=blowUpBO(BO,C,0);
234  kill locaT;
235  int i;
236  for(i=1;i<=size(blow);i++)
237  {
238     def Q=blow[i];
239     setring Q;
240     ideal aS=BO[1];
241     ideal sT=BO[2];
242     ideal eD=BO[4][size(BO[4])];
243     ideal bM=BO[5];
244     kill lastMap;
245     kill thisChart;
246     export(aS);
247     export(sT);
248     export(eD);
249     export(bM);
250     blow[i]=Q;
251     setring S;
252     kill Q;
253  }
254  return(blow);
255}
256example
257{"EXAMPLE:";
258   echo = 2;
259   ring R=0,(x,y),dp;
260   ideal J=x2-y3;
261   ideal C=x,y;
262   list blow=blowUp(J,C);
263   def Q=blow[1];
264   setring Q;
265   aS;
266   sT;
267   eD;
268   bM;
269}
270///////////////////////////////////////////////////////////////////////////////
271
272proc blowUp2(ideal J,ideal C)
273"USAGE:  blowUp2(J,C);
274         J,C = ideals,
275ASSUME:  C  = ideal containing J
276NOTE:    C the ideal of the center of the blowup and J the ideal
277         of the variety
278         Important differences to blowUp:
279         - V(J) itself is blown up, not the ambient space
280         - C is not assumed to be non-singular
281COMPUTE: the blowing up of J in C, the exceptional locus and the blow-up
282         map
283RETURN:  list, say l, of size at most size(C),
284
285         l[i] is the affine ring corresponding to the i-th chart
286         each l[i] contains the ideals
287         - Jnew, ideal of the blownup J
288         - eD, ideal of the new exceptional divisor
289         - bM, ideal corresponding to the blowup map
290
291EXAMPLE: example blowUp2;  shows an example
292"
293{
294//----------------------------------------------------------------------------
295// Initialization and basic sanity checks of the input data
296//----------------------------------------------------------------------------
297  int i,j;
298  def S=basering;
299  def laM=maxideal(1);
300//--- number of generators of C should be as small as possible
301  def mstdC=mstd(C);
302  if(ncols(mstdC[1])<ncols(mstdC[2]))
303  {
304     C=mstdC[1];
305  }
306  else
307  {
308     C=mstdC[2];
309  }
310  C=simplify(interred(C),2);
311//--- the empty set is not a good center ;-)
312  if(deg(lead(C[1]))==0)
313  {
314     ERROR("Your chosen center was the empty set. Exiting.");
315  }
316//--- V(C) should be a subset of V(J)
317  if(size(reduce(J,std(C)))>0)
318  {
319     ERROR("V(J) does not contain V(C). Exiting.");
320  }
321//---------------------------------------------------------------------------
322// To compute the blowing up, we need to consider
323//       (K[x_1,\dots,x_n]/J)[t*C[1],...,t*C[m]]
324// which we want to represent as a quotient of K[x_1,...,x_n,y_1,..,y_m]
325// by an ideal obtained by elimination of the extra variable t.
326// (in the comments n=nvars(S) and m=ncols(C))
327//---------------------------------------------------------------------------
328//--- set up rings for the elimination
329  string r1="ring R=("+charstr(basering)+"),(x(1.."+string(nvars(S));
330  r1=r1+"),y(1.."+string(ncols(C))+")),dp;";
331  string r1t="ring Rt=("+charstr(basering)+"),(x(1.."+string(nvars(S));
332  r1t=r1t+"),y(1.."+string(ncols(C))+"),t),(dp,dp(1));";
333  execute(r1);       // ring for describing the transforms
334  execute(r1t);      // like r1, but with additional variable t
335  def J=fetch(S,J);
336  def C=fetch(S,C);
337//--- we need to eliminate t from J,y(1)-t*C[1],...y(m)-t*C[m]
338  ideal elId=J;
339  for(i=1;i<=ncols(C);i++)
340  {
341     elId=elId,y(i)-t*C[i];
342  }
343  elId=eliminate(elId,t);   // ideal describing the transform of J
344  setring R;                // get rid of t
345  def elId=fetch(Rt,elId);
346  def E=fetch(S,C);         // determine exceptional divisor
347  E=E+elId;
348  def laM0=fetch(S,laM);    // the blowup map
349//-----------------------------------------------------------------------------
350// The result is now represented in an A^n \times P^{m-1}, hence
351// involving n+m variables. For further computations we would like
352// to pass to charts to keep the total number of variables as low as
353// possible
354//----------------------------------------------------------------------------
355  list resList;
356  ideal Jsub;
357  ideal Esub;
358  ideal laM;
359  ideal testId;
360  map phi;
361  list templist;
362  for(i=1;i<=nvars(R)-nvars(S);i++)
363  {
364//--- first pass elId and E on to the i-th chart
365      Jsub=std(subst(elId,y(i),1));
366      if(deg(Jsub[1])==0)
367      {
368//--- transform does not meet this chart ==> ignore it
369         i++;
370         continue;
371      }
372      Esub=std(subst(E,y(i),1));
373//--- now get rid of unnecessary variables
374//--- first by appropriate coordinate changes
375      templist=elimpart(Jsub);
376      if(size(templist[2])>0)
377      {
378         phi=R,templist[5];
379         Jsub=phi(Jsub);
380         Jsub=simplify(interred(Jsub),2);
381         Esub=phi(Esub);
382         Esub=simplify(interred(Esub),2);
383         laM=phi(laM0);
384      }
385      else
386      {
387         laM=laM0;
388      }
389//--- then by dropping non-occuring variables
390      testId=Jsub,Esub,laM;
391      templist=findvars(testId,0);
392      if(size(templist[1])<nvars(R))
393      {
394         for(j=1;j<=size(templist[4]);j++)
395         {
396            laM=subst(laM,var(templist[4][j]),0);
397         }
398         execute("ring Rnew=("+charstr(basering)+"),("+string(templist[1])+"),dp;");
399         ideal Jnew=imap(R,Jsub);
400         ideal eD=imap(R,Esub);
401         ideal bM=imap(R,laM);
402      }
403      else
404      {
405         def Rnew=basering;
406         ideal Jnew=Jsub;
407         ideal eD=Esub;
408         ideal bM=laM;
409      }
410//--- export the relevant data of this ring and add the ring to the list
411      export Jnew;
412      export eD;
413      export bM;
414      resList[size(resList)+1]=Rnew;
415      setring R;
416      kill Rnew;
417  }
418  setring S;
419  return(resList);
420}
421example
422{  "EXAMPLE:";
423   echo = 2;
424   ring r=0,(x,y,z),dp;
425   ideal I=z2-x^3*y^2;
426   ideal C=z,xy;
427   list li=blowUp2(I,C);
428   size(li);                  // number of charts
429   def S1=li[1];
430   setring S1;                // chart 1
431   basering;
432   Jnew;
433   eD;
434   bM;
435   def S2=li[2];
436   setring S2;                // chart 2
437   basering;
438   Jnew;
439   eD;
440   bM;
441}
442///////////////////////////////////////////////////////////////////////////////
443
444proc Center(ideal J,list #)
445"USAGE:  Center(J[,W][,E])
446@*       J,W = ideals
447@*       E   = list
448ASSUME:  J  = ideal containing W ( W = 0 if not specified)
449@*       E  = list of smooth hypersurfaces (e.g. exceptional divisors)
450COMPUTE: the center of the blow-up of J for the resolution algorithm
451         of [Bravo,Encinas,Villamayor]
452RETURN:  ideal, describing the center
453EXAMPLE: example Center;  shows an example
454"
455{
456  ideal W;
457  list E;
458  ideal abb=maxideal(1);
459  intvec v;
460  intvec bvec;
461  intvec w=-1;
462  matrix intE;
463  if(size(#)>0)
464  {
465    if(typeof(#[1])=="ideal")
466    {
467      W=#[1];
468    }
469    if(typeof(#[1])=="list")
470    {
471      E=#[1];
472    }
473    if(size(#)>1)
474    {
475      if((typeof(#[2])=="list") && (size(E)==0))
476      {
477        E=#[2];
478      }
479      if((typeof(#[2])=="ideal") && (size(W)==0))
480      {
481        W=#[2];
482      }
483      if(size(#)==3){bvec=#[3];}
484    }
485  }
486
487  list BO=W,J,bvec,E,abb,v,w,intE,intvec(0);
488  if(defined(invSat)){kill invSat;}
489  list invSat=ideal(0),intvec(0);
490  export(invSat);
491
492  list re=CenterBO(BO);
493  ideal cent=re[1];
494  return(cent);
495}
496example
497{  "EXAMPLE:";
498   echo = 2;
499   ring R=0,(x,y),dp;
500   ideal J=x2-y3;
501   Center(J);
502}
503///////////////////////////////////////////////////////////////////////////////
504
505proc blowUpBO(list BO, ideal C,int e)
506"USAGE:  blowUpBO (BO,C,e);
507@*       BO = basic object, a list: ideal W,
508@*                                  ideal J,
509@*                                  intvec b,
510@*                                  list Ex,
511@*                                  ideal ab,
512@*                                  intvec v,
513@*                                  intvec w,
514@*                                  matrix M
515@*       C  = ideal
516@*       e  = integer (0 usual blowing up, 1 deleting extra charts, 2 deleting
517@*            no charts )
518ASSUME:  R  = basering, a polynomial ring, W an ideal of R,
519@*       J  = ideal containing W,
520@*       C  = ideal containing J
521COMPUTE: the blowing up of BO[1] in C, the exeptional locus, the strict
522         transform of BO[2]
523NOTE:    blowUpBO may be applied to basic objects in the sense of
524@*       [Bravo, Encinas, Villamayor] in the following referred to as BO and
525@*       to presentations in the sense of [Bierstone, Milman] in the following
526@*       referred to as BM.
527RETURN:  a list l of length at most size(C),
528         l[i] is a ring containing an object BO resp. BM:
529@*         BO[1]=BM[1] an ideal, say Wi, defining the ambient space of the
530               i-th chart of the blowing up
531@*         BO[2]=BM[2] an ideal defining the strict transform
532@*         BO[3] intvec, the first integer b such that in the original object
533               (Delta^b(BO[2]))==1
534               the subsequent integers have the same property for Coeff-Objects
535               of BO[2] used when determining the center
536@*         BM[3] intvec, BM[3][i] is the assigned multiplicity of BM[2][i]
537@*         BO[4]=BM[4] the list of exceptional divisors
538@*         BO[5]=BM[5] an ideal defining the map K[W] ----> K[Wi]
539@*         BO[6]=BM[6] an intvec BO[6][j]=1 indicates that <BO[4][j],BO[2]>=1,
540               i.e. the strict transform does not meet the j-th exceptional
541               divisor
542@*         BO[7] intvec,
543               the index of the first blown-up object in the resolution process
544               leading to this object for which the value of b was BO[3]
545               the subsequent ones are the indices for the Coeff-Objects
546               of BO[2] used when determining the center
547@*         BM[7] intvec, BM[7][i] is the index at which the (2i-1)st entry
548               of the invariant first reached its current maximal value
549@*         BO[8]=BM[8] a matrix indicating that BO[4][i] meets BO[4][j] by
550               BO[8][i,j]=1 for i < j
551@*         BO[9] empty
552@*         BM[9] the invariant
553
554EXAMPLE: example blowUpBO;  shows an example
555"
556{
557//---------------------------------------------------------------------------
558// Initialization  and sanity checks
559//---------------------------------------------------------------------------
560   def R0=basering;
561   if(!defined(locaT)){int locaT;}
562   if(locaT){poly pp=@p;}
563   intvec v;
564   int shortC=defined(shortcut);
565   int invS=defined(invSat);
566   int eq,hy;
567   int extra,noDel,keepDiv;
568   if(e==1){extra=1;}
569//---keeps all charts
570   if(e==2){noDel=1;}
571//---this is only for curves and surfaces
572//---keeps all charts with relevant informations on the exceptional divisors
573   if(e==3){keepDiv=1;}
574   if( typeof(attrib(BO[2],"isEqui"))=="int" )
575   {
576      eq=attrib(BO[2],"isEqui");
577   }
578   if( typeof(attrib(BO[2],"isHy"))=="int" )
579   {
580      hy=attrib(BO[2],"isHy");
581   }
582   string newvar;
583   int n=nvars(R0);
584   int i,j,l,m,x,jj,ll,haveCenters,co;
585//---the center should neither be the whole space nor empty
586   if((size(C)==0)||(deg(C[1])==0))
587   {
588      list result=R0;
589      return(result);
590   }
591   if(!defined(debugBlowUp))
592   {
593     int debugBlowUp=0;
594   }
595//---------------------------------------------------------------------------
596// Drop unnecessary variables
597//---------------------------------------------------------------------------
598//---step 1: substitution
599   if(!((keepDiv)||(noDel)))
600   {
601//!!! in case keepDiv and noDel:
602//!!! maybe simplify the situation by an appropriate coordinate change
603//!!! of this kind -- without dropping variables?
604      list L=elimpart(BO[1]);
605      if(size(L[2])!=0)
606      {
607         map psi=R0,L[5];
608         C=psi(C);
609         BO=psi(BO);
610      }
611
612      if(size(BO[1])==0)
613      {
614         ideal LL;
615         for(j=1;j<=size(BO[4]);j++)
616         {
617            LL=LL,BO[4][j];
618         }
619         LL=findvars(LL);
620         L=elimpart(BO[2]);
621         if((size(L[2])!=0)&&(size(std(LL+L[2]))==size(L[2])+size(LL)))
622         {
623            map chi=R0,L[5];
624            C=chi(C);
625            BO=chi(BO);
626         }
627      }
628//---step 2: dropping non-occurring variables
629      int s=size(C);
630      ideal K=BO[1],BO[2],C;
631      for(j=1;j<=size(BO[4]);j++)
632      {
633        K=K,BO[4][j];
634      }
635      list N=findvars(K,0);
636      if(size(N[1])<n)
637      {
638         newvar=string(N[1]);
639         v=N[4];
640         for(j=1;j<=size(v);j++){BO[5]=subst(BO[5],var(v[j]),0);}
641         execute("ring R1=("+charstr(R0)+"),("+newvar+"),dp;");
642         list BO=imap(R0,BO);
643         ideal C=imap(R0,C);
644         n=nvars(R1);
645      }
646      else
647      {
648         def R1=basering;
649      }
650   }
651   else
652   {
653     int s=size(C);
654     def R1=basering;
655   }
656   if(debugBlowUp)
657   {
658     "---> In BlowUp: After dropping unnecessary variables";
659     "BO:";
660     BO;
661     "C:";
662     C;
663   }
664//---------------------------------------------------------------------------
665// Do the actual blow-up
666//---------------------------------------------------------------------------
667//--- control the names of the variables
668   execute("ring R=("+charstr(R0)+"),(x(1..n)),dp;");
669   list BO=fetch(R1,BO);
670   ideal C=fetch(R1,C);
671   list Cmstd=mstd(C);
672   C=Cmstd[2];
673   if(size(Cmstd[1])<=size(Cmstd[2]))
674   {
675      C=Cmstd[1];
676   }
677   else
678   {
679      C=interred(C);
680   }
681   list result;
682//--- the blow-up process
683   ideal W  =BO[1];
684   ideal J  =BO[2];
685   intvec bvec  =BO[3];
686   list Ex  =BO[4];
687   ideal abb=BO[5];
688   intvec wvec=BO[7];
689   ideal laM=maxideal(1);
690   if((typeof(BO[9])=="intmat")||(typeof(BO[9])=="intvec"))
691   {
692      def @invmat=BO[9];
693   }
694   if(size(BO)>9)
695   {
696//--- check whether a previous center had been split into connected components
697      if(size(BO[10])>0)
698      {
699         list knownCenters=BO[10];
700         haveCenters=1;
701       }
702   }
703
704   matrix intE=BO[8];
705   Ex[size(Ex)+1]=var(1);
706         //to have the list depending on R in case BO[4] is empty
707
708   execute("ring S=("+charstr(R)+"),("+varstr(R)+",y(0..s-1)),dp;");
709   list resu;
710   list B;
711   execute("ring T=("+charstr(R)+"),("+varstr(R)+",y(0..s-1),t),dp;");
712   ideal C=imap(R,C);
713   ideal W=imap(R,W);
714   execute("map phi=S,"+varstr(R)+",t*C;")
715   setring S;
716
717//--- the ideal describing the blow-up map
718   ideal abb=imap(R,abb);
719   ideal laM0=imap(R,laM);
720//--- the ideal of the blowing up of the ambient space
721   ideal W=preimage(T,phi,W);
722//--- the ideal of the exceptional locus
723   ideal E=imap(R,C);
724
725   list E1=imap(R,Ex);
726   E1[size(E1)]=E;
727   ideal J=imap(R,J)+W;
728
729   if(haveCenters){list kN=imap(R,knownCenters);}
730
731//---  the strict transform of the exceptional divisors
732   for(j=1;j<size(E1);j++){E1[j]=sat(E1[j]+W,E)[1];}
733//---  the intersection matrix of the exceptional divisors
734   matrix intEold=imap(R,intE);
735   matrix intE[size(E1)][size(E1)];
736   ideal U,Jsub,sLstd;
737
738   for(j=1;j<size(E1);j++)
739   {
740       for(l=j+1;l<=size(E1);l++)
741       {
742          if(deg(E1[j][1])==0)
743          {
744             if(l<size(E1)){intE[j,l]=intEold[j,l];}
745             else             {intE[j,l]=0;}
746          }
747          else
748          {
749             if(deg(E1[l][1])==0)
750             {
751                if(l<size(E1)){intE[j,l]=intEold[j,l];}
752             }
753             else
754             {
755                U= std(E1[l]+E1[j]);
756                if(dim(U)>0){intE[j,l]=1;}
757                else        {intE[j,l]=0;}
758             }
759          }
760       }
761   }
762   if(debugBlowUp)
763   {
764     "----> In BlowUp: After Blowing-up, before Clean-Up";
765     "W:";
766     W;
767     "J:";
768     J;
769   }
770//----------------------------------------------------------------------------
771// generating and cleaning up the different charts
772//----------------------------------------------------------------------------
773   list M;
774   map psi;
775   list E2;
776   ideal K,JJ,laM,LL,MM;
777   n=nvars(S);
778   list N;
779   list Bstd;
780   intvec delCharts,extraCharts;
781   delCharts[s]=0;
782   extraCharts[s]=0;
783   ideal MA=y(0..s-1);
784   list ZRes,ZlaM,ZsLstd;
785   for(i=0;i<=s-1;i++)
786   {
787     if(haveCenters)
788     {
789        B[10]=kN;
790        for(j=1;j<=size(kN);j++)
791        {
792           B[10][j][1]=subst(B[10][j][1],y(i),1);
793        }
794     }
795     B[8]=intE;
796     B[1]=std(subst(W,y(i),1));
797     if(deg(B[1][1])==0)
798     {
799//--- subsets of the empty set are not really interesting!
800       delCharts[i+1]=1;
801       ZRes[i+1]=B;
802       ZlaM[i+1]=laM;
803       i++;
804       continue;
805     }
806     Jsub=subst(J,y(i),1);
807     attrib(Jsub,"isEqui",eq);
808     attrib(Jsub,"isHy",hy);
809     B[2]=Jsub;
810     B[3]=bvec;
811     for(j=1;j<size(E1);j++){E2[j]=subst(E1[j],y(i),1);}
812     E2[size(E1)]=E+B[1];
813     B[4]=E2;
814     M=elimpart(B[1]);
815     B[5]=abb;
816     laM=laM0;
817     psi=S,maxideal(1);
818     if(size(M[2])!=0)
819     {
820        psi=S,M[5];
821        B=psi(B);
822        laM=psi(laM);
823     }
824     Jsub=B[2];
825     B[2]=sat(Jsub,std(psi(E)))[1];
826     if(!defined(MAtmp)){ideal MAtmp=MA;}
827     MAtmp[i+1]=0;
828     JJ=std(B[2]+MAtmp);
829     if(deg(JJ[1])==0)
830     {
831        delCharts[i+1]=1;
832//--- the i-th chart will be marked for deleting because all informations
833//--- are already contained in the union of the remaining charts
834     }
835     else
836     {
837        if((eq)&&(dim(JJ)<dim(std(B[2]))))
838        {
839           extraCharts[i+1]=1;
840//---  compute the singular locus
841           if((B[3][1]<=1)&&(hy))
842//--- B[2] is a smooth hypersurface
843           {
844              sLstd=ideal(1);
845           }
846           else
847           {
848              Bstd=mstd(B[2]);
849              if(n-dim(Bstd[1])>4)
850              {
851//--- in this case the singular locus is too complicated
852                 sLstd=ideal(0);
853              }
854              JJ=Bstd[2];
855              attrib(JJ,"isEqui",eq);
856              B[2]=JJ;
857              sLstd=slocusE(B[2]);
858           }
859           m=0;
860           if(deg(std(sLstd+MAtmp)[1])==0)
861           {
862//--- the singular locus of B[2] is in the union of the remaining charts
863              m=1;
864              for(l=1;l<=size(B[4]);l++)
865              {
866                 if(deg(std(B[2]+B[4][l]+MAtmp)[1])!=0)
867                 {
868//--- the exceptional divisor meets B[2] at the locus of MAtmp
869//--- we continue only if the option extra=1 and we have transversal
870//--- intersection
871                    m=0;
872                    break;
873                  }
874              }
875           }
876           if(m)
877           {
878//--- the i-th chart will be marked for deleting because all informations
879//--- are already contained in the union of the remaining charts
880              delCharts[i+1]=1;
881           }
882        }
883     }
884     if(delCharts[i+1]==0)
885     {
886        MAtmp[i+1]=MA[i+1];
887        ZsLstd[i+1]=sLstd;
888     }
889     ZRes[i+1]=B;
890     ZlaM[i+1]=laM;
891   }
892//---------------------------------------------------------------------------
893// extra = ignore uninteresting charts even if there is a normal
894//         crosssing intersection in it
895//---------------------------------------------------------------------------
896   if(extra)
897   {
898      for(i=0;i<=s-1;i++)
899      {
900         if((delCharts[i+1]==0)&&(extraCharts[i+1]))
901         {
902           MAtmp[i+1]=0;
903           B=ZRes[i+1];
904           sLstd=ZsLstd[i+1];
905           m=0;
906           if(deg(std(sLstd+MAtmp)[1])==0)
907           {
908//--- the singular locus of B[2] is in the union of the remaining charts
909              m=1;
910              for(l=1;l<=size(B[4]);l++)
911              {
912                 if(deg(std(B[2]+B[4][l]+MAtmp)[1])!=0)
913                 {
914//--- the exceptional divisor meets B[2] at the locus of MAtmp
915//--- we continue only if the option extra=1 and we have transversal
916//--- intersection
917                    m=2;
918                    if(!transversalTB(B[2],list(B[4][l]),MAtmp))
919                    {
920                       m=0;break;
921                    }
922                  }
923              }
924           }
925           if(m)
926           {
927              if(m==1)
928              {
929//--- the i-th chart will be marked for deleting because all informations
930//--- are already contained in the union of the remaining charts
931                 delCharts[i+1]=1;
932              }
933              else
934              {
935//--- the option extra=1 and we have transversal intersection
936//--- we delete the chart in case of normal crossings
937                 if(normalCrossB(B[2],B[4],MAtmp))
938                 {
939//--- in case of the option extra
940//--- the i-th chart will be marked for deleting because all informations
941//--- are already contained in the union of the remaining charts
942
943                   delCharts[i+1]=1;
944                 }
945              }
946           }
947           if(delCharts[i+1]==0)
948           {
949              MAtmp[i+1]=MA[i+1];
950           }
951         }
952      }
953      for(i=0;i<=s-1;i++)
954      {
955         if(!delCharts[i+1]){break;}
956      }
957      if(i==s){delCharts[s]=0;}
958   }
959   for(i=0;i<=s-1;i++)
960   {
961     B=ZRes[i+1];
962     laM=ZlaM[i+1];
963     if(noDel){delCharts[i+1]=0;}
964//---keeps chart if the exceptional divisor is not in any other chart
965     if((delCharts[i+1])&&(keepDiv))
966     {
967         for(j=1;j<=size(B[4]);j++)
968         {
969            if(deg(std(B[4][j])[1])>0)
970            {
971               x=0;
972               for(l=0;l<=s-1;l++)
973               {
974                  if((l!=i)&&(!delCharts[l+1])&&(deg(std(ZRes[l+1][4][j])[1])>0))
975                  {
976                     x=1;
977                     break;
978                  }
979               }
980               if(!x)
981               {
982                  delCharts[i+1]=0;
983//!!!evtl. diese Karten markieren und nicht weiter aufblasen???
984                  break;
985               }
986            }
987         }
988     }
989//---keeps charts if the intersection of 2 divisors is not in any other chart
990     if((delCharts[i+1])&&(keepDiv))
991     {
992         for(j=1;j<=size(B[4])-1;j++)
993         {
994            for(l=j+1;l<=size(B[4]);l++)
995            {
996               if(deg(std(B[4][j]+B[4][l])[1])>0)
997               {
998                  x=0;
999                  for(ll=0;ll<=s-1;ll++)
1000                  {
1001                     if((ll!=i)&&(!delCharts[ll+1]))
1002                     {
1003                        if(deg(std(ZRes[ll+1][4][j]+ZRes[ll+1][4][l])[1])>0)
1004                        {
1005                           x=1;
1006                           break;
1007                        }
1008                     }
1009                  }
1010                  if(!x)
1011                  {
1012                     delCharts[i+1]=0;
1013                     break;
1014                  }
1015               }
1016            }
1017            if(!delCharts[i+1]){break;}
1018         }
1019     }
1020//---keeps charts if the intersection of 3 divisors is not in any other chart
1021     if((delCharts[i+1])&&(keepDiv))
1022     {
1023         for(j=1;j<=size(B[4])-2;j++)
1024         {
1025            for(l=j+1;l<=size(B[4])-1;l++)
1026            {
1027               for(ll=l+1;ll<=size(B[4]);ll++)
1028               {
1029                  if(deg(std(B[4][j]+B[4][l]+B[4][ll])[1])>0)
1030                  {
1031                     x=0;
1032                     for(jj=0;jj<=s-1;jj++)
1033                     {
1034                        if((jj!=i)&&(!delCharts[jj+1]))
1035                        {
1036                           if(deg(std(ZRes[jj+1][4][j]
1037                              +ZRes[jj+1][4][l]+ZRes[jj+1][4][ll])[1])>0)
1038                           {
1039                              x=1;
1040                              break;
1041                           }
1042                        }
1043                     }
1044                     if(!x)
1045                     {
1046                       delCharts[i+1]=0;
1047                       break;
1048                     }
1049                  }
1050               }
1051               if(!delCharts[i+1]){break;}
1052            }
1053            if(!delCharts[i+1]){break;}
1054        }
1055     }
1056     if(delCharts[i+1]==0)
1057     {
1058//--- try to decrease the number of variables by substitution
1059        if((!keepDiv)&&(!noDel))
1060        {
1061           list WW=elimpart(B[1]);
1062           map phiW=basering,WW[5];
1063           B=phiW(B);
1064           laM=phiW(laM);
1065           kill WW;
1066           kill phiW;
1067           if(size(B[1])==0)
1068           {
1069              LL=0;
1070              for(j=1;j<=size(B[4]);j++)
1071              {
1072                 MM=std(B[4][j]);
1073                 if(deg(MM[1])>0){LL=LL,MM;}
1074              }
1075              LL=findvars(LL);
1076              M=elimpart(B[2]);
1077              if((size(M[2])!=0)&&(size(std(LL+M[2]))==size(M[2])+size(LL)))
1078              {
1079                psi=S,M[5];
1080                B=psi(B);
1081                laM=psi(laM);
1082              }
1083           }
1084         }
1085//---- interreduce B[1],B[2] and all B[4][j]
1086        B[1]=interred(B[1]);
1087        B[2]=interred(B[2]);
1088        E2=B[4];
1089        for(j=1;j<=size(E2);j++){E2[j]=interred(E2[j]);}
1090        B[4]=E2;
1091        v=0;v[size(E2)]=0;
1092//--- mark those j for which B[4] does not meet B[2]
1093        for(j=1;j<=size(E2);j++)
1094        {
1095            K=E2[j],B[2];
1096            K=std(K);
1097            if(deg(K[1])==0)
1098            {
1099              v[j]=1;
1100            }
1101        }
1102        B[6]=v;
1103        B[7]=wvec;
1104//--- throw away variables which do not occur
1105        K=B[1],B[2],B[5];          //Aenderung!!!
1106        for(j=1;j<=size(B[4]);j++){K=K,B[4][j];}
1107        N=findvars(K,0);
1108        if(size(N[1])<n)
1109        {
1110           newvar=string(N[1]);
1111           v=N[4];
1112           for(j=1;j<=size(v);j++)
1113           {
1114              B[5]=subst(B[5],var(v[j]),0);
1115              laM=subst(laM,var(v[j]),0);
1116           }
1117           execute("ring R2=("+charstr(S)+"),("+newvar+"),dp;");
1118           list BO=imap(S,B);
1119           ideal laM=imap(S,laM);
1120        }
1121        else
1122        {
1123           def R2=basering;
1124           list BO=B;
1125        }
1126        ideal JJ=BO[2];
1127        attrib(JJ,"isEqui",eq);
1128        attrib(JJ,"isHy",hy);
1129        BO[2]=JJ;
1130//--- strict transforms of the known centers
1131        if(haveCenters)
1132        {
1133           ideal tt;
1134           list tList;
1135           for(j=1;j<=size(BO[10]);j++)
1136           {
1137              tt=std(BO[10][j][1]);
1138              if(deg(tt[1])>0)
1139              {
1140                tt=sat(tt,BO[4][size(BO[4])])[1];
1141              }
1142              if((deg(tt[1])>0)&&(deg(std(tt+BO[2]+BO[1])[1])>0))
1143              {
1144                 tList[size(tList)+1]=
1145                            list(tt,BO[10][j][2],BO[10][j][3],BO[10][j][4]);
1146              }
1147           }
1148           BO[10]=tList;
1149           kill tList;
1150        }
1151//--- marking  variables which do not occur in BO[1] and BO[2]
1152//--- and occur in exactly one BO[4][j], which is the hyperplane given by
1153//--- this variable
1154//!!!! not necessarily in exactly one BO[4][j]
1155        list N=findvars(BO[1]+BO[2],0);
1156        if(size(N[1])<nvars(basering))
1157        {
1158           v=N[4];
1159           if(defined(H)){kill H;}
1160           if(defined(EE)){kill EE;}
1161           if(defined(vv)){kill vv;}
1162           list EE;
1163           intvec vv;
1164           ideal H=maxideal(1);
1165           for(j=1;j<=size(v);j++)
1166           {
1167             H[v[j]]=0;
1168           }
1169           H=std(H);
1170           for(l=1;l<=size(BO[4]);l++)
1171           {
1172              if(BO[6][l]==1){l++;continue;}
1173              if(size(ideal(matrix(reduce(BO[4][l],H))-BO[4][l]))==0)
1174              {
1175//!!! need further cleanup:
1176//!!! this part of the code is no longer used since it did not glue properly
1177//                 BO[6][l]=2;
1178                 EE[size(EE)+1]=BO[4][l];
1179                 vv[size(vv)+1]=l;
1180              }
1181           }
1182           if((size(vv)>dim(std(BO[2])))&&(deg(BO[2][1])>0))
1183           {
1184              list BOtemp3=BO;
1185              BOtemp3[4]=EE;
1186              intvec ivtemp3;
1187              ivtemp3[size(BOtemp3[4])]=0;
1188              BOtemp3[6]=ivtemp3;
1189              BOtemp3[7][1]=-1;
1190              list iEtemp3=inters_E(BOtemp3);
1191              if(iEtemp3[2]>=dim(std(BOtemp3[2])))
1192              {
1193                 for(l=2;l<=size(vv);l++)
1194                 {
1195                    BO[6][vv[l]]=0;
1196                 }
1197              }
1198              kill BOtemp3,ivtemp3,iEtemp3;
1199           }
1200        }
1201        list thisChart=ideal(0),i;
1202        export thisChart;
1203
1204//----------------------------------------------------------------------------
1205// export the basic object and append the ring to the list of rings
1206//----------------------------------------------------------------------------
1207        if(debugBlowUp)
1208        {
1209           "----> In BlowUp: Adding a single chart";
1210           "BO:";
1211           BO;
1212        }
1213        if(locaT)
1214        {
1215           map locaPhi=R0,laM;
1216           poly @p=locaPhi(pp);
1217           export(@p);
1218        }
1219        ideal lastMap=laM;
1220        export lastMap;
1221        if(invS){list invSat=imap(R0,invSat);export invSat;}
1222        if(defined(@invmat)){BO[9]=@invmat;}
1223        if(shortC){list shortcut=imap(R0,shortcut);export(shortcut);}
1224        export BO;
1225        result[size(result)+1]=R2;
1226        setring S;
1227        kill R2;
1228     }
1229   }
1230   setring R0;
1231   return(result);
1232}
1233example
1234{"EXAMPLE:";
1235   echo = 2;
1236   ring R=0,(x,y),dp;
1237
1238   ideal W;
1239   ideal J=x2-y3;
1240   intvec b=1;
1241   list E;
1242   ideal abb=maxideal(1);
1243   intvec v;
1244   intvec w=-1;
1245   matrix M;
1246   intvec ma;
1247   list BO=W,J,b,E,abb,v,w,M,ma;
1248
1249   ideal C=CenterBO(BO)[1];
1250
1251   list blow=blowUpBO(BO,C,0);
1252   def Q=blow[1];
1253   setring Q;
1254   BO;
1255}
1256//////////////////////////////////////////////////////////////////////////////
1257static proc slocusE(ideal i)
1258"Internal procedure - no help and no example available
1259"
1260{
1261//--- do slocus in equidimensional case directly -- speed up
1262  if(size(i)==0){return(ideal(1));}
1263  if( typeof(attrib(i,"isEqui"))=="int" )
1264  {
1265     if(attrib(i,"isEqui")==1)
1266     {
1267        ideal j=std(i);
1268        if(deg(j[1])==0){return(j);}
1269        int cod  = nvars(basering)-dim(j);
1270        i        = i+minor(jacob(i),cod);
1271        return(i);
1272     }
1273  }
1274  return(slocus(i));
1275}
1276///////////////////////////////////////////////////////////////////////////////
1277static proc inters_E(list BO)
1278"USAGE:  inters_E(BO);
1279@*       BO = basic object, a list: ideal W,
1280@*                                  ideal J,
1281@*                                  intvec b,
1282@*                                  list Ex,
1283@*                                  ideal ab,
1284@*                                  intvec v,
1285@*                                  intvec w,
1286@*                                  matrix M
1287ASSUME:  R  = basering, a polynomial ring, W an ideal of R,
1288@*       J  = ideal containing W,
1289@*       BO in the setting of case 2 of [Bravo,Encinas,Villamayor]
1290@*             BO[4]=E, BO[4][1..count]=E^-
1291@*             BO[7][1]=count
1292COMPUTE: (W,(P,1),E^+) in the notation of [Bravo,Encinas,Villamayor]
1293RETURN:  a list l ,
1294         l[1]: P = product of ideals I(H_i1)+..+I(H_in) over all
1295                   n-tuples of indices i1..in from 1..count
1296         l[2]: n = maximal number of H_i from E^- meeting J simultaneously
1297         l[3]: maximal number of H_i from E meeting J simultaneously
1298EXAMPLE: internal procedure - no example available
1299"
1300{
1301//---------------------------------------------------------------------
1302// Initialization
1303//---------------------------------------------------------------------
1304  int kk,jj,ii,updated,n,count2,kkdiff;
1305  def rb=basering;
1306  def W=BO[1];
1307  ideal J=BO[1],BO[2];
1308  int nonnormal;
1309  int maxkk=dim(std(J));
1310  int dimJ=maxkk;
1311  ideal test2;
1312  list merklist1,merklist2;
1313  if(size(BO[4])==0)
1314  {
1315    list retlist=BO[2],n;
1316    return(retlist);
1317  }
1318  def E=BO[4];
1319  intvec stoplist=BO[6];
1320//--- fill in all known information about exceptional divisors not meeting
1321//--- current chart
1322  for(ii=1;ii<=size(E);ii++)
1323  {
1324    if(deg(std(E[ii])[1])==0)
1325    {
1326      stoplist[ii]=1;
1327    }
1328  }
1329
1330  int count=BO[7][1];
1331  if(!defined(debug_Inters_E))
1332  {
1333    int debug_Inters_E=0;
1334  }
1335//---------------------------------------------------------------------
1336// we only want to look at E^-, not at all of E
1337//---------------------------------------------------------------------
1338  if (count>-1)
1339  {
1340    if (count>0)
1341    {
1342      list E_new=E[1..count];
1343      count2=size(E);
1344    }
1345    else
1346    {
1347      list E_new;
1348      count2=size(E);
1349    }
1350  }
1351  else
1352  {
1353    list E_new=E;
1354    count=size(E_new);
1355    count2=count;
1356  }
1357//---------------------------------------------------------------------
1358// combinatorics is expensive in an interpreted language,
1359// we leave it to the kernel by translating it into monomial
1360// ideals in a new ring with variables t(i)
1361//---------------------------------------------------------------------
1362  string rstr="ring rcomb=0,(t(1.." + string(size(E)) + ")),dp;";
1363  execute(rstr);
1364  ideal potid,potid2;
1365  list monlist,comblist,merkmon;
1366  for(kk=1;kk<=count;kk++)
1367  {
1368    if(stoplist[kk]==0)
1369    {
1370//**************************************************************************/
1371// it does not make sense to intersect twice by the same E_i
1372// ===> reduce by t(i)^2
1373//**************************************************************************/
1374      potid=potid,t(kk)^2;
1375    }
1376    else
1377    {
1378//**************************************************************************/
1379// it does not make sense to consider E_i with v[i]==1 ===> reduce by t(i)
1380//**************************************************************************/
1381      potid=potid,t(kk);
1382//**************************************************************************/
1383// if stoplist[kk]==2 then J and all E_i automatically intersect E_kk
1384// hence we need not test it, but we have to lower maxkk by one
1385//**************************************************************************/
1386      if(stoplist[kk]==2)
1387      {
1388        maxkk=maxkk-1;
1389        kkdiff++;             // count these for dimension check later on
1390
1391      }
1392    }
1393  }
1394
1395  potid2=std(potid);
1396  if(count2>count)
1397  {
1398    potid=potid,t((count+1)..count2);
1399    for(kk=max(1,count);kk<=count2;kk++)
1400    {
1401      potid2=potid2,t(kk)^2;
1402    }
1403  }
1404  potid=std(potid);
1405  potid2=std(potid2);
1406  for(kk=1;(((kk<=count)||(kk<=maxkk+1))&&(kk<=count2));kk++)
1407  {
1408//-------------------------------------------------------------------------
1409// monlist[kk]=lists of kk entries of E_new, not containing an E_i twice,
1410// not containing an E_i where v[i]==1
1411//-------------------------------------------------------------------------
1412    monlist[kk]=redMax(kk,potid);
1413//*************************************************************************/
1414// in the case of n<=maxkk we also need to know whether n would still be
1415// below this bound if we considered all of E instead of E_new
1416// ===> merkmon contains previously ignored tuples E_i1,..,E_im
1417//*************************************************************************/
1418    if(kk<=maxkk+1)
1419    {
1420      merkmon[kk]=redMax(kk,potid2);
1421      merkmon[kk]=simplify(reduce(merkmon[kk],std(monlist[kk])),2);
1422    }
1423  }
1424  if(debug_Inters_E)
1425  {
1426    "----> In Inters_E: the tuples";
1427    "tuples of E^-:";
1428    monlist;
1429    "the remaining tuples:";
1430    merkmon;
1431  }
1432//-------------------------------------------------------------------------
1433// check whether there is a kk-tuple of E_i intersecting J,
1434// kk running from 1 to count
1435//-------------------------------------------------------------------------
1436  for(kk=1;kk<=count;kk++)
1437  {
1438    if(size(monlist[kk]==0)) break;
1439    kill comblist;
1440    list comblist;
1441//--- transscribe the tuples from monomial notation to intvec notation
1442    for(jj=1;jj<=ncols(monlist[kk]);jj++)
1443    {
1444      comblist[jj]=leadexp(monlist[kk][jj]);
1445    }
1446    setring rb;
1447    updated=0;
1448//------------------------------------------------------------------------
1449// Do the intersections
1450//------------------------------------------------------------------------
1451    for(jj=1;jj<=size(comblist);jj++)
1452    {
1453//--- jj-th tuple from list of tuples of kk E_i
1454      test2=J;
1455      for(ii=1;ii<=count;ii++)
1456      {
1457        if(comblist[jj][ii]==1)
1458        {
1459          test2=test2,E_new[ii];
1460        }
1461      }
1462      test2=std(test2);
1463//--- check whether this intersection is non-empty and store it accordingly
1464      if(deg(test2[1])!=0)
1465      {
1466//--- it is non-empty
1467        if(updated!=0)
1468        {
1469          merklist1[size(merklist1)+1]=comblist[jj];
1470        }
1471        else
1472        {
1473          kill merklist1;
1474          list merklist1;
1475          merklist1[1]=comblist[jj];
1476          updated=1;
1477          n=kk;
1478        }
1479        if(dim(test2)!=maxkk-kk+kkdiff)
1480        {
1481          nonnormal=1;
1482        }
1483      }
1484      else
1485      {
1486//--- it is empty
1487        merklist2[size(merklist2)+1]=jj;
1488      }
1489    }
1490    setring rcomb;
1491    ideal redid;
1492//---------------------------------------------------------------------
1493// update monlist and merkmon by the knowledge what intersections are
1494// empty in the kk-th step
1495//---------------------------------------------------------------------
1496    for(jj=1;jj<=size(merklist2);jj++)
1497    {
1498      redid=redid,monlist[kk][merklist2[jj]];
1499    }
1500    for(jj=kk+1;jj<=count;jj++)
1501    {
1502      monlist[jj]=simplify(reduce(monlist[jj],std(redid)),2);
1503      if(jj<=maxkk+1)
1504      {
1505        merkmon[jj]=simplify(reduce(merkmon[jj],std(redid)),2);
1506      }
1507    }
1508    kill redid;
1509    kill merklist2;
1510    list merklist2;
1511  }
1512  if(debug_Inters_E)
1513  {
1514    "----> In Inters_E: intersections found:";
1515    merklist1;
1516  }
1517//---------------------------------------------------------------------
1518// form the union of the intersections of the appropriate E_i
1519//---------------------------------------------------------------------
1520  setring rb;
1521  ideal center,dummy;
1522  list centlist;
1523  for(kk=1;kk<=size(merklist1);kk++)
1524  {
1525    for(jj=1;jj<=size(merklist1[kk]);jj++)
1526    {
1527      if(merklist1[kk][jj]==1)
1528      {
1529        dummy=dummy,E_new[jj];
1530      }
1531    }
1532    if(size(center)==0)
1533    {
1534      center=dummy;
1535      centlist[1]=dummy;
1536    }
1537    else
1538    {
1539      center=intersect(center,dummy);
1540      centlist[size(centlist)+1]=dummy;
1541    }
1542    dummy=0;
1543  }
1544  if(debug_Inters_E)
1545  {
1546    "----> In Inters_E: intersection of E_i";
1547    "maximal number of E_i encountered in:";
1548    center;
1549    "the components of this locus:";
1550    centlist;
1551    "maximal number of E_i from E^- intersecting simultaneously:",n;
1552    if(nonnormal)
1553    {
1554      "flag nonnormal is set";
1555    }
1556  }
1557  list retlist=center,n;
1558//-------------------------------------------------------------------------
1559// If n<=maxkk, then test if this is the case for all of E not just E_new
1560// using the pairs indicated by merkmon
1561//-------------------------------------------------------------------------
1562  int ntotal=n;
1563  if((n<=maxkk)&&(n<count2)&&(!nonnormal))
1564  {
1565//--- check kk-tuples
1566    for(kk=n+1;kk<=maxkk+1;kk++)
1567    {
1568      setring rcomb;
1569//--- check if there are combinations to be checked
1570      if(kk>size(merkmon))
1571      {
1572        setring rb;
1573        break;
1574      }
1575      if(size(merkmon[kk])!=0)
1576      {
1577        kill comblist;
1578        list comblist;
1579//--- transscribe tuples from  monomial notation to intvec notation
1580        for(jj=1;jj<=size(merkmon[kk]);jj++)
1581        {
1582          comblist[jj]=leadexp(merkmon[kk][jj]);
1583        }
1584        setring rb;
1585//--- check jj-th tuple from the list of kk-tuples
1586        for(jj=1;jj<=size(comblist);jj++)
1587        {
1588          test2=J;
1589          for(ii=1;ii<=nvars(rcomb);ii++)
1590          {
1591            if(comblist[jj][ii]==1)
1592            {
1593              test2=test2,E[ii];
1594            }
1595          }
1596          test2=std(test2);
1597//--- as soon as we found one we can proceed to the subsequent kk
1598          if(deg(test2[1])!=0)
1599          {
1600            ntotal=kk;
1601            if(dim(test2)-kkdiff!=maxkk-kk)
1602            {
1603              nonnormal=2;
1604              break;
1605            }
1606          }
1607        }
1608//--- if we already know that too many E_i intersect simultaneously,
1609//--- we need not proceed any further
1610        if(nonnormal)
1611        {
1612          break;
1613        }
1614      }
1615      else
1616      {
1617        setring rb;
1618        break;
1619      }
1620    }
1621  }
1622//-------------------------------------------------------------------------
1623// update the result accordingly and return it
1624//-------------------------------------------------------------------------
1625  if(maxkk<dimJ)
1626  {
1627    n=n+dimJ-maxkk;
1628    ntotal=ntotal+dimJ-maxkk;
1629  }
1630  retlist[2]=n;
1631  retlist[3]=ntotal;
1632  if(n<=dimJ)
1633  {
1634    retlist[4]=centlist;
1635    retlist[5]=merklist1;
1636    if(nonnormal)
1637    {
1638      retlist[6]=nonnormal;
1639    }
1640  }
1641  return(retlist);
1642}
1643///////////////////////////////////////////////////////////////////////////
1644
1645proc Delta(list BO)
1646"USAGE:  Delta (BO);
1647@*       BO = basic object, a list: ideal W,
1648@*                                  ideal J,
1649@*                                  intvec b,
1650@*                                  list Ex,
1651@*                                  ideal ab,
1652@*                                  intvec v,
1653@*                                  intvec w,
1654@*                                  matrix M
1655ASSUME:  R  = basering, a polynomial ring, W an ideal of R,
1656@*       J  = ideal containing W
1657COMPUTE: Delta-operator applied to J in the notation of
1658         [Bravo,Encinas,Villamayor]
1659RETURN:  ideal
1660EXAMPLE: example Delta;  shows an example
1661"
1662{
1663//---------------------------------------------------------------------------
1664// Initialization  and sanity checks
1665//---------------------------------------------------------------------------
1666   ideal W=BO[1];
1667   ideal J=BO[2];
1668   ideal C=simplify(reduce(J,std(W)),2);
1669   list LC;
1670   int n=nvars(basering);
1671//---------------------------------------------------------------------------
1672// Simple case: W is the empty set
1673//---------------------------------------------------------------------------
1674   if(size(W)==0)
1675   {
1676      C=C,jacob(J);
1677      C=std(C);
1678      return(C);
1679   }
1680//---------------------------------------------------------------------------
1681// General case: W is non-empty
1682// Step 1: Find a minor of the Jacobian of W which is not identically zero
1683//         and look at the complement of the zero-set given by this minor;
1684//         this leads to the system of local parameters
1685// Step 2: Form the derivatives w.r.t. this system of parameters
1686//---------------------------------------------------------------------------
1687//--- Step 1
1688   list re=findMinor(W);
1689   list L;
1690   int ii,i,j,l,k;
1691   J=C;
1692   ideal D=ideal(1);
1693   intvec v,w;
1694   ideal V;
1695   poly m;
1696
1697   for(ii=1;ii<=size(re);ii++)
1698   {
1699      C=J;
1700      L=re[ii];
1701      matrix A=L[1];  //(1/m)*A is the inverse matrix of the Jacobian of W
1702                      //corresponding to m
1703      m=L[2];  //a k- minor of jacob(W), not identically zero
1704                      //k=n-dim(W)
1705      V=L[3];  //the elements of W corresponding to m
1706      v=L[4];  //the indices of variables corresponding to m
1707      w=L[5];  //the indices of variables not corresponding to m
1708
1709//--- Step 2
1710//--- first some initialization depending on results of step 1
1711      k=size(V);
1712      matrix dg[1][k];
1713      matrix df[k][1];
1714//--- derivatives of the generators of J w.r.t. system of parameters
1715      for(i=1;i<=size(J);i++)
1716      {
1717         for(j=1;j<=n-k;j++)
1718         {
1719            for(l=1;l<=k;l++)
1720            {
1721               dg[1,l]=diff(V[l],var(w[j]));
1722               df[l,1]=diff(J[i],var(v[l]));
1723             }
1724             C=C,m*diff(J[i],var(w[j]))-dg*A*df;
1725         }
1726      }
1727//--- everything should live in W, not just in the intersection of
1728//--- D(m) with W
1729      C=C+W;
1730      C=sat(C,m)[1];
1731//--- intersect ideal with previously computed ones to make sure that no
1732//--- components are lost
1733      D=intersect(D,C);
1734      kill dg,df,A;
1735   }
1736//--- return minimal set of generators of the result
1737   list li=mstd(D);
1738   D=li[2];
1739   if(size(li[1])<=size(D)){D=li[1];}
1740   return(D);
1741}
1742example
1743{ "EXAMPLE:";
1744   echo = 2;
1745   ring R=0,(x,y,z),dp;
1746
1747   ideal W=z^2-x;
1748   ideal J=x*y^2+x^3;
1749   intvec b=1;
1750   list E;
1751   ideal abb=maxideal(1);
1752   intvec v;
1753   intvec w=-1;
1754   matrix M;
1755
1756   list BO=W,J,b,E,abb,v,w,M;
1757
1758   Delta(BO);
1759}
1760
1761//////////////////////////////////////////////////////////////////////////////
1762static proc redMax(int k,ideal J)
1763"Internal procedure - no help and no example available
1764"
1765{
1766//--- reduce maxideal(k) by J, more efficient approach
1767    int i;
1768    ideal K=simplify(reduce(maxideal(1),J),2);
1769    for(i=2;i<=k;i++)
1770    {
1771       K=simplify(reduce(K*maxideal(1),J),2);
1772    }
1773    return(K);
1774}
1775
1776//////////////////////////////////////////////////////////////////////////////
1777static proc findMinor(ideal W)
1778"Internal procedure - no help and no example available
1779"
1780{
1781//---------------------------------------------------------------------------
1782// Initialization and sanity checks
1783//---------------------------------------------------------------------------
1784   list L;
1785   intvec v,w;
1786   ideal Wstd=std(W);
1787   int n=nvars(basering);   // total number of columns of Jacobian
1788   int k=n-dim(Wstd);       // size of minors of Jacobian
1789   int a=size(W);           // total number of rows of Jacobian
1790   matrix A[k][k];
1791   list LW=indexSet(a,k);   // set of tuples of k rows
1792   list LV=indexSet(n,k);   // set of tuples of k columns
1793   ideal IW,IV;
1794   int i,j,l,e;
1795   list re;
1796//---------------------------------------------------------------------------
1797// We need to know which minor corresponds to which variable and to which
1798// generator of W - therefore we cannot use the function minor()!
1799//---------------------------------------------------------------------------
1800//--- choose the generators which we want to differentiate
1801   for(i=1;i<=size(LW);i++)
1802   {
1803      IW=0;
1804      for(l=1;l<=a;l++)
1805      {
1806         if(LW[i][l]!=0){IW[size(IW)+1]=W[l];}
1807      }
1808//--- choose the variables by which to differentiate and apply diff
1809      for(j=1;j<=size(LV);j++)
1810      {
1811          IV=0;v=0;w=0;
1812          for(l=1;l<=n;l++)
1813          {
1814             if(LV[j][l]!=0)
1815             {
1816                v[size(v)+1]=l;
1817                IV[size(IV)+1]=var(l);
1818             }
1819             else
1820             {
1821                w[size(w)+1]=l;
1822
1823             }
1824          }
1825          A=diff(IV,IW);      // appropriate submatrix of Jacobian
1826//--- if the minor is non-zero, then it might be the one we need
1827//--- ==> put it in the list of candidates
1828          if(det(A)!=0)
1829          {
1830            v=v[2..size(v)];  // first entry is zero for technical reasons
1831            w=w[2..size(w)];  // first entry is zero for technical reasons
1832            L=inverse_L(A);
1833            L[3]=IW;
1834            L[4]=v;
1835            L[5]=w;
1836            re[size(re)+1]=L;
1837          }
1838      }
1839   }
1840//---------------------------------------------------------------------------
1841// return the result
1842//---------------------------------------------------------------------------
1843   return(re);
1844}
1845
1846/////////////////////////////////////////////////////////////////////////////
1847static proc indexSet(int a, int b)
1848"Internal procedure - no help and no example available
1849"
1850{
1851//---------------------------------------------------------------------------
1852// Find all tuples of size b containing pairwise distict elements from a
1853// list of a elements
1854//---------------------------------------------------------------------------
1855//**************************************************************************/
1856// Combinatorics is expensive in an interpreted language
1857// ==> shift it into the kernel
1858//**************************************************************************/
1859   def R=basering;
1860   list L;
1861   ring S=2,x(1..a),dp;
1862   ideal I=maxideal(b);
1863   int i;
1864   ideal J=x(1)^2;
1865   for(i=2;i<=a;i++){J=J,x(i)^2;}
1866   attrib(J,"isSB",1);
1867   I=reduce(I,J);
1868   I=simplify(I,2);
1869   for(i=1;i<=size(I);i++){L[i]=leadexp(I[i]);}
1870   setring R;
1871   return(L);
1872}
1873
1874/////////////////////////////////////////////////////////////////////////////
1875
1876proc DeltaList(list BO)
1877"USAGE:  DeltaList (BO);
1878@*       BO = basic object, a list: ideal W,
1879@*                                  ideal J,
1880@*                                  intvec b,
1881@*                                  list Ex,
1882@*                                  ideal ab,
1883@*                                  intvec v,
1884@*                                  intvec w,
1885@*                                  matrix M
1886ASSUME:  R  = basering, a polynomial ring, W an ideal of R,
1887@*       J  = ideal containing W
1888COMPUTE: Delta-operator iteratively applied to J in the notation of
1889         [Bravo,Encinas,Villamayor]
1890RETURN:  list l of length ((max w-ord) * b),
1891         l[i+1]=Delta^i(J)
1892EXAMPLE: example DeltaList;  shows an example
1893"
1894{
1895//----------------------------------------------------------------------------
1896// Iteratively apply proc Delta
1897//----------------------------------------------------------------------------
1898   int i;
1899   list L;
1900   ideal C=BO[2];
1901   while(deg(C[1])!=0)
1902   {
1903      L[size(L)+1]=C;
1904      C=Delta(BO);
1905      BO[2]=C;
1906   }
1907   return(L);
1908}
1909example
1910{
1911   "EXAMPLE:";
1912   echo = 2;
1913   ring R=0,(x,y,z),dp;
1914
1915   ideal W=z^2-x;
1916   ideal J=x*y^2+x^3;
1917   intvec b=1;
1918   list E;
1919   ideal abb=maxideal(1);
1920   intvec v;
1921   intvec w=-1;
1922   matrix M;
1923
1924   list BO=W,J,b,E,abb,v,w,M;
1925
1926   DeltaList(BO);
1927}
1928/////////////////////////////////////////////////////////////////////////////
1929proc CenterBM(list BM)
1930"USAGE:  CenterBM(BM);
1931@*       BM = object related to a presentation,
1932                            a list: ideal W,
1933@*                                  ideal J,
1934@*                                  intvec b,
1935@*                                  list Ex,
1936@*                                  ideal ab,
1937@*                                  intvec v,
1938@*                                  intvec w,
1939@*                                  matrix M
1940ASSUME:  R  = basering, a polynomial ring, W an ideal of R,
1941@*       J  = ideal containing W
1942COMPUTE: the center of the next blow-up of BM in the resolution algorithm
1943         of [Bierstone, Milman]
1944RETURN:  list l,
1945         l[1]: ideal describing the center
1946         l[2]: intvec w obtained in the process of determining l[1]
1947         l[3]: intvec b obtained in the process of determining l[1]
1948         l[4]: intmat invmat obtained in the process of determining l[1]
1949EXAMPLE: example CenterBM;  shows an example
1950"
1951{
1952//!!! NOCH NICHT IN BETRIEB
1953ERROR("Not implemented yet");
1954   int i,j;
1955   intmat tmat[2][1]=0,-1;
1956//--- re=center,E^- indices, b vector, n vector
1957   list re=ideal(1),BM[7],BM[3],tmat;
1958   ideal J=BM[2];
1959   if(size(J)==0)
1960   {
1961     re[1]=ideal(0);
1962     return(re);
1963   }
1964//--- find Delta^(b-1)(J)
1965   if(size(reduce(J,std(BM[1])))!=0)
1966   {
1967     list L=DeltaList(BM);
1968   }
1969   else
1970   {
1971     list L;
1972     L[1]=J;
1973   }
1974   if(!defined(debugCenter))
1975   {
1976     int debugCenter;
1977   }
1978   if(debugCenter)
1979   {
1980     "----> In Center: after DeltaList";
1981     "W";
1982     BM[1];
1983     "J";
1984     BM[2];
1985     "The Delta List:";
1986     L;
1987   }
1988   int b=size(L);
1989   if(b==0)
1990   {
1991//--- if J=W, we do not need to do anything
1992//--- returning center=1 marks this chart as completed
1993     return(re);
1994   }
1995//---------------------------------------------------------------------------
1996// check whether max ord is constant
1997//---------------------------------------------------------------------------
1998   if((BM[9][2,1]<0)||(BM[9][1,1]>b))
1999   {
2000//--- we are either at the beginning or the invariant has dropped
2001     intvec tempvec=size(BM[4]);
2002     BM[7]=tempvec;
2003//!!!! nur fuer hyperflaechen!!!!!!!!
2004     tempvec=b;
2005     BM[3]=tempvec;
2006//!!!! Ende !!!!!!!!!!!!
2007     kill tempvec;
2008     BM[9][1,1]=b;
2009     BM[9][2,1]=1;
2010   }
2011//---------------------------------------------------------------------------
2012// prepare for intersection with E_i
2013//---------------------------------------------------------------------------
2014   ideal C=L[b];
2015   re[2]=BM[7];
2016   re[3]=BM[3];
2017   BM[2]=C;
2018      if(debugCenter)
2019   {
2020     "----> In Center: before intersection with E_i:";
2021     "bmax:",b;
2022     "Sing(J,bmax):";
2023      C;
2024     "E:";
2025      BO[4];
2026     "list marking a priori known intersection properties:",BO[6];
2027     "index of last element of E^- in E:",BO[7][1];
2028   }
2029   list E=inters_E(BM);
2030// !!!!!!!!! Drop Redundant fehlt noch!!!!!!!!
2031//---------------------------------------------------------------------------
2032// Check whether it is a single point
2033//---------------------------------------------------------------------------
2034   ideal C1=std(ideal(L[b])+E[1]);
2035   if(dim(C1)==0)
2036   {
2037     if(size(E[4])==1)
2038     {
2039        tmat[1,1]=BM[9][1,1];
2040        tmat[2,1]=BM[9][2,1];
2041        re[4]=tmat;
2042        re[1]=radical(C1);
2043        return(re);
2044     }
2045   }
2046   if(size(BM[9])>2)
2047   {
2048      BM[9][1,2]=E[2];
2049      BM[9][2,2]=1;
2050   }
2051   else
2052   {
2053      intmat tempInt[2][1]=E[2],1;
2054      BM[9]=concatInt(BM[9],tempInt);
2055      kill tempInt;
2056   }
2057   BM[2]=J;
2058   list BM1=dropDim(BM);
2059   list BMlist,hilfList;
2060   ideal hilf;
2061   intvec tempvec;
2062
2063   if(!attrib(BM1[1],"isSB")){BM1[1]=std(BM1[1]);}
2064   for(i=1;i<=size(BM1[4]);i++)
2065   {
2066      hilf=simplify(reduce(BM1[4][i],BM1[1]),2);
2067      if(size(hilf)>1){"Problem with BM1[4]in CenterBM";~;}
2068      hilfList[i]=hilf[1];
2069   }
2070   for(i=1;i<=size(E[4]);i++)
2071   {
2072      BMlist[i]=BM1;
2073      for(j=size(E[5][i]);j>=1;j--)
2074      {
2075         if(E[5][i][j]!=0)
2076         {
2077            BMlist[i][2][size(BMlist[i][2])+1]=hilfList[j];
2078            BMlist[i][3][size(BMlist[i][3])+1]=1;
2079         }
2080         BMlist[i][4]=delete(BMlist[i][4],j);
2081         BMlist[i][6]=deleteInt(BMlist[i][6],j,0);
2082      }
2083      BMlist[i][7]=deleteInt(BMlist[i][7],1,-1);
2084      if(size(BMlist[i][9])>4)
2085      {
2086         intmat tempInt[2][ncols(BMlist[i][9])]=BMlist[i][9];
2087         intmat tempInt2[2][ncols(BMlist[i][9])-2]=
2088                 tempInt[1..2,3..ncols(BMlist[i][9])];
2089         BMlist[i][9]=tempInt2;
2090         kill tempInt,tempInt2;
2091      }
2092      else
2093      {
2094         BMlist[i][9]=tmat;
2095      }
2096   }
2097   kill hilfList;list hilfList;
2098   hilfList[1]=CenterTail(BMlist[1],C);
2099   intmat maxmat=hilfList[1][9];
2100   intvec maxiv=E[4][1];
2101   int pos=1;
2102   for(i=2;i<=size(E[4]);i++)
2103   {
2104      hilfList[i]=CenterTail(BMlist[i],C);
2105      if(invGreater(hilfList[i][4]),maxmat,E[4][i],maxiv)
2106      {
2107         maxmat=hilfList[i][4];
2108         maxiv=E[4][i];
2109         pos=i;
2110      }
2111   }
2112   re[1]=hilfList[pos][1];
2113   intmat tempint=BM[9];
2114   intmat tempint1[2][2]=tempint[1..2,1..2];
2115   re[4]=concatInt(tempint1,maxmat);
2116   re[2]=re[2][1],hilfList[pos][2];
2117   re[3]=b;
2118~;
2119   return(re);
2120}
2121example
2122{  "EXAMPLE:";
2123   echo = 2;
2124   ring R=0,(x,y),dp;
2125
2126   ideal W;
2127   ideal J=x2-y3;
2128   intvec b=1;
2129   list E;
2130   ideal abb=maxideal(1);
2131   intvec v;
2132   intvec w=-1;
2133   matrix M;
2134   intmat invmat[2][1]=0,-1;
2135
2136   list BM=W,J,b,E,abb,v,w,M,invmat;
2137
2138   CenterBM(BM);
2139}
2140
2141/////////////////////////////////////////////////////////////////////////////
2142static proc invGreater(intmat M1, intmat M2, intvec iv1, intvec iv2)
2143{
2144// Auxilliary procedure, BM-algorithm
2145   int i;
2146   for(i=1;i<=min(ncols(M1),ncols(M2));i++)
2147   {
2148      if(M1[2,i]==-1)
2149      {
2150         if(M1[1,i]==0){ERROR("Invariant not set");}
2151         if(M2[2,i]!=-1){return(1);}
2152         if(M2[1,i]==0){ERROR("Invariant not set");}
2153         break;
2154      }
2155      else
2156      {
2157         if(M2[2,i]==-1)
2158         {
2159            if(M2[1,i]==0){ERROR("Invariant not set");}
2160            return(0);
2161         }
2162         if(M1[1,i]*M2[2,i]!= M2[1,i]*M1[2,i])
2163         {
2164            return(M1[1,i]*M2[2,i]> M2[1,i]*M1[2,i]);
2165         }
2166      }
2167   }
2168   return(iv1>iv2);
2169}
2170/////////////////////////////////////////////////////////////////////////////
2171proc CenterTail(list BM, ideal C)
2172{
2173//!!! Auxilliary procedure, BM-algorithm
2174//!!!!!!!!Rueckgabe im Zentrumsformat
2175  int i,j,bmin;
2176  int alpha=lcm(BM[3]);
2177  vector w;
2178  list re;
2179  if(size(BM[2])==0)
2180  {
2181    re[1]=C+BM[1];
2182    intvec tvec;
2183    re[3]=tvec;
2184    intmat tmat[2][1]=-1,-1;
2185    re[4]=tmat;
2186    tvec=size(BM[4]);
2187    re[2]=tvec;
2188    return(re);
2189  }
2190  for(i=1;i<=size(BM[3]);i++)
2191  {
2192     if(BM[2][i]!=0)
2193     {
2194       w[size(w)+1]=BM[2][i]^(alpha/BM[3][i]);
2195     }
2196  }
2197  module M=w;
2198  intvec satex;
2199  list satList;
2200  for(i=1;i<=size(BM[4]);i++)
2201  {
2202     satList=sat(M,BM[4][i]);
2203     satex[i]=satList[2];
2204     M=satList[1];
2205  }
2206//!!!!Hilfsobjekt G bilden!!!!!!!!!!!!!!!!
2207//!!!!Hilfsobjekt H,codim -1 bilden!!!!!!!!!!!!!!!!
2208//!!!! ???? an welcher stelle?????????
2209  list deltaL;
2210  list BMtemp=BM;
2211  for(i=1;i<=nrows(M[1]);i++)
2212  {
2213     BMtemp[2]=ideal(M[1][i])+BMtemp[1];
2214     deltaL[i]=DeltaList(BMtemp);
2215     if(i==1)
2216     {
2217        bmin=size(deltaL[i]);
2218     }
2219     else
2220     {
2221        bmin=min(size(deltaL[i]),bmin);
2222     }
2223  }
2224  if(bmin==0)
2225  {
2226     re[1]=C+BM[1];
2227     intvec tvec;
2228     re[3]=tvec;
2229     if((BM[9][2,1]==-1)||(BM[9][1,1]!=0))
2230     {
2231        tvec=size(BM[4]);
2232        re[2]=tvec;
2233        intmat tmat[2][1]=0,1;
2234        re[4]=tmat;
2235     }
2236     else
2237     {
2238        re[2]=BM[7];
2239        re[4]=BM[9];
2240     }
2241     return(re);
2242  }
2243  ideal Ctemp=ideal(1);
2244  while(deg(Ctemp[1])==0)
2245  {
2246    Ctemp=C;
2247    for(i=1;i<=nrows(M[1]);i++)
2248    {
2249       Ctemp=Ctemp,deltaL[i][bmin];
2250    }
2251    Ctemp=std(Ctemp);
2252    bmin--;
2253    if(bmin==0){ERROR("empty set");}
2254  }
2255  //!!!!!!!!!!!!!Invariante ist bmin/alpha
2256  // naechster Eintrag s_i wie in CenterBM
2257  // dann dropDim  ....
2258}
2259
2260/////////////////////////////////////////////////////////////////////////////
2261static proc deleteInt(intvec v,int i,int ini)
2262{
2263//!!! Should be in kernel of Singular
2264//--- delete i-th entry in intvec v,
2265//--- if necessary reinitializing v with value ini
2266   int s=size(v);
2267   intvec w;
2268   if((i<s)&&(i>1)){w=v[1..i-1],v[i+1..s];}
2269   if(s==1){w=ini;return(w);}
2270   if(i==1){w=v[2..s];}
2271   if(i==s){w=v[1..s-1];}
2272   return(w);
2273}
2274/////////////////////////////////////////////////////////////////////////////
2275static proc concatInt(intmat A, intmat B)
2276{
2277//!!! Should be in kernel of Singular
2278//--- concatenate two intmats
2279  if(nrows(A)!=nrows(B)){ERROR("could not concat, wrong number of rows");}
2280  intmat tempmat[nrows(A)][ncols(A)+ncols(B)];
2281  tempmat[1..nrows(A),1..ncols(A)]=A[1..nrows(A),1..ncols(A)];
2282  tempmat[1..nrows(A),ncols(A)+1..ncols(tempmat)]=B[1..nrows(A),1..ncols(B)];
2283  return(tempmat);
2284}
2285/////////////////////////////////////////////////////////////////////////////
2286proc dropDim(list BM)
2287{
2288ERROR("Not implemented yet");
2289}
2290/////////////////////////////////////////////////////////////////////////////
2291proc CenterBO(list BO,list #)
2292"USAGE:  CenterBO(BO);
2293@*       BO = basic object, a list: ideal W,
2294@*                                  ideal J,
2295@*                                  intvec b,
2296@*                                  list Ex,
2297@*                                  ideal ab,
2298@*                                  intvec v,
2299@*                                  intvec w,
2300@*                                  matrix M
2301ASSUME:  R  = basering, a polynomial ring, W an ideal of R,
2302@*       J  = ideal containing W
2303COMPUTE: the center of the next blow-up of BO in the resolution algorithm
2304         of [Bravo,Encinas,Villamayor]
2305RETURN:  list l,
2306         l[1]: ideal describing the center@*
2307         l[2]: intvec w obtained in the process of determining l[1]@*
2308         l[3]: intvec b obtained in the process of determining l[1]@*
2309         l[4]: intvec inv obtained in the process of determining l[1]
2310EXAMPLE: example CenterBO;  shows an example
2311"
2312{
2313//---------------------------------------------------------------------------
2314// Initialization  and sanity checks
2315//---------------------------------------------------------------------------
2316   int i,bo7save;
2317   intvec tvec;
2318//--- re=center,E^- indices, b vector, n vector
2319   list re=ideal(1),BO[7],BO[3],tvec;
2320   ideal J=BO[2];
2321   if(size(J)==0)
2322   {
2323     re[1]=ideal(0);
2324     return(re);
2325   }
2326//--- find Delta^(b-1)(J)
2327   if(size(reduce(J,std(BO[1])))!=0)
2328   {
2329     list L=DeltaList(BO);
2330   }
2331   else
2332   {
2333     list L;
2334     L[1]=J;
2335   }
2336   if(!defined(debugCenter))
2337   {
2338     int debugCenter;
2339   }
2340   if(debugCenter)
2341   {
2342     "----> In Center: after DeltaList";
2343     "W";
2344     BO[1];
2345     "J";
2346     BO[2];
2347     "The Delta List:";
2348     L;
2349   }
2350   int b=size(L);
2351   if(b==0)
2352   {
2353//--- if J=W, we do not need to do anything
2354//--- returning center=1 marks this chart as completed
2355     return(re);
2356   }
2357//---------------------------------------------------------------------------
2358// check whether max w-ord is constant
2359//---------------------------------------------------------------------------
2360   if(b==BO[3][1])
2361   {
2362//--- max w-ord is constant
2363     if(BO[7][1]==-1)
2364     {
2365//--- first got its value in the previous step ==> initialize BO[7]
2366       tvec[1]=size(BO[4])-1;
2367       for(i=2;i<=size(BO[7]);i++)
2368       {
2369         tvec[i]=BO[7][i];
2370       }
2371       re[2]=tvec;
2372       BO[7]=tvec;
2373     }
2374   }
2375   else
2376   {
2377//--- max w-ord changed ==> reset BO[7], correct BO[3]
2378     tvec[1]=-1;
2379     re[2]=tvec;
2380     BO[7]=tvec;
2381     tvec[1]=b;
2382     BO[3]=tvec;
2383     if(defined(invSat))
2384     {
2385        invSat[2]=intvec(0);
2386     }
2387   }
2388   re[3]=BO[3];
2389//---------------------------------------------------------------------------
2390// reduce from case 2 to case 1 of [Bravo, Encinas, Villamayor]
2391//---------------------------------------------------------------------------
2392   ideal C=L[b];
2393   BO[2]=C;
2394   if(debugCenter)
2395   {
2396     "----> In Center: before intersection with E_i:";
2397     "bmax:",b;
2398     "Sing(J,bmax):";
2399      C;
2400     "E:";
2401      BO[4];
2402     "list marking a priori known intersection properties:",BO[6];
2403     "index of last element of E^- in E:",BO[7][1];
2404   }
2405//--- is intermediate result in iteration already good?
2406//--- return it to calling proc CenterBO
2407   if(size(#)>0)
2408   {
2409     if(#[1]==2)
2410     {
2411       re[1]=C;
2412       kill tvec;
2413       intvec tvec=re[2][1];
2414       re[2]=tvec;
2415       kill tvec;
2416       intvec tvec=re[3][1];
2417       re[3]=tvec;
2418       kill tvec;
2419       intvec tvec;
2420       re[4]=tvec;
2421       return(re);
2422     }
2423   }
2424//--- do the reduction to case 1
2425   list E=inters_E(BO);
2426//--- if J is smooth and not too many E_i intersect simultaneously, let us
2427//--- try to drop redundant components of the candidate for the center
2428   if((b==1)&&(size(E)>3))
2429   {
2430//--- if J is not smooth we do not want to drop any information
2431     if((size(E[4])>0) && (dim(std(slocusE(BO[2])))<0))
2432     {
2433//--- BO[2]==J because b==1
2434//--- DropRedundant is the counterpart to DropCoeff
2435//--- do not leave out one of them separately!!!
2436      E=DropRedundant(BO,E);
2437      if(size(E)==1)
2438       {
2439         kill tvec;
2440         intvec tvec=re[2][1];
2441         re[2]=tvec;
2442         tvec[1]=re[3][1];
2443         re[3]=tvec;
2444         tvec[1]=re[4][1];
2445         re[4]=tvec;
2446         re[1]=E[1];
2447         return(re);
2448       }
2449     }
2450   }
2451//--- set n correctly
2452   if(E[2]<BO[9][1])
2453   {
2454//--- if n dropped, the subsequent Coeff-object will not be the same
2455//--- ===> set BO[3][2] to zero to make sure that no previous data is used
2456      if(defined(tvec)) {kill tvec;}
2457      intvec tvec=BO[7][1],-1;
2458      BO[7]=tvec;
2459      tvec=BO[3][1],0;
2460      BO[3]=tvec;
2461      if(defined(invSat))
2462      {
2463         invSat[2]=intvec(0);
2464      }
2465   }
2466   re[4][1]=E[2];
2467   C=E[1]^b+J;
2468   C=mstd(C)[2];
2469   ideal C1=std(ideal(L[b])+E[1]);
2470   if(debugCenter)
2471   {
2472     "----> In Center: reduction of case 2 to case 1";
2473     "Output of inters_E, after dropping redundant components:";
2474      E;
2475     "result of intersection with E^-, i.e.(E^-)^b+J:";
2476     C;
2477     "candidate for center:";
2478     C1;
2479   }
2480//---------------------------------------------------------------------------
2481// Check whether we have a hypersurface component
2482//---------------------------------------------------------------------------
2483   if(dim(C1)==dim(std(BO[1]))-1)
2484   {
2485     if((size(reduce(J,C1))==0)&&(size(reduce(C1,std(J)))==0))
2486     {
2487//--- C1 equals J and is of codimension 1 in W
2488        re[1]=C1;
2489     }
2490     else
2491     {
2492//--- C1 has a codimension 1 (in W) component
2493        re[1]=std(equiRadical(C1));
2494     }
2495     kill tvec;
2496     intvec tvec=re[2][1];
2497     re[2]=tvec;
2498     tvec[1]=re[3][1];
2499     re[3]=tvec;
2500     tvec[1]=re[4][1];
2501     re[4]=tvec;
2502//--- is the codimension 1 component a good choice or do we need to reset
2503//--- the information from the previous steps
2504     if(transversalT(re[1],BO[4]))
2505     {
2506        if(size(E)>2)
2507        {
2508          if(E[3]>E[2])
2509          {
2510            if(defined(shortcut)){kill shortcut;}
2511            list shortcut=ideal(0),size(BO[4]),BO[7];
2512            export(shortcut);
2513          }
2514        }
2515        return(re);
2516     }
2517
2518     ERROR("reset in Center, please send the example to the authors.");
2519   }
2520//---------------------------------------------------------------------------
2521// Check whether it is a single point
2522//---------------------------------------------------------------------------
2523   if(dim(C1)==0)
2524   {
2525      C1=std(radical(C1));
2526      if(vdim(C1)==1)
2527      {
2528//--- C1 is one point
2529        re[1]=C1;
2530        kill tvec;
2531        intvec tvec=re[2][1];
2532        re[2]=tvec;
2533        kill tvec;
2534        intvec tvec=re[3][1];
2535        re[3]=tvec;
2536        return(re);
2537      }
2538    }
2539//---------------------------------------------------------------------------
2540// Prepare input for forming the Coeff-Ideal
2541//---------------------------------------------------------------------------
2542   BO[2]=C;
2543   if(size(BO[2])>5)
2544   {
2545      BO[2]=mstd(BO[2])[2];
2546   }
2547//--- drop leading entry of BO[3]
2548   tvec=BO[3];
2549   if(size(tvec)>1)
2550   {
2551     tvec=tvec[2..size(tvec)];
2552     BO[3]=tvec;
2553   }
2554   else
2555   {
2556     BO[3][1]=0;
2557   }
2558   tvec=BO[9];
2559   if(size(tvec)>1)
2560   {
2561     tvec=tvec[2..size(tvec)];
2562     BO[9]=tvec;
2563   }
2564   else
2565   {
2566     BO[9][1]=0;
2567   }
2568   bo7save=BO[7][1];          // original value needed for result
2569   if(defined(shortcut))
2570   {
2571     if((bo7save!=shortcut[3][1])&&(size(shortcut[3])!=1))
2572     {
2573       kill shortcut;
2574     }
2575     else
2576     {
2577       shortcut[2]=shortcut[2]-bo7save;
2578       tvec=shortcut[3];
2579       if(size(tvec)>1)
2580       {
2581         tvec=tvec[2..size(tvec)];
2582         shortcut[3]=tvec;
2583       }
2584       else
2585       {
2586         shortcut[3]=intvec(shortcut[2]);
2587       }
2588     }
2589   }
2590   if(BO[7][1]>-1)
2591   {
2592//--- drop E^- and the corresponding information from BO[6]
2593     for(i=1;i<=BO[7][1];i++)
2594     {
2595       BO[4]=delete(BO[4],1);
2596       intvec bla1=BO[6];
2597       BO[6]=intvec(bla1[2..size(bla1)]);
2598       kill bla1;
2599     }
2600//--- drop leading entry of BO[7]
2601     tvec=BO[7];
2602     if(size(tvec)>1)
2603     {
2604       tvec=tvec[2..size(tvec)];
2605       BO[7]=tvec;
2606     }
2607     else
2608     {
2609       BO[7][1]=-1;
2610     }
2611   }
2612   else
2613   {
2614     if(BO[7][1]==-1)
2615     {
2616       list tplist;
2617       BO[4]=tplist;
2618       kill tplist;
2619     }
2620   }
2621   if(debugCenter)
2622   {
2623     "----> In Center: Input to Coeff";
2624     "b:",b;
2625     "BO:";
2626     BO;
2627   }
2628//--- prepare the third entry of the invariant tuple
2629   int invSatSave=invSat[2][1];
2630   tvec=invSat[2];
2631   if(size(tvec)>1)
2632   {
2633     tvec=tvec[2..size(tvec)];
2634     invSat[2]=tvec;
2635   }
2636   else
2637   {
2638     invSat[2][1]=0;
2639   }
2640//---------------------------------------------------------------------------
2641// Form the Coeff-ideal, if possible and useful; otherwise use the previous
2642// candidate for the center
2643//---------------------------------------------------------------------------
2644   list BO1=Coeff(BO,b);
2645   if(debugCenter)
2646   {
2647     "----> In Center: Output of Coeff";
2648     BO1;
2649   }
2650//--- Coeff returns int if something went wrong
2651   if(typeof(BO1[1])=="int")
2652   {
2653      if(BO1[1]==0)
2654      {
2655//--- Coeff ideal was already resolved
2656        re[1]=C1;
2657        return(re);
2658      }
2659      else
2660      {
2661//--- no global hypersurface found
2662        re=CoverCenter(BO,b,BO1[2]);
2663        kill tvec;
2664        intvec tvec=invSatSave;
2665        for(i=1;i<=size(invSat[2]);i++)
2666        {
2667           tvec[i+1]=invSat[2][i];
2668        }
2669        invSat[2]=tvec;
2670        return(re);
2671      }
2672   }
2673   int coeff_invar;
2674   ideal Idropped=1;
2675//--- if b=1 drop redundant components of the Coeff-ideal
2676   if(b==1)
2677   {
2678//--- Counterpart to DropRedundant -- do not leave out one of them separately
2679     Idropped=DropCoeff(BO1);         // blow-up in these components
2680                                      // is unnecessary
2681   }
2682//--- to switch off DropCoeff, set Idropped=1;
2683   BO1[2]=sat(BO1[2],Idropped)[1];
2684   if(deg(BO1[2][1])==0)
2685   {
2686//--- Coeff ideal is trivial
2687     C1=radical(C1);
2688     ideal C2=sat(C1,Idropped)[1];
2689     if(deg(std(C2)[1])!=0)
2690     {
2691        C1=C2;
2692     }
2693//Aenderung: Strategie: nur im Notfall ganze except. Divisoren
2694     if(deg(std(BO1[2])[1])==0)
2695     {
2696       list BOtemp=BO;
2697       int bo17save=BO1[7][1];
2698       BOtemp[7]=0;
2699       BO1=Coeff(BOtemp,b,int(0));
2700       BO1[2]=sat(BO1[2],Idropped)[1];
2701       if(deg(std(BO1[2])[1])==0)
2702       {
2703//--- there is really nothing left to do for the Coeff ideal
2704//--- the whole original BO1[2], i.e. Idropped, is the upcoming center
2705          re[1]=Idropped;
2706          re[2]=intvec(bo7save);
2707          re[3]=intvec(b);
2708          re[4]=intvec(E[2]);
2709          return(re);
2710       }
2711       if(deg(std(slocus(radical(BO1[2])))[1])==0)
2712       {
2713          re[1]=BO1[2];
2714//          re[2]=intvec(bo7save,BO1[7][1]);
2715          re[2]=intvec(bo7save,bo17save);
2716          re[3]=intvec(b,1);
2717          re[4]=intvec(E[2],1);
2718          invSat[2]=intvec(1,0);
2719          return(re);
2720       }
2721//!!! effizienter machen???
2722       list pr=primdecGTZ(BO1[2]);
2723       ideal Itemp1=1;
2724       int aa,bb;
2725       for(aa=1;aa<=size(pr);aa++)
2726       {
2727          if(dim(std(pr[aa][2])) < (dim(std(BO1[1]))-1))
2728          {
2729//--- drop components which are themselves exceptional diviosrs
2730             Itemp1=intersect(Itemp1,pr[aa][1]);
2731          }
2732       }
2733       if(deg(std(Itemp1)[1])!=0)
2734       {
2735//--- treat the remaining components of the weak Coeff ideal
2736          BO1[2]=Itemp1;
2737       }
2738       BO1[7]=BO[7];
2739       for(aa=1;aa<=size(BO1[4]);aa++)
2740       {
2741         if(deg(std(BO1[4][aa])[1])==0){aa++;continue;}
2742         if(defined(satlist)){kill satlist;}
2743         list satlist=sat(BO1[2],BO1[4][aa]+BO1[1]);
2744         if(deg(std(satlist[1])[1])==0)
2745         {
2746            coeff_invar++;
2747            if(satlist[2]!=0)
2748            {
2749              for(bb=1;bb<=satlist[2]-1;bb++)
2750              {
2751                 BO1[2]=quotient(BO1[2],BO1[4][aa]+BO1[1]);
2752              }
2753            }
2754            else
2755            {
2756              ERROR("J of temporary object had unexpected value;
2757                     please send this example to the authors.");
2758            }
2759         }
2760         else
2761         {
2762            BO1[2]=satlist[1];
2763         }
2764       }
2765       if(deg(std(Itemp1)[1])==0)
2766       {
2767          re[1]=BO1[2];
2768          re[2]=intvec(bo7save,BO1[7][1]);
2769          re[3]=intvec(b,1);
2770          re[4]=intvec(E[2],1);
2771          invSat[2]=intvec(1,0);
2772          return(re);
2773       }
2774       kill aa,bb;
2775     }
2776   }
2777   if(invSatSave<coeff_invar)
2778   {
2779     invSatSave=coeff_invar;
2780   }
2781//---------------------------------------------------------------------------
2782// Determine Center of Coeff-ideal and use it as the new center
2783//---------------------------------------------------------------------------
2784   if(!defined(templist))
2785   {
2786     if(size(BO1[2])>5)
2787     {
2788        BO1[2]=mstd(BO1[2])[2];
2789     }
2790     list templist=CenterBO(BO1,2);
2791//--- only a sophisticated guess of a good center computed by
2792//--- leaving center before intersection with the E_i.
2793//--- whether the guess was good, is stored in 'good'.
2794//--- (this variant saves charts in cases like the Whitney umbrella)
2795     list E0,E1;
2796
2797     ideal Cstd=std(radical(templist[1]));
2798     int good=((deg(std(slocusE(Cstd))[1])==0)&&(dim(std(BO1[2]))<=2));
2799//if(defined(satlist)){good=0;}
2800     if(good)
2801     {
2802        for(i=1;i<=size(BO[4]);i++)
2803        {
2804           if((deg(BO[4][i][1])>0)&&(size(reduce(BO[4][i],Cstd))!=0))
2805           {
2806              E0[size(E0)+1]=BO[4][i]+Cstd;
2807              E1[size(E1)+1]=BO[4][i];
2808           }
2809        }
2810        good=transversalT(Cstd,E1);
2811        if(good)
2812        {
2813           good=normalCross(E0);
2814        }
2815     }
2816     if(good)
2817     {
2818        list templist2=CenterBO(BO1,1);
2819        if(dim(std(templist2[1]))!=dim(Cstd))
2820        {
2821           templist[1]=Cstd;
2822           if(defined(shortcut)){kill shortcut;}
2823           list shortcut=ideal(0),size(BO1[4]),templist[2];
2824           export(shortcut);
2825        }
2826        else
2827        {
2828           templist=templist2;
2829        }
2830        kill templist2;
2831     }
2832     else
2833     {
2834//--- sophisticated guess was wrong, follow Villamayor's approach
2835        kill templist;
2836        list templist=CenterBO(BO1,1);
2837     }
2838   }
2839   if((dim(std(templist[1]))==dim(std(BO1[1]))-1)
2840       &&(size(templist[4])==1))
2841   {
2842     if(templist[4][1]==0)
2843     {
2844        for(i=1;i<=size(BO1[4]);i++)
2845        {
2846           if(size(reduce(templist[1],std(BO1[4][i])))==0)
2847           {
2848              templist[4][1]=1;
2849              break;
2850           }
2851        }
2852     }
2853   }
2854//!!! subsequent line should be deleted
2855//if(defined(satlist)){templist[3][1]=BO[3][1];}
2856   if(debugCenter)
2857   {
2858     "----> In Center: Iterated Center returned:";
2859     templist;
2860   }
2861//--------------------------------------------------------------------------
2862// set up the result and return it
2863//--------------------------------------------------------------------------
2864   re[1]=templist[1];
2865   kill tvec;
2866   intvec tvec;
2867   tvec[1]=bo7save;
2868   for(i=1;i<=size(templist[2]);i++)
2869   {
2870     tvec[i+1]=templist[2][i];
2871   }
2872   re[2]=tvec;
2873   if(defined(shortcut))
2874   {
2875      shortcut[2]=shortcut[2]+bo7save;
2876      shortcut[3]=tvec;
2877   }
2878   kill tvec;
2879   intvec tvec;
2880   tvec[1]=invSatSave;
2881   for(i=1;i<=size(invSat[2]);i++)
2882   {
2883      tvec[i+1]=invSat[2][i];
2884   }
2885   invSat[2]=tvec;
2886   kill tvec;
2887   intvec tvec;
2888   tvec[1]=b;
2889   for(i=1;i<=size(templist[3]);i++)
2890   {
2891     tvec[i+1]=templist[3][i];
2892   }
2893   re[3]=tvec;
2894   kill tvec;
2895   intvec tvec;
2896   tvec[1]=E[2];
2897   for(i=1;i<=size(templist[4]);i++)
2898   {
2899     tvec[i+1]=templist[4][i];
2900   }
2901   re[4]=tvec;
2902
2903   return(re);
2904}
2905example
2906{  "EXAMPLE:";
2907   echo = 2;
2908   ring R=0,(x,y),dp;
2909
2910   ideal W;
2911   ideal J=x2-y3;
2912   intvec b=1;
2913   list E;
2914   ideal abb=maxideal(1);
2915   intvec v;
2916   intvec w=-1;
2917   matrix M;
2918
2919   list BO=W,J,b,E,abb,v,w,M,v;
2920
2921   CenterBO(BO);
2922}
2923//////////////////////////////////////////////////////////////////////////////
2924static proc CoverCenter(list BO,int b, ideal Jb)
2925{
2926//----------------------------------------------------------------------------
2927// Initialization
2928//----------------------------------------------------------------------------
2929   def R=basering;
2930   int i,j,k;
2931   intvec merk,merk2,maxv,fvec;
2932   list L,ceList,re;
2933   ceList[1]=ideal(0);
2934   poly @p,@f;
2935   ideal K,dummy;
2936   if(!attrib(BO[2],"isSB"))
2937   {
2938     BO[2]=std(BO[2]);
2939   }
2940   for(i=1;i<=size(Jb);i++)
2941   {
2942      list tempmstd=mstd(slocus(Jb[i]));
2943      if(size(tempmstd[1])>size(tempmstd[2]))
2944      {
2945         dummy=tempmstd[2];
2946      }
2947      else
2948      {
2949         dummy=tempmstd[1];
2950      }
2951      kill tempmstd;
2952      L[i]=dummy;
2953      K=K,dummy;
2954   }
2955   K=simplify(K,2);
2956//---------------------------------------------------------------------------
2957// The intersection of the singular loci of the L[i] is empty.
2958// Find a suitable open covering of the affine chart, such that a global
2959// hypersurface can be found in each open set.
2960//---------------------------------------------------------------------------
2961   matrix M=lift(K,ideal(1));
2962   j=1;
2963   for(i=1;i<=nrows(M);i++)
2964   {
2965      if(M[i,1]!=0)
2966      {
2967         merk[size(merk)+1]=i;
2968         fvec[size(merk)]=j;
2969      }
2970      if((i-k)==size(L[j]))
2971      {
2972        k=i;
2973        j++;
2974      }
2975   }
2976//--------------------------------------------------------------------------
2977// Find a candidate for the center in each open set
2978//--------------------------------------------------------------------------
2979//--- first entry of merk is 0 by construction of merk
2980   for(i=2;i<=size(merk);i++)
2981   {
2982//--- open set is D(@p)
2983       @p=K[merk[i]];
2984//--- hypersurface is V(@f)
2985       @f=Jb[fvec[i]];
2986       execute("ring R1=("+charstr(R)+"),(@y,"+varstr(R)+"),dp;");
2987       poly p=imap(R,@p);
2988       poly f=imap(R,@f);
2989       list @ce;
2990       list BO=imap(R,BO);
2991       BO[1]=BO[1]+ideal(@y*p-1);
2992       BO[2]=BO[2]+ideal(@y*p-1);
2993       for(j=1;j<=size(BO[4]);j++)
2994       {
2995          BO[4][j]=BO[4][j]+ideal(@y*p-1);
2996       }
2997//--- like usual Coeff, but hypersurface is already known
2998       list BO1=SpecialCoeff(BO,b,f);
2999//--- special situation in SpecialCoeff are marked by an error code of
3000//--- type int
3001       if(typeof(BO1[1])=="int")
3002       {
3003         if(BO1[1]==0)
3004         {
3005//--- Coeff ideal was already resolved
3006           @ce[1]=BO[2];
3007           @ce[2]=BO[7];
3008           @ce[3]=BO[3];
3009         }
3010         else
3011         {
3012           if(BO[3]!=0)
3013           {
3014//--- intersections with E do not meet conditions ==> reset
3015              ERROR("reset in Coeff, please send the example to the autors");
3016           }
3017         }
3018       }
3019       else
3020       {
3021//--- now do the recursion as usual
3022         @ce=CenterBO(BO1);
3023       }
3024//---------------------------------------------------------------------------
3025// Go back to the whole affine chart and form a suitable union of the
3026// candidates
3027//---------------------------------------------------------------------------
3028//--- pass from open set to the whole affine chart by taking the closure
3029       @ce[1]=eliminate(@ce[1],@y);
3030       setring R;
3031       ceList[i]=imap(R1,@ce);
3032//--- set up invariant vector and determine maximum value of it
3033       if(size(ceList[i][3])==size(ceList[i][4]))
3034       {
3035          kill merk2,maxv;
3036          intvec merk2,maxv;
3037          for(j=1;j<=size(ceList[i][3]);j++)
3038          {
3039            merk2[2*j-1]=ceList[i][3][j];
3040            merk2[2*j]=ceList[i][4][j];
3041            ceList[i][5]=merk2;
3042            if(maxv<merk2)
3043            {
3044               maxv=merk2;
3045            }
3046          }
3047       }
3048       else
3049       {
3050          ERROR("This situation should not occur, please send the example
3051                 to the authors.");
3052       }
3053       kill R1;
3054   }
3055   kill merk2;
3056   intvec merk2=-2;
3057//--- form the union of the components of the center with maximum invariant
3058   for(i=1;i<=size(ceList);i++)
3059   {
3060     if(size(reduce(ceList[i][1],BO[2]))==0)
3061     {
3062//--- already resolved ==> ignore
3063       i++;
3064       continue;
3065     }
3066     if(ceList[i][5]==maxv)
3067     {
3068        if(merk2!=ceList[i][2])
3069        {
3070//--- E^- not of the same size as before resp. initialization
3071          if(merk2[1]==-2)
3072          {
3073//--- initialization: save size of E^-
3074            merk2=ceList[i][2];
3075            re[1]=ceList[i][1];
3076            re[2]=ceList[i][2];
3077            re[3]=ceList[i][3];
3078            re[4]=ceList[i][4];
3079          }
3080          else
3081          {
3082//--- otherwise ignore
3083            i++;
3084            continue;
3085          }
3086        }
3087        else
3088        {
3089          re[1]=intersect(re[1],ceList[i][1]);
3090        }
3091     }
3092   }
3093//--------------------------------------------------------------------------
3094// Perform last checks and return the result
3095//--------------------------------------------------------------------------
3096   if(size(re)!=4)
3097   {
3098//--- oops: already resolved in all open sets
3099     re[1]=BO[2];
3100     re[2]=-1;
3101     re[3]=0;
3102     re[4]=intvec(0);
3103   }
3104   return(re);
3105}
3106//////////////////////////////////////////////////////////////////////////////
3107static proc SpecialCoeff(list BO,int b,poly f)
3108{
3109//----------------------------------------------------------------------------
3110// Coeff with given hypersurface -- no checks of the hypersurface performed
3111//----------------------------------------------------------------------------
3112   int i,ch;
3113   int e=int(factorial(b));
3114   ideal C;
3115   list L=DeltaList(BO);
3116   int d=size(L);
3117//--- set up ideal
3118   for(i=0;i<b;i++)
3119   {
3120      C=C+L[i+1]^(e/(b-i));
3121   }
3122//--- move to hypersurface V(Z)
3123   ideal Z=f;
3124   C=C,Z;
3125   BO[1]=BO[1]+Z;
3126   BO[2]=C;
3127   for(i=1;i<=size(BO[4]);i++)
3128   {
3129      BO[6][i]=0;             // reset intersection indicator
3130      BO[4][i]=BO[4][i]+Z;    // intersect the E_i
3131      BO[2]=sat(BO[2],BO[4][i]+BO[1])[1];
3132                              // "strict transform" of J w.r.t E, not "total"
3133   }
3134   return(BO);
3135}
3136
3137//////////////////////////////////////////////////////////////////////////////
3138static proc DropCoeff(list BO)
3139"Internal procedure - no help and no example available
3140"
3141{
3142//--- Initialization
3143  int i,j;
3144  list pr=minAssGTZ(BO[2]);
3145  ideal I=BO[2];
3146  ideal Itemp;
3147  ideal Idropped=1;
3148//--- Tests
3149  for(i=1;i<=size(pr);i++)
3150  {
3151     if(i>size(pr))
3152     {
3153//--- the continue statement does not test the loop condition *sigh*
3154        break;
3155     }
3156     if(deg(std(slocus(pr[i]))[1])!=0)
3157     {
3158//--- this component is singular ===> we still need it
3159        i++;
3160        continue;
3161     }
3162     Itemp=sat(I,pr[i])[1];
3163     if(deg(std(Itemp+pr[i])[1])!=0)
3164     {
3165//--- this component is not disjoint from the other ones ===> we still need it
3166        i++;
3167        continue;
3168     }
3169     if(!transversalT(pr[i],BO[4]))
3170     {
3171//--- this component does not meet one of the remaining E_i transversally
3172//--- ===> we still need it
3173        i++;
3174        continue;
3175     }
3176     if(!normalCross(BO[4],pr[i]))
3177     {
3178//--- this component is not normal crossing with the remaining E_i
3179//--- ===> we still need it
3180        i++;
3181        continue;
3182     }
3183     if(defined(EE)){kill EE;}
3184     list EE;
3185     for(j=1;j<=size(BO[4]);j++)
3186     {
3187        EE[j]=BO[4][j]+pr[i];
3188     }
3189     if(!normalCross(EE))
3190     {
3191//--- we do not have a normal crossing situation for this component after all
3192//--- ===> we still need it
3193        i++;
3194        continue;
3195     }
3196     Idropped=intersect(Idropped,pr[i]);
3197     I=Itemp;
3198  }
3199  return(Idropped);
3200}
3201
3202//////////////////////////////////////////////////////////////////////////////
3203static proc DropRedundant(list BO,list E)
3204"Internal procedure - no help and no example available
3205"
3206{
3207//---------------------------------------------------------------------------
3208// Initialization and sanity checks
3209//---------------------------------------------------------------------------
3210  int ii,jj,kkdiff,nonnormal,ok;
3211  ideal testid,dummy;
3212  ideal center;
3213  intvec transverse,dontdrop,zerovec;
3214  transverse[size(BO[4])]=0;
3215  dontdrop[size(E[4])]=0;
3216  zerovec[size(E[4])]=0;
3217  ideal J=BO[2];
3218  int dimJ=dim(std(BO[2]));
3219  list templist;
3220  if(size(E)<5)
3221  {
3222//--- should not occur
3223    return(E);
3224  }
3225  for(ii=1;ii<=BO[7][1];ii++)
3226  {
3227     if(BO[6][ii]==2)
3228     {
3229       kkdiff++;
3230     }
3231  }
3232  int expDim=dimJ-E[2]+kkdiff;
3233  if(size(E)==6)
3234  {
3235    nonnormal=E[6];
3236  }
3237//---------------------------------------------------------------------------
3238// if dimJ were smaller than E[2], we would not have more than 3 entries in
3239//    in the list E
3240// if dimJ is also at least E[3] and nonnormal is 0, we only need to test that
3241// * the intersection is of the expected dimension
3242// * the intersections of the BO[4][i] and J are normal crossing
3243// * the elements of E^+ have no influence (is done below)
3244//---------------------------------------------------------------------------
3245  if((E[3]<=dimJ)&&(!nonnormal))
3246  {
3247    ideal bla=E[1]+BO[2]+BO[1];
3248    bla=radical(bla);
3249    bla=mstd(bla)[2];
3250
3251    if(dim(std(slocusE(bla)))<0)
3252    {
3253      if(transversalT(J,BO[4]))
3254      {
3255        ok=1;
3256        if(E[2]==E[3])
3257        {
3258//--- no further intersection with elements from E^+
3259          for(ii=1;ii<=size(E[4]);ii++)
3260          {
3261            if(dim_slocus(BO[2]+E[4][ii])!=-1)
3262            {
3263              dontdrop[ii]=1;
3264            }
3265          }
3266          if(dontdrop==zerovec)
3267          {
3268            list relist;
3269            relist[1]=std(J);
3270            return(relist);
3271          }
3272        }
3273      }
3274    }
3275  }
3276//---------------------------------------------------------------------------
3277// now check whether the E_i actually occurring in the intersections meet
3278// J transversally (one by one) and mark those elements of E[4] where it is
3279// not the case
3280//---------------------------------------------------------------------------
3281 if(!ok)
3282 {
3283     for(ii=1;ii<=size(E[5]);ii++)
3284     {
3285//--- test the ii-th tuple of E[4] resp. its indices E[5]
3286       for(jj=1;jj<=size(E[5][ii]);jj++)
3287       {
3288//--- if E[5][ii][jj]==1, E_jj is involved in E[4][ii]
3289        if(E[5][ii][jj]==1)
3290         {
3291//--- transversality not yet determined
3292           if(transverse[jj]==0)
3293           {
3294             templist[1]=BO[4][jj];
3295             if(transversalT(BO[2],templist))
3296             {
3297               transverse[jj]=1;
3298             }
3299             else
3300             {
3301               transverse[jj]=-1;
3302               dontdrop[ii]=1;
3303             }
3304           }
3305           else
3306           {
3307//--- already computed transversality
3308             if(transverse[jj]<0)
3309             {
3310               dontdrop[ii]=1;
3311             }
3312           }
3313         }
3314       }
3315     }
3316   }
3317//---------------------------------------------------------------------------
3318// if one of the non-marked tuples from E^- in E[4] has an intersection
3319// of the expected dimension and does not meet any E_i from E^+
3320// - except the ones which are met trivially - , it should be
3321// dropped from the list.
3322// it can also be dropped if an intersection occurs and normal crossing has
3323// been checked.
3324//---------------------------------------------------------------------------
3325  for(ii=1;ii<=size(E[4]);ii++)
3326  {
3327//--- if E[4][ii] does not have transversal intersections, we cannot drop it
3328    if(dontdrop[ii]==1)
3329    {
3330      ii++;
3331      continue;
3332    }
3333//--- testing ii-th tuple from E[4]
3334    testid=BO[1]+BO[2]+E[4][ii];
3335    if(dim(std(testid))!=expDim)
3336    {
3337//--- not expected dimension
3338      dontdrop[ii]=1;
3339      ii++;
3340      continue;
3341    }
3342    testid=mstd(testid)[2];
3343
3344    if(dim(std(slocusE(testid)))>=0)
3345    {
3346//--- not smooth, i.e. more than one component which intersect
3347      dontdrop[ii]=1;
3348      ii++;
3349      continue;
3350    }
3351//--- if E^+ is empty, we are done; otherwise check intersections with E^+
3352    if(BO[7][1]!=-1)
3353    {
3354      if(defined(pluslist)){kill pluslist;}
3355      list pluslist;
3356      for(jj=BO[7][1]+1;jj<=size(BO[4]);jj++)
3357      {
3358        dummy=BO[4][jj]+testid;
3359        dummy=std(dummy);
3360        if(expDim==dim(dummy))
3361        {
3362//--- intersection has wrong dimension
3363          dontdrop[ii]=1;
3364          break;
3365        }
3366        pluslist[jj-BO[7][1]]=BO[4][jj]+testid;
3367      }
3368      if(dontdrop[ii]==1)
3369      {
3370        ii++;
3371        continue;
3372      }
3373      if(!normalCross(pluslist))
3374      {
3375//--- unfortunately, it is not normal crossing
3376        dontdrop[ii]=1;
3377      }
3378    }
3379  }
3380//---------------------------------------------------------------------------
3381// The returned list should look like the truncated output of inters_E
3382//---------------------------------------------------------------------------
3383  list retlist;
3384  for(ii=1;ii<=size(E[4]);ii++)
3385  {
3386    if(dontdrop[ii]==1)
3387    {
3388      if(size(center)>0)
3389      {
3390        center=intersect(center,E[4][ii]);
3391      }
3392      else
3393      {
3394        center=E[4][ii];
3395      }
3396    }
3397  }
3398  retlist[1]=center;
3399  retlist[2]=E[2];
3400  retlist[3]=E[3];
3401  return(retlist);
3402}
3403//////////////////////////////////////////////////////////////////////////////
3404static proc transversalT(ideal J, list E,list #)
3405"Internal procedure - no help and no example available
3406"
3407{
3408//----------------------------------------------------------------------------
3409// check whether J and each element of the list E meet transversally
3410//----------------------------------------------------------------------------
3411   def R=basering;
3412   if(size(#)>0)
3413   {
3414     ideal pp=#[1];
3415   }
3416   int i;
3417   ideal T,M;
3418   ideal Jstd=std(J);
3419   ideal Tstd;
3420   int d=nvars(basering)-dim(Jstd)+1;   // d=n-dim(V(J) \cap hypersurface)
3421   for(i=1;i<=size(E);i++)
3422   {
3423      if(size(reduce(E[i],Jstd))==0)
3424      {
3425//--- V(J) is contained in E[i]
3426        return(0);
3427      }
3428      T=J,E[i];
3429      Tstd=std(T);
3430      d=nvars(basering)-dim(Tstd);
3431      if(deg(Tstd[1])!=0)
3432      {
3433//--- intersection is non-empty
3434//!!! abgeklemmt, da es doch in der Praxis vorkommt und korrekt sein kann!!!
3435//!!! wenn ueberhaupt dann -1 zurueckgeben!!!
3436//         if((d>=4)&&(size(T)>=10)){return(0);}
3437         M=minor(jacob(T),d,Tstd)+T;
3438         M=std(M);
3439         if(deg(M[1])>0)
3440         {
3441//--- intersection is not transversal
3442           if(size(#)==0)
3443           {
3444              return(0);
3445           }
3446           M=std(radical(M));
3447           if(size(reduce(pp,M))>0){return(0);}
3448         }
3449      }
3450   }
3451//--- passed all tests
3452   return(1);
3453}
3454///////////////////////////////////////////////////////////////////////////////
3455static proc transversalTB(ideal J, list E,ideal V)
3456"Internal procedure - no help and no example available
3457"
3458{
3459//----------------------------------------------------------------------------
3460// check whether J and each element of the list E meet transversally
3461//----------------------------------------------------------------------------
3462   def R=basering;
3463
3464   int i;
3465   ideal T,M;
3466   ideal Jstd=std(J);
3467   ideal Tstd;
3468   int d=nvars(basering)-dim(Jstd)+1;   // d=n-dim(V(J) \cap hypersurface)
3469   for(i=1;i<=size(E);i++)
3470   {
3471      if(size(reduce(E[i],Jstd))==0)
3472      {
3473//--- V(J) is contained in E[i]
3474        return(0);
3475      }
3476      T=J,E[i];
3477      Tstd=std(T);
3478      d=nvars(basering)-dim(Tstd);
3479      if(deg(Tstd[1])!=0)
3480      {
3481//--- intersection is non-empty
3482         if((d>=4)&&(size(T)>=10)){return(0);}
3483         M=minor(jacob(T),d,Tstd)+T;
3484         M=std(M+V);
3485         if(deg(M[1])>0)
3486         {
3487            return(0);
3488         }
3489      }
3490   }
3491//--- passed all tests
3492   return(1);
3493}
3494///////////////////////////////////////////////////////////////////////////////
3495static proc powerI(ideal I,int n,int m)
3496{
3497//--- compute (n!/m)-th power of I, more efficient variant
3498   int i;
3499   int mon=1;
3500   for(i=1;i<=size(I);i++)
3501   {
3502     if(size(I[i])>1){mon=0;break;}
3503   }
3504   if(mon)
3505   {
3506      if(size(reduce(I,std(radical(I[1]))))<size(I)-1){mon=0;}
3507   }
3508   if((mon)&&(size(I)>3))
3509   {
3510      int e=int(factorial(n))/m;
3511      ideal J=1;
3512      poly p=I[1];
3513      I=I[2..size(I)];
3514      ideal K=p^e;
3515      for(i=1;i<=e;i++)
3516      {
3517         J=interred(J*I);
3518         K=K,(p^(e-i))*J;
3519      }
3520      return(K);
3521   }
3522   for(i=n;i>1;i--)
3523   {
3524      if(i!=m)
3525      {
3526         I=I^i;
3527      }
3528   }
3529   return(I);
3530}
3531
3532///////////////////////////////////////////////////////////////////////////////
3533static proc Coeff(list BO, int b, list #)
3534"USAGE:  Coeff (BO);
3535@*       BO = basic object, a list: ideal W,
3536@*                                  ideal J,
3537@*                                  intvec b (already truncated for Coeff),
3538@*                                  list Ex  (already truncated for Coeff),
3539@*                                  ideal ab,
3540@*                                  intvec v,
3541@*                                  intvec w (already truncated for Coeff),
3542@*                                  matrix M
3543@*       b  = integer indication bmax(BO)
3544ASSUME:  R  = basering, a polynomial ring, W an ideal of R,
3545@*       J  = ideal containing W
3546COMPUTE: Coeff-Ideal of BO as defined in [Bravo,Encinas,Villamayor]
3547RETURN:  basic object of the Coeff-Ideal
3548EXAMPLE: example Coeff;    shows an example
3549"
3550{
3551//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3552//!!! TASK: lower dimension by more than one in a single step if possible !!!
3553//!!!       (improve bookkeeping of invariants in Coeff and Center)       !!!
3554//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3555//---------------------------------------------------------------------------
3556// Initialization  and sanity checks
3557//---------------------------------------------------------------------------
3558   int i,k,dummy,errtype;
3559   int ma=size(BO[4]);
3560   intvec merk;
3561   if(!defined(debugCoeff))
3562   {
3563     int debugCoeff;
3564   }
3565   ideal C;
3566   list L;
3567   if(size(#)!=0)
3568   {
3569      if(typeof(#[1])=="ideal")
3570      {
3571        L=#;
3572      }
3573      else
3574      {
3575        ma=#[1];
3576        L=DeltaList(BO);
3577      }
3578   }
3579   else
3580   {
3581      L=DeltaList(BO);
3582   }
3583
3584   if(debugCoeff)
3585   {
3586     "----> In Coeff: result of DeltaList:";
3587     L;
3588   }
3589   int d=size(L);                   // bmax of BO
3590   if((debugCoeff)&&(d!=b))
3591   {
3592      "!!!!!! Length of DeltaList does not equal second argument !!!!!!";
3593      "!!!!!! BO might not have been ord ~ 1 or wrong b          !!!!!!";
3594   }
3595   if(b>=6){return(0);}              // b is too big
3596   int e=int(factorial(b));          // b of Coeff-Ideal
3597   if(e==0)
3598   {
3599     ERROR( "// integer size too small for forming b! .");
3600   }
3601   if(b==0)
3602   {
3603     ERROR( "// second argument to Coeff should never be zero."  );
3604   }
3605//----------------------------------------------------------------------------
3606// Form the Coeff-Ideal
3607// Step 1: choose hypersurface
3608// Step 2: sum over correct powers of Delta^i(BO[2])
3609// Step 3: do the intersection
3610//----------------------------------------------------------------------------
3611//--- Step 1
3612   ideal Z;
3613   poly p;
3614   for(i=1;i<=ncols(L[d]);i++)
3615   {
3616//--- Look for smooth hypersurface in generators of Delta^(bmax-1)(BO[2])
3617     dummy=goodChoice(BO,L[d][i]);
3618     if(!dummy)
3619     {
3620       Z= L[d][i];
3621       break;
3622     }
3623     else
3624     {
3625       if(dummy>1)
3626       {
3627          merk[size(merk)+1]=i;
3628       }
3629       if(dummy>errtype)
3630       {
3631         errtype=dummy;
3632       }
3633     }
3634   }
3635   if(size(Z)==0)
3636   {
3637//--- no suitable element in generators of Delta^(bmax-1)(BO[2])
3638//--- try random linear combination
3639      for(k=1;k<=10;k++)
3640      {
3641         for(i=2;i<=size(merk);i++)
3642         {
3643            p=p+random(-100,100)*L[d][merk[i]];
3644         }
3645         dummy=goodChoice(BO,p);
3646         if(!dummy)
3647         {
3648            Z=p;
3649            break;
3650         }
3651         else
3652         {
3653            p=0;
3654         }
3655      }
3656      if(dummy)
3657      {
3658        for(i=1;i<=size(L[d]);i++)
3659        {
3660           p=p+random(-100,100)*L[d][i];
3661        }
3662        dummy=goodChoice(BO,p);
3663        if(!dummy)
3664        {
3665//--- found a suitable one
3666           Z=p;
3667        }
3668      }
3669      if(dummy)
3670      {
3671//--- did not find a suitable random linear combination either
3672         if(dummy>errtype)
3673         {
3674           errtype=dummy;
3675         }
3676         list retlist=errtype,L[d];
3677         return(retlist);
3678      }
3679   }
3680   if(debugCoeff)
3681   {
3682     "----> In Coeff: Chosen hypersurface";
3683     Z;
3684   }
3685//--- Step 2
3686   C=Z;
3687   for(i=0;i<b;i++)
3688   {
3689      C=C,powerI(simplify(reduce(L[i+1],std(Z)),2),b,b-i);
3690   }
3691   C=interred(C);
3692
3693   if(debugCoeff)
3694   {
3695     "----> In Coeff: J before saturation";
3696     C;
3697   }
3698
3699//--- Step 3
3700   BO[1]=BO[1]+Z;
3701   BO[2]=C;
3702   for(i=1;i<=size(BO[4]);i++)
3703   {
3704      BO[6][i]=0;             // reset intersection indicator
3705      BO[4][i]=BO[4][i]+Z;    // intersect the E_i
3706      if(i<=ma)
3707      {
3708        BO[2]=sat(BO[2],BO[4][i]+BO[1])[1];
3709                       // "strict transform" of J w.r.t E, not "total"
3710      }
3711   }
3712   if(debugCoeff)
3713   {
3714     "----> In Coeff:";
3715     " J after saturation:";
3716     BO[2];
3717   }
3718   return(BO);
3719}
3720example
3721{"EXAMPLE:";
3722   echo = 2;
3723   ring R=0,(x,y,z),dp;
3724
3725   ideal W;
3726   ideal J=z^2+x^2*y^2;
3727   intvec b=0;
3728   list E;
3729   ideal abb=maxideal(1);
3730   intvec v;
3731   intvec w=-1;
3732   matrix M;
3733
3734   list BO=W,J,b,E,abb,v,w,M;
3735
3736   Coeff(BO,2);
3737}
3738//////////////////////////////////////////////////////////////////////////////
3739static proc goodChoice(list BO, poly p)
3740"Internal procedure - no help and no example available
3741"
3742{
3743//---------------------------------------------------------------------------
3744// test whether new W is smooth
3745//---------------------------------------------------------------------------
3746   ideal W=BO[1]+ideal(p);
3747   if(size(reduce(p,std(BO[1])))==0)
3748   {
3749//--- p is already in BO[1], i.e. does not define a hypersurface in W
3750     return(1);
3751   }
3752   if(dim(std(slocusE(W)))>=0)
3753//   if(dim(timeStd(slocusE(W),20))>=0)
3754   {
3755//--- new W would not be smooth
3756     return(1);
3757   }
3758   if(size(BO[4])==0)
3759   {
3760//--- E is empty, no further tests necessary
3761     return(0);
3762   }
3763//--------------------------------------------------------------------------
3764// test whether the hypersurface meets the E_i transversally
3765//--------------------------------------------------------------------------
3766   list E=BO[4];
3767   int i,d;
3768   ideal T=W;
3769   ideal Tstd=std(T);
3770   d=nvars(basering)-dim(Tstd)+1;
3771   ideal M;
3772   for(i=1;i<=size(E);i++)
3773   {
3774      T=W,E[i];
3775      M=minor(jacob(T),d,Tstd)+T;
3776      M=std(M);
3777      if(deg(M[1])>0)
3778      {
3779//--- intersection not transversal
3780        return(2);
3781      }
3782   }
3783//--------------------------------------------------------------------------
3784// test whether the new E_i have normal crossings
3785//--------------------------------------------------------------------------
3786   for(i=1;i<=size(E);i++)
3787   {
3788      E[i]=E[i],p;
3789   }
3790   if(normalCross(E))
3791   {
3792     return(0);
3793   }
3794   else
3795   {
3796     return(2);
3797   }
3798}
3799//////////////////////////////////////////////////////////////////////////////
3800
3801proc presentTree(list L)
3802"USAGE:  presentTree(L);
3803         L=list, output of resolve
3804RETURN:  nothing, only pretty printing of the output data of resolve()
3805EXAMPLE: none
3806"
3807{
3808  def r=basering;
3809  int i,j,k;
3810  if(size(L[2])==1)
3811  {
3812     "The object was already resolved or the list L does not";
3813     "have required input format. There is just one chart in";
3814     "the tree.";
3815     return();
3816  }
3817  for(i=1;i<=size(L[1]);i++)
3818  {
3819    "                        ";
3820    "/////////////////////////// Final Chart",i,"/////////////////////////";
3821    def s=L[1][i];
3822    setring s;
3823    "======================== History of this chart ======================";
3824    for(j=2;j<=ncols(path);j++)
3825    {
3826      "                  ";
3827      "Blow Up",j-1,":";
3828      "     Center determined in L[2]["+string(path[1,j])+"],";
3829      "     Passing to chart ",path[2,j]," in resulting blow up.";
3830    }
3831    "                 ";
3832    "======================== Data of this chart ========================";
3833    showBO(BO);
3834    setring r;
3835    kill s;
3836    pause();
3837  }
3838  "////////////////////////////////////////////////////////////////////";
3839  "For identification of exceptional divisors please use the tools";
3840  "provided by reszeta.lib, e.g. collectDiv.";
3841  "For viewing an illustration of the tree of charts please use the";
3842  "procedure ResTree from resgraph.lib.";
3843  "////////////////////////////////////////////////////////////////////";
3844  return();
3845}
3846//////////////////////////////////////////////////////////////////////////////
3847
3848proc showBO(list BO)
3849"USAGE:  showBO(BO);
3850@*        BO=basic object, a list: ideal W,
3851@*                                  ideal J,
3852@*                                  intvec b (already truncated for Coeff),
3853@*                                  list Ex  (already truncated for Coeff),
3854@*                                  ideal ab,
3855@*                                  intvec v,
3856@*                                  intvec w (already truncated for Coeff),
3857@*                                  matrix M
3858RETURN:   nothing, only pretty printing
3859EXAMPLE:  none
3860"
3861{
3862  "                       ";
3863  "==== Ambient Space: ";BO[1];"      ";
3864  "==== Ideal of Variety: ";BO[2];"      ";
3865  int i;
3866  list M;
3867  for(i=1;i<=size(BO[4]);i++)
3868  {
3869    M[i]=ideal(BO[4][i]);
3870  }
3871  "==== Exceptional Divisors: ";print(M);"   ";
3872  "==== Images of variables of original ring:";BO[5];"   ";
3873}
3874//////////////////////////////////////////////////////////////////////////////
3875static proc max(int i,int j)
3876"Internal procedure - no help and no example available
3877"
3878{
3879//--- why is there no proc for max in general.lib?
3880  if(i>j){return(i);}
3881  return(j);
3882}
3883//////////////////////////////////////////////////////////////////////////////
3884static proc min(int i,int j)
3885"Internal procedure - no help and no example available
3886"
3887{
3888//--- why is there no proc for max in general.lib?
3889  if(i<j){return(i);}
3890  return(j);
3891}
3892//////////////////////////////////////////////////////////////////////////////
3893////////////////////////    main procedure    ////////////////////////////////
3894//////////////////////////////////////////////////////////////////////////////
3895proc resolve(ideal J, list #)
3896"USAGE:  resolve (J); or resolve (J,i[,k]);
3897@*       J ideal
3898@*       i,k int
3899COMPUTE: a resolution of J,
3900@*       if i > 0 debugging is turned on according to the following switches:
3901@*       j1: value 0 or 1; turn off or on correctness checks in all steps
3902@*       j2: value 0 or 2; turn off or on debugCenter
3903@*       j3: value 0 or 4; turn off or on debugBlowUp
3904@*       j4: value 0 or 8; turn off or on debugCoeff
3905@*       j5: value 0 or 16:turn off or on debugging of Intersection with E^-
3906@*       j6: value 0 or 32:turn off or on stop after pass throught the loop
3907@*       i=j1+j2+j3+j4+j5+j6
3908RETURN:  a list l of 2 lists of rings
3909         l[1][i] is a ring containing a basic object BO, the result of the
3910         resolution.
3911         l[2] contains all rings which occured during the resolution process
3912NOTE:    result may be viewed in a human readable form using presentTree()
3913EXAMPLE: example resolve;  shows an example
3914"
3915{
3916//----------------------------------------------------------------------------
3917// Initialization and sanity checks
3918//----------------------------------------------------------------------------
3919   def R=basering;
3920   list allRings;
3921   allRings[1]=R;
3922   list endRings;
3923   module path=[0,-1];
3924   ideal W;
3925   list E;
3926   ideal abb=maxideal(1);
3927   intvec v;
3928   intvec bvec;
3929   intvec w=-1;
3930   matrix intE;
3931   int extra,bm;
3932   if(defined(BO)){kill BO;}
3933   if(defined(cent)){kill cent;}
3934
3935   ideal Jrad=equiRadical(J);
3936   if(size(reduce(Jrad,std(J)))!=0)
3937   {
3938      "WARNING! The input is not reduced or not equidimensional!";
3939      "We will continue with the reduced top-dimensional part of input";
3940      J=Jrad;
3941   }
3942
3943   int i,j,debu,loca,locaT,ftemp,debugResolve,smooth;
3944//--- switches for local and for debugging may occur in any order
3945   i=size(#);
3946   extra=3;
3947   for(j=1;j<=i;j++)
3948   {
3949     if(typeof(#[j])=="int")
3950     {
3951       debugResolve=#[j];
3952//--- debu: debug switch for resolve, smallest bit in debugResolve
3953       debu=debugResolve mod 2;
3954     }
3955     else
3956     {
3957        if(#[j]=="M")
3958        {
3959           bm=1;
3960           ERROR("Not implemented yet");
3961        }
3962        if(#[j]=="E"){extra=0;}
3963        if(#[j]=="A"){extra=2;}
3964        if(#[j]=="K"){extra=3;}
3965        if(#[j]=="L"){loca=1;}
3966     }
3967   }
3968   if(loca)
3969   {
3970      list qs=minAssGTZ(J);
3971      ideal K=ideal(1);
3972      for(j=1;j<=size(qs);j++)
3973      {
3974         if(size(reduce(qs[j],std(maxideal(1))))==0)
3975         {
3976            K=intersect(K,qs[j]);
3977         }
3978      }
3979      J=K;
3980      list qr=minAssGTZ(slocus(J));
3981      K=ideal(1);
3982      for(j=1;j<=size(qr);j++)
3983      {
3984         if(size(reduce(qr[j],std(maxideal(1))))!=0)
3985         {
3986            K=intersect(K,qr[j]);
3987            smooth++;
3988         }
3989         else
3990         {
3991            if(dim(std(qr[j]))>0){loca=0;}
3992//---- test for isolated singularity at 0
3993         }
3994      }
3995      K=std(K);
3996//---- if deg(K[1])==0 the point 0 is on all components of the singular
3997//---- locus and we can work globally
3998      if(smooth==size(qr)){smooth=-1;}
3999//---- the point 0 is not on the singular locus
4000      if((deg(K[1])>0)&&(smooth>=0)&&(!loca))
4001      {
4002         locaT=1;
4003         poly @p;
4004         for(j=1;j<=size(K);j++)
4005         {
4006            if(jet(K[j],0)!=0)
4007            {
4008               @p=K[j];
4009               break;
4010            }
4011         }
4012         export(@p);
4013      }
4014      if((loca)&&(!smooth)){loca=0;}
4015//---- the case that 0 is isolated singularity and the only singular point
4016   }
4017   export(locaT);
4018//---In case of option "L" the following holds
4019//---loca=0 and locaT=0  we perform the global case
4020//---loca !=0: 0 is isolated singular point, but there are other singularities
4021//---locaT!=0: O is singular point, but not isolated, and there is a componente//---          of the singular locus not containing 0
4022
4023//--- if necessary, set the corresponding debugFlags
4024   if(defined(debugResolve))
4025   {
4026//--- 2nd bit from the right
4027     int debugCenter=(debugResolve div 2) mod 2;
4028     export debugCenter;
4029//--- 3rd bit from the right
4030     int debugBlowUp=(debugResolve div 4) mod 2;
4031     export debugBlowUp;
4032//--- 4th bit from the right
4033     int debugCoeff=(debugResolve div 8) mod 2;
4034     export debugCoeff;
4035//--- 5th bit from the right
4036     int debug_Inters_E=(debugResolve div 16) mod 2;
4037     export debug_Inters_E;
4038//--- 6th bit from the right
4039     int praes_stop=(debugResolve div 32) mod 2;
4040   }
4041//--- set the correct attributes to J for speed ups
4042   if( typeof(attrib(J,"isEqui"))!="int" )
4043   {
4044      if(size(J)==1)
4045      {
4046         attrib(J,"isEqui",1);
4047      }
4048      else
4049      {
4050         attrib(J,"isEqui",0);
4051      }
4052   }
4053   if(size(J)==1)
4054   {
4055      attrib(J,"isHy",1);
4056   }
4057   else
4058   {
4059      attrib(J,"isHy",0);
4060   }
4061//--- create the BO
4062   list BO=W,J,bvec,E,abb,v,w,intE;
4063   if(defined(invSat)){kill invSat;}
4064   list invSat=ideal(0),intvec(0);
4065   export(invSat);
4066   if(bm)
4067   {
4068     intmat invmat[2][1]=0,-1;
4069     BO[9]=invmat;
4070   }
4071   else
4072   {
4073     BO[9]=intvec(0);
4074   }
4075   export BO;
4076   list tmpList;
4077   int blo;
4078   int k,Ecount,tmpPtr;
4079   i=0;
4080   if(smooth==-1)
4081   {
4082      endRings[1]=R;
4083      list result=endRings,allRings;
4084      if(debu)
4085      {
4086         "============= result will be tested ==========";
4087         "                                              ";
4088         "the number of charts obtained:",size(endRings);
4089         "=============     result is o.k.    ==========";
4090      }
4091      kill debugCenter,debugBlowUp,debugCoeff,debug_Inters_E;
4092      return(result);
4093   }
4094//-----------------------------------------------------------------------------
4095// While there are rings to be considered, determine center and blow up
4096//-----------------------------------------------------------------------------
4097   while(i<size(allRings))
4098   {
4099      i++;
4100      def S=allRings[i];
4101      setring S;
4102      list pr;
4103      ideal Jstd=std(BO[2]);
4104//-----------------------------------------------------------------------------
4105// Determine Center
4106//-----------------------------------------------------------------------------
4107      if(i==1)
4108      {
4109         list deltaL=DeltaList(BO);
4110         ideal sL=radical(deltaL[size(deltaL)]);
4111         if((deg(std(slocus(sL))[1])==0)&&(size(minAssGTZ(sL))==1))
4112         {
4113            list @ce=sL,intvec(-1),intvec(0),intvec(0);
4114            ideal cent=@ce[1];
4115         }
4116      }
4117//--- before computing a center, check whether we have a stored one
4118      if(size(BO)>9)
4119      {
4120         while(size(BO[10])>0)
4121         {
4122            list @ce=BO[10][1];
4123//--- check of the center
4124           // @ce=correctC(BO,@ce,bm);
4125//--- use stored center
4126            BO[10]=delete(BO[10],1);
4127            if(size(@ce[1])==0)
4128            {
4129//--- stored center was not ok
4130              continue;
4131            }
4132            tmpPtr=0;
4133            for(Ecount=1;Ecount <= size(@ce[2]); Ecount++)
4134            {
4135              if(@ce[2][Ecount]>-1)
4136              {
4137                tmpPtr=tmpPtr+@ce[2][Ecount];
4138              }
4139              else
4140              {
4141                @ce[2][Ecount]=size(BO[4])-tmpPtr-1;
4142                for(int cnthlp=1;cnthlp<=size(BO[10]);cnthlp++)
4143                {
4144                   BO[10][cnthlp][2][Ecount]=@ce[2][Ecount];
4145                }
4146                kill cnthlp;
4147                break;
4148              }
4149            }
4150            if(Ecount<size(@ce[2]))
4151            {
4152              for(tmpPtr=Ecount+1;tmpPtr<=size(@ce[2]);tmpPtr++)
4153              {
4154                @ce[2][tmpPtr]=0;
4155                for(int cnthlp=1;cnthlp<=size(BO[10]);cnthlp++)
4156                {
4157                   BO[10][cnthlp][2][tmpPtr]=@ce[2][tmpPtr];
4158                }
4159                kill cnthlp;
4160              }
4161            }
4162            break;
4163         }
4164         if(defined(@ce))
4165         {
4166           if(size(@ce[1])==0)
4167           {
4168              kill @ce;
4169           }
4170           else
4171           {
4172              ideal cent=@ce[1];
4173           }
4174         }
4175         if(size(BO[10])==0)
4176         {
4177//--- previously had stored centers, all have been used; we need to clean up
4178            BO=delete(BO,10);
4179         }
4180      }
4181      if((loca)&&(i==1))
4182      {
4183//--- local case: initial step is blow-up in origin
4184         if(defined(@ce)){kill @ce;}
4185         if(defined(cent)){kill cent;}
4186         if(size(reduce(slocusE(BO[2]),std(maxideal(1))))==0)
4187         {
4188            list @ce=maxideal(1),intvec(-1),intvec(0),intvec(0);
4189         }
4190         else
4191         {
4192            list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4193         }
4194         ideal cent=@ce[1];
4195      }
4196      if(((loca)||(locaT))&&(i!=1))
4197      {
4198         int JmeetsE;
4199         for(j=1;j<=size(BO[4]);j++)
4200         {
4201            if(deg(std(BO[2]+BO[4][j])[1])!=0)
4202            {
4203               JmeetsE=1;
4204               break;
4205            }
4206         }
4207         if(!JmeetsE)
4208         {
4209            list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4210            ideal cent=@ce[1];
4211         }
4212         kill JmeetsE;
4213      }
4214      if((locaT)&&(!defined(@ce)))
4215      {
4216          if(@p!=1)
4217          {
4218             list tr=minAssGTZ(slocusE(BO[2]));
4219             ideal L=ideal(1);
4220             for(j=1;j<=size(tr);j++)
4221             {
4222                if(size(reduce(ideal(@p),std(tr[j])))==0)
4223                {
4224                   L=intersect(L,tr[j]);
4225                }
4226             }
4227             L=std(L);
4228             if(deg(L[1])==0)
4229             {
4230                @p=1;
4231             }
4232             else
4233             {
4234               ideal fac=factorize(@p,1);
4235               if(size(fac)==1)
4236               {
4237                 @p=fac[1];
4238               }
4239               else
4240               {
4241                  for(j=1;j<=size(fac);j++)
4242                  {
4243                     if(reduce(fac[j],L)==0)
4244                     {
4245                        @p=fac[j];
4246                        break;
4247                     }
4248                  }
4249               }
4250             }
4251             kill tr,L;
4252          }
4253          execute("ring R1=("+charstr(S)+"),(@z,"+varstr(S)+"),dp;");
4254          poly p=imap(S,@p);
4255          list BO=imap(S,BO);
4256          list invSat=imap(S,invSat);
4257          export(invSat);
4258          ideal te=std(BO[2]);
4259          BO[1]=BO[1]+ideal(@z*p-1);
4260          BO[2]=BO[2]+ideal(@z*p-1);
4261          for(j=1;j<=size(BO[4]);j++)
4262          {
4263             BO[4][j]=BO[4][j]+ideal(@z*p-1);
4264          }
4265//--- for computation of center: drop components not meeting the Ei
4266          def BO2=BO;
4267          list qs=minAssGTZ(BO2[2]);
4268          ideal K=ideal(1);
4269          for(j=1;j<=size(qs);j++)
4270          {
4271             if(CompMeetsE(qs[j],BO2[4]))
4272             {
4273                K=intersect(K,qs[j]);
4274             }
4275          }
4276          BO2[2]=K;
4277//--- check whether we are done
4278          if(deg(std(BO2[2])[1])==0)
4279          {
4280             list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4281          }
4282          if(!defined(@ce))
4283          {
4284             if(bm)
4285             {
4286                list @ce=CenterBM(BO2);
4287             }
4288             else
4289             {
4290                list @ce=CenterBO(BO2);
4291             }
4292          }
4293//--- if computation of center returned BO2[2], we are done
4294//--- ==> set @ce to BO[2], because later checks work with BO instead of BO2
4295          if((size(reduce(@ce[1],std(BO2[2])))==0)&&
4296             (size(reduce(BO2[2],std(@ce[1])))==0))
4297          {
4298             @ce[1]=BO[2];
4299          }
4300          if(size(specialReduce(@ce[1],te,p))==0)
4301          {
4302             BO=imap(S,BO);
4303             @ce[1]=BO[2];
4304          }
4305          else
4306          {
4307             //@ce=correctC(BO,@ce,bm);
4308             @ce[1]=eliminate(@ce[1],@z);
4309          }
4310          setring S;
4311          list @ce=imap(R1,@ce);
4312          kill R1;
4313
4314         if((size(reduce(BO[2],std(@ce[1])))==0)
4315             &&(size(reduce(@ce[1],Jstd))==0))
4316         {
4317//--- J and center coincide
4318            pr[1]=@ce[1];
4319            ideal cent=@ce[1];
4320         }
4321         else
4322         {
4323//--- decompose center and use first component
4324            pr=minAssGTZ(@ce[1]);
4325               if(size(reduce(@p,std(pr[1])))==0){"Achtung";~;}
4326               if(deg(std(slocus(pr[1]))[1])>0){"singulaer";~;}
4327            ideal cent=pr[1];
4328         }
4329         if(size(pr)>1)
4330         {
4331//--- store the other components
4332            for(k=2;k<=size(pr);k++)
4333            {
4334               if(size(reduce(@p,std(pr[k])))==0){"Achtung";~;}
4335               if(deg(std(slocus(pr[k]))[1])>0){"singulaer";~;}
4336               if(size(reduce(@p,std(pr[k])))!=0)
4337               {
4338                  tmpList[size(tmpList)+1]=list(pr[k],@ce[2],@ce[3],@ce[4]);
4339               }
4340            }
4341            BO[10]=tmpList;
4342            kill tmpList;
4343            list tmpList;
4344         }
4345      }
4346      if(!defined(@ce))
4347      {
4348//--- no previously determined center, we need to compute one
4349         if(loca)
4350         {
4351//--- local case: center should be inside exceptional locus
4352            ideal Ex=ideal(1);
4353            k=0;
4354            for(j=1;j<=size(BO[4]);j++)
4355            {
4356               if(deg(BO[4][j][1])!=0)
4357               {
4358                 Ex=Ex*BO[4][j];   //----!!!!hier evtl. Durchschnitt???
4359                 k++;
4360               }
4361            }
4362//--- for computation of center: drop components not meeting the Ei
4363            list BOloc=BO;
4364            list qs=minAssGTZ(BOloc[2]);
4365            ideal K=ideal(1);
4366            for(j=1;j<=size(qs);j++)
4367            {
4368              if(CompMeetsE(qs[j],BOloc[4]))
4369              {
4370                 K=intersect(K,qs[j]);
4371              }
4372            }
4373            BOloc[2]=K;
4374//--- check whether we are done
4375            if(deg(std(BOloc[2])[1])==0)
4376            {
4377               list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4378            }
4379            if(!defined(@ce))
4380            {
4381              if(BO[3][1]!=0)
4382              {
4383                 BOloc[2]=BO[2]+Ex^((BO[3][1] div k)+1);//!!!!Vereinfachen???
4384              }
4385              else
4386              {
4387                 BOloc[2]=BO[2]+Ex^((size(DeltaList(BO)) div k)+1);
4388              }
4389              if(bm)
4390              {
4391                 list @ce=CenterBM(BOloc);
4392              }
4393              else
4394              {
4395                 list @ce=CenterBO(BOloc);
4396              }
4397              if(size(reduce(Ex,std(@ce[1])))!=0)
4398              {
4399                 list tempPr=minAssGTZ(@ce[1]);
4400                 for(k=size(tempPr);k>=1;k--)
4401                 {
4402                    if(size(reduce(Ex,std(tempPr[k])))!=0)
4403                    {
4404                      tempPr=delete(tempPr,k);
4405                    }
4406                 }
4407                 @ce[1]=1;
4408                 for(k=1;k<=size(tempPr);k++)
4409                 {
4410                    @ce[1]=intersect(@ce[1],tempPr[k]);
4411                 }
4412                 if(deg(std(@ce[1])[1])==0)
4413                 {
4414                    @ce[1]=BO[2];
4415                 }
4416              }
4417            }
4418//--- test whether we are done
4419            if(size(reduce(slocusE(BO[2]),std(@ce[1])))!=0)
4420            {
4421               if(transversalT(BO[2],BO[4]))
4422               {
4423                  if(defined(E)){kill E;}
4424                  list E=BO[4];
4425                  for(j=1;j<=size(E);j++){if(deg(E[j][1])>0){E[j]=E[j]+BO[2];}}
4426                  if(normalCross(E))
4427                  {
4428                     @ce[1]=BO[2];
4429                  }
4430                  kill E;
4431               }
4432            }
4433         }
4434         else
4435         {
4436//--- non-local
4437            if(bm)
4438            {
4439               list @ce=CenterBM(BO);
4440            }
4441            else
4442            {
4443               list @ce=CenterBO(BO);
4444            }
4445//--- check of the center
4446            //@ce=correctC(BO,@ce,bm);
4447            if((size(@ce[1])==0)&&(size(@ce[4])<(size(@ce[3])-1)))
4448            {
4449               intvec xxx=@ce[3];
4450               xxx=xxx[1..size(@ce[4])];
4451               @ce[3]=xxx;
4452               xxx=@ce[2];
4453               xxx=xxx[1..size(@ce[4])];
4454               @ce[2]=xxx;
4455               kill xxx;
4456            }
4457         }
4458         if((size(reduce(BO[2],std(@ce[1])))==0)
4459             &&(size(reduce(@ce[1],Jstd))==0))
4460         {
4461//--- J and center coincide
4462            pr[1]=@ce[1];
4463            ideal cent=@ce[1];
4464         }
4465         else
4466         {
4467//--- decompose center and use first component
4468            pr=minAssGTZ(@ce[1]);
4469            ideal cent=pr[1];
4470         }
4471         if(size(pr)>1)
4472         {
4473//--- store the other components
4474            for(k=2;k<=size(pr);k++)
4475            {
4476               tmpList[k-1]=list(pr[k],@ce[2],@ce[3],@ce[4]);
4477            }
4478            BO[10]=tmpList;
4479            kill tmpList;
4480            list tmpList;
4481         }
4482      }
4483//--- do not forget to update BO[7] and BO[3]
4484      export cent;
4485      BO[7]=@ce[2];
4486      BO[3]=@ce[3];
4487      if((loca||locaT)&&(size(@ce)<4)){@ce[4]=0;} //Provisorium !!!
4488      if((size(@ce[4])<size(@ce[2])-1)||(size(@ce[4])<size(@ce[3])-1))
4489      {
4490        if((deg(std(@ce[1])[1])==0)&&(deg(std(BO[2])[1])==0))
4491        {
4492           intvec nullvec;
4493           nullvec[size(@ce[2])-1]=0;
4494           @ce[4]=nullvec;
4495           kill nullvec;
4496        }
4497        else
4498        {
4499           "ERROR:@ce[4] hat falsche Laenge - nicht-trivialer Fall";
4500           ~;
4501        }
4502      }
4503      if((typeof(@ce[4])=="intvec") || (typeof(@ce[4])=="intmat"))
4504      {
4505        BO[9]=@ce[4];
4506      }
4507//---------------------------------------------------------------------------
4508// various checks and debug output
4509//---------------------------------------------------------------------------
4510      if((debu) || (praes_stop))
4511      {
4512//--- Show BO of this step
4513         "++++++++++++++ Overview of Current Chart +++++++++++++++++++++++";
4514         "Current number of final charts:",size(endRings);
4515         "Total number of charts currently in chart-tree:",size(allRings);
4516         "Index of the current chart in chart-tree:",i;
4517         "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++";
4518         showBO(BO);
4519         "-------------------------- Upcoming Center ---------------------";
4520         interred(cent);
4521         "----------------------------------------------------------------";
4522      }
4523      if(praes_stop)
4524      {
4525         ~;
4526      }
4527      if(debu)
4528      {
4529//--- various checks, see output for comments
4530         if(size(BO[1])>0)
4531         {
4532            if(deg(BO[1][1])==0)
4533            {
4534               "!!!  W is empty  !!!";
4535               path;
4536               setring R;
4537               kill S;
4538               list result=endRings,allRings;
4539               return(result);
4540            }
4541            if(deg(std(slocusE(BO[1]))[1])>0)
4542            {
4543               "!!!  W not smooth  !!!";
4544               path;
4545               setring R;
4546               kill S;
4547               list result=endRings,allRings;
4548               return(result);
4549            }
4550         }
4551         if((!loca)&&(!locaT))
4552         {
4553            if(deg(std(slocusE(cent+BO[1]))[1])>0)
4554            {
4555               "!!!  Center not smooth  !!!";
4556               path;
4557               std(cent+BO[1]);
4558               ~;
4559               setring R;
4560               kill S;
4561               list result=endRings,allRings;
4562               return(result);
4563            }
4564         }
4565         for(j=1;j<=size(BO[4]);j++)
4566         {
4567            if(deg(BO[4][j][1])>0)
4568            {
4569               if(deg(std(slocusE(BO[4][j]+BO[1]))[1])>0)
4570               {
4571                  "!!!  exceptional divisor is not smooth  !!!";
4572                  path;
4573                  setring R;
4574                  kill S;
4575                  list result=endRings,allRings;
4576                  return(result);
4577               }
4578            }
4579         }
4580         if((!loca)&&(!locaT))
4581         {
4582            if((norC(BO,cent))&&(size(reduce(cent,Jstd))!=0))
4583            {
4584               "!!!  this chart is already finished  !!!";
4585               cent=BO[2];
4586               ~;
4587            }
4588         }
4589      }
4590//----------------------------------------------------------------------------
4591// Do the blow up
4592//----------------------------------------------------------------------------
4593//!!!! Change this as soon as there is time!!!
4594//!!!! quick and dirty bug fix for old shortcut which has not yet been killed
4595if((dim(std(cent))==0)&&defined(shortcut)) {kill shortcut;}
4596//!!! end of bugfix
4597      if(size(reduce(cent,Jstd))!=0)
4598      {
4599//--- center does not equal J
4600         tmpList=blowUpBO(BO,cent,extra);
4601         if((debu)&&(!loca)&&(!locaT))
4602         {
4603//--- test it, if debu is set
4604            if(!testBlowUp(BO,cent,tmpList,i,extra))
4605            {
4606               "non-redundant chart has been killed!";
4607               ~;
4608            }
4609         }
4610//--- extend the list of all rings
4611         allRings[size(allRings)+1..size(allRings)+size(tmpList)]=
4612                         tmpList[1..size(tmpList)];
4613         for(j=1;j<=size(tmpList);j++)
4614         {
4615            def Q=allRings[size(allRings)-j+1];
4616            setring Q;
4617            def path=imap(S,path);
4618            path=path,[i,size(tmpList)-j+1];
4619            export path;
4620            setring S;
4621            kill Q;
4622         }
4623         kill tmpList;
4624         list tmpList;
4625      }
4626      else
4627      {
4628//--- center equals J
4629           k=0;
4630           for(j=1;j<=size(BO[6]);j++)
4631           {
4632              if(BO[6][j]!=1)
4633              {
4634//--- there is an E_i which meets J in this chart
4635                k=1;
4636                break;
4637              }
4638           }
4639           if((k)||(extra==2))
4640           {
4641//--- chart finished, non-redundant
4642              endRings[size(endRings)+1]=S;
4643           }
4644      }
4645      kill pr;
4646      setring R;
4647      kill S;
4648   }
4649//---------------------------------------------------------------------------
4650// set up the result, test it (if debu is set) and return it
4651//---------------------------------------------------------------------------
4652   list result=endRings,allRings;
4653   if(debu)
4654   {
4655      "============= result will be tested ==========";
4656      "                                              ";
4657      "the number of charts obtained:",size(endRings);
4658      if(locaT){loca=2;}
4659      int tes=testRes(J,endRings,loca);
4660      if(tes)
4661      {
4662         "=============     result is o.k.    ==========";
4663      }
4664      else
4665      {
4666         "============     result is wrong    =========="; ~;
4667      }
4668   }
4669   kill debugCenter,debugBlowUp,debugCoeff,debug_Inters_E;
4670   if(locaT){kill @p;}
4671   kill locaT;
4672   return(result);
4673}
4674example
4675{"EXAMPLE:";
4676   echo = 2;
4677   ring R=0,(x,y,z),dp;
4678   ideal J=x3+y5+yz2+xy4;
4679   list L=resolve(J,0);
4680   def Q=L[1][7];
4681   setring Q;
4682   showBO(BO);
4683}
4684//////////////////////////////////////////////////////////////////////////
4685//static
4686proc CompMeetsE(ideal J, list E)
4687"Internal procedure - no help and no example available
4688"
4689{
4690   int i;
4691   for(i=1;i<=size(E);i++)
4692   {
4693      if(deg(std(E[i])[1])!=0)
4694      {
4695         if(deg(std(J+E[i])[1])!=0)
4696         {
4697            return(1);
4698         }
4699      }
4700   }
4701   return(0);
4702}
4703
4704//========================================================================
4705//--------------  procedures for testing the result ----------------------
4706//                (not yet commented)
4707//========================================================================
4708
4709//////////////////////////////////////////////////////////////////////////
4710static proc testRes(ideal J,list L,int loca)
4711"Internal procedure - no help and no example available
4712"
4713{
4714   int loc;
4715   if(defined(locaT)){loc=locaT;}
4716   if(loc){loca=0;}
4717   def R=basering;
4718   ideal M=maxideal(1);
4719   int i,j,tr;
4720   for(i=1;i<=size(L);i++)
4721   {
4722      def Q=L[i];
4723      setring Q;
4724      ideal J=BO[2];
4725      list E=BO[4];
4726      map phi=R,BO[5];
4727      ideal K=phi(J)+BO[1];
4728      ideal stTK=std(K);
4729
4730      if(loca)
4731      {
4732         ideal M=phi(M)+BO[1];
4733         ideal stTM=std(M);
4734      }
4735      for(j=1;j<=size(E);j++)
4736      {
4737         if(deg(E[j][1])>0)
4738         {
4739            stTK=sat(stTK,E[j])[1];
4740         }
4741         if(loca)
4742         {
4743            stTM=sat(stTM,E[j])[1];
4744         }
4745      }
4746      ideal sL=slocusE(J);
4747      if(loca){sL=sL+stTM;}
4748      ideal sLstd=std(sL);
4749      if(deg(sLstd[1])>0)
4750      {
4751         if(!loc)
4752         {
4753            "J is not smooth";i;
4754            setring R;
4755            return(0);
4756         }
4757         if(size(reduce(@p,std(radical(sLstd))))>0)
4758         {
4759            "J is not smooth";i;
4760            setring R;
4761            return(0);
4762         }
4763      }
4764      if(!((size(reduce(J,std(stTK)))==0)
4765           &&(size(reduce(stTK,std(J)))==0)))
4766      {
4767         "map is wrong";i;
4768         setring R;
4769         return(0);
4770      }
4771      if(loc){tr=transversalT(J,E,@p);}
4772      else{tr=transversalT(J,E);}
4773      if(!tr)
4774      {
4775         "E not transversal with J";i;
4776         setring R;
4777         return(0);
4778      }
4779      if(!normalCross(E))
4780      {
4781         "E not normal crossings";i;
4782         setring R;
4783         return(0);
4784      }
4785      for(j=1;j<=size(E);j++)
4786      {
4787         if(deg(E[j][1])>0){E[j]=E[j]+J;}
4788      }
4789      if(!normalCross(E))
4790      {
4791         "E not normal crossings with J";i;
4792         setring R;
4793         return(0);
4794      }
4795      kill J,E,phi,K,stTK;
4796      if(loca){kill M,stTM;}
4797      setring R;
4798      kill Q;
4799   }
4800   return(1);
4801}
4802//////////////////////////////////////////////////////////////////////////////
4803static proc testBlowUp(list BO,ideal cent,list tmpList, int j, int extra)
4804{
4805   def R=basering;
4806   int n=nvars(basering);
4807   int i;
4808   if((extra!=3)&&(extra!=2))
4809   {
4810      ideal K=BO[1],BO[2],cent;
4811      for(i=1;i<=size(BO[4]);i++)
4812      {
4813        K=K,BO[4][i];
4814      }
4815      list N=findvars(K,0);
4816      //list N=findvars(BO[2],0);
4817      if(size(N[1])<n)
4818      {
4819         string newvar=string(N[1]);
4820         execute("ring R1=("+charstr(R)+"),("+newvar+"),dp;");
4821         list BO=imap(R,BO);
4822         ideal cent=imap(R,cent);
4823         n=nvars(R1);
4824      }
4825      else
4826      {
4827         def R1=basering;
4828      }
4829   }
4830   else
4831   {
4832      def R1=basering;
4833   }
4834
4835   i=0;
4836   ideal T=cent;
4837   ideal TW;
4838   for(i=1;i<=size(tmpList);i++)
4839   {
4840      def Q=tmpList[i];
4841      setring Q;
4842      map phi=R1,lastMap;
4843      ideal TE=radical(slocusE(BO[2]));
4844      setring R1;
4845      TW=preimage(Q,phi,TE);
4846      T=intersect(T,TW);
4847      kill Q;
4848   }
4849   ideal sL=intersect(slocusE(BO[2]),cent);
4850   if(size(reduce(sL,std(radical(T))))>0){setring R;return(0);}
4851   if(size(reduce(T,std(radical(sL))))>0){setring R;return(0);}
4852   setring R;
4853   return(1);
4854}
4855//////////////////////////////////////////////////////////////////////////////
4856static proc normalCross(list E,list #)
4857"Internal procedure - no help and no example available
4858"
4859{
4860   int loc;
4861   if((defined(locaT))&&(defined(@p)))
4862   {
4863      loc=1;
4864      ideal pp=@p;
4865   }
4866   int i,d,j;
4867   int n=nvars(basering);
4868   list E1,E2;
4869   ideal K,M,Estd,cent;
4870   intvec v,w;
4871   if(size(#)>0){cent=#[1];}
4872
4873   for(i=1;i<=size(E);i++)
4874   {
4875      Estd=std(E[i]);
4876      if(deg(Estd[1])>0)
4877      {
4878         E1[size(E1)+1]=Estd;
4879      }
4880   }
4881   E=E1;
4882   for(i=1;i<=size(E);i++)
4883   {
4884      v=i;
4885      E1[i]=list(E[i],v);
4886   }
4887   list ll;
4888   int re=1;
4889   int ok;
4890   while(size(E1)>0)
4891   {
4892      K=E1[1][1];
4893      v=E1[1][2];
4894      attrib(K,"isSB",1);
4895      E1=delete(E1,1);
4896      d=n-dim(K);
4897      M=minor(jacob(K),d)+K;
4898      if(deg(std(M)[1])>0)
4899      {
4900         if(size(#)>0)
4901         {
4902            if(size(reduce(M,std(cent)))>0)
4903            {
4904               ll[size(ll)+1]=std(M);
4905            }
4906            else
4907            {
4908               ok=1;
4909            }
4910         }
4911         if(!loc)
4912         {
4913            re=0;
4914         }
4915         else
4916         {
4917            if(size(reduce(pp,std(radical(M))))>0){re=0;}
4918         }
4919      }
4920      for(i=1;i<=size(E);i++)
4921      {
4922         for(j=1;j<=size(v);j++){if(v[j]==i){break;}}
4923         if(j<=size(v)){if(v[j]==i){i++;continue;}}
4924         Estd=std(K+E[i]);
4925         w=v;
4926         if(deg(Estd[1])==0){i++;continue;}
4927         if(d==n-dim(Estd))
4928         {
4929            if(size(#)>0)
4930            {
4931               if(size(reduce(Estd,std(cent)))>0)
4932               {
4933                  ll[size(ll)+1]=Estd;
4934               }
4935               else
4936               {
4937                  ok=1;
4938               }
4939            }
4940            if(!loc)
4941            {
4942               re=0;
4943            }
4944            else
4945            {
4946               if(size(reduce(pp,std(radical(M))))>0){re=0;}
4947            }
4948         }
4949         w[size(w)+1]=i;
4950         E2[size(E2)+1]=list(Estd,w);
4951      }
4952      if(size(E2)>0)
4953      {
4954         if(size(E1)>0)
4955         {
4956            E1[size(E1)+1..size(E1)+size(E2)]=E2[1..size(E2)];
4957         }
4958         else
4959         {
4960            E1=E2;
4961         }
4962      }
4963      kill E2;
4964      list E2;
4965   }
4966/*
4967   if((!ok)&&(!re)&&(size(#)==1))
4968   {
4969
4970      "the center is wrong";
4971      "it could be one of the following list";
4972      ll;
4973       ~;
4974   }
4975*/
4976   if((!ok)&&(!re)&&(size(#)==2))
4977   {
4978      return(2);   //for Center correction
4979   }
4980   return(re);
4981}
4982//////////////////////////////////////////////////////////////////////////////
4983static proc normalCrossB(ideal J,list E,ideal V)
4984"Internal procedure - no help and no example available
4985"
4986{
4987   int i,d,j;
4988   int n=nvars(basering);
4989   list E1,E2;
4990   ideal K,M,Estd;
4991   intvec v,w;
4992
4993   for(i=1;i<=size(E);i++)
4994   {
4995      Estd=std(E[i]+J);
4996      if(deg(Estd[1])>0)
4997      {
4998         E1[size(E1)+1]=Estd;
4999      }
5000   }
5001   E=E1;
5002   for(i=1;i<=size(E);i++)
5003   {
5004      v=i;
5005      E1[i]=list(E[i],v);
5006   }
5007   list ll;
5008   int re=1;
5009
5010   while((size(E1)>0)&&(re==1))
5011   {
5012      K=E1[1][1];
5013      v=E1[1][2];
5014      attrib(K,"isSB",1);
5015      E1=delete(E1,1);
5016      d=n-dim(K);
5017      M=minor(jacob(K),d)+K;
5018      if(deg(std(M+V)[1])>0)
5019      {
5020         re=0;
5021         break;
5022      }
5023      for(i=1;i<=size(E);i++)
5024      {
5025         for(j=1;j<=size(v);j++){if(v[j]==i){break;}}
5026         if(j<=size(v)){if(v[j]==i){i++;continue;}}
5027         Estd=std(K+E[i]);
5028         w=v;
5029         if(deg(Estd[1])==0){i++;continue;}
5030         if(d==n-dim(Estd))
5031         {
5032            if(deg(std(Estd+V)[1])>0)
5033            {
5034               re=0;
5035               break;
5036            }
5037         }
5038         w[size(w)+1]=i;
5039         E2[size(E2)+1]=list(Estd,w);
5040      }
5041      if(size(E2)>0)
5042      {
5043         if(size(E1)>0)
5044         {
5045            E1[size(E1)+1..size(E1)+size(E2)]=E2[1..size(E2)];
5046         }
5047         else
5048         {
5049            E1=E2;
5050         }
5051      }
5052      kill E2;
5053      list E2;
5054   }
5055   return(re);
5056}
5057//////////////////////////////////////////////////////////////////////////////
5058static proc norC(list BO,ideal cent)
5059"Internal procedure - no help and no example available
5060"
5061{
5062   int j;
5063   list  E=BO[4];
5064   ideal N=BO[2];
5065   if(BO[3][1]>1){return(0);}
5066   if(deg(std(slocusE(BO[2]))[1])>0){return(0);}
5067   if(!transversalT(N,E)){return(0);}
5068   for(j=1;j<=size(E);j++){if(deg(E[j][1])>0){E[j]=E[j]+N;}}
5069   if(!normalCross(E,cent)){return(0);}
5070   return(1);
5071}
5072//////////////////////////////////////////////////////////////////////////////
5073static proc specialReduce(ideal I,ideal J,poly p)
5074{
5075   matrix M;
5076   int i,j;
5077   for(i=1;i<=ncols(I);i++)
5078   {
5079      M=coeffs(I[i],@z);
5080      I[i]=0;
5081      for(j=1;j<=nrows(M);j++)
5082      {
5083         I[i]=I[i]+M[j,1]*p^(nrows(M)-j);
5084      }
5085      I[i]=reduce(I[i],J);
5086   }
5087   return(I);
5088}
Note: See TracBrowser for help on using the repository browser.