source: git/Singular/LIB/resjung.lib @ 1b2216

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