source: git/Singular/LIB/classifyci.lib @ 5b6807

spielwiese
Last change on this file since 5b6807 was 5b6807, checked in by Hans Schoenemann <hannes@…>, 10 years ago
add: classifyci.lib
  • Property mode set to 100644
File size: 25.8 KB
Line 
1////////////////////////////////////////////////////////////////////////////////
2version=="version classifyci.lib 4.0.0.0 Jun_2013 ";
3category="Singularities";
4
5info="
6LIBRARY:  classifyci.lib     Isolated complete intersection singularities in characteristic  0
7AUTHORS:  Gerhard Pfister     pfister@mathematik.uni-kl.de
8          Deeba Afzal         deebafzal@gmail.com
9
10
11OVERVIEW:
12   A library for classifying isolated complete intersection singularities for the base field of characteristic  0
13   and for computing weierstrass semigroup of the space curve.Isolated complete intersection singularities were
14   classified by M.Giusti [1] for the base field of characteristic 0. Algorithm for the semigroup of a space
15   curve singularity is given in [2].
16
17REFERENCES:
18[1] Giusti,M:Classification des singularities isolees simples d'intersections completes,
19             C,R.Acad.Sci.Paris Ser.A-B 284(1977),167-169.
20[2] Castellanos,A.,Castellanos,J.,2005:Algorithm for the semigroup of a space curve singularity.
21             Semigroup Forum 70,44-66.
22PROCEDURES:
23 classifyicis(I);  Isolated simple complete intersection singularities for the base field of charateristic  0
24 Semigroup(I);     Weierstrass semigroup of the space curve given by equations
25";
26LIB "classify.lib";
27LIB "classify_aeq.lib";
28LIB "poly.lib";
29LIB "curvepar.lib";
30LIB "algebra.lib";
31
32////////////////////////////////////////////////////////////////////////////////
33proc classifyicis(ideal I)
34"USAGE: classifyicis(I);  I ideal
35ASSUME: I is given by two generators
36PURPOSE:Check whether the ideal defines a complete intersection singularity or not
37RETURN: String  type in the classification of Giusti,M
38@           or The given singularity is not simple
39EXAMPLE: example classifyicis;  shows an example
40"
41{
42   def R=basering;
43   def SS=changeord(list(list("ds",nvars(R))),R);
44   setring SS;
45   ideal I=imap(R,I);
46   string re;
47   if(char(basering)==0)
48   {
49       re=ICIS12(I);
50   }
51    if(char(basering)!=0)
52   {
53          re="The characteristic of basering should be 0";
54   }
55   setring R;
56   return(re);
57}
58example
59{
60  "EXAMPLE:"; echo=2;
61   ring R=0,(x,y,z),ds;
62   ideal I=x2+yz,xy+z4;
63   classifyicis(I);
64}
65////////////////////////////////////////////////////////////////////////////////
66static proc ICIS12(ideal I)
67{
68  int n=nvars(basering);
69  if(n==2)
70  {
71     return(zerodim_ICIS(I));
72  }
73  if(n==3)
74  {
75     return(onedim_ICIS(I));
76  }
77  if(n>=4)
78  {
79     return("The given singularity is not simple");
80  }
81}
82////////////////////////////////////////////////////////////////////////////////
83static proc zerodim_ICIS(ideal I)
84"USAGE: zerodim_ICIS(l);  I is an ideal
85ASSUME: I is given by two generators
86PURPOSE: Check whether the ideal defines a complete intersection singularity of dimension zero or not
87RETURN: String  type in the classification of Giusti,M, of the 0-dimensional complete inetersection
88@*       or The given singularity is not simple
89EXAMPLE: example zerodim_ICIS; shows an example
90"
91{
92  def R=basering;
93  poly g,h,r;
94  ideal J;
95  int a,b,c,d;
96  map phi;
97  list L;
98  string re;
99  // g=g[0]+g[1]   where ord(g[1])>=3   ,g[0] can be zero
100  // h=h[0]+h[1]    where ord(h[1])>=3    ,h[0] can be zero
101  d=vdim(std(I));
102  if(d==-1)
103  {
104    return("The given singularity is not simple");
105  }
106  a=ord(I[1]);
107  b=ord(I[2]);
108  if((a>=3)&&(b>=3))
109  {                                             //start case1
110    return("The given singularity is not simple");
111  }                                           // end case 1
112  if((a==2&&b>=3)||(a>=3&&b==2))                // start case 2
113  {
114    if(a==2)
115    {
116      g=I[1];
117      h=I[2];
118    }
119    if(b==2)
120    {
121      g=I[2];
122      h=I[1];
123    }
124    L=factorize(jet(g,2));
125    if(size(L[1])==3)
126    {
127      re=findwhichF(g,h,L);
128      return(re);
129    }                                     // end size(L[1]=3)
130    if((size(L[1])==2)&&(L[2][2]==2))
131    {                                     // start (size(L[1])==2)&&(L[2][2]==2)
132      // case (x2,h);
133      r=L[1][2];
134      if(size(r)==2)                //  ax+by  goes to x
135      {
136        matrix M3=coef(r,var(1));
137        M3=subst(M3,var(2),1);
138        matrix A[2][2]=M3[2,1],M3[2,2],0,1;
139        matrix B=inverse(A);
140        phi=R,B[1,1]*var(1)+B[1,2]*var(2),B[2,1]*var(1)+B[2,2]*var(2);
141        g=phi(g);
142        h=phi(h);
143        J=(g,h);
144      }            // end size(r)=2
145      if(size(r)==1)                 // jet(g,2)=ax2 or  by2   goes to x2
146      {
147        if(leadmonom(r)==var(1))
148        {
149          phi=R,var(1)/leadcoef(r),var(2);
150          g=phi(g);
151          h=phi(h);
152        }
153        if(leadmonom(r)==var(2))
154        {
155          phi=R,var(2)/leadcoef(r),var(1);
156          g=phi(g);
157          h=phi(h);
158        }
159        J=(g,h);
160      }                // end size(r)=1
161      c=milnor(g);
162      if((d>=7)&&(c==2))
163      {
164        //"I-series";
165        if(d mod(2)==0)
166        {
167          return("I_"+string(d-1)+":(x2+y3,y"+string(d div 2)+")");
168        }
169        if(d mod(2)!=0)
170        {
171          return("I_"+string(d-1)+":(x2+y3,xy"+string((d-3) div 2)+")");
172        }
173      }
174      if(d==6)
175      {
176        return("G_5:(x2,y3)");
177      }
178      ring R1=0,(var(2),var(1)),ds;
179      setring R1;
180      ideal J=imap(R,J);
181      poly h1=reduce(J[2],std(J[1]),d);
182      poly h2=leadmonom(h1);
183      int ss=deg(h2)-1;
184      if((h2==var(1)^4)&&((c>=3)||(c==-1)))
185      {
186        setring R;
187        return("G_7:(x2,y4)");
188      }
189      if((h2==var(2)*var(1)^2)&&((c>=3)||(c==-1)))
190      {
191        setring R;
192        return("H_"+string(d-1)+":(x2+y"+string(d-4)+",xy2)");
193      }
194      setring R;
195      return("The given singularity is not simple");
196    }                                          // end (size(L[1])==2)&&(L[2][2]==2)
197    if((size(L[1])==2)&&(L[2][2]==1))
198    {
199      def S=factorExt(g);
200      setring S;
201      poly h=imap(R,h);                     // poly g=imap(R,g);  we need not S has already g
202      re= findwhichF(g,h,L);
203      setring R;
204      return(re);
205    }
206  }                                            // end case 2
207  if((a==2)&&(b==2))                  // start case 3
208  {
209    g=I[1];
210    h=I[2];
211    poly Q=testDiv(jet(h,2),jet(g,2));
212    if(Q!=0)
213    {
214      I=(g,h-Q*g);
215      return(zerodim_ICIS(I));
216    }
217    if(Q==0)
218    {
219      L=factorize(jet(g,2));
220      if(size(L[1])==3)
221      {
222        re=findwhichF(g,h,L);
223        return(re);
224      }
225      if((size(L[1])==2)&&(L[2][2]==1))
226      {
227        def S=factorExt(g);
228        setring S;
229        poly h=imap(R,h);
230        re=findwhichF(g,h,L);
231        setring R;
232        return(re);
233      }
234      L=factorize(jet(h,2));
235      if(size(L[1])==3)
236      {
237        re=findwhichF(h,g,L);
238        return(re);
239      }
240      if((size(L[1])==2)&&(L[2][2]==1))
241      {
242        def S=factorExt(h);
243        setring S;
244        poly g=imap(R,g);
245        re=findwhichF(h,g,L);
246        setring R;
247        return(re);
248      }
249      else
250      {    // there exist a s.t g[0]+ah[0] has two different factors.
251        int e=Finda(g,h);
252        I=(g+e*h,h);
253        re=zerodim_ICIS(I);
254        return(re);
255      }
256    }
257  }                                              // end case 3
258}                                                //  proc
259example
260{
261  "EXAMPLE:"; echo=2;
262   ring R=0,(x,y),ds;
263   ideal I=x2+8xy+16y2+y3,xy7+4y8;
264   zerodim_ICIS(I);
265}
266//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
267static proc Finda(poly g,poly h)
268"USAGE:  Finda(g,h);  g,h are polynomials
269PURPOSE: Find a such that jet(g,2)+a*jet(h,2) has two different factors
270RETURN: integer a
271{
272  // find  a s.t jet(g,2)+a*jet(h,2) has two different factors.
273  int o;
274  list L=factorize(jet(h,2));
275  if(L[2][2]==1){return(0);}
276  poly r= jet(g,2);
277  list T=factorize(r);
278  while(T[2][2]!=1)
279  {
280    o++;
281    r= r+jet(h,2);
282    T=factorize(r);
283  }
284  return(o);
285}
286/*
287ring R=0,(x,y),ds;
288poly g=x2+xy3;
289poly h=y2+x3+y7;
290finda(g,h);
291*/
292//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
293static proc testDiv(poly f,poly g)
294"USAGE:  testDiv(f,g);  f,g are polynomials
295ASSUME: I is given by two generators
296PURPOSE: Check whether  f divides g or not.
297RETURN: poly h(quotient)  if f divides g
298 @*              0 if f does not divide g
299{
300  poly h=f/g;
301  if(f-h*g==0)
302  {
303    return(h);
304  }
305  return(0);
306}
307/*
308   ring R=0,(x,y,ds;
309   poly f=x2+y2;
310   poly g=3x2+3y2
311   testDiv(f,g);
312*/
313//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
314 static proc findwhichF(poly g,poly h,list L)
315"USAGE:  findwhichF(g,h,L);  g,h are polynomials,L list of factors of jet(g,2)
316RETURN: string  type F^n,p_n+p-1 in the classification of Giusti,M
317{
318  // return("F^n,p_n+p-1")
319  def R=basering;
320  matrix M,N;
321  ideal J;
322  map si;
323  list T;
324  string rem;
325  //  in each case we want to transform jet(g,2) which has two factors to xy
326// and then find std and know about the type of  "F^n,p_n+p-1"
327  if((size(L[1][2])==2)&&(size(L[1][3])==2))
328  {
329    M=coef(L[1][2],var(1));
330    N=coef(L[1][3],var(1));
331    matrix A[2][2]=M[2,1],M[2,2],N[2,1],N[2,2];
332    A=subst(A,var(1),1,var(2),1);
333    matrix B=inverse(A);
334    si=R,B[1,1]*var(1)+B[1,2]*var(2),B[2,1]*var(1)+B[2,2]*var(2);
335    g=si(g);
336    h=si(h);
337    J=(g,h);
338    J=std(J);
339    // if g and h both of order 2 no problem in that case because deg(J[2])=2 so lead(J[2])=x2,y2,xy does not matter
340    T[1]=deg(lead(J[2]));
341    T[2]=deg(lead(J[3]))-1;
342    rem="F^"+string(T[1])+","+string(T[2])+"_"+string(T[1]+T[2]-1)+":(xy,x"+string(T[1])+"+y"+string(T[2])+")";
343    return(rem);
344  }
345  if((size(L[1][2])==2)&&(size(L[1][3])==1))
346  {
347    // two cases 1- jet(g,2)=(ax+by)*cx,  2-  jet(g,2)=(ax+by)*cy
348    if(leadmonom(L[1][3])==var(1))
349    {
350      M=coef(L[1][2],var(1));
351      M=subst(M,var(2),1);
352      si=R,var(1),-M[2,1]/M[2,2]*var(1)+var(2);
353      g=si(g);
354      h=si(h);
355      J=(g,h);
356      J=std(J);
357      T[1]=deg(lead(J[2]));
358      T[2]=deg(lead(J[3]))-1;
359      rem="F^"+string(T[1])+","+string(T[2])+"_"+string(T[1]+T[2]-1)+":(xy,x"+string(T[1])+"+y"+string(T[2])+")";
360      return(rem);
361    }
362    if(leadmonom(L[1][3])==var(2))
363    {
364      M=coef(L[1][2],var(1));
365      M=subst(M,var(2),1);
366      matrix A[2][2]=M[2,1],M[2,2],0,1;
367      matrix  B=inverse(A);
368      si=R,B[1,1]*var(1)+B[1,2]*var(2),B[2,1]*var(1)+B[2,2]*var(2);
369      g=si(g);
370      h=si(h);
371      J=(g,h);
372      J=std(J);
373      T[1]=deg(lead(J[2]));
374      T[2]=deg(lead(J[3]))-1;
375      rem="F^"+string(T[1])+","+string(T[2])+"_"+string(T[1]+T[2]-1)+":(xy,x"+string(T[1])+"+y"+string(T[2])+")";
376      return(rem);
377    }
378  }
379  else
380  {
381    J=(g,h);
382    J=std(J);
383    T[1]=deg(lead(J[2]));
384    T[2]=deg(lead(J[3]))-1;
385    rem="F^"+string(T[1])+","+string(T[2])+"_"+string(T[1]+T[2]-1)+":(xy,x"+string(T[1])+"+y"+string(T[2])+")";
386    return(rem);
387  }
388}
389/*
390 ring R=0,(x,y),ds;
391poly g=xy+y2+y4;
392poly h=x4+4x3y+6x2y2+4xy3+y4+y7+xy7;
393L=factorize(jet(g,2),2);
394findwhichF(g,h,L);
395*/
396////////////////////////////////////////////////////////////////////////////////
397static proc factorExt(poly g)
398"USAGE:  procExt(g); jet(g,2) is an irreducible polynomial
399PURPOSE: Find the field extension in which jet(g,2) has two different factors
400RETURN: ring S in which jet(g,2) has factors
401{
402  def R=basering;
403  g=simplify(g,1);
404  poly f=jet(g,2);
405  list L=factorize(f);
406  if(L[2][2]==1)
407  {
408    ring S=(0,t),(var(1),var(2)),ds;
409    poly f=fetch(R,f);
410    poly g=fetch(R,g);
411    minpoly=t2+leadcoef(f/(var(1)*var(2)))*t+leadcoef(f/var(2)^2);
412    list L=factorize(f);
413    export L;
414    export g;
415  }
416  else
417  {
418    def S=R;
419    export L;
420  }
421  setring R;
422  return(S);
423}
424/*
425ring R=0,(x,y),ds;
426poly g=x2=y2+x4+xy11;
427factorExt(g);
428*/
429////////////////////////////////////////////////////////////////////////////////
430static proc onedim_ICIS(ideal I)
431"USAGE:  onedim_ICIS(l);  I is an ideal
432ASSUME: I is given by two generators
433PURPOSE: Check whether the ideal defines a complete intersection singularity of dimension 1 or not
434RETURN: String  type in the classification of Giusti,M, of 1-dimesnional complete inetersection singualrity
435@*       or The given singularity is not simple
436EXAMPLE: example onedim_ICIS; shows an example
437"
438{
439  int m,t,r;
440  poly g1,g2,f1,f2;
441  string rem;
442  list A,B;
443  f1=I[1];
444  f2=I[2];
445  // I=nf_icis(I);
446  m=genericmilnor(I);
447  t=tjurina(I);
448  if(m==-1)
449  {
450    return("The given singularity is not simple");
451  }
452  if(m!=t)                // in ICIS milnor=tjurina
453  {
454    return("The given singularity is not simple");
455  }
456  g1=jet(f1,2);
457  g2=jet(f2,2);
458  if((ord(g1)==1)||(ord(g2)==1)){return(arnoldsimple(I,m));}
459  if(g1==0)
460  {
461    return("The given singularity is not simple");
462  }
463  if(g2==0)
464  {
465    return("The given singularity is not simple");
466  }
467  rem=typejet2(g1,g2);
468  if(rem=="type1")
469  {
470    return("S_5:(x2+y2+z2,yz)");
471  }
472  if(rem=="type2")
473  {
474    return("S_"+string(m)+":(x2+y2+z"+string(m-3)+",yz)");
475  }
476  if(rem=="type3")
477  {
478    if(m==7)
479    {
480      return("T_7:(x2+y3+z3,yz)");
481    }
482    if(m==8)
483    {
484      return("T_8:(x2+y3+z4,yz)");
485    }
486    if(m==9)
487    {
488      B=Semigroup(I);
489      B=changeType(B);
490      A=list(list(2,5),list(2,3));
491      if(compLL(A,B))
492      {
493        return("T_9:(x2+y3+z5,yz)");
494      }
495    }
496    return("The given singularity is not simple");
497  }
498  if(rem=="type4")
499  {
500    if(m==7)
501    {
502      return("U_7:(x2+yz,xy+z3)");
503    }
504    if(m==8)
505    {
506      return("U_8:(x2+yz+z3,xy)");
507    }
508    if(m==9)
509    {
510      return("U_9:(x2+yz,xy+z4)");
511    }
512    return("The given singularity is not simple");
513  }
514  if(rem=="type5")
515  {
516    if(m==8)
517    {
518      return("W_8:(x2+z3,y2+xz)");
519    }
520    if(m==9)
521    {
522      return("W_9:(x2+yz2,y2+xz)");
523    }
524    return("The given singularity is not simple");
525  }
526  if(rem=="type6")
527  {
528    if(m==9)
529    {
530      return("Z_9:(x2+z3,y2+z3)");
531    }
532    if(m==10)
533    {
534      return("Z_10:(x2+yz2,y2+z3)");
535    }
536    return("The given singularity is not simple");
537  }
538  if(rem=="not simple")
539  {
540    return("The given singularity is not simple");
541  }
542}
543example
544{
545  "EXAMPLE:"; echo=2;
546   ring R=0,(x,y,z),ds;
547   ideal I=x2+8xy+16y2+2xz+8yz+z2+yz2+9z3,y2+xz+22yz+82z2;
548   onedim_ICIS(I);
549}
550////////////////////////////////////////////////////////////////////////////////
551static proc arnoldsimple(ideal I,int m)
552"USAGE:  arnoldsimple(I,m);  I is an ideal, m is an integer greater or equal to milnor number of the ideal I
553ASSUME: I is given by two generators and one of the generator of the ideal I is of order 1
554PURPOSE: check whether the ideal defines a  hypersurface simple complete intersection singularity or not
555RETURN: string  type in the classification of Arnold,
556@*       or The given singularity is not simple
557EXAMPLE: example arnoldsimple; shows an example
558"
559{
560 //if one generator is of order 1 we reduce case to hypersurface case
561  def R=basering;
562  if(ord(I[1])==1)
563  {
564    poly g=specialNF(I[2],I[1],m);
565    ring S=0,(var(2),var(3)),ds;
566    setring S;
567    poly g=imap(R,g);
568    string dd=complexSingType(g);
569    int e=modality(g);
570    setring R;
571    if(e==0)
572    {
573      return(dd);
574    }
575    if(e!=0)
576    {
577      return( "The given singularity is not simple");
578    }
579  }
580  if(ord(I[2])==1)
581  {
582    I=I[2],I[1];
583    return(arnoldsimple(I,m));
584  }
585}
586example
587{
588  "EXAMPLE:"; echo=2;
589   ring R=0,(x,y,z),ds;
590   ideal I=x+y2,y3+z4+xy11;
591   arnoldsimple(I,6);
592}
593////////////////////////////////////////////////////////////////////////////////
594static proc specialNF(poly g,poly f,int m)
595"USAGE:  specialNF(g,f);  g,f are polynomials, m is an integer greater or equal to milnor number of the ideal I=(g,h)
596ASSUME: f is of order 1 and f=x+higher
597RETURN: poly g not involving x (Using implicit fn thm)
598{
599  poly k;
600  list T=linearpart(g,f);
601  f=T[1];
602  g=T[2];
603  poly h=var(1)-f;
604  while(1)
605  {
606    g=subst(g,var(1),h);
607    k=jet(g,m);
608    if(diff(k,var(1))==0)
609    {
610      break;
611    }
612  }
613  return(k);
614}
615/*
616ring R=0,(x,y,z),ds;
617poly f=x+y2;
618poly g=y3+z4+xy11;
619speciaNF(g,f,6);
620*/
621////////////////////////////////////////////////////////////////////////////////
622static proc linearpart(poly g,poly f)          //   f=linear part+higher,g   output list T,T[1]=x+higher term,T[2]=g
623"USAGE:  specialNF(g,f);  g,f are polynomials
624ASSUME: f=lineat part+higher that is f is of order 1
625RETURN: list T, T[1]=x+higher and T[2]=g'
626{
627  def R=basering;
628  poly i,j,k;
629  list T;
630  i=diff(jet(f,1),var(1));
631  j=diff(jet(f,1),var(2));
632  k=diff(jet(f,1),var(3));
633  if(i!=0)
634  {
635    ideal M=maxideal(1);
636    M[1]=(var(1)-((j*var(2)+k*var(3))/leadcoef(f)))/leadcoef(f);
637    map phi=R,M;
638    f=phi(f);
639    g=phi(g);
640  }
641  if(i==0)
642  {
643    if(j!=0)
644    {
645      map phi=R,var(2),var(1),var(3);
646      f=phi(f);
647      g=phi(g);
648      return(linearpart(g,f));
649    }
650    if(k!=0)
651    {
652      map phi=R,var(3),var(2),var(1);
653      f=phi(f);
654      g=phi(g);
655      return(linearpart(g,f));
656    }
657  }
658  T[1]=f;
659  T[2]=g;
660  return(T);
661}
662/*
663ring R=0,(x,y,z),ds;
664poly f=x+2y+y2;
665poly g=y3+z4+zy11;
666lineatpart(g,f);
667*/
668////////////////////////////////////////////////////////////////////////////////
669static  proc typejet2(poly g1,poly g2)
670"USAGE:  typejet2(g1,g2);  g1,g2 are polynomials
671ASSUME: g1,g2 are homogenous polynomials of degree 2
672PURPOSE: Check whether (g1,g2) is a quadratic form in the list of Guisti or not
673RETURN: string  type for the quadratic forms appearing in Guist's list
674@*       or not simple
675{
676  def R=basering;
677  ideal I=(g1,g2);
678  def S=absPrimdecGTZ(I);
679  setring S;
680  list L,T;
681  int e,i,j;
682  intvec a,a1;
683  L=primary_decomp;
684  for(j=1;j<=size(L);j++)
685  {
686    if(dim(std(L[j][1]))!=2)
687    {
688      return("not simple");
689    }
690  }
691  T=absolute_primes;
692  for(i=1;i<=size(T);i++)
693  {
694    e=e+T[i][2];
695  }
696  if(e==4)
697  {
698    setring R;
699    return("type1");
700  }
701  if(e==3)
702  {
703    setring R;
704    return("type2");
705  }
706  if(e==2)
707  {
708    ideal J=std(L[1][1]);
709    ideal J1=std(L[2][1]);
710    a=hilbPoly(J);
711    a1=hilbPoly(J1);
712    if((a[2]==2)&&(a1[2]==2))
713    {
714      setring R;
715      return("type3");
716    }
717    if(((a[2]==3)&&(a1[2]==1))||((a[2]==1)&&(a1[2]==3)))                      //||(a[2]==1)&&(a1[2]==3))
718    {
719      setring R;
720      return("type4");
721    }
722  }
723  if(e==1)
724  {
725    setring R;             // I lies in R and zero in S1
726    ideal JJ=radical(I);
727    JJ=JJ^3;
728    ideal JJJ=reduce(JJ,std(I));
729    if(size(JJJ)==0)
730    {
731      return("type6");
732    }
733    if(size(JJJ)!=0)
734    {
735      return("type5");
736    }
737  }
738  setring R;
739  return("not simple");
740}
741/*
742ring R=0,(x,y,z),ds;
743poly g1=x2+yz;
744poly g2=xy;
745typejet2(g1,g2);
746*/
747////////////////////////////////////////////////////////////////////////////////
748static proc compL(list L,list M)
749{
750  int l,m,i,j;
751  l=size(L);
752  m=size(M);
753  if(l!=m)
754  {return(0);}
755  for(i=1;i<=m;i++)
756  {
757    if(L[i]!=M[i])
758    {return(0);}
759  }
760  return(1);
761}
762////////////////////////////////////////////////////////////////////////////////
763static proc compLL(list L,list M)
764"USAGE:  compLL(L,M);  L, M are lists
765PURPOSE: Check whether the lists are equal or not
766RETURN: 1 if both lists are equal upto a permutation
767 @*              0 if both are not equal
768{
769  int l,m,i,j,s;
770  l=size(L);
771  m=size(M);
772  if(l!=m)
773  {return(0);}
774  for(i=1;i<=m;i++)
775  {
776    for(j=1;j<=m;j++)
777    {
778      if(compL(L[i],M[j]))
779      {
780        s++;
781        break;
782      }
783    }
784  }
785  if(s==m)
786  {return(1);}
787  if(s!=m)
788  {
789    return(0);
790  }
791}
792/*
793ring R=0,(x,y,z),ds;
794list L=(1),(2,3),(2,5);
795list T=(2,3),(1),(2,5);
796compLL(L,T);
797*/
798////////////////////////////////////////////////////////////////////////////////
799static proc changeType(list L)
800"USAGE:  changeType(L);  L is a list of intvectors
801PURPOSE: Change the list of intvectors to the list of lists
802RETURN: List of lists
803{
804  int i,j;
805  list T;
806  for(i=1;i<=size(L);i++)
807  {
808    list S;
809    for(j=1;j<=size(L[i]);j++)
810    {
811      S[j]=L[i][j];
812    }
813    T[size(T)+1]=S;
814    kill S;
815  }
816  return(T);
817}
818/*
819ring R=0,(x,y,z),ds;
820list B=(4,6,7);
821changeType(B);
822*/
823////////////////////////////////////////////////////////////////////////////////
824static proc genericmilnor(ideal I)
825"USAGE:  genericmilnor(l);  I is an ideal
826PURPOSE: Computes the milnor number of generic linear combination of the ideal I
827RETURN: Milnor number of I if it is finite
828@*       or -1 if it is not finite
829{
830  int m=milnor(I);
831  int i,a,b;
832  if(m>=0)
833  {
834    return(m);
835  }
836  def R=basering;
837  def R1=changechar(32003,R);
838  setring R1;
839  ideal I;
840  while(i<10)
841  {
842    i++;
843    a=random(-100,100);
844    b=random(-100,100);
845    while(a==0)
846    {
847      a=random(-100,100);
848    }
849    while(b==0)
850    {
851      b=random(-100,100);
852    }
853    I=imap(R,I);
854    I[1]=a*I[1]+b*I[2];
855    m=milnor(I);
856    if(m>=0)
857    {
858      setring R;
859      return(m);
860    }
861  }
862  setring R;
863  return(-1);
864}
865/*
866   ring R=0,(x,y,z),ds;
867   ideal I=x2+z3,y2+z3;
868   genericmilnor(I);
869*/
870////////////////////////////////////////////////////////////////////////////////
871static proc Semigroup(ideal I)
872"USAGE:  Semigroup(l);  I is an ideal
873PURPOSE: Computes the semigroup of the ideal I corresponding to each branch
874RETURN: list of semigroup of ideal I corresponding to each branch
875EXAMPLE: Semigroup; shows an example
876"
877{
878  list L=facstd(I);
879  list RE,JE,PE;
880  if(size(L)==1)
881  {
882     RE=CurveRes(L[1]);
883     RE=semi_group(RE);
884     return(RE);
885  }
886  ideal J,K;
887  list T,T1,T2,T3,T4,T5,H;
888  int i,j,l;
889  for(i=1;i<=size(L);i++)
890  {
891    RE=CurveRes(radical(L[i])) ;
892    T1[i]=semi_group(RE);
893    for(j=i+1;j<=size(L);j++)
894    {
895      JE=CurveRes(radical(L[j]));
896      T2[j]=semi_group(JE);
897      J=L[i]+L[j];
898      if(dim(std(J))!=1)
899      {
900        break;
901      }
902      K=slocus(J);
903      if(K[1]==1)
904      { T3=1;}
905      else
906      {
907        PE=CurveRes(radical(J));
908        T3=semi_group(PE);
909      }
910      T4=commonpartlists(T1[i],T3);
911      T5=commonpartlists(T2[j],T3);
912      if(compLL(T4,T5))
913      {
914        T1[i]=del(T1[i],T4);
915      }
916    }
917    for(l=1;l<=size(T1[i]);l++)
918    {
919      H[size(H)+1]=T1[i][l];
920    }
921  }
922  return(H);
923}
924example
925{
926  "EXAMPLE:"; echo=2;
927   ring R=0,(x,y,z),ds;
928   ideal I=x2+y3+z5,yz;
929   classifyicis(I);
930}
931////////////////////////////////////////////////////////////////////////////////
932static proc del(list L,list M)
933"USAGE:  del(L,M);  L and M are two lists
934PURPOSE: Delete common part of list M from List L
935RETURN: list L
936{
937  int i,j;
938  for(i=1;i<=size(M);i++)
939  {
940    for(j=1;j<=size(L);j++)
941    {
942      if(compL(L[j],M[i]))
943      {L=delete(L,j);}
944    }
945  }
946  return(L);
947}
948////////////////////////////////////////////////////////////////////////////////
949static proc commonpartlists(list L,list M)
950"USAGE:  commonpart(L,M);  L and M are two lists
951PURPOSE: Computes the intersetion of two list
952RETURN: list T
953{
954   list T;
955   int i,m,l,j,k;
956   m=size(M);
957   l=size(L);
958   if(l>=m)
959   {
960     for(i=1;i<=m;i++)
961     {
962       for(j=1;j<=l;j++)
963       {
964         if(compLL(M[i],L[j]))
965         {
966           T[k+1]=M[i];
967           k++;
968         }
969       }
970     }
971   }
972   return(T);
973}
974////////////////////////////////////////////////////////////////////////////////
975static proc semi_group(list H)
976"USAGE:  semi_group(H);  H list
977COMPUTE:Weierstrass semigroup of space curve C,which is given by an ideal
978RETURN: list A  , which gives generators set of the Weierstrass semigroup corresponding to each irreducible component of C
979{
980   int i,d,u,v,w,k;
981   int j=1;
982   int e=1;
983   def R=basering;
984
985   list A;
986   string mpo;
987   list LL;
988   for(k=1;k<=size(H);k++)
989   {
990     LL=CurveParam(H[k]);
991     def S=LL[1];
992     setring S;
993     list TT;
994     for(i=1;i<=size(Param);i++)
995     {
996       d=deg(Param[i][2]);
997       TT=Param[i];
998       mpo=string(Param[i][2]);
999       ring S1=(0,a),(t),ds;
1000       setring S1;
1001       execute("minpoly="+mpo+";");
1002       list TT=imap(S,TT);
1003       list T;
1004       ideal J1;
1005       for(u=1;u<=size(TT[1]);u++)
1006       {
1007         J1[u]=TT[1][u];
1008       }
1009       J1=simplify(J1,2);
1010       J1=sagbiM(J1);
1011       w=ConductorBound(J1);
1012       J1=lead(J1);
1013       list TTT;
1014       for(v=1;v<=size(J1);v++)
1015       {
1016         TTT[v]=J1[v];
1017       }
1018       for(j=1;j<=d;j++)
1019       {
1020         T=WSemigroup(TTT,w);
1021         A[e]=T[1];       //  intersted only in semigroup
1022         e++;
1023       }
1024       setring S;
1025       kill S1;
1026       kill T;
1027     }
1028     setring R;
1029     kill S;
1030   }
1031   return(A);
1032}
1033//==============================Examples======================================
1034/*
1035//=========Examples of Isolated simple complete intersection singularities======
1036ring R=0,(x,y),ds;
1037ideal M=maxideal(1);
1038//======================
1039ideal I=x2+y3,xy11;
1040M[1]=x;
1041M[2]=x+3y+xy;
1042map phi=R,M;
1043I=phi(I);
1044classifyicis(I);
1045//======================
1046ideal I=xy,x5+y4;
1047M[1]=x+4y;
1048M[2]=y;
1049map phi=R,M;
1050I=phi(I);
1051classifyicis(I);
1052//======================
1053ideal I=x2,y4;
1054M[1]=x+xy2;
1055M[2]=x+y+x2+y2;
1056map phi=R,M;
1057I=phi(I);
1058classifyicis(I)
1059//===========================================
1060ideal I=x2+y11,x2y3+xy4;
1061classifyicis(I);
1062//======================
1063ring S=0,(u,v),dp;
1064ideal N=maxideal(1);
1065//======================
1066ideal J=u2+v7,uv2;
1067N[1]=u+3v+uv+u3v;
1068N[2]=v;
1069map si=S,N;
1070J=si(J);
1071classifyicis(J);
1072//======================
1073ideal J=u2+v2+uv5+v11,uv4+v5;
1074classifyicis(I);
1075//===========================================
1076ring R=0,(x,y,z),ds;
1077ideal M=maxideal(1);
1078//======================
1079ideal I=x2+y3+z5,yz;
1080classifyicis(I);
1081//======================
1082ideal I=x2+z3,y2+z3;
1083classificis(I);
1084//======================
1085ideal I=x2+yz+z3,xy;
1086M[3]=x+4y+3z+x2y;
1087map phi=R,M;
1088I=phi(I);
1089classifyicis(I);
1090//======================
1091ideal I=x2+y3+z6,yz+xy3;
1092classifyicis(I);
1093//============================================
1094ideal I=x2+z3,y2+xz;
1095M[2]=x+3y;
1096map phi=R,M;
1097I=phi(I);
1098classifyicis(I);
1099//============================================
1100ring S=0,(u,v,w),ds;
1101ideal M=maxideal(1);
1102ideal I=u2+vw+w3,uv;
1103M[1]=u+3v+3vw+w2;
1104 map phi=S,M;
1105 I=phi(I);
1106classifyicis(I);
1107//==========Examples of Semigroup of the space curves====================
1108ring R=0,(x,y,z),ds;
1109ideal I=xy+z3,xz+z2y2+y6;
1110Semigroup(I);
1111//======================
1112ideal I=xy,xz+z3+z2y3+y11;
1113Semigroup(I);
1114//======================
1115ideal I=xy+z4,xz+y6+yz2;
1116Semigroup(I);
1117//======================
1118ideal I=xy+z2,x2+z2y+5y4;
1119Semigroup(I);
1120//======================
1121ideal I=x2+yz2,y2+z3;
1122Semigroup(I);
1123//======================
1124ideal I=x2+yz,xy+z4;
1125Semigroup(I);
1126//======================
1127ideal I=xy,xz+z3+z2y5+2y15;
1128Semigroup(I);
1129//======================
1130ideal I=xy,xz+z3+zy9;
1131Semigroup(I);
1132//======================
1133*/
Note: See TracBrowser for help on using the repository browser.