source: git/Singular/LIB/resjung.lib @ e31d29

spielwiese
Last change on this file since e31d29 was e31d29, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13455 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 26.9 KB
Line 
1version="$Id$";
2category="Algebraic Geometry";
3info="
4LIBRARY:  resjung.lib      Resolution of surface singularities (Desingularization)
5                           by Jung's Algorithm
6AUTHORS:  Philipp Renner,   philipp_renner@web.de
7@*        Anne Fruehbis-Krueger, anne@math.uni-hannover.de
8OVERVIEW: This library implements resolution of singularities by Jung's algorithm,
9@*        which is only applicable to surfaces and persues the following strategy:
10@*        1) project surface to the plane
11@*        2) resolve singularities of the branch locus
12@*        3) pull-back the original surface by this resolution morphism
13@*        4) normalize the resulting surface so that the remaining singularities
14@*           are of Hirzebruch-Jung type
15@*        5) resolve Hirzebruch-Jung singularities explicitly
16@*        Currently, the Hirzebruch-Jung singularities are resolved by calling
17@*        the procedure resolve from the library resolve.lib, because this is not
18@*        overly expensive and the original last step of Jung's algorithm is not
19@*        implemented yet.
20REFERENCES:
21[1] Jung, H.: Darstellung der Funktionen eines algebraischen Koerpers zweier unabhaengigen
22Veraenderlichen x,y in der Umgebung x=a, y= b, Journal fuer Reine und Angewandte Mathematik
23133,289-314 (1908)
24@*  (the origin of this method)
25@*[2] J.Kollar: Lectures on Resolution of Singularities, Princeton University Press (2007)
26@*  (contains large overview over various known methods for curves and surfaces as well as
27   a detailed description of the approach in the general case)
28
29PROCEDURES:
30 jungresolve(J,i)  computes a resolution of the surface given by the
31                   ideal J using Jungs Method
32
33 clocus(J)         computes the critical locus of the projection of V(J)
34                   onto the coordinate plane of the last two coordinates
35 embR(C)           computes a strong embedded resolution of the plane curve V(C)
36 jungnormal(J,i)   computes intermediate step in Jung's algorithm
37                   such that all singularities are of Hirzebruch-Jung type
38";
39
40LIB "resolve.lib";
41LIB "mregular.lib";
42LIB "sing.lib";
43LIB "normal.lib";
44LIB "primdec.lib";
45
46//////////////////////////////////////////////////////////////////////////////
47//----------------------------------------------------------------------------
48//Critical locus for the Weierstrass map induced by the noether normalization
49//----------------------------------------------------------------------------
50proc clocus(ideal id)
51"USAGE:  clocus(ideal J);
52@*       J = ideal
53ASSUME:  J  = two dimensional ideal in noether position with respect
54         to the last two variables
55RETURN:  A ring containing the ideal Clocus respresenting the critical
56         locus of the projection V(J)-->C^2 onto the coordinate plane
57         of the last two coordinates
58EXAMPLE: example clocus;  shows an example
59"
60{
61  def A = basering;
62  int n = nvars(A);
63  list l = equidim(id);
64  int i,j;
65  int k = size(l);
66  ideal LastTwo = var(n-1),var(n);
67  ideal lowdim = 1;   //the components of id with dimension smaller 2
68  if(k>1){
69    for(j=1;j<k;j++){
70      lowdim = intersect(lowdim,l[j]);
71    }
72  }
73  //lowdim = radical(lowdim);  // affects performance
74  ideal I = l[size(l)];
75  poly product=1;
76  kill l;
77  for(i=1; i < n-1; i++){
78  //elimination of all variables exept var(i),var(n-1),var(n)
79    intvec v;
80    for(j=1; j < n-1; j++){
81      if(j<>i){
82        v[j]=1;
83      }
84      else{
85        v[j]=0;
86      }
87    }
88    v[size(v)+1]=0;
89    v[size(v)+1]=0;
90    if(defined(ringl)) {kill ringl;}
91    list ringl = ringlist(A);
92    list l;
93    l[1] = "a";
94    l[2] = v;
95    list ll = insert(ringl[3],l);
96    ringl[3]=ll;
97    kill l,ll;
98    def R = ring(ringl); //now x_j > x_i > x_n-1 > x_n forall j != i,n-1,n
99    setring R;
100    ideal J = groebner(fetch(A,I));  //this eliminates the variables
101    setring A;
102    ideal J = fetch(R,J);
103    attrib(J,"isPrincipal",0);
104    if(size(J)==1){
105      attrib(J,"isPrincipal",1);
106    }
107    int index = 1;
108    if(attrib(J,"isPrincipal")==0){
109      setring R;
110      for(j = 1;j<=size(J);j++){   //determines the monic polynomial
111                                   //in var(i) with coefficents in C2
112        if(defined(w)) {kill w;}
113        intvec w = leadexp(J[j]);
114        attrib(w,"isMonic",1);
115        for(k = 1;k<=size(w);k++){
116          if(w[k] <> 0 && k <> i){
117            attrib(w,"isMonic",0);
118            break;
119          }
120        }
121        if(attrib(w,"isMonic")==1){
122          index = j;
123          break;
124        }
125        kill w;
126      }
127      setring A;
128    }
129    product = product*resultant(J[index],diff(J[index],var(i)),var(i));
130                       //Product of the discriminants, which lies in C2
131    kill index,J,v;
132  }
133  ring C2 = 0,(var(n-1),var(n)),dp;
134  setring C2;
135  ideal Clocus = imap(A,product);      //the critical locus is contained in this
136  ideal I = preimage(A,LastTwo,lowdim);
137  Clocus = radical(intersect(Clocus,I));  //radical is necessary since the
138                                          //resultant is not reduced in general
139  export(Clocus);
140  return(C2);
141}
142example
143{"EXAMPLE:";
144   echo = 2;
145   ring R=0,(x,y,z),dp;
146   ideal I=x2-y3z3;
147   list li=clocus(I);
148   def S=li[1];
149   setring S;
150   Clocus;
151}
152///////////////////////////////////////////////////////////////////////////////
153//-----------------------------------------------------------------------------
154// Build the fibre product of the embedded resolution and
155// the coordinate ring of the variety
156//-----------------------------------------------------------------------------
157static proc buildFP(list embR,ideal NoetherN, map phi){
158  def A = basering;
159  int i,j,k;
160  list fibreP;
161  int n = nvars(A);
162  for(i=1;i<=size(embR);i++){
163    if(defined(R)) {kill R;}
164    def R = embR[i];
165    setring R;
166    list temp = ringlist(A);
167// create data for the new ring
168// e.g. if A=K[x_1,..,x_n] and R=K[y_1,..,y_m], K[x_1,..,x_n-2,y_1,..,y_m]
169    for(j = 1; j<= nvars(R);j++){
170       string st = string(var(j));
171       temp[2][n-2+j] = st;
172       kill st;
173    }
174    temp[4]    = BO[1];
175    ideal J = BO[5];        //ideal of the resolution map
176    export(J);
177    int m = size(J);
178    def R2 = ring(temp);
179    kill temp;
180    setring R2;
181    ideal Temp=0;           //defines map from R to R2 which is the inclusion
182    for(k=n-1;k<n-1+nvars(R);k++){
183      Temp = Temp + ideal(var(k));
184    }
185    map f = R,Temp;
186    kill Temp;
187    ideal FibPMI = ideal(0);       //defines the map from A to R2
188    for(k=1;k<=nvars(A)-m;k++){
189      FibPMI=FibPMI+var(k);
190    }
191    FibPMI= FibPMI+ideal(f(J));
192    map FibMap = A,FibPMI;
193    kill f,FibPMI;
194    ideal TotalT = groebner(FibMap(NoetherN));
195    ideal QIdeal = TotalT;
196    export(QIdeal);
197    ideal FibPMap = ideal(FibMap(phi));
198    ideal BMap = FibPMap;
199    export(BMap);
200    fibreP[i] = R2;
201    setring R;
202    kill J,R,R2,m;
203  }
204  return(fibreP);
205};
206
207///////////////////////////////////////////////////////////////////////////////
208//-----------------------------------------------------------------------------
209// embedded resolution for curves -- optimized for our situation
210//-----------------------------------------------------------------------------
211proc embR(ideal C)
212"USAGE:  embR(ideal C);
213@*       C = ideal
214ASSUME:  C  = ideal of plane curve
215RETURN:  a list l of rings
216         l[i] is a ring containing a basic object BO, the result of the
217         resolution.
218THEORY:  Given a plain curve C, an embedded resolution of this curve
219@*         by means of a sequence of blow ups at singular points of C is
220@*         computed such that the resulting total transform consists
221@*         of non-singular components and only has strict normal crossings.
222@*         The result of each blow up is represented by means of
223@*         affine charts and hence also the final result is represented in
224@*         charts, which are collected in a list l of rings. Each ring
225@*         in this list corresponds to one chart. The collection of
226@*         data describing the ambient space, the strict transform, the
227@*         exceptional divisors and the combination of the blow ups
228@*         leading to this chart is contained in the list BO which resides
229@*         in this ring.
230NOTE:      - The best way to look at the data contained in BO is the
231@*           procedure showBO from the library resolve.lib.
232@*         - The algorithm does not touch normal crossings of V(C).
233EXAMPLE: example embR; shows an example
234"
235{
236  ideal J = 1;
237  attrib(J,"iswholeRing",1);
238  list primdec = equidim(C);
239  if(size(primdec)==2){
240// zero dimensional components of the discrimiant curve
241// are smooth an cross normally so they can be ignored
242// in the resolution process
243    ideal Lowdim = radical(primdec[1]);
244  }
245  else{
246    J=radical(C);
247  }
248  kill primdec;
249  list l;
250  list BO = createBO(J,l);
251  kill J,l;
252  list result = resolve2(BO);
253  if(defined(Lowdim)){
254    for(int i = 1;i<=size(result);i++){
255// had zero dimensional components which are now added to the end result
256      if(defined(R)) {kill R;}
257      def R = result[i];
258      setring R;
259      map f = R2,BO[5];
260      BO[2]=BO[2]*f(Lowdim);
261      kill R,f;
262    }
263  }
264  return(result);
265}
266example
267{"EXAMPLE:";
268  echo=2;
269  ring R=0,(x,y),dp;
270  ideal C=x2-y3;
271  list li=embR(C);
272  def S=li[1];
273  setring S;
274  showBO(BO);
275}
276///////////////////////////////////////////////////////////////////////////////
277static proc resolve2(list BO){
278// computes an embedded resolution for the basic object BO
279// and returns a list of rings with BO -- specifically optimized
280// to our situation
281  def H = basering;
282  int i,j,k;
283  setring H;  attrib(BO[2],"smoothC",0);
284  export(BO);
285  list result;
286  result[1]=H;
287  attrib(result[1],"isResolved",0);    // has only simple normal crossings
288  attrib(result[1],"smoothC",0);       // has smooth components
289  int safety=0;                        // number of runs restricted to 30
290  while(1){
291    int count2 = 0;         // counts the number of smooth charts
292    int p = size(result);
293    for(j = 1;j<=p;j++){
294      if(attrib(result[j],"isResolved")==0){
295        if(defined(R)){kill R;}
296        def R = result[j];
297        setring R;
298        if(attrib(result[j],"smoothC")==0){
299// has possibly singular components so choose a singular point and blow up
300          list primdecPC = primdecGTZ(BO[2]);
301          attrib(result[j],"smoothC",1);
302          for(i = 1;i<=size(primdecPC);i++){
303            ideal Sl = groebner(slocus(primdecPC[i][2]));
304            if(deg(NF(1,Sl))!=-1){
305              list primdecSL = primdecGTZ(Sl);
306              for(int h =1;h<=size(primdecSL);h++){
307                attrib(primdecSL[h],"isRational",1);
308              }
309              kill h;
310              if(!defined(index)){int index = 1;}
311              if(defined(blowup)){kill blowup;}
312              list blowup = blowUpBO(BO,primdecSL[index][2],3);
313// if it has only a rational singularity, blow it up,
314// else choose some arbitary singular point
315              if(attrib(primdecSL[1],"isRational")==0){
316// if we blew up a non rational singularity, the exeptional divisors
317// are reduzible, so we need to separate them
318                for(k=1;k<=size(blowup);k++){
319                  def R2=blowup[k];
320                  setring R2;
321                  list L;
322                  for(int l = 1;l<=size(BO[4]);l++){
323                    list primdecED=primdecGTZ(BO[4][l]);
324                    L = L + primdecED;
325                    kill primdecED;
326                  }
327                  kill l;
328                  BO[4] = L;
329                  blowup[k]=R2;
330                  kill L,R2;
331                }
332              }
333              kill primdecSL;
334              list hlp;
335              for(k = 1;k<j;k++){
336                hlp[k]=result[k];
337                attrib(hlp[k],"isResolved",attrib(result[k],"isResolved"));
338                attrib(hlp[k],"smoothC",attrib(result[k],"smoothC"));
339              }
340              for(k =1;k<=size(blowup);k++){
341                hlp[size(hlp)+1]=blowup[k];
342                attrib(hlp[size(hlp)],"isResolved",0);
343                attrib(hlp[size(hlp)],"smoothC",0);
344              }
345              for(k = j+1;k<=size(result);k++){
346                hlp[size(hlp)+1]=result[k];
347                attrib(hlp[size(hlp)],"isResolved",attrib(result[k],"isResolved"));
348                attrib(hlp[size(hlp)],"smoothC",attrib(result[k],"smoothC"));
349              }
350              result = hlp;
351              kill hlp;
352              i=size(primdecPC);
353            }
354            else{
355              attrib(result[j],"smoothC",1);
356            }
357            kill Sl;
358          }
359          kill primdecPC;
360          j=p;
361          break;
362        }
363        else{
364// if it has smooth components, determine all the intersection points
365// and check whether they are snc or not
366          int count = 0;
367          ideal Collect = BO[2];
368          for(i = 1;i<=size(BO[4]);i++){
369            Collect = Collect*BO[4][i];
370          }
371          list primdecSL = primdecGTZ(slocus(Collect));
372          for(k = 1;k<=size(primdecSL);k++){
373             attrib(primdecSL[k],"isRational",1);
374          }
375          if(defined(blowup)){kill blowup;}
376          list blowup = blowUpBO(BO,primdecSL[1][2],3);
377          if(attrib(primdecSL[1],"isRational")==0){
378            for(k=1;k<=size(blowup);k++){
379              def R2=blowup[k];
380              setring R2;
381              list L;
382              for(int l = 1;l<=size(BO[4]);l++){
383                list primdecED=primdecGTZ(BO[4][l]);
384                L = L + primdecED;
385                kill primdecED;
386              }
387              kill l;
388              BO[4] = L;
389              blowup[k]=R2;
390              kill L,R2;
391            }
392          }
393          kill Collect;
394          for(i=1;i<=size(primdecSL);i++){
395            list L = BO[4];
396            L[size(L)+1]=BO[2];
397            for(int l = 1;l<=size(L);l++){
398              if(L[l][1]==1){L=delete(L,l);}
399            }
400            kill l;
401            if(normalCrossing(ideal(0),L,primdecSL[i][2])==0){
402              if(defined(blowup)){kill blowup;}
403              list blowup = blowUpBO(BO,primdecSL[i][2],3);
404              list hlp;
405              for(k = 1;k<j;k++){
406                hlp[k]=result[k];
407                attrib(hlp[k],"isResolved",attrib(result[k],"isResolved"));
408                attrib(hlp[k],"smoothC",attrib(result[k],"smoothC"));
409              }
410              for(k =1;k<=size(blowup);k++){
411                hlp[size(hlp)+1]=blowup[k];
412                attrib(hlp[size(hlp)],"isResolved",0);
413                attrib(hlp[size(hlp)],"smoothC",1);
414              }
415              for(k = j+1;k<=size(result);k++){
416                hlp[size(hlp)+1]=result[k];
417                attrib(hlp[size(hlp)],"isResolved",attrib(result[k],"isResolved"));
418                attrib(hlp[size(hlp)],"smoothC",attrib(result[k],"smoothC"));
419              }
420              result = hlp;
421              kill hlp;
422              j = p;
423              break;
424            }
425            else{
426              count++;
427            }
428            kill L;
429          }
430          if(count == size(primdecSL)){
431            attrib(result[j],"isResolved",1);
432          }
433          kill count,primdecSL;
434        }
435        kill R;
436      }
437      else{
438        count2++;
439      }
440    }
441    if(count2==size(result)){
442      break;
443    }
444    kill count2,p;
445    safety++;
446  }
447  return(result);
448};
449///////////////////////////////////////////////////////////////////////////////
450static proc jungfib(ideal id, int noeth)
451"USAGE:  jungfib(ideal J, int i);
452@*       J = ideal
453@*       i = int
454ASSUME:  J  = two dimensional ideal
455RETURN:  a list l of rings
456         l[i] is a ring containing two Ideals: QIdeal and BMap.
457         BMap defines a birational morphism from V(QIdeal)-->V(J), such that
458         V(QIdeal) has only quasi-ordinary singularities.
459         If i!=0 then it's assumed that J is in noether position with respect
460         to the last two variables.
461         If i=0 the algorithm computes a coordinate change such that J is in
462         noether position.
463EXAMPLE: none, as it is a static procedure
464"
465{
466  int i;
467  if(!defined(noeth)){
468    int noeth = 0;
469  }
470  if(noeth <> 1){print("//WARNING: Noether normalization can make the algorithm unstable");}
471  ideal I = std(radical(id));
472  def A = basering;
473  int n = nvars(A);
474  if(deg(NF(1,groebner(slocus(id)))) == -1){
475    list result;
476    ideal QIdeal = I;
477    ideal BMap = maxideal(1);
478    export(QIdeal);
479    export(BMap);
480    result[1] = A;
481    return(result);
482  }
483  if(char(A) <> 0){ERROR("only works for characterisitc 0");}  //dummy check
484  if(dim(I)<> 2){ERROR("dimension is unequal 2");}  //dummy check
485
486// Noether Normalization
487  if(noeth == 0){
488    if(n==3){
489      int pos = NoetherP_test(I);
490      if(pos ==0){
491        if(size(I) == 1){
492          ideal noethpos = var(1)+var(3),var(2)+var(3),var(3);
493        }
494        else{
495          ideal noethpos = NoetherPosition(I);
496        }
497        map phi = A,noethpos;
498        kill noethpos;
499      }
500      else{
501        ideal NoetherPos = var(pos);
502        for(i = 1;i<=3;i++){
503          if(i<>pos){
504            NoetherPos = NoetherPos + var(i);
505          }
506        }
507        map phi = A,NoetherPos;
508        kill i,pos,NoetherPos;
509      }
510    }
511    else{
512      map phi = A,NoetherPosition(I);
513    }
514    ideal NoetherN = ideal(phi(I));
515//image of id under the NoetherN coordinate change
516  }
517  else{
518    ideal NoetherN = I;
519    map phi = A,maxideal(1);
520  }
521  kill I;
522//Critical Locus
523  def C2 = clocus(NoetherN);
524  setring C2;
525
526//dim of critical locus is 0 then the normalization is an resolution
527  if(dim(Clocus) == 0){
528    setring A;
529    list nor = normal(NoetherN);
530    list result;
531    for(i = 1;i<=size(nor[1]);i++){
532      def R = nor[1][i];
533      setring R;
534      ideal QIdeal = norid;
535      ideal BMap = BMap;
536      export(QIdeal);
537      export(BMap);
538      result[size(result)+1] = R;
539      kill R;
540    }
541    print("This is a resolution.");
542    return(result);
543  }
544
545// dim of critical locus is 1, so compute embedded resolution of the discriminant curve
546  list embRe = embR(Clocus);
547
548// build the fibreproduct
549  setring A;
550  list fibreP = buildFP(embRe,NoetherN,phi);
551// a list of lists, where fibreP[i] contains the information
552// concerning  the i-th chart of the fibrepoduct
553// fibreP[i] is the ring; QIdeal the quotientideal; BMap is the map from A
554  return(fibreP);
555}
556///////////////////////////////////////////////////////////////////////////////
557
558proc jungnormal(ideal id,int noeth)
559"USAGE:  jungnormal(ideal J, int i);
560@*       J = ideal
561@*       i = int
562ASSUME:  J  = two dimensional ideal
563RETURN:  a list l of rings
564         l[k] is a ring containing two Ideals: QIdeal and BMap.
565         BMap defines a birational morphism from V(QIdeal)-->V(J), such that
566         V(QIdeal) has only singularities of Hizebuch-Jung type.
567         If i!=0 then it's assumed that J is in noether position with respect
568         to the last two variables.
569         If i=0 the algorithm computes a coordinate change such that J is in
570         noether position.
571EXAMPLE: example jungnormal; shows an example
572"
573{
574  def A = basering;
575  list fibreP = jungfib(id,noeth);
576  list result;
577  int i,j;
578  for(i =1;i<=size(fibreP);i++){
579    def R1 = fibreP[i];
580    setring R1;
581    map f1 = A,BMap;
582    def nor = normal(QIdeal);
583    for(j = 1;j<=size(nor[1]);j++){
584      def R2 = nor[1][j];
585      setring R2;
586      map f2 = R1,normap;
587      ideal BMap = ideal(f2(f1));
588      ideal QIdeal = norid;
589      export(BMap);
590      export(QIdeal);
591      result[size(result)+1] = R2;
592      setring R1;
593      kill R2;
594    }
595    kill R1;
596  }
597  return(result);
598}
599example
600{"EXAMPLE:";
601   echo = 2;
602   ring R=0,(x,y,z),dp;
603   ideal J=x2+y3z3;
604   list li=jungnormal(J,1);
605   li;
606   def S=li[1];
607   setring S;
608   QIdeal;
609   BMap;
610}
611///////////////////////////////////////////////////////////////////////////////
612
613proc jungresolve(ideal id,int noeth)
614"USAGE:  jungresolve(ideal J, int i);
615@*       J = ideal
616@*       i = int
617ASSUME:  J  = two dimensional ideal
618RETURN:  a list l of rings
619         l[k] is a ring containing two Ideals: QIdeal and BMap.
620         BMap defines a birational morphism from V(QIdeal)-->V(J), such that
621         V(QIdeal) is smooth. For this the algorithm computes first
622         a representation of V(J) with Hirzebruch-Jung singularities
623         and then it currently uses Villamayor's algorithm to resolve
624         these singularities.
625         If i!=0 then it's assumed that J is in noether position with respect
626         to the last two variables.
627         If i=0 the algorithm computes a coordinate change such that J is in
628         noether position.
629EXAMPLE: example jungresolve; shows an example
630"
631{
632  def A = basering;
633  list result;
634  int i,j;
635  list nor = jungnormal(id,noeth);
636  for(i = 1;i<=size(nor);i++){
637    if(defined(R)){kill R;}
638    def R3 = nor[i];
639    setring R3;
640    def R = changeord("dp");
641    setring R;
642    ideal QIdeal = imap(R3,QIdeal);
643    ideal BMap = imap(R3,BMap);
644    map f = A,BMap;
645    if(QIdeal <> 0){
646      list res = resolve(QIdeal);
647      for(j =1;j<=size(res[1]);j++){
648        def R2 = res[1][j];
649        setring R2;
650        if(defined(QIdeal)){kill QIdeal;}
651        if(defined(BMap)){kill BMap;}
652        if(BO[1]<>0){ideal QIdeal = BO[1]+BO[2];}
653        else{ideal QIdeal = BO[2];}
654        map g = R,BO[5];
655        ideal BMap = ideal(g(f));
656        export(QIdeal);
657        export(BMap);
658        result[size(result)+1] = R2;
659        kill R2;
660      }
661      kill res;
662    }
663    else{
664      result[size(result)+1] = nor[i];
665    }
666    kill R,R3;
667  }
668  return(result);
669}
670example
671{"EXAMPLE:";
672   echo = 2;
673   ring R=0,(x,y,z),dp;
674   ideal J=x2+y3z3+y2z5;
675   list li=jungresolve(J,1);
676   li;
677   def S=li[1];
678   setring S;
679   QIdeal;
680   BMap;
681}
682///////////////////////////////////////////////////////////////////////////////
683static proc NoetherP_test(ideal id){
684  def A = basering;
685  list ringA=ringlist(A);
686  int i,j,k;
687  int index = 0;
688  if(size(id)==1 && nvars(A)){  // test if V(id) = C[x,y,z]/<f>
689    list L;
690    intvec v = 1,1,1;
691    L[1] = "lp";
692    L[2] = v;
693    kill v;
694    poly f = id[1];
695    for(i = 1;i<=3;i++){
696      setring A;
697      list l = ringA;  //change ordering to lp and var(i)>var(j) j!=i
698      list vari = ringA[2];
699      string h = vari[1];
700      vari[1] = vari[i];
701      vari[i] = h;
702      l[2] = vari;
703      kill h,vari;
704      l[3][1] = L;
705      def R = ring(l);
706      kill l;
707      setring R;
708      ideal I = imap(A,id);
709      if(defined(v)){kill v;}
710      intvec v = leadexp(I[1]);
711      attrib(v,"isMonic",1);
712      for(k = 2;k<=3;k++){
713// checks whether f is monic in var(i)
714        if(v[k] <> 0 || v[1] == 0){
715          attrib(v,"isMonic",0);
716          j++;
717          break;
718        }
719      }
720      if(attrib(v,"isMonic")==1){
721        index = i;
722        return(index);
723      }
724      kill R;
725    }
726    if(j == 3){
727      return(0);
728    }
729  }
730  else{
731// not yet a test for more variables
732    return(index);
733  }
734}
735/////////////////////////////////////////////////////////////////////////
736//// copied from resolve.lib, deleting parts of procedures which are  ///
737//// not necessary in this setting                                    ///
738/////////////////////////////////////////////////////////////////////////
739static proc normalCrossing(ideal J,list E,ideal V)
740"Internal procedure - no help and no example available
741"
742{
743   int i,d,j;
744   int n=nvars(basering);
745   list E1,E2;
746   ideal K,M,Estd;
747   intvec v,w;
748
749   for(i=1;i<=size(E);i++)
750   {
751      Estd=std(E[i]+J);
752      if(deg(Estd[1])>0)
753      {
754         E1[size(E1)+1]=Estd;
755      }
756   }
757   E=E1;
758   for(i=1;i<=size(E);i++)
759   {
760      v=i;
761      E1[i]=list(E[i],v);
762   }
763   list ll;
764   int re=1;
765
766   while((size(E1)>0)&&(re==1))
767   {
768      K=E1[1][1];
769      v=E1[1][2];
770      attrib(K,"isSB",1);
771      E1=delete(E1,1);
772      d=n-dim(K);
773      M=minor(jacob(K),d)+K;
774      if(deg(std(M+V)[1])>0)
775      {
776         re=0;
777         break;
778      }
779      for(i=1;i<=size(E);i++)
780      {
781         for(j=1;j<=size(v);j++){if(v[j]==i){break;}}
782         if(j<=size(v)){if(v[j]==i){i++;continue;}}
783         Estd=std(K+E[i]);
784         w=v;
785         if(deg(Estd[1])==0){i++;continue;}
786         if(d==n-dim(Estd))
787         {
788            if(deg(std(Estd+V)[1])>0)
789            {
790               re=0;
791               break;
792            }
793         }
794         w[size(w)+1]=i;
795         E2[size(E2)+1]=list(Estd,w);
796      }
797      if(size(E2)>0)
798      {
799         if(size(E1)>0)
800         {
801            E1[size(E1)+1..size(E1)+size(E2)]=E2[1..size(E2)];
802         }
803         else
804         {
805            E1=E2;
806         }
807      }
808      kill E2;
809      list E2;
810   }
811   return(re);
812}
813//////////////////////////////////////////////////////////////////////////////
814
815static proc zariski(ideal id){
816  def A = basering;
817  ideal QIdeal = std(id);
818  ideal BMap= maxideal(1);
819  export(BMap);
820  int i,j;
821  export(QIdeal);
822  if(dim(QIdeal)<>2){
823     print("wrong dimension");
824     list result;
825     return(result);
826  }
827  list result;
828  result[1]= A;
829  attrib(result[1],"isSmooth",0);
830  while(1){
831    if(defined(count)){kill count;}
832    int count =0;
833    for(i = 1;i<= size(result);i++){
834      attrib(result[i],"isSmooth",1);
835      def R = result[i];
836      setring R;
837      if(!defined(Slocus)){
838        if(QIdeal[1] == 0){
839          ideal Slocus = 1;
840        }
841        else{
842          ideal Slocus = std(slocus(QIdeal));
843        }
844      }
845      if(NF(1,Slocus)<>0){
846        attrib(result[i],"isSmooth",0);
847      }
848      else{
849        count++;
850      }
851      kill R;
852    }
853    if(count == size(result)){return(result);}
854    count = 0;
855    list hlp;
856    for(i = 1;i<=size(result);i++){
857      if(attrib(result[i],"isSmooth")==1){
858        hlp[size(hlp)+1] = result[i];
859        attrib(hlp[size(hlp)],"isSmooth",attrib(result[i],"isSmooth"));
860        count++;
861        i++;
862        continue;
863      }
864      def R = result[i];
865      setring R;
866      //print(N1);
867      list nor = normal(QIdeal);
868      //print(N2);
869      for(j = 1;j<= size(nor[1]);j++){
870        def R3 = nor[1][j];
871        setring R3;
872        def R2 = changeord("dp");
873        setring R2;
874        ideal norid = imap(R3,norid);
875        ideal normap = imap(R3,normap);
876        kill R3;
877        map f = R,normap;
878        if(defined(BMap)){kill BMap;}
879        if(defined(QIdeal)){kill QIdeal;}
880        ideal BMap = f(BMap);
881        ideal QIdeal = norid;
882        export(BMap);
883        export(QIdeal);
884        if(QIdeal[1]<> 0){
885          ideal Slocus = slocus(QIdeal);
886          Slocus = std(radical(Slocus));
887        }
888        else{
889          ideal Slocus = 1;
890        }
891        if(NF(1,Slocus)<> 0){
892          list blowup = blowUp(norid,Slocus);
893          for(int k = 1;k<=size(blowup);k++){
894            def R3 = blowup[k];
895            setring R3;
896            ideal QIdeal = sT + aS;
897            map f = R2,bM;
898            ideal BMap = f(BMap);
899            export(BMap);
900            export(QIdeal);
901            hlp[size(hlp)+1] = R3;
902            attrib(hlp[size(hlp)],"isSmooth",0);
903            kill R3;
904          }
905          kill k,blowup;
906        }
907        else{
908          hlp[size(hlp)+1] = R2;
909          attrib(hlp[size(hlp)],"isSmooth",1);
910        }
911        kill R2;
912      }
913      kill R,j,nor;
914    }
915    kill i;
916    if(count == size(result)){
917      return(result);
918    }
919    else{
920      result = hlp;
921      attrib(result,"isSmooth",0);
922      kill hlp;
923    }
924  }
925}
926//////////////////////////////////////////////////////////////////////////////
927// End of copied part and edited part from resolve.lib
928//////////////////////////////////////////////////////////////////////////////
929
Note: See TracBrowser for help on using the repository browser.