source: git/Singular/LIB/classifyceq.lib @ 380a17b

spielwiese
Last change on this file since 380a17b was 380a17b, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: new version numbers for libs
  • Property mode set to 100644
File size: 61.1 KB
Line 
1/////////////////////////////////////////////////////////////////////////////////
2version="version classifyceq.lib 4.0.0.0 Jun_2013 ";
3category="Singularities";
4
5info="
6LIBRARY:  classifyCeq.lib       simple hypersurface singularities in characteristic p > 0
7AUTHORS:  Deeba Afzal           deebafzal@gmail.com
8          Faira Kanwal Janjua   fairakanwaljanjua@gmail.com
9
10OVERVIEW:
11   A library for classifying the simple singularities with respect to contact equivalence in charateristic p > 0 .
12   Simple hypersurface singularities in charateristic p > 0 were classified by Greuel and
13   Kroening [1] with respect to contact equivalence. The classifier we use has been proposed in [2].
14
15REFERENCES:
16[1] Greuel, G.-M.; Kroening, H.: Simple singularities in positive characteristic.
17    Math.Z. 203, 339-354 (1990).
18[2] Afzal,D.;Binyamin,M.A.;Janjua,F.K.: On the classification of simple singularities in positive characteristic.
19
20PROCEDURES:
21  classifyCeq(f);     simple hypersurface singularities in charateristic p > 0
22";
23LIB "sing.lib";
24LIB "classify.lib";
25LIB "primdec.lib";
26LIB "ring.lib";
27/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
28proc classifyCeq(poly f)
29"USAGE:  classifyCeq(f);  f poly
30RETURN: string including the Tjurina number of f and its type in the classification of Greuel and Kroening
31EXAMPLE: example classifyCeq;  shows an example
32"
33{
34   def R=basering;
35   def S=changeord(list(list("ds",1:nvars(basering))));
36   setring S;
37   poly f=imap(R,f);
38   string re;
39   if(char(basering)==0)
40   {
41     re=(string(classify(f)));
42   }
43   if(char(basering)!=2)
44   {
45     re=classifyCeq1(f);
46   }
47   if(char(basering)==2)
48   {
49     re=classifyCeq2(f);
50   }
51   setring R;
52   return(re);
53}
54example
55{
56  "EXAMPLE:"; echo=2;
57   ring R=3,(x,y,z,u,v,w),ds;
58   classifyCeq(-x2+xy+y2+xz-yz-z2+w2+u3+v4);
59}
60///////////////////////////////////////////////////////////  char p > 2 //////////////////////////////////////////
61static proc blowupone(poly f)
62{
63   //=== input f smooth or isolated simple singularity at zero
64   //=== output var(1) or poly with isolated singularity at zero
65   def R=basering;
66   def S=changeord(list(list("ds",1:nvars(basering))));
67   setring S;
68   int n=nvars(basering);
69   int i;
70   poly f=imap(R,f);
71   if(deg(lead(f))<=1){setring R;return(var(1));}
72   def T=changeord(list(list("lp",1:nvars(basering))));
73   setring T;
74   map phi;
75   ideal mphi, sing;
76   poly p,q;
77   //===========  blow up =======================================================
78   for(i=1;i<=n;i++)
79   {
80      mphi=var(i)*maxideal(1);
81      mphi[i]=var(i);
82      phi=S,mphi;
83       p=phi(f);
84      q=p/var(i);
85      while(size(p)==size(q))
86      {
87         p=q;
88         q=q/var(i);
89      }
90      //===============  p is the strict transform var(i) exceptional divisor ====
91       //=============== analysis of singularities ================================
92      sing=jacob(p),p,var(i);
93      sing=radical(sing);
94      option(redSB);
95      sing=std(sing);
96      sing=simplify(sing,1);
97       if(size(sing)>1)
98      {
99         if(size(sing)!=n){ERROR("not simple");}
100         ideal mpsi=std(maxideal(1));
101         for(i=1;i<=n;i++){mpsi[i]=var(i)-sing[n-i+1][2];}
102         map psi=T,mpsi;
103         p=psi(p);
104         setring R;
105         poly p=imap(T,p);
106         return(p);
107      }
108    }
109   setring R;
110   return(var(1));
111}
112////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
113 static proc classifyCeq1(poly f)
114//====Classify Simple hypersurface singularities when charateristic p!=2
115{
116      // char!=2
117      //====input poly f
118      //====output  The function defines ......not an isolated singularity or  not a simple singularity  and tjurina
119      //          of the singularity or a simple singularity ,tjurina of the singularity and type of the singularity.
120      def R=basering;
121      int d;
122      int m=tjurina(f);
123
124           if(m==-1)
125           {
126                 return("The given function defines not an isolated singularity" );
127
128           }
129           poly g;
130           int c;
131           c=corank(f);
132           if(c<=1)
133           {
134
135               if(c==0)
136               {
137                  return("The given function defines a simple singularity.
138The tjurina number is "+string(m)+".
139A[1]:"+string(var(1))+"^2+"+string(var(2))+"^2");
140               }
141
142               if(c==1)
143
144              {
145                  poly w=blowupone(f);
146                  int n=tjurina(w);
147
148                 int b=m-n;
149                 if(b<=2)
150                 {
151                   return( "The given function defines a simple singularity.
152The tjurina number is "+string(m)+".
153A["+string(m)+"]:"+string(var(1))+"^2+"+string(var(2))+"^"+string(m+1)+"");
154                 }
155                 else
156                 {
157
158                    return("The given function defines a simple singularity.
159The tjurina number is "+string(m)+".
160A["+string(m-1)+"]:"+string(var(1))+"^2+"+string(var(2))+"^"+string(m)+"");
161                 }
162
163   }
164}
165           if(c==2)
166           {
167                 g=jet(f,3);
168                 if(g==0)
169                 {
170                    return( "The given function defines not a simple singularity.
171The tjurina number is "+string(m)+".");
172                 }
173                 if(g!=0)
174                 {
175
176                       def S=GDsplitting(f);
177                       setring S;
178
179
180                      poly g=jet(f,3);
181                         if(g==0)
182                         {
183                             setring(R);
184                            return("The given function defines not a simple singularity.
185The tjurina number is "+string(m)+".");
186                          }
187                      list L=factorize(g);
188                      if(size(L[1])==2)
189                      {
190                              ideal M=var(1),var(2)^2;
191                              ideal N=std(M^3);
192                               poly h=reduce(f,N);
193                                if(h==0)
194                                {
195                                   return(" The given function defines not a simple singularity
196The tjurina number is "+string(m)+".");
197                                }
198                                if(h!=0)
199                                {
200                                       if((char(R)!=3)&&(char(R)!=5))
201                                        {
202                                             setring R;
203                                             if(m==6)
204                                             {
205                                                 return("The given function defines a simple singularity.
206The tjurina number is "+string(m)+".
207E^0[6]:"+string(var(1))+"^3+"+string(var(2))+"^4");
208                                             }
209                                             if(m==7)
210                                             {
211                                                 return("The given function defines a simple singularity.
212The tjurina number is "+string(m)+".
213E^0[7]:"+string(var(1))+"^3+"+string(var(1))+"*"+string(var(2))+"^3");
214                                             }
215                                             if(m==8)
216                                             {
217                                                 return("The given function defines a simple singularity.
218The tjurina number is "+string(m)+".
219E^0[8]:"+string(var(1))+"^3+"+string(var(2))+"^5");
220                                             }
221                                        }
222                                        if(char(R)==5)
223                                        {
224                                            setring R;
225                                            if(m==6)
226                                            {
227                                               return("The given function defines a simple singularity.
228The tjurina number is "+string(m)+".
229E^0[6]:"+string(var(1))+"^3+"+string(var(2))+"^4");
230                                            }
231                                            if(m==7)
232                                            {
233                                               return("The given function defines a simple singularity.
234The tjurina number is "+string(m)+".
235E^0[7]:"+string(var(1))+"^3+"+string(var(1))+"*"+string(var(2))+"^3");
236                                            }
237                                            if(m==10)
238                                            {
239                                                 return("The given function defines a simple singularity.
240The tjurina number is "+string(m)+".
241E^0[8]:"+string(var(1))+"^3+"+string(var(2))+"^5");
242                                             }
243                                             if(m==8)
244                                             {
245                                                return("The given function defines a simple singularity.
246The tjurina number is "+string(m)+".
247E^1[8]:"+string(var(1))+"^3+"+string(var(2))+"^5+"+string(var(1))+"*"+string(var(2))+"^4");
248                                             }
249                                        }
250                                        if(char(R)==3)
251                                        {
252                                             poly p=blowupone(f);
253                                             int e=(std(jacob(p)+ideal(p))==1);
254                                             setring R;
255                                             if((m==7)&&e)
256                                             {
257                                                 return("The given function defines a simple singularity.
258The tjurina number is "+string(m)+".
259E^1[6]:"+string(var(1))+"^3+"+string(var(2))+"^4+"+string(var(1))+"^2*"+string(var(2))+"^2");
260                                             }
261                                             if((m==7)&&!e)
262                                             {
263                                                 return("The given function defines a simple singularity.
264The tjurina number is "+string(m)+".
265E^1[7]:"+string(var(1))+"^3+"+string(var(1))+"*"+string(var(2))+"^3+"+string(var(1))+"^2*"+string(var(2))+"^2");
266                                             }
267                                             if(m==8)
268                                             {
269
270                                                return("The given function defines a simple singularity.
271The tjurina number is "+string(m)+".
272E^2[8]:"+string(var(1))+"^3+"+string(var(2))+"^5+"+string(var(1))+"^2*"+string(var(2))+"^2");
273                                             }
274                                             if(m==10)
275                                             {
276                                                 return("The given function defines a simple singularity.
277The tjurina number is "+string(m)+".
278E^1[8]:"+string(var(1))+"^3+"+string(var(2))+"^5+"+string(var(1))+"^2*"+string(var(2))+"^3");
279                                             }
280                                            if((m==9)&&e)
281                                            {
282
283                                                return("The given function defines a simple singularity.
284The tjurina number is "+string(m)+".
285E^0[6]:"+string(var(1))+"^3+"+string(var(2))+"^4");
286                                            }
287                                            if((m==9)&&!e)
288                                            {
289
290                                                 return("The given function defines a simple singularity.
291The tjurina number is "+string(m)+".
292E^0[7]:"+string(var(1))+"^3+"+string(var(1))+"*"+string(var(2))+"^3");
293                                            }
294                                           if(m==12)
295                                           {
296
297                                                return("The given function defines a simple singularity.
298The tjurina number is "+string(m)+".
299E^0[8]:"+string(var(1))+"^3+"+string(var(2))+"^5");
300                                           }
301
302                                       }
303                                       else
304                                       {
305                                             return( "The given function defines not a simple singularity");
306                                       }
307
308                                   }
309
310                   }
311
312                   if(size(L[1])==3)
313                   {
314                           setring R;
315                           return("The given function defines a simple singularity.
316The tjurina number is "+string(m)+".
317D["+string(m)+"]:"+string(var(1))+"^2*"+string(var(2))+"+"+string(var(2))+"^"+string(m-1)+"");
318                   }
319                    if(size(L[1])==4)
320                   {
321                                         setring R;
322                                        return("The given function defines a simple singularity.
323The tjurina number is "+string(m)+".
324D[4]:"+string(var(1))+"^2*"+string(var(2))+"+"+string(var(1))+"^3");
325                   }
326          }
327     }
328  if(c>=3)
329  {
330     return( "The given function defines not a simple singularity.
331The tjurina number is "+string(m)+".");
332  }
333
334}                    // ends classifyCeq1
335 example
336{
337"EXAMPLE:"; echo=2;
338  ring R=5,(x,y),ds;
339  classifyCeq1(x2y+y22);
340}
341//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
342static proc GDsplitting(poly f)
343{      // the result of the splitting lemma (myMorsesplitting) in a special ring of 2 variables
344       // input a poly f of order >=2
345       // output  ring S=char(basering),(x,y),ds;  and a poly f in <x,y>^3 in S with the following properties:
346       // f + a sum of squares of further variables is contact equivalentt to the input polynomial
347       def R=basering;
348       intvec t;
349       int b,i;
350       b=tjurina(f);
351       f=myMorsesplitting(f,b);
352       t=forv(findVar(f));
353       ring S=char(basering),(x,y),ds;
354       ideal M;
355       M[t[1]]=var(1);
356       M[t[2]]=var(2);
357       i=1;
358       if((i!=t[1])&&(i!=t[2]))
359       {
360         M[i]=0;
361       }
362       map phi=R,M;
363       poly f=phi(f);
364       export(f);
365       setring(R);
366       return(S);
367}
368////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
369static proc myMorsesplitting(poly f,int b)
370{
371   // splitting lemma
372   // input  poly f=f(x(1),...,x(n))      jet(f,2)!=0
373   // output a polynomial in 2 variables of order >=3 such that a sum of squares of the remaining variables plus this
374   // polynomial is contact equivalent to the input polynomial
375   intvec w;
376   f=simplify(jet(f,b),1);
377   while(jet(f,2)!=0)
378   {
379      w=findlead(f);
380      if(w[1]==w[2])
381      {
382             f=splitting_one(f,w[1],b);
383      }
384      else
385      {
386             f=splitting_two(f,w[1],w[2],b);
387      }
388    }
389    return(f);
390}
391example
392{
393 "EXAMPLE:"; echo=2;
394  ring R=3,(x,y,z,u,v),ds;
395  poly f=x2+x3z+u2+v2+z3+y11;
396  myMorsesplitting(f,30);
397}
398////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
399static proc findlead(poly f)
400
401{
402  //  input poly f=x(i)^2+h      or   poly f=x(i)*x(j)+h  , h of order >=2
403  //  output intvec w   w[1]=i,w[2]=i   or  w[1]=i ,w[2]=j
404  intvec v=leadexp(f);
405  int i,k,n;
406  intvec w;
407  n=nvars(basering);
408  for(i=1;i<=n;i++)
409  {
410      if(v[i]==2)
411      {
412          w[1]=i;
413          w[2]=i;
414          break;
415      }
416     if(v[i]==1)
417      {
418          k++;
419          w[k]=i;
420      }
421  }
422  return(w);
423}
424example
425{
426 "EXAMPLE:"; echo=2;
427  ring R=3,(x,y,z,u,v),ds;
428  poly f=x2+x3z+u2+v2+z3+y11;
429  findlead(f);
430}
431////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
432static proc splitting_one(poly p, int i, int b)
433{
434// assumes that p=x_i^2+h, no x_i^2 in h, h of order >=2
435// returns q(x_1,??x_i-1,x_i+1,...,x_n) such that x_i^2 +q is right equivalent to p mod
436// <x_1,??,x_n>^b
437
438   if(b<2){b=2;}
439   def R=basering;
440   ideal M=maxideal(1);
441   map phi;
442   poly q=jet((p-var(i)^2)/var(i),b);
443   while(q!=0)
444   {
445      p=quickSubst(p,var(i)-1/2*q,i,b);
446      q=jet((p-var(i)^2)/var(i),b);
447   }
448   return(simplify(jet(p,b)-var(i)^2,1));                           //make the leading coefficient 1
449}
450example
451{
452 "EXAMPLE:"; echo=2;
453 ring R=3,(x,y,z,u,v),ds;
454 poly f=x2+y3+z4+xy3+u2+v2;
455 splitting_one(f,1,18);
456}
457////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
458static proc splitting_two(poly p, int i, int j, int b)
459{
460// assumes that p=x_i*x_j+h, no x_i^2 in h, h of order >=2
461// returns q(x_1,???,x_i-1,x_i+1,...,x_j-1,x_j+1,...,x_n) such that x_i*x_j +q is right equivalent to p mod
462// <x_1,??,x_n>^b
463
464   if(b<2){b=2;}
465   def R=basering;
466   ideal M=maxideal(1);
467   map phi;
468   poly q=jet((p-var(i)*var(j))/var(i),b);
469   while(q!=0)
470   {
471      p=quickSubst(p,var(j)-q,j,b);
472      q=jet((p-var(i)*var(j))/var(i),b);
473   }
474   return(simplify(substitute(jet(p,b),var(j),0),1));              //make the leading coefficient 1
475}
476example
477{
478 "EXAMPLE:"; echo=2;
479ring R=5,(u,v,w,s,t),ds;
480poly f=uv+t2+u4s+u11+v7+s9+w8;
481 splitting_two(f,1,2,128);
482}
483////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
484static proc quickSubst(poly h, poly r, int i, int b)
485{
486   //=== assume h, r is in Q[x_1,...,x_n], computes jet(h(x_i=r),b)
487   h=jet(h,b);
488   r=jet(r,b);
489   matrix M=coef(h,var(i));
490   poly q = 1;
491   int j,k,d;
492   intvec v;
493   d=deg(M[1,1]);
494   v[d+1]=1;
495   for(k = 2; k <= ncols(M); k++)
496   {
497        v[deg(M[1,k])+1]=1;
498   }
499   h=0;
500   for(k=1;k<=d+1;k++)
501   {
502      if(v[k]==1)
503      {
504         h=h+jet(q*M[2,ncols(M)-j],b);
505         j++;
506      }
507      q=jet(q*r,b);
508   }
509   return(h);
510}
511////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
512static proc findVar(poly p)
513{
514   //  input poly f
515   //  output intvec v   v[i]=0 if f is not depending on var(i) and v[j]=1 if f is depending on var(j)
516    intvec v;
517    int i,n;
518     n=nvars(basering);
519    for(i=1;i<=n;i++)
520    {
521           if(subst(p,var(i),0)==p)
522           {
523               v[i]=0;
524           }
525           else
526           {
527                v[i]=1;
528           }
529    }
530    return(v);
531}
532example
533{
534 "EXAMPLE:"; echo=2;
535  ring R=3,(x,y,z,u,v),ds;
536  poly f=y3+u5;
537  findlead(f);
538}
539////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
540static proc forv(intvec v)
541{
542       //returns the places of v which are different from zero
543       intvec w;
544       int i,j;
545       j=1;
546       for(i=1;i<=nvars(basering);i++)
547       {
548            if(v[i]!=0)
549            {
550               w[j]=i;
551                j++;
552            }
553        }
554        return(w);
555}
556
557/////////////////////////////////////////// char p=2 ////////////////////////////////////////////////////////////////////////////
558static proc classifyCeq2(poly p)
559//====Classification of hypersurface singularities in Characteristic p=2.
560{
561       //    input    poly p
562       //    output   The normal form to which f is contact equvalent or the function is not simple.
563       def R=basering;
564       int t=tjurina(p);
565       list T;
566       if(t==-1)
567       {
568          return("The given function defines not an isolated singularity");
569       }
570        int b=t+1;
571       p=SPILPRO(p,b);
572       list L=findVAR(p);
573       if(L[2]>=4)
574       {
575
576           return("The given function defines not a simple singularity.
577The Tjurina Number is "+string(t)+". ");
578       }
579       if(L[2]==1)
580       {
581          def S=redvar2(p);
582           setring S;
583
584          string a="The given function defines an isolated Singularity.
585The Tjurina number is "+string(t)+".
586A_"+string(leadexp(p)[1]-1)+":xy+z"+string(leadexp(p)[1])+".";
587          setring R;
588          return(a);
589       }
590
591       if(jet(p,2)==0)
592       {
593           if(L[2]==3)
594           {
595              return("The given function defines not a simple singularity.
596The Tjurina Number is "+string(t)+". ");
597
598           }
599           def S=redvar2(p);
600           setring S;
601           string a=curCeq2(p);
602           setring R;
603           return(a);
604       }
605       else
606       {
607         if(L[2]==3)
608         {
609             int B=t div 2+1;
610             p=splitting_SQUA(p,B);
611             def S=redvar2(p);
612             setring S;
613             string a=surCeq2(p);
614             setring R;
615             return(a);
616         }
617           p=splitting_SQUA(p,b);
618           def S=redvar2(p);
619           setring S;
620           string a=curCeq2(p);
621           setring R;
622           return(a);
623      }
624}
625example
626{
627   "EXAMPLE:"; echo=2;
628    ring r=2,(x,y,z,w,t,v),ds;
629    poly f=xy+zw+t5+v3;
630    classifyCeq2(f);
631}
632//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
633static proc SPILPRO(poly p,int B)
634{
635    //Splitting lemma in characteristic 2
636    //input a polynomial f and a bound B
637    //output a polynomial q, the result of the splitting lemma applied to f up to the order B
638  int i,j;
639  def R=basering;
640  int n=nvars(R);
641  poly q=splittingLchar2(p);
642  list T=FindPRO(q);
643
644  while(size(T)!=0)
645  {
646     i=T[1]; j=T[2];
647     q=splitting_two2(q,i,j,B);
648
649     T=FindPRO(q);
650  }
651  return(q);
652}
653////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
654static proc surCeq2(poly f)
655//====Classification of Surface Singularaties in case when p=2
656{
657        //===input a funation defined in the ring whose Characteristic is 2.
658        //===output is a Normal Form to which Given function is contact equivalent or function is not simple.
659        def R=basering;
660        int n=nvars(R);
661        list L;
662        int k=tjurina(f);
663        if(k==-1)
664        {
665           return("The given function defines not an isolated singularity.");
666        }
667        return(surEsing(f));
668}
669example
670{
671  "EXAMPLE:"; echo=2;
672   ring R=2,(x,y,z),Ds;
673   surCeq(xy+y2+yz+x2y+xy2+xyz+z21+xz21);
674
675}
676/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
677static proc surEsing(poly f)
678//====Classifcation of Ek singularities.
679{
680     //===input: A function which is an isolated singularity.
681     //===output: Normal form to which given function is contact equivalent.
682     list M;
683     def R=basering;
684     f=formf(f);
685     poly p=jet(f,2);
686     int k=tjurina(f);
687
688     list L=factorize(p);
689     poly g=subst(f,leadmonom(L[1][2]),0);
690     poly j=jet(g,3)-jet(g,2);
691     poly C=f-g;
692     C=C/leadmonom(L[1][2]);
693     poly h=subst(C,leadmonom(L[1][2]),0);
694     h=jet(h,2);
695     poly h1=j+h*leadmonom(L[1][2]);
696
697     if(h!=0)
698     {
699            list P=factorize(h1);
700            list N=factorize(j);
701
702            if((size(P[1])==2)&&(size(P[2])==2)&&(P[2][2]==1)&&(N[2][2]==3))
703            {
704                if(k==6)
705                {
706
707                   return(" The given function defines an isolated Singularity.
708The Tjurina number is "+string(tjurina(f))+".
709E^1[6]:z2+x3+y2z+xyz.");
710
711                }
712                if(k==8)
713                 {
714                    return("The given function defines an isolated Singularity.
715The Tjurina number is "+string(tjurina(f))+".
716E^0[6]:z2+x3+y2z.");
717
718                }
719            }
720            if(((size(P[1])==4)&&(size(P[2])==4))||((size(N[1])==4)&&(size(N[2])==4)))
721            {
722                if(k==8)
723                {
724                   return("The given function defines an isolated Singularity.
725The Tjurina number is "+string(tjurina(f))+".
726D[4]:z2+x2y+xy2.");
727                }
728                if(k==6)
729                {
730                   return("The given function defines an isolated Singularity.
731The Tjurina number is "+string(tjurina(f))+".
732D^1[4]:z2+x2y+xy2+xyz.");
733
734                 }
735            }
736            if((size(P[1])<=3)&&(size(P[2])<=3))
737            {
738
739                    if(k==6)
740                    {
741                        return("The given function defines an isolated Singularity.
742The Tjurina number is "+string(tjurina(f))+".
743D^1[5]:z2+x2y+y2z+xyz.");
744                    }
745                    if(k==8)
746                    {
747
748                           if(size(lengthBL(f))==4)
749                           {
750                                 return("The given function defines an isolated Singularity.
751The Tjurina number is "+string(tjurina(f))+".
752E^3[7]:z2+x3+xy3+xyz.");
753
754                           }
755                            else
756                           {
757                                if(size(lengthBL(f))==5)
758                                {
759
760                                    return("The given function defines an isolated Singularity.
761The Tjurina number is "+string(tjurina(f))+".
762E^4[8]:z2+x3+y5+xyz.");
763                                }
764                                else
765                                {
766                                    return("The given function defines an isolated Singularity.
767The Tjurina number is "+string(tjurina(f))+".
768D^0[5]:z2+x2y+y2z.");
769
770                                 }
771                            }
772
773                       }
774
775               }
776               list Q=factorize(j);
777               if((size(Q[1])==3)&&(size(Q[2])==3))
778               {
779
780                      return(surDsing(f));
781               }
782               else
783               {
784                        if(k==10)
785                        {
786                             if(size(lengthBL(f))==4)
787                             {
788                                return("The given function defines an isolated Singularity.
789The Tjurina number is "+string(tjurina(f))+".
790E^2[7]:z2+x3+xy3+y3z.");
791
792                             }
793                             if(size(lengthBL(f))==5)
794                             {
795                                return("The given function defines an isolated Singularity.
796The Tjurina number is "+string(tjurina(f))+".
797E^3[8]:z2+x3+y5+y3z.");
798
799                             }
800                        }
801                        if(k==12)
802                        {
803
804                             if(size(lengthBL(f))==4)
805                             {
806                                   return("The given function defines an isolated Singularity.
807The Tjurina number is "+string(tjurina(f))+".
808E^1[7]:z2+x3+xy3+xy3z.");
809
810                              }
811                              if(size(lengthBL(f))==5)
812                              {
813                                   return("The given function defines an isolated Singularity.
814The Tjurina number is "+string(tjurina(f))+".
815E^2[8]:z2+x3+y5+xy2z.");
816
817                              }
818                         }
819                         if(k==14)
820                         {
821
822                              if(size(lengthBL(f))==4)
823                              {
824                                     retrun("The given function defines an isolated Singularity.
825The Tjurina number is "+string(tjurina(f))+".
826E^0[7]:z2+x3+xy3.");
827
828                              }
829                              if(size(lengthBL(f))==5)
830                              {
831                                    return("The given function defines an isolated Singularity.
832The Tjurina number is "+string(tjurina(f))+".
833E^1[8]:z2+x3+y5+xy3z.");
834                              }
835                          }
836                          if(k==16)
837                          {
838                               return("The given function defines an isolated Singularity.
839The Tjurina number is "+string(tjurina(f))+".
840E^0[8]:z2+x3+y5.");
841                          }
842                  }
843     }
844     if(h==0)
845     {
846            list P=factorize(j);
847            if((size(P[1])==2)&&(size(P[2])==2)&&((P[2][2])==1))
848            {
849               if(k==6)
850               {
851                    return("The given function defines an isolated Singularity.
852The Tjurina number is "+string(tjurina(f))+".
853E^1[6]:z2+x3+y2z+xyz.");
854               }
855               if(k==8)
856               {
857                    return("The given function defines an isolated Singularity.
858The Tjurina number is "+string(tjurina(f))+".
859E^0[6]:z2+x3+y2z.");
860
861               }
862            }
863            if((size(P[1])==4)&&(size(P[2])==4))
864            {
865               if(k==8)
866               {
867                   return("The given function defines an isolated Singularity.
868The Tjurina number is "+string(tjurina(f))+".
869D^0[4]:z2+x2y+xy2.");
870
871               }
872               if(k==6)
873               {
874                   return("The given function defines an isolated Singularity.
875The Tjurina number is "+string(tjurina(f))+".
876D^1[4]:z2+x2y+xy2+xyz.");
877
878               }
879            }
880
881            if((size(P[1])==3)&&(size(P[2])==3))
882            {
883                if(((P[2][2])==1)&&((P[2][3])==1))
884                {
885                      if(k==6)
886                      {
887                           return("The given function defines an isolated Singularity.
888The Tjurina number is "+string(tjurina(f))+".
889D^1[5]:z2+x2y+y2z+xyz.");
890                      }
891                      if(k==8)
892                      {
893                          if(size(lengthBL(f))==4)
894                          {
895                              return("The given function defines an isolated Singularity.
896The Tjurina number is "+string(tjurina(f))+".
897E^3[7]:z2+x3+xy3+xyz.");
898                          }
899                          else
900                          {
901                                if(size(lengthBL(f))==5)
902                                {
903                                    return("The given function defines an isolated Singularity.
904The Tjurina number is "+string(tjurina(f))+".
905E^4[8]:z2+x3+y5+xyz.");
906
907                                 }
908                                 else
909                                 {
910                                     return("The given function defines an isolated Singularity.
911The Tjurina number is "+string(tjurina(f))+".
912D^0[5]:z2+x2y+y2z.");
913
914                                 }
915                            }
916                        }
917                   }
918
919              }
920
921              if((size(P[1])==2)&&((P[2][2])==3))
922              {
923                         if(k==10)
924                         {
925                              if(size(lengthBL(f))==4)
926                              {
927                                     return("The given function defines an isolated Singularity.
928The Tjurina number is "+string(tjurina(f))+".
929E^2[7]:z2+x3+xy3+y3z.");
930
931                              }
932                              if(size(lengthBL(f))==5)
933                              {
934                                     return("The given function defines an isolated Singularity.
935The Tjurina number is "+string(tjurina(f))+".
936E^3[8]:z2+x3+y5+y3z.");
937
938                               }
939                          }
940                          if(k==12)
941                          {
942
943                               if(size(lengthBL(f))==4)
944                               {
945                                     return("The given function defines an isolated Singularity.
946The Tjurina number is "+string(tjurina(f))+".
947E^1[7]:z2+x3+xy3+xy3z.");
948
949                                }
950                                if(size(lengthBL(f))==5)
951                                {
952                                     return("The given function defines an isolated Singularity.
953The Tjurina number is "+string(tjurina(f))+".
954E^2[8]:z2+x3+y5+xy2z.");
955
956                                 }
957                           }
958                           if(k==14)
959                           {
960
961                                if(size(lengthBL(f))==4)
962                                {
963                                     return("The given function defines an isolated Singularity.
964The Tjurina number is "+string(tjurina(f))+".
965E^0[7]:z2+x3+xy3.");
966
967                                }
968                                if(size(lengthBL(f))==5)
969                                {
970                                     return("The given function defines an isolated Singularity.
971The Tjurina number is "+string(tjurina(f))+".
972E^1[8]:z2+x3+y5+xy3z.");
973
974                                }
975                           }
976                           if(k==16)
977                           {
978                                return("The given function defines an isolated Singularity.
979The Tjurina number is "+string(tjurina(f))+".
980E^0[8]:z2+x3+y5.");
981
982                           }
983               }
984               return(surDsing(f));
985      }
986}
987
988///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
989static proc surDsing(poly f)
990//====Classification of Dk singularaties========================
991
992{
993      //===input: A function which is an isolated singularity.
994      //===output: Normal form to which given function is contact equivalent.
995
996     def R=basering;
997     int n=nvars(R);
998     f=formf(f);
999     int k=tjurina(f);
1000
1001     list M=factorize(jet(f,2));
1002
1003     poly g=subst(f,leadmonom(M[1][2]),0);
1004     poly j=jet(g,3);
1005     if(j==0)
1006     {
1007
1008       return("The given function defines not a simple singularity.
1009The Tjurina Number is "+string(k)+".");
1010
1011     }
1012     list P=factorize(j);
1013     if((size(P[1])==4)&&(size(P[2])==4))
1014     {
1015           if(k==8)
1016           {
1017
1018              return("The given function defines an isolated Singularity.
1019The Tjurina number is "+string(tjurina(f))+".
1020D^0 [4]:z2+x2y+xy2.");
1021
1022           }
1023           if(k==6)
1024           {
1025               return("The given function defines an isolated Singularity.
1026The Tjurina number is "+string(tjurina(f))+".
1027D^1 [4]:z2+x2y+xy2+xyz.");
1028
1029           }
1030     }
1031     if((size(P[1])==3)&&(size(P[2])==3))
1032     {
1033
1034         poly q=BlowUpO(f);
1035
1036         if((tjurina(f)-tjurina(q))==4)
1037         {
1038               list Q=whichSUR(f);
1039
1040               j=subst(Q[2],leadmonom(M[1][2]),0);
1041               poly j1=jet(j,3);
1042               list L=factorize(j1);
1043
1044               if((size(L[1])==4)&&(size(L[2])==4))
1045               {
1046                     return("The given function defines an isolated Singularity.
1047The Tjurina number is "+string(tjurina(f))+".
1048D^0 ["+string(tjurina(f) div 2)+"]:z2+x2y+xy"+string(tjurina(f) div 4)+". ");
1049
1050               }
1051               if((size(L[1])==3)&&(size(L[2])==3))
1052               {
1053                    return("The given function defines an isolated Singularity.
1054The Tjurina number is "+string(tjurina(f))+".
1055D^0 ["+string((tjurina(f) div 2)+1)+"]:z2+x2y+y"+string(tjurina(f) div 4)+"z.");
1056
1057               }
1058          }
1059          if((tjurina(f)-tjurina(q))==2)
1060          {
1061
1062               list Q=whichSUR(f);
1063               if(Q[1]==8){Q[2]=BlowUpO(Q[2]);}
1064               j=subst(Q[2],leadmonom(M[1][2]),0);
1065               poly j1=jet(j,3);
1066
1067               list L=factorize(j1);
1068               if((size(L[1])==4)&&(size(L[2])==4))
1069               {
1070                  return("The given function defines an isolated Singularity.
1071The Tjurina number is "+string(tjurina(f))+".
1072D^"+string(findSUR(f))+" ["+string((tjurina(f)+2*findSUR(f)) div 2)+"]:z2+x2y+xy^"+string((tjurina(f)+2*findSUR(f)) div 4)+"+ xy"+string(((tjurina(f)+2*findSUR(f)) div 4)-findSUR(f))+" z.");
1073
1074               }
1075               if((size(L[1])==3)&&(size(L[2])==3))
1076               {
1077               return("The given function defines an isolated Singularity.
1078The Tjurina number is "+string(tjurina(f))+".
1079D^"+string(findRSUR(f))+" ["+string(((tjurina(f)+2*findRSUR(f)) div 2)+1)+"]:z2+x2y+y"+string((tjurina(f)+2*findRSUR(f)) div 4)+"z+xy"+string(((tjurina(f)+2*findRSUR(f)) div 4)-findRSUR(f))+"z.");
1080
1081               }
1082          }
1083     }
1084}
1085/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1086static proc formf(poly f)
1087//=====Transform jet(f,2) in to the polynomial var(3)^2.
1088//=====This is the case when We are in Dk and Ek surface singularaties.
1089{
1090
1091        poly j=jet(f,2);
1092        list L=factorize(j);
1093        def r=basering;
1094        if(leadmonom(L[1][2])==var(1))
1095        {
1096            f=subst(f,var(1),L[1][2]);
1097            map phi=r,var(3),var(2),var(1);
1098            return(phi(f));
1099
1100        }
1101        if(leadmonom(L[1][2])==var(2))
1102        {
1103            f=subst(f,var(2),L[1][2]);
1104             map phi=r,var(1),var(3),var(2);
1105                return(phi(f));
1106
1107        }
1108        if(leadmonom(L[1][2])==var(3))
1109        {
1110             return(f);
1111        }
1112}
1113///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1114static proc whichSUR(poly f)
1115//====This procedure is required to separate the Surface Case D_2m from D_2m+1 as discribes in [2].
1116{
1117   int d=tjurina(f);
1118   list L;
1119   while(1)
1120   {
1121       f=BlowUpO(f);
1122       d=tjurina(f);
1123
1124       if((d==6)||(d==8))
1125       {
1126         L[1]=d;
1127         L[2]=f;
1128
1129         return(L);
1130       }
1131   }
1132   return(d);
1133}
1134///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1135static proc findSUR(poly f)//D2m,r
1136//====Number of blowup required in order either to get the difference equal to 4 or the Tjurina number is less than 6.
1137{
1138
1139     int a, r,b;
1140     while(1)
1141     {
1142         r++;
1143         a=tjurina(f);
1144         if(a<=1)
1145         {
1146            return(r);
1147         }
1148         f=BlowUpO(f);
1149         b=tjurina(f);
1150         if((a-b)==4)
1151         {
1152              if(b==2)
1153              {
1154                return(r);
1155              }
1156              return(r-1);
1157         }
1158         if(b<6)
1159         {
1160             return(r-1);
1161         }
1162
1163      }
1164}
1165/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1166static proc findRSUR(poly f)//D2m+1,r
1167//====Calculate the number of blow ups required by a polynomial in order to get the difference 4.
1168{
1169
1170    int a, r,b;
1171    while(1)
1172    {
1173       r++;
1174       a=tjurina(f);
1175
1176       if(a==2)
1177       {
1178           return(r);
1179       }
1180
1181       f=BlowUpO(f);
1182       b=tjurina(f);
1183
1184       if((a-b)==4)
1185       {
1186           return(r-1);
1187       }
1188       if(b<6)
1189       {
1190            return(r);
1191       }
1192    }
1193}
1194/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1195static proc BlowUpO(poly f)
1196//====Gives the procedure of Blowing up ingeneral
1197{
1198   //=== input f smooth or isolated singularity at zero
1199   //=== output var(1) or poly with isolated singularity at zero, the transformation
1200   //=== of the singularity of the blowing up to zero
1201    def R=basering;
1202   def S=changeord(list(list("ds",1:nvars(basering))));
1203   setring S;
1204   int n=nvars(basering);
1205   int i,t,c,d,e,j,k;
1206   poly f=imap(R,f);
1207   if(deg(lead(f))<=1){setring R;return(var(1));}
1208   poly p;
1209   def T=changeord(list(list("lp",1:nvars(basering))));
1210   setring T;
1211   list L;
1212   map phi,psi;
1213   ideal mphi, mpsi,sing;
1214   poly p,q,m,l;
1215   poly h=var(1);
1216   //===========  blow up=======================================================
1217   for(i=1;i<=n;i++)
1218   {
1219      mphi=var(i)*maxideal(1);
1220      mphi[i]=var(i);
1221      phi=S,mphi;
1222       p=phi(f);
1223      q=p/var(i);
1224      while(size(p)==size(q))
1225      {
1226         p=q;
1227         q=q/var(i);
1228      }
1229      //===============  p is the strict transform var(i) exceptional divisor ====
1230       //=============== analysis of singularities ================================
1231      sing=jacob(p),p,var(i);
1232      sing=radical(sing);
1233      option(redSB);
1234      sing=std(sing);
1235      if(dim(sing)>0){ERROR("not simple");}
1236      sing=std(simplify(sing,1));
1237      if(dim(sing)==0)
1238      {
1239        if(vdim(sing)==1)
1240            {
1241                  mpsi=std(maxideal(1));
1242                  for(k=1;k<=n;k++){mpsi[k]=var(k)-sing[n-k+1][2];}
1243                  psi=T,mpsi;
1244                  p=psi(p);
1245            }
1246            else
1247            {     //this can only happen in case of a D-singularity
1248                   L=minAssGTZ(sing);
1249                   d=0;
1250                   for(j=1;j<=size(L);j++)
1251                   {
1252                       sing=std(L[j]);
1253                       sing=std(simplify(sing,1));
1254                       if(vdim(sing)!=1){ERROR("something is wrong in blowUpO");}
1255                       mpsi=std(maxideal(1));
1256                       for(k=1;k<=n;k++){mpsi[k]=var(k)-sing[n-k+1][2];}
1257                       psi=T,mpsi;
1258                       m=psi(p);
1259                       setring S;
1260                       p=imap(T,m);
1261                       e=tjurina(p);
1262                       setring T;
1263                       if(e>d)
1264                       {
1265                            d=e;
1266                            l=m;
1267                       }
1268                    }
1269                   p=l;
1270            }
1271            setring S;
1272            p=imap(T,p);
1273            c=tjurina(p);
1274            setring T;
1275            if(c>t)
1276            {
1277                 t=c;
1278                 h=p;
1279            }
1280         }
1281       }
1282     setring R;
1283     poly h=imap(T,h);
1284     return(h);
1285}
1286////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1287static proc lengthBL(poly f)
1288//====Return the list of Tjurina numbers of each blowing up in the resolution before it becomes smooth.
1289{
1290    list L;
1291    int i=1;
1292    int d=tjurina(f);
1293    while(d>=2)
1294    {
1295        f=BlowUpO(f);
1296         L[i]=d;
1297        d=tjurina(f);
1298
1299         i=i+1;
1300    }
1301    return(L);
1302}
1303/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1304static proc curCeq2(poly f)
1305//====Classification of Curves in Case when p=2.
1306{
1307
1308     int c;
1309     int k;
1310     number m,a,b;
1311     ideal I;
1312     int d;
1313     k=tjurina(f);
1314     if(k==-1)
1315     {
1316           return("The given function defines not an isolated singularity");
1317
1318     }
1319     if(deg(lead(f))==2)
1320     {
1321          poly f1=jet(f,2);
1322          list T=factorize(f1);
1323          if(size(T[1])==3)
1324          {
1325              return("The given function defines an isolated Singularity.
1326The Tjurina number is "+string(tjurina(f))+".
1327A1:x2+xy.");
1328
1329          }
1330          if(size(T[1])==2)
1331          {
1332                poly g=BlowUpO(f);
1333                if((tjurina(f)-tjurina(g))==4)
1334                {
1335                      return("The given function defines an isolated Singularity.
1336The Tjurina number is "+string(tjurina(f))+".
1337A["+string(k div 2)+"]:x2+y "+string(tjurina(f) div 2+1)+". where r=0. ");
1338                }
1339                else
1340                {
1341                     d=whichtru(f);
1342                    if(d==1)
1343                    {
1344                         if((k mod 4)==0)
1345                         {
1346
1347                            return("The given function defines an isolated Singularity.
1348The Tjurina number is "+string(tjurina(f))+".
1349A["+string(2*(tjurina(f) div 2)-1)+"]:x2+xy"+string(tjurina(f) div 2)+".");
1350                         }
1351                         else
1352                         {
1353                            return("The given function defines an isolated Singularity.
1354The Tjurina number is "+string(tjurina(f))+".
1355A["+string((2*(tjurina(f)+1) div 2)-1)+"]:x2+xy"+string((tjurina(f)+1) div 2)+".");
1356
1357                         }
1358                    }
1359                    if(d==0)
1360                    {
1361                        if((findR(f) mod 2)==0)
1362                        {
1363                           return("The given function defines an isolated Singularity.
1364The Tjurina number is "+string(tjurina(f))+".
1365A^"+string(findR(f))+" ["+string(2*((tjurina(f)+2*findR(f)) div 4))+"]:x2+y"+string(((tjurina(f)+2*findR(f)) div 2)+1)+"+xy"+string(((tjurina(f)+2*findR(f)) div 2)-findR(f))+".");
1366                        }
1367                        if((findR(f) mod 2)!=0)
1368                        {
1369                           return("The given function defines an isolated Singularity.
1370The Tjurina number is "+string(tjurina(f))+".
1371A^"+string(findR(f))+" ["+string(2*(((tjurina(f)+1)+2*findR(f)) div 4))+"]:x2+y"+string((((tjurina(f)+1)+2*findR(f)) div 2)+1)+"+xy"+string((((tjurina(f)+1)+2*findR(f)) div 2)-findR(f))+".");
1372
1373                        }
1374                     }
1375
1376                 }
1377
1378           }
1379
1380       }
1381       if(deg(lead(f))==3)
1382       {
1383           if(jet(f,3)==0)
1384           {return("The given function defines not a simple Singularity.
1385             The Tjurina number is "+string(tjurina(f))+".")}
1386           poly f1=jet(f,3);
1387           list L=factorize(f1);
1388           if(size(L[1])==4)
1389           {
1390              return("The given function defines an isolated Singularity.
1391The Tjurina number is "+string(tjurina(f))+".
1392D[4]:x2y+xy2.");
1393
1394            }
1395            if(size(L[1])==3)
1396            {
1397                poly g=BlowUpO(f);
1398                if((tjurina(f)-tjurina(g))==8)
1399                {
1400                     return("The given function defines an isolated Singularity.
1401The Tjurina number is "+string(tjurina(f))+".
1402D^0["+string((tjurina(f) div 2)+1)+"]:x2y+y"+string(tjurina(f) div 2)+". where m="+string(tjurina(f)/4)+".");
1403
1404                 }
1405                 else
1406                 {
1407                      if(whichtru(f)==1)
1408                      {
1409                         return("The given function defines an isolated Singularity.
1410The Tjurina number is "+string(tjurina(f))+".
1411D["+string(tjurina(f))+"]:x2y+xy"+string(tjurina(f) div 2)+". ");
1412                       }
1413                       if(whichtru(f)==0)
1414                       {
1415                             if((findRD(f) mod 2)==0)
1416                             {
1417                            return("The given function defines an isolated Singularity.
1418The Tjurina number is "+string(tjurina(f))+".
1419D^"+string(findRD(f))+"["+string(2*((tjurina(f) div 2+findRD(f)) div 2)+1)+"]:x2y+y"+string((tjurina(f) div 2+findRD(f)))+"+xy"+string((tjurina(f) div 2+findRD(f))-findRD(f))+" where even r="+string(findRD(f))+".");
1420                              }
1421                              else
1422                              {
1423                               return("The given function defines an isolated Singularity.
1424The Tjurina number is "+string(tjurina(f))+".
1425D^"+string(findRD(f))+"["+string(2*(((tjurina(f)+1) div 2+findRD(f)) div 2)+1)+"]:x2y+y"+string(((tjurina(f)+1) div 2+findRD(f)))+"+xy"+string(((tjurina(f)+1) div 2+findRD(f))-findRD(f))+",odd r="+string(findRD(f))+".");
1426                               }
1427                         }
1428                 }
1429             }
1430             a=leadcoef(f);
1431             b=leadcoef(f-lead(f));
1432             f=subst(f,var(1),1/a*(var(1)-a*b*var(2)));
1433             I=var(1),var(2)^2;
1434             I=std(I^3);
1435             if(reduce(f,I)!=0)
1436             {
1437                if(k==6)
1438                 {
1439                     return("The given function defines an isolated Singularity.
1440The Tjurina number is "+string(tjurina(f))+".
1441E^1[6]:x3+y4+xy3.");
1442
1443                 }
1444                 if(k==7)
1445                 {
1446                      return("The given function defines an isolated Singularity.
1447The Tjurina number is "+string(tjurina(f))+".
1448E[7]:x3+xy3.");
1449                  }
1450                  if(k==8)
1451                  {
1452                       poly t=BlowUpO(f);
1453                       if(tjurina(t)==0)
1454                       {
1455                           return("The given function defines an isolated Singularity.
1456The Tjurina number is "+string(tjurina(f))+".
1457E^0[6]:x3+y4.");
1458                       }
1459                        else
1460                       {
1461                           return("The given function defines an isolated Singularity.
1462The Tjurina number is "+string(tjurina(f))+".
1463E[8]:x3+y5.");
1464                        }
1465                   }
1466               }
1467               else
1468               {return("The given function defines not a simple singularity.
1469                        The Tjurina number is "+string(tjurina(f))+".");}
1470      }
1471
1472      if(deg(lead(f))>3)
1473      {
1474
1475          return("The given function defines not a simple singularity.
1476                  The Tjurina number is "+string(tjurina(f))+".");
1477
1478      }
1479}
1480example
1481{
1482   "EXAMPLE:"; echo=2;
1483   ring R=2,(x,y),Ds;
1484   curCeq2(x3+x2y+xy2+y3+xy3+y4);
1485
1486}
1487///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1488static proc whichtru(poly f)
1489//====Gives the information that we do successive blowups and the last singularaty before becoming smmoth has Tjurina
1490//====number 1 if not than returns 0.(returns either a singularaty or smooth surve)============
1491{
1492   int d=tjurina(f);
1493   while(1)
1494   {
1495       f=BlowUpO(f);
1496       d=tjurina(f);
1497       if(d==1)
1498       {
1499         return(d);
1500       }
1501       if(d==0)
1502       {
1503          return(d);
1504       }
1505   }
1506  return(d);
1507}
1508///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1509static proc findR(poly f)
1510//====Find the number of blowups required by polynomial such that
1511//====difference between two consecutive blowups become 4
1512//====this procedure find the value of r in case of Ak singularaties as in [1].
1513{
1514    int a, r,b;
1515    while(1)
1516    {
1517        r++;
1518        a=tjurina(f);
1519        if(a<=1)
1520        {
1521            return(-1);
1522        }
1523        f=BlowUpO(f);
1524
1525        b=tjurina(f);
1526
1527        if((a-b)==4)
1528        {
1529            return(r-1);
1530        }
1531     }
1532}
1533////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1534static proc findRD(poly p)
1535//=====Find the number of blowups required by polynomial such
1536//=====that difference between two consecutive blowups become 4
1537//=====this procedure find the value of r in case of Dk singularaties as in [1].
1538{
1539
1540   int k=tjurina(p);
1541   poly q=BlowUpO(p);
1542   int t=tjurina(q);
1543   int r;
1544
1545   if((k mod 4)==0)
1546   {
1547      r=findR(BlowUpO(p))+2;
1548      return(r);
1549   }
1550   if((k mod 4 !=0)&&(t mod 4==0))
1551   {
1552     if(t!=0)
1553     {
1554        r=findR(BlowUpO(p))+1;
1555       return(r);
1556     }
1557   }
1558   if((k mod 2)==0)
1559   {
1560      r=findR(BlowUpO(p))+2;
1561      return(r);
1562   }
1563
1564}
1565//////////////////////////////////////////////////////////////////////////////////////////////////////////
1566static proc splitting_SQUA(poly p, int b)
1567 {
1568         // assumes that p=x_i^2+h, h of order >=3
1569         // returns q(x_1,??,x_n) such that x_i^2 +q is right equivalent to p mod
1570          // <x_1,??,x_n>^b and x_i is only linear in q
1571            def R=basering; int n=nvars(R);
1572            ideal M=maxideal(1);
1573            int j;
1574             map phi;
1575             p=jet(p,b);
1576             p=simplify(p,1);
1577             list T=FindSQUA(p);
1578             int i=T[1];
1579             poly q=p/var(i)^2;
1580             while(q!=1)
1581             {
1582                for(j=1;j<=n;j++)
1583                {
1584                   M[j]=var(j)*q;
1585                }
1586                phi=R,M;
1587                p=phi(p);
1588                p=p*inverseUnit(q,b);
1589                p=jet(p,b);
1590                q=p/var(i)^2;
1591             }
1592             return(p);
1593 }
1594///////////////////////////////////////////////////////////////////////////////////////////////////////////////
1595static proc splitting_two2(poly p, int i, int j, int b)
1596{
1597     // assumes that p=x_i*x_j+h, no x_i^2, no in h, h of order >=2
1598     // returns q(x_1,??x_i-1,x_i+1,...,x_j-1,x_j+1,...,x_n) such that x_i*x_j +q is right equivalent to p mod
1599     // <x_1,??,x_n>^b
1600      if(b<2){b=2;}
1601      def R=basering;
1602      poly q=jet((p-var(i)*var(j))/var(i),b);
1603      while(q!=0)
1604      {
1605         p=quickSubst2(p,var(j)-q,j,b);
1606         q=jet((p-var(i)*var(j))/var(i),b);
1607      }
1608      return(simplify(substitute(jet(p,b),var(j),0),1));              //make the leading coefficient 1
1609
1610}
1611////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1612static proc quickSubst2(poly h, poly r, int i, int b)
1613{
1614    //=== assume h, r is in Q[x_1,...,x_n], computes jet(h(x_i=r),b)
1615     h=jet(h,b);
1616     r=jet(r,b);
1617     matrix M=coef(h,var(i));
1618     poly q = 1;
1619     int j,k,d;
1620     intvec v;
1621     d=deg(M[1,1]);
1622     v[d+1]=1;
1623     for(k = 2; k <= ncols(M); k++)
1624     {
1625         v[deg(M[1,k])+1]=1;
1626     }
1627     h=0;
1628     for(k=1;k<=d+1;k++)
1629     {
1630        if(v[k]==1)
1631        {
1632            h=h+jet(q*M[2,ncols(M)-j],b);
1633            j++;
1634        }
1635        q=jet(q*r,b);
1636      }
1637   return(h);
1638}
1639////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1640static proc FindPRO(poly f)
1641{
1642        //input the polynomial
1643      // output is a list T where T[1]=the first variable which appear as product in the quadraic form.
1644      //                          T[2]=the second variable which appear as product with the var(T[1]) in the quadraic form.
1645  def R=basering;
1646  int n=nvars(R);
1647  int i,j,k;
1648  list T;
1649  for(i=1;i<=n;i++)
1650  {
1651     for(k=i+1;k<=n;k++)
1652     {
1653         for(j=1;j<=size(f);j++)
1654         {
1655             if(leadmonom(f[j])==var(i)*var(k))
1656             {
1657
1658                T[1]=i;
1659                T[2]=k;
1660                return(T);
1661             }
1662         }
1663     }
1664
1665  }
1666   return(T);
1667}
1668///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1669static proc FindSQUA(poly f)
1670{
1671      //input the polynomial
1672      // output is a list T where T[1]=the first variable which appears as a square in the quadraic part of f.
1673      //                          T[2]= the number of squares appearing in the quadraic part of f.
1674  def R=basering;
1675  int n=nvars(R);
1676  int i,j,k;
1677  list T;
1678  for(i=1;i<=n;i++)
1679  {
1680       for(j=1;j<=size(f);j++)
1681       {
1682           if(leadmonom(f[j])==var(i)^2)
1683           {
1684                k++;
1685                T[1]=i;
1686                T[2]=k;
1687                return(T);
1688            }
1689       }
1690   }
1691   return(T);
1692}
1693/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1694static proc findVAR(poly f)
1695{
1696     //input: a poly f
1697     //output a list L=v,k v an intvec, v[i]=0, if depends on var(i), v[i]=1 else
1698     //k is the number of variables occurring in f
1699     intvec v;
1700     int i,k;int n=nvars(basering);
1701     list L;
1702     for(i=1;i<=n;i++)
1703     {
1704        if(f==subst(f,var(i),0))
1705        {
1706           v[i]=1;
1707        }
1708        else
1709        {
1710           v[i]=0;
1711           k++;
1712        }
1713     }
1714     L=v,k;
1715     return(L);
1716}
1717//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1718static proc redvar2(poly p)
1719{
1720       //input  a polynomial depending on x_i_1,...,x_i_k
1721       //output a new ring with k variables containing the polynomial
1722   def R=basering;
1723   int n=nvars(R);
1724   list N=findVAR(p);
1725   int m=N[2];
1726   intvec v=N[1];
1727   int l;
1728   int i;
1729   if(n==m)
1730   {
1731
1732        def S=R;
1733        if(defined(h))
1734        {
1735           kill h;
1736        }
1737        poly h;
1738        export h;
1739
1740   }
1741   else
1742   {
1743      def S=defring("2",m,"u","ds");
1744      setring S;
1745      ideal M;
1746      for(i=1;i<=n;i++)
1747      {
1748           if(v[i]==1)
1749           {
1750              M[i]=0;
1751           }
1752           else
1753           {
1754              l++;
1755              M[i]=var(l);
1756           }
1757       }
1758       map phi=R,M;
1759       poly p=phi(p);
1760       export p;
1761       setring R;
1762    }
1763    return(S);
1764}
1765////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1766static proc inverseUnit(poly q,int b)
1767{
1768   //input a polynomial q with non-zero constant part
1769   //output the inverse of q up to order b as a power series
1770   number c=leadcoef(q);
1771   q=q/c;
1772   int i;
1773   poly u=1;
1774   poly a=q-1;
1775   poly s=-a;
1776   for(i=1;i<=b;i++)
1777   {
1778      u=u+s;
1779      s=s*a;
1780   }
1781   return(u/c);
1782}
1783/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1784static proc splittingLchar2(poly f)
1785{
1786   //the quadratic part of the splitting lemma in characteristic 2
1787   //input: poly f
1788   //output: a polynomial right equivalent to the input polynomial such that the quadratic part is in normal form
1789   poly p=jet(f,2);
1790   if(p==0){return(f);}
1791   if(!homog(p)){return(jet(f,1));}
1792   def R=basering;
1793   def S=changeord(list(list("lp",1:nvars(basering))));
1794   setring S;
1795   ideal ma=maxideal(1);
1796   number a;
1797   int c=ringlist(S)[1][1];
1798   poly f=imap(R,f);
1799   poly p=imap(R,p);
1800   poly h,q;
1801   intvec v=leadexp(p);
1802   int i,j,k;
1803   intvec w;
1804   for(k=1;k<=size(v);k++)
1805   {
1806      if(v[k]!=0){w[size(w)+1]=k;}
1807   }
1808   i=w[2];
1809   if(size(w)==3){j=w[3];}
1810   if(j>0)
1811   {
1812      a=1/leadcoef(p);
1813      ma[j]=a*(var(j)+(p-lead(p))/var(i));
1814      map phi=S,ma;
1815      p=phi(p);
1816      f=phi(f);
1817      q=(p-lead(p))/var(j);
1818      ma=maxideal(1);
1819      ma[i]=var(i)+q;
1820      phi=S,ma;
1821      f=phi(f);
1822      h=splittingLchar2(f-var(i)*var(j));
1823      if(jet(h,2)==0)
1824     {
1825         setring R;
1826         f=imap(S,f);
1827         return(f);
1828      }
1829      v=leadexp(h);
1830      for(k=1;k<=size(v);k++)
1831      {
1832         if(v[k]!=0){break;}
1833      }
1834      if(v[k]==2)
1835      {
1836         ma=maxideal(1);
1837         ma[i]=var(j);
1838         ma[j]=var(k);
1839         ma[k]=var(i);
1840         phi=S,ma;
1841         h=phi(h);
1842         f=h+var(j)*var(k);
1843      }
1844      else
1845      {
1846         f=h+var(i)*var(j);
1847      }
1848      setring R;
1849      f=imap(S,f);
1850      return(f);
1851   }
1852
1853   a= leadcoef(p)^(c div 2);
1854   ma[i]=1/a*var(i);
1855   map phi=S,ma;
1856   p=phi(p);
1857   f=phi(f);
1858   q=p/var(i);
1859   if(size(q)==1)
1860  {
1861      if(p==var(i)^2)
1862      {
1863         setring R;
1864         f=imap(S,f);
1865         return(f);
1866      }
1867      h=splittingLchar2(f-var(i)^2);
1868      p=jet(h,2);
1869      v=leadexp(p);
1870      h=h+var(i)^2;
1871
1872      for(k=1;k<=size(v);k++)
1873      {
1874         if(v[k]!=0){j=k;break;}
1875      }
1876      if(v[j]==2)
1877      {
1878         ma=maxideal(1);
1879         ma[i]=var(i)+var(j);
1880         phi=S,ma;
1881         h=phi(h);
1882         h=splittingLchar2(h);
1883      }
1884      setring R;
1885      f=imap(S,h);
1886      return(f);
1887   }
1888   else
1889   {
1890      q=q-lead(q);
1891      v=leadexp(q);
1892      for(k=1;k<=size(v);k++)
1893      {
1894         if(v[k]!=0){j=k;break;}
1895      }
1896      ma=maxideal(1);
1897      ma[j]=1/leadcoef(q)*(var(j)+q-lead(q));
1898      phi=S,ma;
1899      f=phi(f);
1900      ma=maxideal(1);
1901      ma[j]=var(i);
1902      ma[i]=var(j);
1903      phi=S,ma;
1904      f=phi(f);
1905      p=jet(f,2);
1906      if(leadmonom(p)==var(i)^2)
1907      {
1908                                            f=phi(f);
1909         p=jet(f,2);
1910         ma=maxideal(1);
1911         ma[i]=p/var(j)-var(j);
1912         phi=S,ma;
1913         f=phi(f);
1914         h=splittingLchar2(f-var(i)^2-var(i)*var(j)-var(j)^2);
1915         p=jet(h,2);
1916
1917         f=var(i)^2+var(i)*var(j)+var(j)^2+h;
1918         if(p!=0)
1919         {
1920            v=leadexp(p);
1921            for(k=1;k<=size(v);k++)
1922            {
1923               if(v[k]!=0){break;}
1924            }
1925            if(v[k]==2)
1926            {
1927               ma=maxideal(1);
1928               ma[k]=var(i)+var(k);
1929               phi=S,ma;
1930               f=phi(f);
1931               setring R;
1932               f=imap(S,f);
1933               return(splittingLchar2(f));
1934            }
1935            setring R;
1936            f=imap(S,f);
1937            return(f);
1938         }
1939         else                     //we return (x_i)^2+x_i*x_j+(x_j)^2 since we need a field extension
1940         {                        //to transform it to x_i*x_j
1941            setring R;
1942            f=imap(S,f);
1943            return(f);
1944         }
1945      }
1946      setring R;
1947      f=imap(S,f);
1948      return(splittingLchar2(f));
1949   }
1950}
1951///////////////////////////////////////////////////////////////////////////////////////////////////////////
1952/*
1953===============================   Examples for characteristic 2 ==========================================
1954ring r=2,(x,y,z,w,t,v),ds;
1955poly p=xy+zw+t5+v3;
1956classifyCeq(p);
1957
1958map phi=r,v,t,w,z,y,x;
1959poly c=1+x+y+z;
1960poly q=c*p;
1961q=phi(q);
1962q=zw+tv+x3+zw2+zwt+zwv+wtv+t2v+tv2+x3w+x3t+x3v+y5+y5w+y5t+y5v;
1963classifyCeq(q);
1964
1965poly p=ztv+x3+zw2+zwt+zwv+wtv+tv2+tv2+x3w+x3t+x3v+xy3+xy3w+xy3v;
1966classifyCeq(p);
1967
1968poly p=zw+tv+x3+zw2+zwt+wtv+t2v+tv2+x3w+x3t+x3v+y4+y4w+y4t+y4v;
1969classifyCeq(p);//E^0[6]
1970
1971poly p=zw+tv+x3+zw2+zwt+zwv+wtv+t2v+tv2+x3w+x3t+x3v+xy3+y4+xy3w+xy3t+xy3v+y4w+y4t+y4v;
1972classifyCeq(p);//E^1[6]
1973
1974poly p=y2+zw+tv+y2w+y2t+y2v+zw2+zwt+zwv+wtv+t2v+tv2+x11+x11w+x11t+x11v;
1975classifyCeq(p);//A[10] r=0;
1976
1977poly p=zw+tv+xy2+zw2+zwv+zwt+wtv+t2v+tv2+xy2w+xy2t+xy2v+x14y+x14yw+x14yt+x14yv+x26+x26w+x26t+x26v;
1978classifyCeq(p);
1979
1980ring r=2,(x,y,z,t,w,u,v),Ds;
1981ideal I=z+t,x+t+w,v,y,u,t,w+z;
1982det(jacob(I));
1983map phi=r,I;
1984poly a=1+x+u+v;
1985poly p=xy+zt+wu+v19;
1986poly q=a*p;
1987poly j=phi(q);
1988classifyCeq(j);
1989classifyCeq(p);
1990
1991poly p=x2+yz+tw+u3+v2x;
1992classifyCeq(p);
1993
1994poly q=a*p;
1995poly j=phi(q);
1996classifyCeq(j);
1997
1998poly p=xy+zw+t2+u2v+v19u;
1999classifyCeq(p);
2000
2001poly q=a*p;
2002poly j=phi(q);
2003classifyCeq(j);
2004
2005poly p=xy+zw+t2+u2v+v5t+uvt;
2006classifyCeq(p);
2007
2008poly q=a*p;
2009poly j=phi(q);
2010classifyCeq(j);
2011
2012poly p=xy+zw+t2+u2v+v5u+uvt;
2013classifyCeq(p);
2014
2015poly q=a*p;
2016poly j=phi(q);
2017classifyCeq(j);
2018
2019ring r=2,(x,y,z,t,w,u,v,e),Ds;
2020ideal I=x+y+z+t+w,y+z+t+w,e,v,u,t,w,z;
2021map si=r,I;
2022
2023poly p=xy+zt+wu+v2+e81;
2024poly j=si(p);
2025classifyCeq(j);
2026
2027poly p=xy+zt+wu+v2+e81+e^(80-12)*v;
2028classifyCeq(p);
2029poly j=si(p);
2030classifyCeq(j);
2031
2032poly p=xy+zt+wu+v2e+e60+e31v;
2033classifyCeq(p);
2034
2035poly p=xy+zt+wu+v2e+e60+e59v;
2036poly j=si(p);
2037classifyCeq(j);
2038===============================   Examples for characteristic > 2 ========================================
2039ring R=3,(x,y,z,u,v,w,s,t),ds;
2040ideal M=maxideal(1);
2041M[3]=x-y+z+u+w-s+t;
2042poly f=u3+v4+u2v2+x2+y2+z2+s2+t2+w2;
2043classifyCeq(f);
2044
2045map phi=R,M;
2046f=phi(f);
2047classifyCeq(f);
2048
2049f=u2+y32+s2+t2+v2+x2+z2+w2;
2050classifyCeq(f);
2051
2052map phi=R,M;
2053f=phi(f);
2054classifyCeq(f);
2055
2056
2057f=v2w+w61+x2+y2+z2+t2+z2+s2+u2;
2058classifyCeq(f);
2059
2060f=phi(f);
2061classifyCeq(f);
2062
2063f=u2+v2+w2+s2+t2-x2y+xy2+z6+y11+xy11+x14+z53;  // f is of corank=3
2064classifyCeq(f);
2065
2066f=x3+y4+xy5+s2+t2;
2067classifyCeq(f);
2068
2069f=u3+v5+u2v2+x2+y2+z2+t2+s2+w2;
2070classifyCeq(f);
2071
2072ideal M=maxideal(1);
2073M[2]=x+y+2s+t;
2074map phi=R,M;
2075f=phi(f);
2076classifyCeq(f);
2077
2078ring R=3,(u,v,w,s,t),ds;
2079poly f=w2s+s21+u2+v2+t2;
2080classifyCeq(f);
2081
2082ideal M=maxideal(1);
2083M[3]=u+2v+w+2t+2s;
2084map phi=R,M;
2085f=phi(f);
2086classifyCeq(f);
2087
2088ring R=5,(x,y,z,u,v,w,s,t),dp;
2089poly f=v2w+w61+x2+y2+z2+t2+z2+s2+u2;
2090classifyCeq(f);
2091
2092*/
Note: See TracBrowser for help on using the repository browser.