source: git/Singular/LIB/resjung.lib @ 5895e8

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