source: git/Singular/LIB/realclassify.lib @ 7d24a2

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