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

spielwiese
Last change on this file since 6942859 was 6942859, checked in by Andreas Steenpass <steenpass@…>, 10 years ago
chg: update realclassify.lib (cherry picked from commit 9eb2f474dffc7af6a2fe0ffdd8d65a80f37027c8) Signed-off-by: Andreas Steenpass <steenpass@mathematik.uni-kl.de>
  • Property mode set to 100644
File size: 18.5 KB
Line 
1////////////////////////////////////////////////////////////////////////////
2version="version realclassify.lib 4.0.0.0 Jun_2013 "; // $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 "elim.lib";
31LIB "primdec.lib";
32LIB "classify.lib";
33LIB "rootsur.lib";
34LIB "rootsmr.lib";
35LIB "atkins.lib";
36LIB "solve.lib";
37///////////////////////////////////////////////////////////////////////////////
38proc realclassify(poly f, list #)
39"
40USAGE:    realclassify(f[, format]); f poly, format string
41RETURN:   A list containing (in this order)
42          @* - the type of the singularity as a string,
43          @* - the normal form,
44          @* - the corank, the Milnor number, the inertia index and
45               a bound for the determinacy as integers.
46          @* The normal form involves parameters for singularities of modality
47             greater than 0. The actual value of the parameters is not computed
48             in most of the cases. If the value of the parameter is unknown,
49             the normal form is given as a string with an \"a\" as the
50             parameter. Otherwise, it is given as a polynomial.
51          @* An optional string @code{format} can be provided. Its default
52             value is \"short\" in which case the return value is the list
53             described above. If set to \"nice\", a string is added at the end
54             of this list, containing the result in a more readable form.
55NOTE:     The classification is done over the real numbers, so in contrast to
56          classify.lib, the signs of coefficients of monomials where even
57          exponents occur matter.
58          @* The ground field must be Q (the rational numbers). No field
59             extensions of any kind nor floating point numbers are allowed.
60          @* The monomial order must be local.
61          @* The input polynomial must be contained in maxideal(2) and must be
62             an isolated singularity of modality 0 or 1. The Milnor number is
63             checked for being finite.
64SEE ALSO: classify
65KEYWORDS: Classification of singularities
66EXAMPLE:  example realclassify; shows an example"
67{
68  /* auxiliary variables */
69  int i, j;
70
71  /* name for the basering */
72  def br = basering;
73
74  /* read optional parameters */
75  int printcomments;
76  if(size(#) > 0)
77  {
78    if(size(#) > 1 || typeof(#[1]) != "string")
79    {
80      ERROR("Wrong optional parameters.");
81    }
82    if(#[1] != "short" && #[1] != "nice")
83    {
84      ERROR("Wrong optional parameters.");
85    }
86    if(#[1] == "nice")
87    {
88      printcomments = 1;
89    }
90  }
91
92  /* error check */
93  if(charstr(br) != "0")
94  {
95    ERROR("The ground field must be Q (the rational numbers).");
96  }
97  int n = nvars(br);
98  for(i = 1; i <= n; i++)
99  {
100    if(var(i) > 1)
101    {
102      ERROR("The monomial order must be local.");
103    }
104  }
105  if(jet(f, 1) != 0)
106  {
107    ERROR("The input polynomial must be contained in maxideal(2).");
108  }
109
110  /* compute Milnor number before continuing the error check */
111  int mu = milnornumber(f);
112
113  /* continue error check */
114  if(mu < 1)
115  {
116    ERROR("The Milnor number of the input polynomial must be"+newline
117      +"positive and finite.");
118  }
119
120  /* call classify before continuing the error check */
121  list dataFromClassify = prepRealclassify(f);
122  int m = dataFromClassify[1];                // the modality of f
123  string complextype = dataFromClassify[2];   // the complex type of f
124
125  /* continue error check */
126  if(m > 1)
127  {
128    ERROR("The input polynomial must be a singularity of modality 0 or 1.");
129  }
130
131  /* apply splitting lemma */
132  list morse = realmorsesplit(f, mu);
133  int cr = morse[1];
134  int lambda = morse[2];
135  int d = morse[3];
136  poly rf = morse[4];
137
138  /* determine the type */
139  string typeofsing;
140  poly nf;
141  poly monparam;   // the monomial whose coefficient is the parameter
142                   // in the modality 1 cases, 0 otherwise
143  string morecomments = newline;
144
145  if(cr == 0)   // case A[1]
146  {
147    typeofsing, nf = caseA1(rf, lambda, n);
148  }
149  if(cr == 1)   // case A[k], k > 1
150  {
151    typeofsing, nf = caseAk(rf, n);
152  }
153  if(cr == 2)
154  {
155    if(complextype[1,2] == "D[")   // case D[k]
156    {
157      if(mu == 4)   // case D[4]
158      {
159        typeofsing, nf = caseD4(rf);
160      }
161      else   // case D[k], k > 4
162      {
163        typeofsing, nf = caseDk(rf, mu);
164      }
165    }
166    if(complextype == "E[6]")   // case E[6]
167    {
168      typeofsing, nf = caseE6(rf);
169    }
170    if(complextype == "E[7]")   // case E[7]
171    {
172      typeofsing, nf = caseE7();
173    }
174    if(complextype == "E[8]")   // case E[8]
175    {
176      typeofsing, nf = caseE8();
177    }
178    if(typeofsing == "")
179    {
180      ERROR("This case is not yet implemented.");
181    }
182  }
183  if(cr > 2)
184  {
185    ERROR("This case is not yet implemented.");
186  }
187
188  /* add the non-corank variables to the normal forms */
189  nf = addnondegeneratevariables(nf, lambda, cr);
190
191  /* write normal form as a string in the cases with modality greater than 0 */
192  if(monparam != 0)
193  {
194    poly nf_tmp = nf;
195    kill nf;
196    def nf = modality1NF(nf_tmp, monparam);
197  }
198
199  /* write comments */
200  if(printcomments)
201  {
202    string comments = newline;
203    comments = comments+"Type of singularity: "   +typeofsing    +newline
204                       +"Normal form:         "   +string(nf)    +newline
205                       +"Corank:              "   +string(cr)    +newline
206                       +"Milnor number:       "   +string(mu)    +newline
207                       +"Inertia index:       "   +string(lambda)+newline
208                       +"Determinacy:         <= "+string(d)     +newline;
209    if(morecomments != newline)
210    {
211      comments = comments+morecomments;
212    }
213  }
214
215  /* return results */
216  if(printcomments)
217  {
218    return(list(typeofsing, nf, cr, mu, lambda, d, comments));
219  }
220  else
221  {
222    return(list(typeofsing, nf, cr, mu, lambda, d));
223  }
224}
225example
226{
227  "EXAMPLE:";
228  echo = 2;
229  ring r = 0, (x,y,z), ds;
230  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
231  realclassify(f, "nice");
232}
233
234///////////////////////////////////////////////////////////////////////////////
235static proc caseA1(poly rf, int lambda, int n)
236{
237  string typeofsing = "A[1]";
238  poly nf = 0;
239  return(typeofsing, nf);
240}
241
242///////////////////////////////////////////////////////////////////////////////
243static proc caseAk(poly rf, int n)
244{
245  /* preliminaries */
246  string typeofsing;
247  poly nf;
248
249  int k = deg(lead(rf), 1:n)-1;
250  if(k%2 == 0)
251  {
252    nf = var(1)^(k+1);
253    typeofsing = "A["+string(k)+"]";
254  }
255  else
256  {
257    if(leadcoef(rf) > 0)
258    {
259      nf = var(1)^(k+1);
260      typeofsing = "A["+string(k)+"]+";
261    }
262    else
263    {
264      nf = -var(1)^(k+1);
265      typeofsing = "A["+string(k)+"]-";
266    }
267  }
268  return(typeofsing, nf);
269}
270
271///////////////////////////////////////////////////////////////////////////////
272static proc caseD4(poly rf)
273{
274  /* preliminaries */
275  string typeofsing;
276  poly nf;
277  def br = basering;
278  map phi;
279
280  rf = jet(rf, 3);
281  number s1 = number(rf/(var(1)^3));
282  number s2 = number(rf/(var(2)^3));
283  if(s2 == 0 && s1 != 0)
284  {
285    phi = br, var(2), var(1);
286    rf = phi(rf);
287  }
288  if(s1 == 0 && s2 == 0)
289  {
290    number t1 = number(rf/(var(1)^2*var(2)));
291    number t2 = number(rf/(var(2)^2*var(1)));
292    if(t1+t2 == 0)
293    {
294      phi = br, var(1)+2*var(2), var(2);
295      rf = phi(rf);
296    }
297    else
298    {
299      phi = br, var(1)+var(2), var(2);
300      rf = phi(rf);
301    }
302  }
303  ring R = 0, y, dp;
304  map phi = br, 1, y;
305  poly rf = phi(rf);
306  int k = nrroots(rf);
307  setring(br);
308  if(k == 3)
309  {
310    nf = var(1)^2*var(2)-var(2)^3;
311    typeofsing = "D[4]-";
312  }
313  else
314  {
315    nf = var(1)^2*var(2)+var(2)^3;
316    typeofsing = "D[4]+";
317  }
318  return(typeofsing, nf);
319}
320
321///////////////////////////////////////////////////////////////////////////////
322static proc caseDk(poly rf, int mu)
323{
324  /* preliminaries */
325  string typeofsing;
326  poly nf;
327  def br = basering;
328  map phi;
329
330  rf = jet(rf, mu-1);
331  list factorization = factorize(jet(rf, 3));
332  list factors = factorization[1][2];
333  if(factorization[2][2] == 2)
334  {
335    factors = insert(factors, factorization[1][3], 1);
336  }
337  else
338  {
339    factors = insert(factors, factorization[1][3]);
340  }
341  factors[2] = factorization[1][1]*factors[2];
342  matrix T[2][2] = factors[1]/var(1), factors[1]/var(2),
343         factors[2]/var(1), factors[2]/var(2);
344  phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
345  rf = phi(rf);
346  rf = jet(rf, mu-1);
347  poly g;
348  int i;
349  for(i = 4; i < mu; i++)
350  {
351    g = jet(rf, i) - var(1)^2*var(2);
352    if(g != 0)
353    {
354      phi = br, var(1)-(g/(var(1)*var(2)))/2,
355          var(2)-(g/var(1)^i)*var(1)^(i-2);
356      rf = phi(rf);
357      rf = jet(rf, mu-1);
358    }
359  }
360  number a = number(rf/var(2)^(mu-1));
361  if(a > 0)
362  {
363    typeofsing = "D["+string(mu)+"]+";
364    nf = var(1)^2*var(2)+var(2)^(mu-1);
365  }
366  else
367  {
368    typeofsing = "D["+string(mu)+"]-";
369    nf = var(1)^2*var(2)-var(2)^(mu-1);
370  }
371  return(typeofsing, nf);
372}
373
374///////////////////////////////////////////////////////////////////////////////
375static proc caseE6(poly rf)
376{
377  /* preliminaries */
378  string typeofsing;
379  poly nf;
380  def br = basering;
381  map phi;
382
383  poly g = jet(rf,3);
384  number s = number(g/(var(1)^3));
385  if(s == 0)
386  {
387    phi = br, var(2), var(1);
388    rf = phi(rf);
389    g = jet(rf,3);
390  }
391  list Factors = factorize(g);
392  poly g1 = Factors[1][2];
393  phi = br, (var(1)-(g1/var(2))*var(2))/(g1/var(1)), var(2);
394  rf = phi(rf);
395  rf = jet(rf,4);
396  number w = number(rf/(var(2)^4));
397  if(w > 0)
398  {
399    typeofsing = "E[6]+";
400    nf = var(1)^3+var(2)^4;
401  }
402  else
403  {
404    typeofsing = "E[6]-";
405    nf = var(1)^3-var(2)^4;
406  }
407  return(typeofsing, nf);
408}
409
410///////////////////////////////////////////////////////////////////////////////
411static proc caseE7()
412{
413  string typeofsing = "E[7]";
414  poly nf = var(1)^3+var(1)*var(2)^3;
415  return(typeofsing, nf);
416}
417
418///////////////////////////////////////////////////////////////////////////////
419static proc caseE8()
420{
421  string typeofsing = "E[8]";
422  poly nf = var(1)^3+var(2)^5;
423  return(typeofsing, nf);
424}
425
426///////////////////////////////////////////////////////////////////////////////
427/*
428  print the normal form as a string for the modality 1 cases.
429  The first argument is the normalform with parameter = 1,
430  the second argument is the monomial whose coefficient is the parameter.
431*/
432static proc modality1NF(poly nf, poly monparam)
433{
434  def br = basering;
435  list lbr = ringlist(br);
436  ring r = (0,a), x, dp;
437  list lr = ringlist(r);
438  setring(br);
439  list lr = fetch(r, lr);
440  lbr[1] = lr[1];
441  def s = ring(lbr);
442  setring(s);
443  poly nf = fetch(br, nf);
444  poly monparam = fetch(br, monparam);
445  nf = nf+(a-1)*monparam;
446  string result = string(nf);
447  setring(br);
448  return(result);
449}
450
451///////////////////////////////////////////////////////////////////////////////
452/*
453  add squares of the non-degenerate variables (i.e. var(cr+1), ...,
454  var(nvars(basering)) for corank cr) to the normalform nf,
455  with signs according to the inertia index lambda
456*/
457static proc addnondegeneratevariables(poly nf, int lambda, int cr)
458{
459  int n = nvars(basering);
460  int i;
461  for(i = cr+1; i <= n-lambda; i++)
462  {
463    nf = nf+var(i)^2;
464  }
465  for(i = n-lambda+1; i <= n ; i++)
466  {
467    nf = nf-var(i)^2;
468  }
469  return(nf);
470}
471
472///////////////////////////////////////////////////////////////////////////////
473proc realmorsesplit(poly f, list #)
474"
475USAGE:    realmorsesplit(f[, mu]); f poly, mu int
476RETURN:   a list consisting of the corank of f, the inertia index, an upper
477          bound for the determinacy, the residual form of f and
478          the transformation
479NOTE:     The characteristic of the basering must be zero, the monomial order
480          must be local, f must be contained in maxideal(2) and the Milnor
481          number of f must be finite.
482          @* The Milnor number of f can be provided as an optional parameter in
483             order to avoid that it is computed again.
484SEE ALSO: morsesplit
485KEYWORDS: Morse lemma; Splitting lemma
486EXAMPLE:  example morsesplit; shows an example"
487{
488  /* auxiliary variables */
489  int i, j;
490
491  /* error check */
492  if(char(basering) != 0)
493  {
494    ERROR("The characteristic must be zero.");
495  }
496  int n = nvars(basering);
497  for(i = 1; i <= n; i++)
498  {
499    if(var(i) > 1)
500    {
501      ERROR("The monomial order must be local.");
502    }
503  }
504  if(jet(f, 1) != 0)
505  {
506    ERROR("The input polynomial must be contained in maxideal(2).");
507  }
508
509  /* get Milnor number before continuing error check */
510  int mu;
511  if(size(#) > 0)   // read optional parameter
512  {
513    if(size(#) > 1 || typeof(#[1]) != "int")
514    {
515      ERROR("Wrong optional parameters.");
516    }
517    else
518    {
519      mu = #[1];
520    }
521  }
522  else              // compute Milnor number
523  {
524    mu = milnornumber(f);
525  }
526
527  /* continue error check */
528  if(mu < 0)
529  {
530    ERROR("The Milnor number of the input polynomial must be"+newline
531      +"non-negative and finite.");
532  }
533
534  /* preliminary stuff */
535  list S;
536  int k = determinacy(f, mu);
537  f = jet(f, k);
538  def br = basering;
539  map Phi = br, maxideal(1);
540  map phi;
541  poly a, p, r;
542
543  /* treat the variables one by one */
544  for(i = 1; i <= n; i++)
545  {
546    if(jet(f, 2)/var(i) == 0)
547    {
548      S = insert(S, i);
549    }
550    else
551    {
552      f, a, p, r = rewriteformorsesplit(f, k, i);
553      if(jet(a, 0) == 0)
554      {
555        for(j = i+1; j <= n; j++)
556        {
557          if(jet(f, 2)/(var(i)*var(j)) != 0)
558          {
559            break;
560          }
561        }
562        phi = br, maxideal(1);
563        phi[j] = var(j)+var(i);
564        Phi = phi(Phi);
565        f = phi(f);
566      }
567      f, a, p, r = rewriteformorsesplit(f, k, i);
568      while(p != 0)
569      {
570        phi = br, maxideal(1);
571        phi[i] = var(i)-p/(2*jet(a, 0));
572        Phi = phi(Phi);
573        f = phi(f);
574        f, a, p, r = rewriteformorsesplit(f, k, i);
575      }
576    }
577  }
578
579  /* sort variables according to corank */
580  int cr = size(S);
581  phi = br, 0:n;
582  j = 1;
583  for(i = size(S); i > 0; i--)
584  {
585    phi[S[i]] = var(j);
586    j++;
587  }
588  for(i = 1; i <= n; i++)
589  {
590    if(phi[i] == 0)
591    {
592      phi[i] = var(j);
593      j++;
594    }
595  }
596  Phi = phi(Phi);
597  f = phi(f);
598
599  /* compute the inertia index lambda */
600  int lambda;
601  list negCoeff, posCoeff;
602  number ai;
603  poly f2 = jet(f, 2);
604  for(i = 1; i <= n; i++)
605  {
606    ai = number(f2/var(i)^2);
607    if(ai < 0)
608    {
609      lambda++;
610      negCoeff = insert(negCoeff, i);
611    }
612    if(ai > 0)
613    {
614      posCoeff = insert(posCoeff, i);
615    }
616  }
617
618  /* sort variables according to lambda */
619  phi = br, maxideal(1);
620  j = cr+1;
621  for(i = size(negCoeff); i > 0; i--)
622  {
623    phi[negCoeff[i]] = var(j);
624    j++;
625  }
626  for(i = size(posCoeff); i > 0; i--)
627  {
628    phi[posCoeff[i]] = var(j);
629    j++;
630  }
631  Phi = phi(Phi);
632  f = phi(f);
633
634  /* compute residual form */
635  phi = br, maxideal(1);
636  for(i = size(S)+1; i <= n; i++)
637  {
638    phi[i] = 0;
639  }
640  f = phi(f);
641
642  return(list(cr, lambda, k, f, Phi));
643}
644example
645{
646  "EXAMPLE:";
647  echo = 2;
648  ring r = 0, (x,y,z), ds;
649  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
650  realmorsesplit(f);
651}
652
653///////////////////////////////////////////////////////////////////////////////
654/*
655   - apply jet(f, k)
656   - rewrite f as f = a*var(i)^2+p*var(i)+r with
657     var(i)-free p and r
658*/
659static proc rewriteformorsesplit(poly f, int k, int i)
660{
661  f = jet(f, k);
662  matrix C = coeffs(f, var(i));
663  poly r = C[1,1];
664  poly p = C[2,1];
665  poly a = (f-r-p*var(i))/var(i)^2;
666  return(f, a, p, r);
667}
668
669///////////////////////////////////////////////////////////////////////////////
670proc milnornumber(poly f)
671"
672USAGE:    milnornumber(f); f poly
673RETURN:   Milnor number of f, or -1 if the Milnor number is not finite
674KEYWORDS: Milnor number
675NOTE:     The monomial order must be local.
676EXAMPLE:  example milnornumber; shows an example"
677{
678  /* error check */
679  int i;
680  for(i = nvars(basering); i > 0; i--)
681  {
682    if(var(i) > 1)
683    {
684      ERROR("The monomial order must be local.");
685    }
686  }
687
688  return(vdim(std(jacob(f))));
689}
690example
691{
692  "EXAMPLE:";
693  echo = 2;
694  ring r = 0, (x,y), ds;
695  poly f = x3+y4;
696  milnornumber(f);
697}
698
699///////////////////////////////////////////////////////////////////////////////
700proc determinacy(poly f, list #)
701"
702USAGE:    determinacy(f[, mu]); f poly, mu int
703RETURN:   an upper bound for the determinacy of f
704NOTE:     The characteristic of the basering must be zero, the monomial order
705          must be local, f must be contained in maxideal(1) and the Milnor
706          number of f must be finite.
707          @* The Milnor number of f can be provided as an optional parameter in
708             order to avoid that it is computed again.
709SEE ALSO: milnornumber, highcorner
710KEYWORDS: Determinacy
711EXAMPLE:  example determinacy; shows an example"
712{
713  /* auxiliary variables */
714  int i;
715
716  /* error check */
717  if(char(basering) != 0)
718  {
719    ERROR("The characteristic must be zero.");
720  }
721  int n = nvars(basering);
722  for(i = 1; i <= n; i++)
723  {
724    if(var(i) > 1)
725    {
726      ERROR("The monomial order must be local.");
727    }
728  }
729  if(jet(f, 0) != 0)
730  {
731    ERROR("The input polynomial must be contained in maxideal(1).");
732  }
733
734  /* get Milnor number before continuing error check */
735  int mu;
736  if(size(#) > 0)   // read optional parameter
737  {
738    if(size(#) > 1 || typeof(#[1]) != "int")
739    {
740      ERROR("Wrong optional parameters.");
741    }
742    else
743    {
744      mu = #[1];
745    }
746  }
747  else              // compute Milnor number
748  {
749    mu = milnornumber(f);
750  }
751
752  /* continue error check */
753  if(mu < 0)
754  {
755    ERROR("The Milnor number of the input polynomial must be"+newline
756      +"non-negative and finite.");
757  }
758
759  int k;   // an upper bound for the determinacy,
760           // we use several methods:
761
762  /* Milnor number */
763  k = mu+1;
764  f = jet(f, k);
765
766  /* highest corner */
767  int hc;
768  for(i = 0; i < 3; i++)
769  {
770    f = jet(f, k);
771    hc = deg(highcorner(std(maxideal(i)*jacob(f))));
772    hc = hc+2-i;
773    if(hc < k)
774    {
775      k = hc;
776    }
777  }
778
779  return(k);
780}
781example
782{
783  "EXAMPLE:";
784  echo = 2;
785  ring r = 0, (x,y), ds;
786  poly f = x3+xy3;
787  determinacy(f);
788}
789
Note: See TracBrowser for help on using the repository browser.