source: git/Singular/LIB/realclassify.lib @ e1b841

spielwiese
Last change on this file since e1b841 was e1b841, checked in by Andreas Steenpass <steenpass@…>, 12 years ago
update documentation for parallel.lib and realclassify.lib
  • Property mode set to 100644
File size: 30.1 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Singularities";
4info="
5LIBRARY:  realclassify.lib   Classification of real singularities
6AUTHOR:   Magdaleen Marais,  magdaleen@aims.ac.za
7          Andreas Steenpass, steenpass@mathematik.uni-kl.de
8
9OVERVIEW:
10   A library for classifying isolated hypersurface singularities over the reals
11   w.r.t. right equivalence, based on the determinator of singularities by
12   V.I. Arnold. This library is based on classify.lib by Kai Krueger, but
13   handles the real case, while classify.lib does the complex classification.
14
15REFERENCES:
16Arnold, Varchenko, Gusein-Zade: Singularities of Differentiable Maps.
17Vol. 1: The classification of critical points caustics and wave fronts.
18Birkh\"auser, Boston 1985
19
20Greuel, Lossen, Shustin: Introduction to singularities and deformations.
21Springer, Berlin 2007
22
23PROCEDURES:
24 realclassify(f);    real classification of singularities of modality 0 and 1
25 realmorsesplit(f);  splitting lemma in the real case
26 milnornumber(f);    Milnor number
27 determinacy(f);     an upper bound for the determinacy
28";
29
30LIB "classify.lib";
31LIB "rootsur.lib";
32LIB "atkins.lib";
33LIB "solve.lib";
34///////////////////////////////////////////////////////////////////////////////
35proc realclassify(poly f, list #)
36"
37USAGE:    realclassify(f[, printcomments]); f poly, printcomments int
38RETURN:   A list containing (in this order)
39          @* - the type of the singularity as a string,
40          @* - the normal form,
41          @* - the corank, the Milnor number, the inertia index and
42               a bound for the determinacy as integers.
43          @* The normal form involves parameters for singularities of modality
44             greater than 0. The actual value of the parameters is not computed
45             in most of the cases. If the value of the parameter is unknown,
46             the normal form is given as a string with an \"a\" as the
47             parameter. Otherwise, it is given as a polynomial.
48          @* An optional integer printcomments can be provided. If its value
49             is non-zero, a string will be added at the end of the returned
50             list, containing the result in more readable form and in some
51             cases also more comments on how to interpret the result. The
52             default is zero.
53NOTE:     The classification is done over the real numbers, so in contrast to
54          classify.lib, the signs of coefficients of monomials where even
55          exponents occur matter.
56          @* The ground field must be Q (the rational numbers). No field
57             extensions of any kind nor floating point numbers are allowed.
58          @* The monomial order must be local.
59          @* The input polynomial must be contained in maxideal(2) and must be
60             an isolated singularity of modality 0 or 1. The Milnor number is
61             checked for being finite.
62SEE ALSO: classify
63KEYWORDS: Classification of singularities
64EXAMPLE:  example realclassify; shows an example"
65{
66  /* auxiliary variables */
67  int i, j;
68
69  /* name for the basering */
70  def br = basering;
71
72  /* read optional parameters */
73  int printcomments;
74  if(size(#) > 0)
75  {
76    if(size(#) > 1 || typeof(#[1]) != "int")
77    {
78      ERROR("Wrong optional parameters.");
79    }
80    printcomments = #[1];
81  }
82
83  /* error check */
84  if(charstr(br) != "0")
85  {
86    ERROR("The ground field must be Q (the rational numbers).");
87  }
88  int n = nvars(br);
89  for(i = 1; i <= n; i++)
90  {
91    if(var(i) > 1)
92    {
93      ERROR("The monomial order must be local.");
94    }
95  }
96  if(jet(f, 1) != 0)
97  {
98    ERROR("The input polynomial must be contained in maxideal(2).");
99  }
100
101  /* compute Milnor number before continuing the error check */
102  int mu = milnornumber(f);
103
104  /* continue error check */
105  if(mu < 1)
106  {
107    ERROR("The Milnor number of the input polynomial must be"+newline
108      +"positive and finite.");
109  }
110
111  /* call classify before continuing the error check */
112  list dataFromClassify = prepRealclassify(f);
113  int m = dataFromClassify[1];                // the modality of f
114  string complextype = dataFromClassify[2];   // the complex type of f
115
116  /* continue error check */
117  if(m > 1)
118  {
119    ERROR("The input polynomial must be a singularity of modality 0 or 1.");
120  }
121
122  /* apply splitting lemma */
123  list morse = realmorsesplit(f, mu);
124  int cr = morse[1];
125  int lambda = morse[2];
126  int d = morse[3];
127  poly rf = morse[4];
128
129  /* determine the type */
130  string typeofsing, typeofsing_alternative;
131  poly nf, nf_alternative;
132  poly monparam;   // the monomial whose coefficient is the parameter
133                   // in the modality 1 cases, 0 otherwise
134  string morecomments = newline;
135  map phi;
136  if(cr == 0)   // case A[1]
137  {
138    if(lambda < n)
139    {
140      typeofsing = "A[1]+";
141      typeofsing_alternative = "A[1]-";
142    }
143    else
144    {
145      typeofsing = "A[1]-";
146      typeofsing_alternative = "A[1]+";
147    }
148  }
149  if(cr == 1)   // case A[k], k > 1
150  {
151    int k = deg(lead(rf), 1:n)-1;
152    if(k%2 == 0)
153    {
154      nf = var(1)^(k+1);
155      nf_alternative = nf;
156      typeofsing = "A["+string(k)+"]";
157      typeofsing_alternative = typeofsing;
158    }
159    else
160    {
161      if(leadcoef(rf) > 0)
162      {
163        nf = var(1)^(k+1);
164        typeofsing = "A["+string(k)+"]+";
165        typeofsing_alternative = "A["+string(k)+"]-";
166      }
167      else
168      {
169        nf = -var(1)^(k+1);
170        typeofsing = "A["+string(k)+"]-";
171        typeofsing_alternative = "A["+string(k)+"]+";
172      }
173      nf_alternative = -nf;
174    }
175  }
176  if(cr == 2)
177  {
178    if(complextype[1,2] == "D[")   // case D[k]
179    {
180      if(mu == 4)   // case D[4]
181      {
182        rf = jet(rf, 3);
183        number s1 = number(rf/(var(1)^3));
184        number s2 = number(rf/(var(2)^3));
185        if(s2 == 0 && s1 != 0)
186        {
187          phi = br, var(2), var(1);
188          rf = phi(rf);
189        }
190        if(s1 == 0 && s2 == 0)
191        {
192          number t1 = number(rf/(var(1)^2*var(2)));
193          number t2 = number(rf/(var(2)^2*var(1)));
194          if(t1+t2 == 0)
195          {
196            phi = br, var(1), 2*var(2);
197            rf = phi(rf);
198          }
199          phi = br, var(1)+var(2), var(2);
200          rf = phi(rf);
201        }
202        ring R = 0, y, dp;
203        map phi = br, 1, y;
204        poly rf = phi(rf);
205        int k = nrroots(rf);
206        setring(br);
207        if(k == 3)
208        {
209          nf = var(1)^2*var(2)-var(2)^3;
210          typeofsing = "D[4]-";
211        }
212        else
213        {
214          nf = var(1)^2*var(2)+var(2)^3;
215          typeofsing = "D[4]+";
216        }
217        typeofsing_alternative = typeofsing;
218        nf_alternative = nf;
219      }
220      else   // case D[k], k > 4
221      {
222        rf = jet(rf, d);
223        list factorization = factorize(jet(rf, 3));
224        list factors = factorization[1][2];
225        if(factorization[2][2] == 2)
226        {
227          factors = insert(factors, factorization[1][3], 1);
228        }
229        else
230        {
231          factors = insert(factors, factorization[1][3]);
232        }
233        factors[2] = factorization[1][1]*factors[2];
234        matrix T[2][2] = factors[1]/var(1), factors[1]/var(2),
235                         factors[2]/var(1), factors[2]/var(2);
236        phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
237        rf = phi(rf);
238        rf = jet(rf, d);
239        poly g;
240        for(i = 4; i < mu; i++)
241        {
242          g = jet(rf, i) - var(1)^2*var(2);
243          if(g != 0)
244          {
245            phi = br, var(1)-(g/(var(1)*var(2)))/2,
246                      var(2)-(g/var(1)^i)*var(1)^(i-2);
247            rf = phi(rf);
248            rf = jet(rf, d);
249          }
250        }
251        number a = number(rf/var(2)^(mu-1));
252        if(a > 0)
253        {
254          typeofsing = "D["+string(mu)+"]+";
255          nf = var(1)^2*var(2)+var(2)^(mu-1);
256          if(mu%2 == 0)
257          {
258            nf_alternative = nf;
259            typeofsing_alternative = typeofsing;
260          }
261          else
262          {
263            nf_alternative = var(1)^2*var(2)-var(2)^(mu-1);
264            typeofsing_alternative = "D["+string(mu)+"]-";
265          }
266        }
267        else
268        {
269          typeofsing = "D["+string(mu)+"]-";
270          nf = var(1)^2*var(2)-var(2)^(mu-1);
271          if(mu%2 == 0)
272          {
273            nf_alternative = nf;
274            typeofsing_alternative = typeofsing;
275          }
276          else
277          {
278            nf_alternative = var(1)^2*var(2)+var(2)^(mu-1);
279            typeofsing_alternative = "D["+string(mu)+"]+";
280          }
281        }
282      }
283    }
284    if(complextype == "E[6]")  // case E[6] ;
285    {
286      poly g = jet(rf,3);
287      number s = number(g/(var(1)^3));
288      if(s == 0)
289      {
290        phi = br, var(2), var(1);
291        rf = phi(rf);
292        g = jet(rf,3);
293      }
294      list Factors = factorize(g);
295      poly g1 = Factors[1][2];
296      phi = br, (var(1)-(g1/var(2))*var(2))/(g1/var(1)), var(2);
297      rf = phi(rf);
298      rf = jet(rf,4);
299      number w = number(rf/(var(2)^4));
300      if(w > 0)
301      {
302        typeofsing = "E[6]+";
303        nf = var(1)^3+var(2)^4;
304        typeofsing_alternative = "E[6]-";
305        nf_alternative = var(1)^3-var(2)^4;
306      }
307      else
308      {
309        typeofsing = "E[6]-";
310        nf = var(1)^3-var(2)^4;
311        typeofsing_alternative = "E[6]+";
312        nf_alternative = var(1)^3+var(2)^4;
313      }
314    }
315    if(complextype == "E[7]")   // case E[7]
316    {
317      typeofsing = "E[7]";
318      nf = var(1)^3+var(1)*var(2)^3;
319      typeofsing_alternative = typeofsing;
320      nf_alternative = nf;
321    }
322    if(complextype == "E[8]")   // case E[8]
323    {
324      typeofsing = "E[8]";
325      nf = var(1)^3+var(2)^5;
326      typeofsing_alternative = typeofsing;
327      nf_alternative = nf;
328    }
329    if(complextype == "J[2,0]")  // case J[10]
330    {
331      int signforJ10;
332      poly g = jet(rf,3);
333      number s = number(g/(var(1)^3));
334      if (s == 0)
335      {
336        phi = br, var(2), var(1);
337        rf = phi(rf);
338        g = jet(rf,3);
339      }
340      list Factors = factorize(g);
341      poly g1 = Factors[1][2];
342      phi = br, (var(1)-(g1/var(2))*var(2))/(g1/var(1)), var(2);
343      rf = phi(rf);
344      poly rf3 = jet(rf,3);
345      number w0 = number(rf3/(var(1)^3));
346      if(w0 < 0)
347      {
348        phi = br, -var(1), var(2);
349        rf = phi(rf);
350      }
351      rf3 = jet(rf,3);
352      poly rf4 = jet(rf,4);
353      poly rf5 = jet(rf,5);
354      poly rf6 = jet(rf,6);
355      poly b1 = 3*(rf3/(var(1)^3))*var(1)^2+2*(rf4/(var(1)^2*var(2)^2))
356                *var(1)+(rf5/(var(1)*var(2)^4));
357      ring R = 0,x,dp;
358      map phi = br, x, 1;
359      poly b11 = phi(b1);
360      int r = nrroots(b11);
361      if( r == 0 || r == 1)
362      {
363        setring(br);
364        signforJ10 = 1;
365      }
366      else
367      {
368        setring(br);
369        poly c1 = (rf3/(var(1)^3))*var(1)^3+(rf4/(var(1)^2*var(2)^2))*
370                  var(1)^2+(rf5/(var(1)*var(2)^4))*
371                  var(1)+rf6/(var(2)^6);
372        list Factors = factorize(c1);
373        if( size(Factors[1])>2)
374        {
375          if( deg(Factors[1][2]) == 1)
376          {
377            poly g1 = Factors[1][2];
378          }
379          else
380          {
381            poly g1 = Factors[1][3];
382          }
383          phi = br, ((g1/var(1))*var(1)-g1)/(g1/var(1)),
384                  var(2);
385          number b = number(phi(b1));
386          if(b > 0)
387          {
388            signforJ10 = 1;
389          }
390          else
391          {
392            signforJ10 = -1;
393          }
394        }
395        else
396        {
397          ring R = (complex,40,40),x,lp;
398          phi = br, x, 1;
399          poly c1 = phi(c1);
400          poly b1 = phi(b1);
401          list L = laguerre_solve(c1,30);
402          list LL;
403          for(i = 1; i <= size(L); i++)
404          {
405            if(impart(L[i]) == 0)
406            {
407              LL = insert(LL,L[i]);
408            }
409          }
410          number r1 = LL[1];
411          for(j = 1; j <= size(LL); j++)
412          {
413            r1 = round(r1*10000000000)/10000000000;
414            number c1min = number(subst(c1,x,r1-0.0000000001));
415            number c1plus = number(subst(c1,x,r1+0.0000000001));
416            number b1min = number(subst(b1,x,r1-0.00000000001));
417            number b1plus = number(subst(b1,x,r1+0.00000000001));
418            if(c1min*c1plus<0)
419            {
420              int c = -1;
421            }
422            if(c1min*c1plus>0)
423            {
424              int c = 1;
425            }
426            if(b1min>0 && b1plus>0)
427            {
428              int b = 1;
429            }
430            if(b1min<0 && b1plus<0)
431            {
432              int b = -1;
433            }
434            if(b1min*b1plus<=0)
435            {
436              int b = 0;
437            }
438            if( c < 0 && b != 0)
439            {
440              r1 = LL[j];
441              break;
442            }
443          }
444          setring(br);
445          if (c == -1 && b == 1)
446          {
447            signforJ10 = 1;
448          }
449          if (c == -1 && b == -1)
450          {
451            signforJ10 = -1;
452          }
453          if (c == 1 || b == 0)
454          {
455            ERROR("Ask Arnold the normal form.)");
456          }
457        }
458      }
459      if(signforJ10 == 1)
460      {
461        typeofsing = "J[10]+";
462      }
463      else
464      {
465        typeofsing = "J[10]-";
466      }
467      nf = var(1)^3+var(1)^2*var(2)^2+signforJ10*var(1)*var(2)^4;
468      monparam = var(1)^2*var(2)^2;
469      typeofsing_alternative = typeofsing;
470      nf_alternative = nf;
471    }
472    if(complextype[1,3] == "X[1") //case X[1,k]
473    {
474      if(mu > 9)
475      {
476        rf = jet(rf,4);
477
478        number s1 = number(rf/(var(1)^4));
479        number s2 = number(rf/(var(2)^4));
480        if(s2 != 0 && s1 == 0)
481        {
482          phi = br, var(2), var(1);
483          rf = phi(rf);
484        }
485        if(s2 == 0 && s1 == 0)
486        {
487          number t1 = number(rf/(var(1)^3*var(2)));
488          number t2 = number(rf/(var(1)^2*var(2)^2));
489          number t3 = number(rf/(var(1)*var(2)^3));
490          if(t1+t2+t3 == 0)
491          {
492            if(2*t1+4*t2+8*t3 != 0)
493            {
494              phi = br, var(1), 2*var(2);
495              rf = phi(rf);
496            }
497            else
498            {
499              phi = br, var(1), 3*var(2);
500              rf = phi(rf);
501            }
502          }
503          phi = br, var(1), var(1)+var(2);
504          rf = phi(rf);
505        }
506        ring R = 0, x, dp;
507        map phi = br, var(1), 1;
508        int k = nrroots(phi(rf));
509        setring(br);
510        if(k == 1)
511        {
512          number w = number(rf/(var(1)^4));
513          if(w > 0)
514          {
515            typeofsing = "X["+string(mu)+"]++";
516            nf = var(1)^4+var(1)^2*var(2)^2+var(2)^(4+(mu-9));
517            typeofsing_alternative = "X["+string(mu)+"]--";
518            nf_alternative = -var(1)^4-var(1)^2*var(2)^2+var(2)^(4+(mu-9));
519          }
520          else
521          {
522            typeofsing = "X["+string(mu)+"]--";
523            nf = -var(1)^4-var(1)^2*var(2)^2+var(2)^(4+(mu-9));
524            typeofsing_alternative = "X["+string(mu)+"]++";
525            nf_alternative = var(1)^4+var(1)^2*var(2)^2+var(2)^(4+(mu-9));
526          }
527        }
528        if(k == 3)
529        {
530          list Factors = factorize(rf);
531          for(i = 2; i <= size(Factors[1]); i++)
532          {
533            if(Factors[2][i] == 2)
534            {
535              poly g1 = Factors[1][i];
536              break;
537            }
538          }
539          map phi = br, (var(1)-(g1/(var(2))*var(2)))/(g1/var(1)), var(2);
540          rf = phi(rf);
541          number w = number(rf/(var(1)^2*var(2)^2));
542          if(w > 0)
543          {
544            typeofsing = "X["+string(mu)+"]-+";
545            nf = -var(1)^4+var(1)^2*var(2)^2+var(2)^(4+(mu-9));
546            typeofsing_alternative = "X["+string(mu)+"]+-";
547            nf_alternative = var(1)^4-var(1)^2*var(2)^2+var(2)^(4+(mu-9));
548          }
549          else
550          {
551            typeofsing = "X["+string(mu)+"]+-";
552            nf = var(1)^4-var(1)^2*var(2)^2+var(2)^(4+(mu-9));
553            typeofsing_alternative = "X["+string(mu)+"]-+";
554            nf_alternative = -var(1)^4+var(1)^2*var(2)^2+var(2)^(4+(mu-9));
555          }
556        }
557        monparam = var(2)^(4+(mu-9));
558      }
559    }
560    if(complextype == "E[12]")   // case E[12]
561    {
562      typeofsing = "E[12]";
563      nf = var(1)^3+var(2)^7+var(1)*var(2)^5;
564      monparam = var(1)*var(2)^5;
565      typeofsing_alternative = typeofsing;
566      nf_alternative = nf;
567    }
568    if(complextype == "E[13]")   // case E[13]
569    {
570      typeofsing = "E[13]";
571      nf = var(1)^3+var(1)*var(2)^5+var(2)^8;
572      monparam = var(2)^8;
573      typeofsing_alternative = typeofsing;
574      nf_alternative = nf;
575    }
576    if(complextype == "E[14]") //case E[14]
577    {
578      poly g = jet(rf,3);
579      number s = number(g/(var(1)^3));
580      if(s == 0)
581      {
582        phi = br, var(2), var(1);
583        rf = phi(rf);
584        g = jet(rf,3);
585        s = number(g/(var(1)^3));
586      }
587      rf = rf/s;
588      list Factors = factorize(g);
589      poly g1 = Factors[1][2];
590      phi = br, (var(1)-(g1/var(2))*var(2))/(g1/var(1)), var(2);
591      rf = phi(rf);
592      g = jet(rf,3);
593      number w0 = number(g/(var(1)^3));
594      phi = br, var(1)-((jet(rf,4)-(w0*var(1)^3))/(3*var(1)^2)), var(2);
595      rf = phi(rf);
596      phi = br, var(1)-((jet(rf,5)-(w0*var(1)^3))/(3*var(1)^2)), var(2);
597      rf = phi(rf);
598      rf = s*rf;
599      rf = jet(rf,8);
600      number w = number(rf/(var(2)^8));
601      if(w > 0)
602      {
603        typeofsing = "E[14]+";
604        nf = var(1)^3+var(2)^8+var(1)*var(2)^6;
605        typeofsing_alternative = "E[14]-";
606        nf_alternative = var(1)^3-var(2)^8+var(1)*var(2)^6;
607      }
608      if(w < 0)
609      {
610        typeofsing = "E[14]-";
611        nf = var(1)^3-var(2)^8+var(1)*var(2)^6;
612        typeofsing_alternative = "E[14]+";
613        nf_alternative = var(1)^3+var(2)^8+var(1)*var(2)^6;
614      }
615      monparam = var(1)*var(2)^6;
616    }
617    if(complextype == "Z[11]")   // case Z[11]
618    {
619      typeofsing = "Z[11]";
620      nf = var(1)^3*var(2)+var(2)^5+var(1)*var(2)^4;
621      monparam = var(1)*var(2)^4;
622      typeofsing_alternative = typeofsing;
623      nf_alternative = nf;
624    }
625    if(complextype == "Z[12]")   // case Z[12]
626    {
627      typeofsing = "Z[12]";
628      nf = var(1)^3*var(2)+var(1)*var(2)^4+var(1)^2*var(2)^3;
629      monparam = var(1)^2*var(2)^3;
630      typeofsing_alternative = typeofsing;
631      nf_alternative = nf;
632    }
633    if(complextype == "Z[13]")
634    {
635      poly g = jet(rf,4);
636      number s = number(g/var(1)^3*var(2));
637      if(s == 0)
638      {
639        phi = br, var(2), var(1);
640        rf = phi(rf);
641        g = jet(rf,4);
642      }
643      list Factors = factorize(g);
644      if(Factors[2][2] == 3)
645      {
646        poly g1 = Factors[1][2];
647      }
648      else
649      {
650        poly g1 = Factors[1][3];
651      }
652      phi = br, var(1)-(g1/var(2))*var(2), var(2);
653      rf = phi(rf);
654      rf = jet(rf,6);
655      number w = number(rf/var(2)^6);
656      if(w > 0)
657      {
658        typeofsing = "Z[13]+";
659        nf = var(1)^3*var(2)+var(2)^6+var(1)*var(2)^5;
660        typeofsing_alternative = "Z[13]-";
661        nf_alternative = var(1)^3*var(2)-var(2)^6+var(1)*var(2)^5;
662      }
663      else
664      {
665        typeofsing = "Z[13]-";
666        nf = var(1)^3*var(2)-var(2)^6+var(1)*var(2)^5;
667        typeofsing_alternative = "Z[13]+";
668        nf_alternative = var(1)^3*var(2)+var(2)^6+var(1)*var(2)^5;
669      }
670      monparam = var(1)*var(2)^5;
671    }
672    if(complextype == "W[12]") //case W[12]
673    {
674      poly g = jet(rf, 4);
675      number s = number(g/(var(1)^4));
676      if(s == 0)
677      {
678        s = number(g/(var(2)^4));
679        phi = br, var(2), var(1);   // maybe we'll need this transformation
680        rf = phi(rf);               // later
681      }
682      if(s > 0)
683      {
684        typeofsing = "W[12]+";
685        nf = var(1)^4+var(2)^5+var(1)^2*var(2)^3;
686        typeofsing_alternative = "W[12]-";
687        nf_alternative = -var(1)^4+var(2)^5+var(1)^2*var(2)^3;
688      }
689      else
690      {
691        typeofsing = "W[12]-";
692        nf = -var(1)^4+var(2)^5+var(1)^2*var(2)^3;
693        typeofsing_alternative = "W[12]+";
694        nf_alternative = var(1)^4+var(2)^5+var(1)^2*var(2)^3;
695      }
696      monparam = var(1)^2*var(2)^3;
697    }
698    if(complextype == "W[13]")   //case W[13]
699    {
700      poly g = jet(rf, 4);
701      number s = number(g/(var(1)^4));
702      if(s == 0)
703      {
704        s = number(g/(var(2)^4));
705        phi = br, var(2), var(1);   // maybe we'll need this transformation
706        rf = phi(rf);               // later
707      }
708      if(s > 0)
709      {
710        typeofsing = "W[13]+";
711        nf = var(1)^4+var(1)*var(2)^4+var(2)^6;
712        typeofsing_alternative = "W[13]-";
713        nf_alternative = -var(1)^4+var(1)*var(2)^4+var(2)^6;
714      }
715      else
716      {
717        typeofsing = "W[13]-";
718        nf = -var(1)^4+var(1)*var(2)^4+var(2)^6;
719        typeofsing_alternative = "W[13]+";
720        nf_alternative = var(1)^4+var(1)*var(2)^4+var(2)^6;
721      }
722      monparam = var(2)^6;
723    }
724    if(typeofsing == "")
725    {
726      ERROR("This case is not yet implemented.");
727    }
728  }
729  if(cr > 2)
730  {
731    ERROR("This case is not yet implemented.");
732  }
733
734  /* add the non-corank variables to the normal forms */
735  nf = addnondegeneratevariables(nf, lambda, cr);
736  nf_alternative = addnondegeneratevariables(nf_alternative, n-cr-lambda, cr);
737
738  /* write normal form as a string in the cases with modality greater than 0 */
739  if(monparam != 0)
740  {
741    poly nf_tmp = nf;
742    poly nf_alternative_tmp = nf_alternative;
743    def nf = modality1NF(nf_tmp, monparam);
744    def nf_alternative = modality1NF(nf_alternative_tmp, monparam);
745  }
746
747  /* write comments */
748  if(printcomments)
749  {
750    string comments = newline;
751    comments = comments+"Type of singularity: "   +typeofsing    +newline
752                       +"Normal form:         "   +string(nf)    +newline
753                       +"Corank:              "   +string(cr)    +newline
754                       +"Milnor number:       "   +string(mu)    +newline
755                       +"Inertia index:       "   +string(lambda)+newline
756                       +"Determinacy:         <= "+string(d)     +newline;
757    if(typeofsing != typeofsing_alternative || nf != nf_alternative
758       || lambda != n-cr-lambda)
759    {
760      comments = comments+newline
761        +"Note: By multiplying the input polynomial with -1,"+newline
762        +"      it can also be regarded as of the following case:"+newline
763        +"Type of singularity: "+typeofsing_alternative+newline
764        +"Normal form:         "+string(nf_alternative)+newline
765        +"Inertia index:       "+string(n-cr-lambda)   +newline;
766    }
767    if(morecomments != newline)
768    {
769      comments = comments+morecomments;
770    }
771  }
772
773  /* return results */
774  if(printcomments)
775  {
776    return(list(typeofsing, nf, cr, mu, lambda, d, comments));
777  }
778  else
779  {
780    return(list(typeofsing, nf, cr, mu, lambda, d));
781  }
782}
783example
784{
785  "EXAMPLE:";
786  echo = 2;
787  ring r = 0, (x,y,z), ds;
788  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
789  realclassify(f, 1);
790}
791
792///////////////////////////////////////////////////////////////////////////////
793/*
794  print the normal form as a string for the modality 1 cases.
795  The first argument is the normalform with parameter = 1,
796  the second argument is the monomial whose coefficient is the parameter.
797*/
798static proc modality1NF(poly nf, poly monparam)
799{
800  def br = basering;
801  list lbr = ringlist(br);
802  ring r = (0,a), x, dp;
803  list lr = ringlist(r);
804  setring(br);
805  list lr = fetch(r, lr);
806  lbr[1] = lr[1];
807  def s = ring(lbr);
808  setring(s);
809  poly nf = fetch(br, nf);
810  poly monparam = fetch(br, monparam);
811  nf = nf+(a-1)*monparam;
812  string result = string(nf);
813  setring(br);
814  return(result);
815}
816
817///////////////////////////////////////////////////////////////////////////////
818/*
819  add squares of the non-degenerate variables (i.e. var(cr+1), ...,
820  var(nvars(basering)) for corank cr) to the normalform nf,
821  with signs according to the inertia index lambda
822*/
823static proc addnondegeneratevariables(poly nf, int lambda, int cr)
824{
825  int n = nvars(basering);
826  int i;
827  for(i = cr+1; i <= n-lambda; i++)
828  {
829    nf = nf+var(i)^2;
830  }
831  for(i = n-lambda+1; i <= n ; i++)
832  {
833    nf = nf-var(i)^2;
834  }
835  return(nf);
836}
837
838///////////////////////////////////////////////////////////////////////////////
839proc realmorsesplit(poly f, list #)
840"
841USAGE:    realmorsesplit(f[, mu]); f poly, mu int
842RETURN:   a list consisting of the corank of f, the inertia index, an upper
843          bound for the determinacy, the residual form of f and
844          the transformation
845NOTE:     The characteristic of the basering must be zero, the monomial order
846          must be local, f must be contained in maxideal(2) and the Milnor
847          number of f must be finite.
848          @* The Milnor number of f can be provided as an optional parameter in
849             order to avoid that it is computed again.
850SEE ALSO: morsesplit
851KEYWORDS: Morse lemma; Splitting lemma
852EXAMPLE:  example morsesplit; shows an example"
853{
854  /* auxiliary variables */
855  int i, j;
856
857  /* error check */
858  if(char(basering) != 0)
859  {
860    ERROR("The characteristic must be zero.");
861  }
862  int n = nvars(basering);
863  for(i = 1; i <= n; i++)
864  {
865    if(var(i) > 1)
866    {
867      ERROR("The monomial order must be local.");
868    }
869  }
870  if(jet(f, 1) != 0)
871  {
872    ERROR("The input polynomial must be contained in maxideal(2).");
873  }
874
875  /* get Milnor number before continuing error check */
876  int mu;
877  if(size(#) > 0)   // read optional parameter
878  {
879    if(size(#) > 1 || typeof(#[1]) != "int")
880    {
881      ERROR("Wrong optional parameters.");
882    }
883    else
884    {
885      mu = #[1];
886    }
887  }
888  else              // compute Milnor number
889  {
890    mu = milnornumber(f);
891  }
892
893  /* continue error check */
894  if(mu < 0)
895  {
896    ERROR("The Milnor number of the input polynomial must be"+newline
897      +"non-negative and finite.");
898  }
899
900  /* preliminary stuff */
901  list S;
902  int k = determinacy(f, mu);
903  f = jet(f, k);
904  def br = basering;
905  map Phi = br, maxideal(1);
906  map phi;
907  poly a, p, r;
908
909  /* treat the variables one by one */
910  for(i = 1; i <= n; i++)
911  {
912    if(jet(f, 2)/var(i) == 0)
913    {
914      S = insert(S, i);
915    }
916    else
917    {
918      f, a, p, r = rewriteformorsesplit(f, k, i);
919      if(jet(a, 0) == 0)
920      {
921        for(j = i+1; j <= n; j++)
922        {
923          if(jet(f, 2)/(var(i)*var(j)) != 0)
924          {
925            break;
926          }
927        }
928        phi = br, maxideal(1);
929        phi[j] = var(j)+var(i);
930        Phi = phi(Phi);
931        f = phi(f);
932      }
933      f, a, p, r = rewriteformorsesplit(f, k, i);
934      while(p != 0)
935      {
936        phi = br, maxideal(1);
937        phi[i] = var(i)-p/(2*jet(a, 0));
938        Phi = phi(Phi);
939        f = phi(f);
940        f, a, p, r = rewriteformorsesplit(f, k, i);
941      }
942    }
943  }
944
945  /* sort variables according to corank */
946  int cr = size(S);
947  phi = br, 0:n;
948  j = 1;
949  for(i = size(S); i > 0; i--)
950  {
951    phi[S[i]] = var(j);
952    j++;
953  }
954  for(i = 1; i <= n; i++)
955  {
956    if(phi[i] == 0)
957    {
958      phi[i] = var(j);
959      j++;
960    }
961  }
962  Phi = phi(Phi);
963  f = phi(f);
964
965  /* compute the inertia index lambda */
966  int lambda;
967  list negCoeff, posCoeff;
968  number ai;
969  poly f2 = jet(f, 2);
970  for(i = 1; i <= n; i++)
971  {
972    ai = number(f2/var(i)^2);
973    if(ai < 0)
974    {
975      lambda++;
976      negCoeff = insert(negCoeff, i);
977    }
978    if(ai > 0)
979    {
980      posCoeff = insert(posCoeff, i);
981    }
982  }
983
984  /* sort variables according to lambda */
985  phi = br, maxideal(1);
986  j = cr+1;
987  for(i = size(negCoeff); i > 0; i--)
988  {
989    phi[negCoeff[i]] = var(j);
990    j++;
991  }
992  for(i = size(posCoeff); i > 0; i--)
993  {
994    phi[posCoeff[i]] = var(j);
995    j++;
996  }
997  Phi = phi(Phi);
998  f = phi(f);
999
1000  /* compute residual form */
1001  phi = br, maxideal(1);
1002  for(i = size(S)+1; i <= n; i++)
1003  {
1004    phi[i] = 0;
1005  }
1006  f = phi(f);
1007
1008  return(list(cr, lambda, k, f, Phi));
1009}
1010example
1011{
1012  "EXAMPLE:";
1013  echo = 2;
1014  ring r = 0, (x,y,z), ds;
1015  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
1016  realmorsesplit(f);
1017}
1018
1019///////////////////////////////////////////////////////////////////////////////
1020/*
1021   - apply jet(f, k)
1022   - rewrite f as f = a*var(i)^2+p*var(i)+r with
1023     var(i)-free p and r
1024*/
1025static proc rewriteformorsesplit(poly f, int k, int i)
1026{
1027  f = jet(f, k);
1028  matrix C = coeffs(f, var(i));
1029  poly r = C[1,1];
1030  poly p = C[2,1];
1031  poly a = (f-r-p*var(i))/var(i)^2;
1032  return(f, a, p, r);
1033}
1034
1035///////////////////////////////////////////////////////////////////////////////
1036proc milnornumber(poly f)
1037"
1038USAGE:    milnornumber(f); f poly
1039RETURN:   Milnor number of f, or -1 if the Milnor number is not finite
1040KEYWORDS: Milnor number
1041NOTE:     The monomial order must be local.
1042EXAMPLE:  example milnornumber; shows an example"
1043{
1044  /* error check */
1045  int i;
1046  for(i = nvars(basering); i > 0; i--)
1047  {
1048    if(var(i) > 1)
1049    {
1050      ERROR("The monomial order must be local.");
1051    }
1052  }
1053
1054  return(vdim(std(jacob(f))));
1055}
1056example
1057{
1058  "EXAMPLE:";
1059  echo = 2;
1060  ring r = 0, (x,y), ds;
1061  poly f = x3+y4;
1062  milnornumber(f);
1063}
1064
1065///////////////////////////////////////////////////////////////////////////////
1066proc determinacy(poly f, list #)
1067"
1068USAGE:    determinacy(f[, mu]); f poly, mu int
1069RETURN:   an upper bound for the determinacy of f
1070NOTE:     The characteristic of the basering must be zero, the monomial order
1071          must be local, f must be contained in maxideal(1) and the Milnor
1072          number of f must be finite.
1073          @* The Milnor number of f can be provided as an optional parameter in
1074             order to avoid that it is computed again.
1075SEE ALSO: milnornumber, highcorner
1076KEYWORDS: Determinacy
1077EXAMPLE:  example determinacy; shows an example"
1078{
1079  /* auxiliary variables */
1080  int i;
1081
1082  /* error check */
1083  if(char(basering) != 0)
1084  {
1085    ERROR("The characteristic must be zero.");
1086  }
1087  int n = nvars(basering);
1088  for(i = 1; i <= n; i++)
1089  {
1090    if(var(i) > 1)
1091    {
1092      ERROR("The monomial order must be local.");
1093    }
1094  }
1095  if(jet(f, 0) != 0)
1096  {
1097    ERROR("The input polynomial must be contained in maxideal(1).");
1098  }
1099
1100  /* get Milnor number before continuing error check */
1101  int mu;
1102  if(size(#) > 0)   // read optional parameter
1103  {
1104    if(size(#) > 1 || typeof(#[1]) != "int")
1105    {
1106      ERROR("Wrong optional parameters.");
1107    }
1108    else
1109    {
1110      mu = #[1];
1111    }
1112  }
1113  else              // compute Milnor number
1114  {
1115    mu = milnornumber(f);
1116  }
1117
1118  /* continue error check */
1119  if(mu < 0)
1120  {
1121    ERROR("The Milnor number of the input polynomial must be"+newline
1122      +"non-negative and finite.");
1123  }
1124
1125  int k;   // an upper bound for the determinacy,
1126           // we use several methods:
1127
1128  /* Milnor number */
1129  k = mu+1;
1130  f = jet(f, k);
1131
1132  /* highest corner */
1133  int hc;
1134  for(i = 0; i < 3; i++)
1135  {
1136    f = jet(f, k);
1137    hc = deg(highcorner(std(maxideal(i)*jacob(f))));
1138    hc = hc+2-i;
1139    if(hc < k)
1140    {
1141      k = hc;
1142    }
1143  }
1144
1145  return(k);
1146}
1147example
1148{
1149  "EXAMPLE:";
1150  echo = 2;
1151  ring r = 0, (x,y), ds;
1152  poly f = x3+xy3;
1153  determinacy(f);
1154}
Note: See TracBrowser for help on using the repository browser.