source: git/Singular/LIB/classify.lib @ 0dd77c2

spielwiese
Last change on this file since 0dd77c2 was 341696, checked in by Hans Schönemann <hannes@…>, 14 years ago
Adding Id property to all files git-svn-id: file:///usr/local/Singular/svn/trunk@12231 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 89.1 KB
Line 
1// KK,GMG last modified: 17.12.00
2///////////////////////////////////////////////////////////////////////////////
3version  = "$Id$";
4category="Singularities";
5info="
6LIBRARY:  classify.lib  Arnold Classifier of Singularities
7AUTHOR:   Kai Krueger, krueger@mathematik.uni-kl.de
8
9OVERVIEW:
10   A library for classifying isolated hypersurface singularities w.r.t. right
11   equivalence, based on the determinator of singularities by V.I. Arnold.
12
13PROCEDURES:
14 basicinvariants(f); computes Milnor number, determinacy-bound and corank of
15 classify(f);        normal form of polynomial f determined with Arnold's method
16 corank(f);          computes the corank of f (i.e. of the Hessian of f)
17 Hcode(v);           coding of intvec v acoording to the number repetitions
18 init_debug([n]);    print trace and debugging information depending on int n
19 internalfunctions();display names of internal procedures of this library
20 milnorcode(f[,e]);  Hilbert polynomial of [e-th] Milnor algebra coded with Hcode
21 morsesplit(f);      residual part of f after applying the splitting lemma
22 quickclass(f)       normal form of f determined by invariants (milnorcode)
23 singularity(s,[]);  normal form of singularity given by its name s and index
24 A_L(s/f);           shortcut for quickclass(f) or normalform(s)
25 normalform(s);      normal form of singularity given by its name s
26 debug_log(lev,[]);  print trace and debugging information w.r.t level>@DeBug
27 swap(a,b);          swaps the arguments
28           (parameters in square brackets [] are optional)
29";
30
31LIB "inout.lib";
32LIB "elim.lib";
33LIB "sing.lib";
34LIB "makedbm.lib";
35
36///////////////////////////////////////////////////////////////////////////////
37proc classify_init
38{
39  string s;
40  link l="DBM:r NFlist";
41  s = read(l,"VERSION");
42  if (s == "" ) {
43    if (printlevel)
44    {
45      "(classify.lib): Need to create database...";
46    }
47    create_sing_dbm();
48  }
49  close(l);
50  l="DBM:r NFlist";
51  s = read(l,"VERSION");
52  //"(classify.lib): Creation done. Current version:", s;
53}
54///////////////////////////////////////////////////////////////////////////////
55
56proc classify (poly f_in)
57"USAGE:    classify(f);  f=poly
58COMPUTE:  normal form and singularity type of f with respect to right
59          equivalence, as given in the book \"Singularities of differentiable
60          maps, Volume I\" by V.I. Arnold, S.M. Gusein-Zade, A.N. Varchenko
61RETURN:   normal form of f, of type poly
62REMARK:   This version of classify is only beta. Please send bugs and
63          comments to: \"Kai Krueger\" <krueger@mathematik.uni-kl.de> @*
64          Be sure to have at least @sc{Singular} version 1.0.1. Updates can be
65          found at: @*
66          URL=http://www.mathematik.uni-kl.de/~krueger/Singular/
67NOTE:     type init_debug(n); (0 <= n <= 10) in order to get intermediate
68          information, higher values of n give more information.
69          The proc creates several global objects with names all starting
70          with @, hence there should be no name conflicts.
71EXAMPLE:  example classify; shows an example"
72{
73//---------------------------- initialisation ---------------------------------
74  poly   f_out;
75  int    n, i, corank, show_nf;
76  string s2;
77  list   v;
78  def ring_ext = basering;
79
80  // Namespaces:
81  if(defined(A_Z)==0) { proc A_Z = general_lib::A_Z; export A_Z; }
82  init_debug();                    // initialize trace/debug mode
83  if(checkring()) { return(f_in); }
84  n = nvars(basering);
85  show_nf  = 1;                // return normal form if set to '1'
86
87  // define new ring
88  ring ring_top=char(basering),(x(1..n)),(c,ds);
89
90  map conv_ext2top=ring_ext,maxideal(1);
91
92  if(defined(@ringdisplay)!=0) { kill @ringdisplay; }
93  string @ringdisplay = "ring_ext";
94  export @ringdisplay;
95
96  v = Klassifiziere(conv_ext2top(f_in));
97  if(typeof(v[1])=="poly") {
98     poly f_out = v[1];
99     s2 = v[2];          // s2: Typ des Polynoms f z.b: E[18]
100     corank = v[3];
101  }
102  else {
103     s2="NoClass";
104  }
105
106//---------------- collect results and create return-value --------------------
107  if( s2=="error!" || s2=="NoClass") {
108      setring ring_ext;
109      return(f_in);
110  }
111
112  if(show_nf==1) {
113    poly f_nf = normalform(s2);
114    for(i=corank+1;i<=n;i=i+1) { f_nf = f_nf + x(i)^2; }
115    debug_log(2, "Normal form NF(f)=", f_nf);
116  }
117  if(corank>1) { for(i=corank+1;i<=n;i=i+1) { f_out = f_out + x(i)^2; } }
118  setring ring_ext;
119  map conv_top2ext=ring_top,maxideal(1);
120
121  if(show_nf == 1) { return(conv_top2ext(f_nf)); }
122  else { return(conv_top2ext(f_out)); }
123}
124example
125{"EXAMPLE"; echo=2;
126   ring r=0,(x,y,z),ds;
127   poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
128   classify(f);
129   init_debug(3);
130   classify(f);
131}
132
133///////////////////////////////////////////////////////////////////////////////
134static proc Klassifiziere (poly f)
135{
136//--------------------------- initialisation ----------------------------------
137  string s1;
138  int    n, cnt, corank_f, K, Mu;
139  list   v, cstn;
140  map    PhiG;
141  def ring_top = basering;
142
143  n = nvars(basering);    // Zahl der Variablen des aktuellen Rings.
144  PhiG = ring_top, maxideal(1);
145  cstn[4] = PhiG;
146  if( defined(@ringdisplay) == 0)
147  {
148    string @ringdisplay;               // Define always 'ringdisplay' to be
149    export @ringdisplay;               // able to run 'Show(f)'
150  }
151  @ringdisplay = "setring RingB;";
152  if(defined(RingB)!=0) { kill RingB; }
153  execute ("ring RingB="+charstr(basering)+",("+A_Z("x", n)+"),(c,ds);");
154  export RingB;
155  setring ring_top;
156
157//---------------------- compute basciinvariants ------------------------------
158  if(jet(f,0) != 0 )
159  {
160    cstn[1] = corank(f); cstn[2]=-1; cstn[3]=1;
161    return(printresult(1, f, "a unit", cstn, -1));
162  }
163
164  debug_log(1, "Computing Basicinvariants of f ...");
165  K, Mu, corank_f = basicinvariants(f);
166  debug_log(0, "About the singularity :");
167  debug_log(0, "          Milnor number(f)   = "+string(Mu));
168  debug_log(0, "          Corank(f)          = "+string(corank_f));
169  debug_log(0, "          Determinacy       <= "+string(K));
170  cstn[1] = corank_f;
171  cstn[2] = Mu;
172  cstn[3] = K;
173
174  if( Mu == 0)
175  {
176    cstn[1]=1; cstn[3]=1;
177    return(printresult(1, f, "A[0]", cstn, 0));
178  }
179
180  if(Mu<0)
181  {
182    debug_log(0, "The Milnor number of the function is infinite.");
183    debug_log(0, "The singularity is not in Arnolds list.");
184    return(printresult(1, 1, "error!", cstn, -1));
185  }
186
187  f = jet(f, K);
188  v = HKclass(milnorcode(f));
189  if(v[2]>0) { debug_log(0, "Guessing type via Milnorcode: ", v[1]);}
190  else
191  {
192    debug_log(0, "Hilbert polynomial not recognised. Milnor code = ",
193              milnorcode(f));
194  }
195  debug_log(0, "");
196  debug_log(0, "Computing normal form ...");
197
198//------------ step 1, classification according to corank(f) ------------------
199  if(corank_f == 0)
200  {
201    return(printresult(2, f, "A["+string(Mu)+"]", cstn, 0));
202  }
203  if(corank_f == 1)
204  {
205    return(printresult(2, f, "A["+string(Mu)+"]", cstn, 0));
206  }
207  cstn[4] = 0;
208  if(corank_f == 2) { return(Funktion1bis(f, cstn)); }
209  if(corank_f == 3) { return(Funktion1bis(f, cstn)); }
210  return(printresult(105, f, "NoClass", cstn, -1));
211}
212
213///////////////////////////////////////////////////////////////////////////////
214static proc Funktion1bis (poly f, list cstn)
215{
216//---------------------------- initialisation ---------------------------------
217   def ring_top=basering;
218   poly   g;
219   int    n, corank, K;
220   map    conv, PhiG;
221   string s1;
222   list   v_res, v_class, v, iv;
223
224   corank = cstn[1];
225   K = cstn[3];
226   n = nvars(basering);
227
228//-------------------- apply morsesplit if n>corank ---------------------------
229   if( n > corank) {
230     debug_log(0,
231      "I have to apply the splitting lemma. This will take some time....:-)");
232     v_res = Morse(f, K, corank, 0);
233     g = v_res[1];
234     PhiG = v_res[2];
235
236     conv = ReOrder(g);
237     g = conv(g);
238     PhiG = conv(PhiG);
239
240     if(defined(RingB) != 0 ) { kill RingB; }
241     ring ring_rest=char(basering),(x(1..corank)),(c,ds);
242
243     map MapReduce=ring_top,maxideal(1);
244     poly G = MapReduce(g);
245     map PhiG=ring_top,maxideal(1);// Konstruiere Id auf r
246     PhiG = MapReduce(PhiG);
247
248     execute("ring RingB="+charstr(basering)+",("+A_Z("x",corank)+"),(c,ds);");
249     export RingB;
250     setring ring_rest;
251   }
252   else {
253     ring ring_rest=char(basering),(x(1..corank)),(c,ds);
254     map  PhiG=ring_top,maxideal(1);
255     poly G = PhiG(f);
256   }
257
258//--------------------- step 1 of Arnold now finished -------------------------
259    map phi = ring_rest,maxideal(1);
260    cstn[4] = phi;
261    if(corank == 2) { v_class = Funktion3(G, cstn); }
262    else {
263      if(corank == 3) { v_class = Funktion50(G, cstn); }
264      else { v_class = printresult(1, f, "error!", cstn, -1); }
265    }
266//-------------------------- classification done ------------------------------
267    if(typeof(v_class[1])!="poly") {
268       return(v);
269    }
270    poly f_result = v_class[1];
271    v[2] = v_class[2];
272    v[3] = v_class[3];
273    map Phi = v_class[4];
274    PhiG = Phi(PhiG);
275    setring ring_top;
276    if(defined(conv)!=0) { kill conv; }
277    map conv=ring_rest,maxideal(1);
278    v[1] = conv(f_result);
279    return(v);
280}
281
282///////////////////////////////////////////////////////////////////////////////
283static proc Funktion3 (poly f, list cstn)
284{
285//---------------------------- initialisation ---------------------------------
286    poly f3 = jet(f, 3);
287    ideal Jf;
288    int Dim, Mult, Mu;
289    debug_log(1, "Step 3");
290
291    if( f3 == 0 ) { return(Funktion13(f, cstn)); }
292
293    // f3 ~ x3 , x2y+y3 , x2y
294    Jf = std(jacob(f3));
295    Dim = dim(Jf);
296    if(Dim == 0) { return(printresult(4, f, "D[4]", cstn, 0)); }
297
298    Mult = mult(Jf);
299    Mu = cstn[2];
300    if(Dim == 1) {
301      if( Mult == 1) { return(printresult(5,f,"D["+string(Mu)+"]", cstn, 0)); }
302      if( Mult == 2) { return(Funktion6(f, cstn));}         // series E,J
303      debug_log(0, "dimension 1 und deg != 1, 2 => error, ",
304                        "this should never occur");
305      return(printresult(3, f, "error!", cstn, -1));
306      // Should never reach this line
307    }
308    // Should never reach this line
309    return(printresult(3, f, "error!", cstn, -1));
310}
311
312///////////////////////////////////////////////////////////////////////////////
313static proc Funktion6 (poly f, list cstn)
314{ // Arnold's steps 6-12
315//---------------------------- initialisation ---------------------------------
316    poly   f3, fk;
317    ideal  JetId, Jf;
318    int    k, Dim, Mult, n, Mu;
319    map    PhiG, Phi;
320    intvec RFlg;
321    list   v;
322
323    def ring_top=basering;
324    f3   = jet(f, 3);            // 3-Jet von f
325    n    = nvars(basering);
326    Mu   = cstn[2];
327    PhiG = cstn[4];
328    k    = 1;
329    debug_log(1, "   Step 6");
330
331    RFlg = GetRf(f, n);
332    v    = Faktorisiere(f, f3, 3, 1, RFlg);
333    f    = v[1];
334    Phi  = v[2];
335    PhiG = Phi(PhiG);
336    cstn[4] = PhiG;
337
338//---------------------------- begin of loop ----------------------------------
339    while( (6*k) <= Mu ) {
340      JetId = x(1)^3+x(2)^(3*k);
341      fk    = jet(f, 3*k, weight(JetId));
342      //--------------------------- step 6(k) ---------------------------------
343      if( fk == Coeff(fk,x(1), x(1)^3)*x(1)^3 ) {
344        JetId = x(1)^3+x(2)^(3*k+1);           // check jet x3,y3k+1  : E[6k]
345        fk    = jet(f, 3*(3*k+1), weight(JetId));
346        if( Coeff(fk,x(2),x(2)^(3*k+1)) != 0 ) {
347          return(printresult(7, f, "E["+string(6*k)+"]", cstn, k-1));
348        }
349
350        JetId = x(1)^3+x(1)*x(2)^(2*k+1);      // check jet x3,xy2k+1 : E[6k+1]
351        fk    = jet(f, 3*(2*k+1), weight(JetId));
352        if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(2*k+1)) != 0 ) {
353          return(printresult(8, f,"E["+string(6*k+1)+"]", cstn, k-1));
354        }
355
356        JetId = x(1)^3+x(2)^(3*k+2);           // check jet x3,y3k+1  : E[6k+2]
357        fk    = jet(f, 3*(3*k+2), weight(JetId));
358        if( Coeff(fk,x(2),x(2)^(3*k+2)) != 0 ) {
359          return(printresult(9, f,"E["+string(6*k+2)+"]", cstn, k-1));
360        }
361
362        //------------------------- step 10(k) --------------------------------
363        k++;
364        JetId = x(1)^3+x(2)^(3*k);
365        fk    = jet(f, 3*k, weight(JetId));
366        Jf    = std(jacob(fk));
367        Dim   = dim(Jf);
368
369        if(Dim==0) { return(printresult(11,f,"J["+string(k)+",0]",cstn,k-1)); }
370        Mult = mult(Jf);
371        if( Dim ==1  && Mult==1) {
372          return(printresult(12,f,"J["+string(k)+","+string(Mu - 6*k +2)+"]",
373                 cstn, k-1));
374        }
375        if( Dim == 1  && Mult == 2) {
376          if(Coeff(fk, x(2), x(2)^(3*k)) != 0) {
377            v    = Faktorisiere(f, fk, 3, k, RFlg);
378            f    = v[1];
379            Phi  = v[2];
380            PhiG = Phi(PhiG);
381            cstn[4] = PhiG;
382          }
383        }
384      }
385    }
386    // Should never reach this line
387    return(printresult(6, f, "error!", cstn, -1));
388}
389
390///////////////////////////////////////////////////////////////////////////////
391static proc Funktion13 (poly f, list cstn)
392{
393//---------------------------- initialisation ---------------------------------
394    poly f4;
395    ideal Jf;
396    int Dim, Mult, Mu;
397
398    debug_log(1, "   Step 13");
399    Mu = cstn[2];
400    f4 = jet(f, 4);
401    if( f4 == 0 ) { return(Funktion47(f, cstn)); }
402
403    // f4 ~ x4+ax2y2+y4, x4+x2y2, x2y2, x3y, x4
404    Jf  = std(jacob(f4));
405    Dim = dim(Jf);
406
407    if(Dim==0) { return(printresult(14,f,"X[9] = X[1,0] = T[2,4,4]",cstn,1)); }
408    Mult = mult(Jf);
409    if( Dim == 1) {
410      if( Mult == 1 ) {
411        return(printresult(15, f,
412              "X[1,"+string(Mu-9)+"] = T[2,4,"+string(Mu-5)+"]", cstn, 1));
413      }
414      if( Mult == 2 ) {
415        Jf = Jf, jacob(Jf);
416        Jf = std(Jf);
417        Dim = dim(Jf);
418        if(Dim==0){return(printresult(16,f,"Y[1,p,q] = T[2,4+p,4+q]",cstn,1));}
419        if( Dim == 1 ) { return(Funktion17(f, cstn)); }
420      }
421      if( Mult == 3 ) { return(Funktion25(f, cstn)); }
422    }
423    // Should never reach this line
424    return(printresult(13, f, "error!", cstn, -1));
425}
426
427///////////////////////////////////////////////////////////////////////////////
428static proc Funktion17 (poly f, list cstn)
429{ // Analog zu Fumktion 6, Kombination 17-24
430//---------------------------- initialisation ---------------------------------
431    poly  fk, ft;
432    ideal JetId, Jf;
433    int   p, Dim, Mult, Mu;
434    list  v;
435    map   PhiG, Phi;
436
437    def ring_top=basering;
438    debug_log(1, "      Step 17");
439    Mu   = cstn[2];
440    PhiG = cstn[4];
441    fk   = jet(f, 4);
442    p    = 1;
443
444//---------------------------- begin of loop ----------------------------------
445    while( 3*p<= Mu) {
446      debug_log(1, "         Step 18("+string(p)+")");
447      v   = Isomorphie_s17(f, fk, p, 1);
448      f, Phi  = v[1..2];
449      PhiG    = Phi(PhiG);
450      cstn[4] = PhiG;
451
452      if ( p>1) {
453        JetId = x(1)^3*x(2) + x(2)^(3*p);
454        fk = jet(f, 3*p, weight(JetId));
455      }
456      //--------------------------- step 18(p) --------------------------------
457      JetId = x(1)^3*x(2) + x(2)^(3*p+2);     // check jet x3y,y3k+2  : Z[6p+5]
458      fk = jet(f, 3*(3*p+2), weight(JetId));
459      if( Coeff(fk, x(2), x(2)^(3*p+2)) != 0) {
460        return(printresult(19,f, "Z["+string(6*p+5)+"]", cstn, p));
461      }
462
463      JetId = x(1)^3*x(2)+x(1)*x(2)^(2*p+2);  // check jet x3y,xy2k+2 : Z[6p+6]
464      fk = jet(f, 2*(3*p+2)+1, weight(JetId));
465      if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(2*p+2)) != 0) {
466        return(printresult(20, f, "Z["+string(6*p+6)+"]", cstn, p));
467      }
468
469      JetId = x(1)^3*x(2) + x(2)^(3*p+3);     // check jet x3y,y3k+3  : Z[6p+7]
470      fk = jet(f, 3*(3*p+3), weight(JetId));
471      if( Coeff(fk, x(2), x(2)^(3*p+3)) != 0) {
472        return(printresult(21, f, "Z["+string(6*p+7)+"]", cstn, p));
473      }
474
475      //---------------------------- step 22 ----------------------------------
476      p = p+1;
477      JetId = x(1)^3*x(2) + x(2)^(3*p+1);
478      fk   = jet(f, 3*p+1, weight(JetId));
479      ft   = Teile(fk, x(2));
480      Jf   = std(jacob(ft));
481      Dim  = dim(Jf);
482      Mult = mult(Jf);
483      if(Dim==0) { return(printresult(23,f,"Z["+string(p-1)+",0]", cstn, p)); }
484      if( Mult == 1 ) {
485         return(printresult(24, f, "Z["+string(p-1)+","+string(Mu-3-6*p)+"]",
486                cstn, p));
487      }
488    }
489    // Should never reach this line
490    return(printresult(17, f, "error!", cstn, -1));
491}
492
493///////////////////////////////////////////////////////////////////////////////
494static proc Funktion25 (poly f, list cstn)
495{ // Analog zu Fumktion 6, Kombination 25-46
496//---------------------------- initialisation ---------------------------------
497    poly   fk, ft;
498    ideal  JetId, Jf;
499    int    k, Dim, Mult, Mu;
500    map    PhiG, Phi;
501    intvec RFlg;
502    list   v;
503    def ring_top=basering;
504
505    debug_log(1, "      Step 25");
506    Mu   = cstn[2];
507    PhiG = cstn[4];
508    fk   = jet(f, 4);
509    k    = 1;
510    RFlg = GetRf(f, 2);
511
512//---------------------------- begin of loop ----------------------------------
513    while (k<Mu) {
514      v    =  Faktorisiere(f, fk, 4 , k, RFlg);
515      f, Phi  = v[1..2];
516      PhiG    = Phi(PhiG);
517      cstn[4] = PhiG;
518
519      //--------------------------- step 26(k) --------------------------------
520      JetId = x(1)^4 + x(2)^(4*k+1);          // check jet x4,y4k+1  : W[12k]
521      fk    = jet(f, 4*(4*k+1), weight(JetId));
522      if( Coeff(fk, x(2), x(2)^(4*k+1)) != 0) {
523        return(printresult(27, f, "W["+string(12*k)+"]", cstn, 3*k-2));
524      }
525
526      JetId = x(1)^4 + x(1)*x(2)^(3*k+1);     // check jet x4,xy3k+1 : W[12k+1]
527      fk    = jet(f, 4*(3*k+1), weight(JetId));
528      if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(3*k+1)) != 0) {
529        return(printresult(28, f, "W["+string(12*k+1)+"]", cstn, 3*k-2));
530      }
531
532      //--------------------------- step 29(k) --------------------------------
533      JetId = x(1)^4 + x(2)^(4*k+2);
534      fk    = jet(f, 2*(4*k+2), weight(JetId));
535      if( Coeff(fk, x(2), x(2)^(4*k+2)) != 0) {
536        Jf  = std(jacob(fk));
537        Dim = dim(Jf);
538        if(Dim==0) {return(printresult(30,f,"W["+string(k)+",0]",cstn,3*k-1));}
539        if(Dim==1) {
540           return(printresult(32, f,
541                  "W#["+string(k)+","+string(Mu-12*k-3)+"]", cstn, 3*k-1));
542        }
543        return(printresult(29, f, "error!", cstn, -1));
544      }
545      else {
546        // x^4 oder x^2(x^2+x(2)^2k+1)
547        ft  = Teile(fk, x(1)^2);
548        Jf  = std(jacob(ft));
549        Dim = dim(Jf);
550        if( Dim == 0 ) {
551           return(printresult(31, f, "W["+string(k)+","+string(Mu-12*k-3)+"]",
552                  cstn, 3*k-1));
553        }
554        if( Dim != 1 ) { return(printresult(29, f, "error!", cstn, -1)); }
555
556        //-------------------------- step 33(k) -------------------------------
557        JetId = x(1)^4 + x(1)*x(2)^(3*k+2);   // check jet x4,xy3k+2 : W[12k+5]
558        fk    = jet(f, 4*(3*k+2), weight(JetId));
559        if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(3*k+2)) != 0) {
560          return(printresult(34, f,"W["+string(12*k+5)+"]", cstn, 3*k-1));
561        }
562
563        JetId = x(1)^4 + x(2)^(4*k+3);        // check jet x4,y4k+3  : W[12k+6]
564        fk    = jet(f, 4*(4*k+3), weight(JetId));
565        if( Coeff(fk, x(2), x(2)^(4*k+3)) != 0){
566          return(printresult(35, f,"W["+string(12*k+6)+"]", cstn, 3*k-1));
567        }
568
569        //-------------------------- step 36(k) -------------------------------
570        k = k+1;
571        JetId = x(1)^4 + x(2)^(4*k);
572        fk    = jet(f, (4*k), weight(JetId));
573        Jf    = std(jacob(fk));
574        Dim   = dim(Jf);
575        Mult  = mult(Jf);
576        if(Dim==0) {return(printresult(37,f,"X["+string(k)+",0]",cstn,3*k-1));}
577        if(Dim==1) {
578          if(Mult==1) {
579             return(printresult(38, f,"X["+string(k)+","+string(Mu-12*k+3)+"]",
580                    cstn, 3*k-1));
581          }
582          if(Mult==2) {
583            ft  = Teile(fk, x(1)^2);
584            Jf  = std(jacob(ft));
585            Dim = dim(Jf);
586            if( Dim == 0) { return(Funktion40(f, cstn, k)); }
587            if( Dim == 1) {
588               return(printresult(39, f, "Y["+string(k)+",r,s]", cstn,3*k-1));
589            }
590          }
591          if(Mult!=3) {
592            return(printresult(36, f, "error!", cstn, -1)); }
593        }
594        else { return(printresult(36, f, "error!", cstn, -1)); }
595      }
596    }  // Ende der While-Schleife
597    // Should never reach this line
598    return(printresult(25, f, "error!", cstn, -1));
599}
600
601///////////////////////////////////////////////////////////////////////////////
602static proc Funktion40 (poly f, list cstn, int k)
603{
604//---------------------------- initialisation ---------------------------------
605    int    r, kr, rr, sr, oldDebug;
606    poly   fk, f2, a, b, c;
607    ideal  JetId, Jfsyz;
608    string Typ, RestRing, s1;
609    list   v1, v2;
610    def ring_top=basering;
611
612    debug_log(1, "         Step 40" );
613
614//------------------------------ compute f2 -----------------------------------
615    JetId = x(1)^4 + x(2)^(4*k);
616    fk  = jet(f, (4*k), weight(JetId));
617    f2  = -fk / x(1)^3;
618    Jfsyz = f - fk, x(1)^3, f2;
619    matrix Mat = matrix(syz(Jfsyz));
620    a = Mat[2,1] / Mat[1,1] - Mat[2,2];
621    b = - Mat[3,1] / Mat[1,1] + Mat[3,2];
622    ring tmp_ring=char(basering), (x(1),x(2)),(c,ds);
623    map map_top2tmp=ring_top,maxideal(1);
624    oldDebug = @DeBug;
625    init_debug(-1);
626//------------------------------ classify f2 ----------------------------------
627    v1=Klassifiziere(map_top2tmp(b));
628    init_debug(oldDebug);
629    Typ = v1[2];
630    v2 = DecodeNormalFormString(Typ);
631    Typ, kr, rr, sr = v2[1..4];
632    r   = kr-k;
633    setring ring_top;
634    if( Typ == "E[6k]" ) {
635       return(printresult(42, f, "Z["+string(k)+","+string(12*k+6*r-1)+"]",
636              cstn, 3*k+r-2));
637    }
638    if( Typ == "E[6k+1]" ) {
639       return(printresult(43, f, "Z["+string(k)+","+string(12*k+6*r)+"]",
640              cstn, 3*k+r-2));
641    }
642    if( Typ == "E[6k+2]" ) {
643       return(printresult(44, f, "Z["+string(k)+","+string(12*k+6*r+1)+"]",
644              cstn, 3*k+r-2));
645    }
646    if( Typ == "J[k,0]" ) {
647       return(printresult(45, f, "Z["+string(k)+","+string(r)+",0]",
648              cstn, 3*k+r-2));
649    }
650    if( Typ == "J[k,r]" ) {
651       return(printresult(46,f,"Z["+string(k)+","+string(r)+","+string(rr)+"]",
652              cstn, 3*k+r-2));
653    }
654    // Should never reach this line
655    return(printresult(40, f, "error!", cstn, -1));
656}
657
658///////////////////////////////////////////////////////////////////////////////
659static proc Funktion50 (poly f, list cstn)
660{
661//---------------------------- initialisation ---------------------------------
662    poly  f3;
663    ideal Jf, Jf1, Jf2;
664    int   Dim, Mult, Mu;
665    debug_log(1, "Step 50");
666
667    f3 = jet(f, 3);
668    if( f3 == 0 ) { return(printresult(104, f, "NoClass", cstn, -1)); }
669
670    // f3 ~
671    Jf1  = jacob(f3);
672    Jf   = std(Jf1);
673    Dim  = dim(Jf);
674    Mult = mult(Jf);
675    Mu   = cstn[2];
676
677    if(Dim == 0) {
678       return(printresult(51, f, "P[8] = T[3,3,3]", cstn, 1));
679    } // x3 + y3 + z3 + axyz
680    if(Dim == 1) {
681      if (Mult == 1) {
682        return(printresult(52, f,"P["+string(Mu)+"] = T[3,3,"+string(Mu-5)+"]",
683               cstn, 1));
684      } // x3 + y3 + xyz
685      if(Mult == 2) {
686        Jf2 = wedge(jacob(Jf1),3-Dim), Jf1;
687        Jf2 = std(Jf2);
688        Dim = dim(Jf2);
689        if (Dim==0) { return(printresult(54,f,"R[p,q] = T[3,p,q]", cstn, 1)); }
690        if (Dim==1) { return(Funktion58(f, cstn)); }  // x3 + yz2
691      }
692      if(Mult == 3) {
693        Jf2 = wedge(jacob(Jf1),3-Dim), Jf1;
694        Jf2 = std(Jf2);
695        Dim = dim(Jf2);
696        if(Dim == 0) { return(printresult(56, f, "T[p,q,r]", cstn, 1)); }
697        if(Dim == 1) { return(Funktion66(f, cstn)); }   // x2z + yz2
698      }
699      if(Mult == 4) { return(Funktion82(f, cstn)); }    // x3 + xz2
700    }
701    if(Dim == 2) {
702      if(Mult == 1) { return(Funktion97(f, cstn)); }    // x2y
703      if(Mult == 2) { return(printresult(103,f,"NoClass", cstn, -1));}
704    }
705
706    // Should never reach this line
707    return(printresult(50, f, "error!", cstn, -1));
708}
709
710///////////////////////////////////////////////////////////////////////////////
711static proc Funktion58 (poly fin, list cstn)
712{
713//---------------------------- initialisation ---------------------------------
714    poly  f, f3, a, b, phi, b1, b2, b3, fa, fb, fc;
715    ideal B, Jf3, J1, J2, S, Jfsyz;
716    int   kx, ky, kz;
717    string tp;
718    matrix M[2][3];
719    matrix C[2][3], D;
720    list v;
721    map  PhiG, VERT;
722    def ring_top=basering;
723    debug_log(1, "   Step 58");
724
725    f    = fin;
726    f3   = jet(f, 3);
727    PhiG = cstn[4];
728    tp   = "Nix";
729    kx   = 1;     // Koordinate x
730    ky   = 2;     // Koordinate y
731    kz   = 3;     // Koordinate z
732    B    = maxideal(1);     // ideal fuer Abbildungen
733    Jf3  = jacob(f3);
734    S    = sat(Jf3, maxideal(1))[1];
735    J1   = diff(S[1], x(kx)), diff(S[1], x(ky)), diff(S[1], x(kz)),
736           diff(S[2], x(kx)), diff(S[2], x(ky)), diff(S[2], x(kz));
737    M    = J1;
738    J2   = minor(M, 2), S;
739
740    //------------------ determine coordinate named 'x' -----------------------
741    S  = sat(J2, maxideal(1))[1];
742    J1 = Coeff(S[1], x(1), x(1)), Coeff(S[1], x(2), x(2)),
743          Coeff(S[1], x(3), x(3)), Coeff(S[2], x(1), x(1)),
744          Coeff(S[2], x(2), x(2)), Coeff(S[2], x(3), x(3));
745    C  = J1;
746    D  = syz(C);
747    kill C;
748
749    b1 = D[1,1];
750    b2 = D[2,1];
751    b3 = D[3,1];
752
753    debug_log(6, "f3,s1=", Show(f3));
754    if( b1 != 0) {
755      VERT=ring_top,-1*b1*x(1), -1*b2*x(1)+x(2), -1*b3*x(1) + x(3);
756      kx=1; ky=2; kz=3;
757    }
758    else {
759      if( b2 != 0) {
760        VERT=ring_top, x(1) + -1*b1*x(2), -1*b2*x(2), -1*b3*x(2) + x(3);
761        kx=2; ky=1; kz=3;
762      }
763      else {
764        if( b3 != 0) {
765          VERT=ring_top,x(1) + -1*b1*x(3), x(2) + -1*b2*x(3), -1*b3*x(3);
766          kx=3; ky=1; kz=2;
767        }
768      }
769    }
770    f       = VERT(f);
771    PhiG    = VERT(PhiG);
772    cstn[4] = PhiG;
773    debug_log(6, VERT);
774    f3 = jet(f,3);
775    debug_log(6, "f3,s2=", Show(f3));
776
777    //---------------- compute f_2 such that j3f = xf_2+f_3 -------------------
778    debug_log(6, "1) x=", kx, "  y=", ky, "  z=", kz );
779    matrix C = Coeffs(f3, x(kx));
780    fb=C[2,1];  // Coeff von x^1
781    fc=C[1,1];  // Coeff von x^0
782    if(diff(fb, x(ky)) != 0) {
783      Jfsyz = fb, diff(fb, x(ky));
784      matrix Mat = matrix(syz(Jfsyz));
785      B = maxideal(1);     // setze Abbildungs-ideal
786      if( nrows(Mat) == 2) {
787        poly Relation = -2 * Mat[2,1] / Mat[1,1];
788        b = Coeff(Relation, x(kz), x(kz));
789        B[rvar(x(ky))] = x(ky)-b*x(kz);
790      }
791      else {
792        Jfsyz = fb, diff(fb, x(kz));
793        Mat = matrix(syz(Jfsyz));
794        poly Relation = -2 * Mat[2,1];
795        a = Coeff(Relation, x(ky), x(ky));
796        B[rvar(x(kz))] = x(kz)-a*x(kz);
797        ky, kz = swap(ky, kz);
798      }
799      VERT = ring_top, B;
800      f    = VERT(f);
801      PhiG = VERT(PhiG);
802      cstn[4] = PhiG;
803      f3 = jet(f,3);
804      kill Mat;
805    }
806    else { ky,kz = swap(ky,kz); }
807
808    //------- compute tschirnhaus for 'z' and get f3=f_1(x,y,z)y^2+z^3 --------
809    C    = Coeffs(f3, x(kx));
810    fb   = C[2,1];  // Coeff von x^1
811    fc   = C[1,1];  // Coeff von x^0
812    v    = tschirnhaus(fc, x(kz));
813    fc, VERT = v[1..2];
814    f    = VERT(f);
815    PhiG = VERT(PhiG);
816    cstn[4] = PhiG;
817    f3 = jet(f,3);
818
819    //------------------- compute f_1 and get f3=xy^2+z^3 ---------------------
820    fb = (f3 - 1*(Coeffs(f3, x(kz))[4,1])*x(kz)^3)/(x(ky)^2);
821    fc=(x(kx)-1*(Coeffs(fb,x(ky))[2,1])*x(ky)-1*(Coeffs(fb,x(kz))[2,1])*x(kz));
822    fa = Coeffs(fb, x(kx))[2,1];
823    if ( fa != 0 ) {
824      B = maxideal(1);
825      B[rvar(x(kx))] = fc / fa;
826      VERT = ring_top, B;
827      f    = VERT(f);
828      PhiG = VERT(PhiG);
829      cstn[4] = PhiG;
830      f3   = jet(f,3);
831    }
832
833    //--------------------- permutation of x,y,z  -----------------------------
834    if(Coeffs(f3, x(1))[4,1]!=0) {
835      kx=1;
836      if(Coeffs(f3, x(2))[3,1]==0) { ky=2; kz=3; }
837      else { ky=3; kz=2; }
838    }
839    else {
840      if(Coeffs(f3, x(2))[4,1]!=0) {
841        kx=2;
842        if(Coeffs(f3, x(3))[3,1]==0) { ky=3; kz=1; }
843        else { ky=1; kz=3; }
844      }
845      else {
846        kx=3;
847        if(Coeffs(f3, x(1))[3,1]==0) { ky=1; kz=2; }
848        else { ky=2; kz=1; }
849      }
850    }
851    VERT = ring_top, x(kx), x(ky), x(kz);
852    f    = VERT(f);
853    PhiG = VERT(PhiG);
854    cstn[4] = PhiG;
855    f3   = jet(f,3);
856    return(Funktion59(f, cstn));
857}
858
859///////////////////////////////////////////////////////////////////////////////
860static proc Funktion59 (poly f, list cstn)
861{
862//---------------------------- initialisation ---------------------------------
863    poly   phi, fr, fk, alpha, beta, f_tmp;
864    ideal JetId;
865    int    p, Dim, Mult, Mu;
866    string tp;
867    list   v;
868    map    PhiG, Phi;
869    intvec w, RFlg;
870    def ring_top=basering;
871    debug_log(1, "      Step 59");
872
873    Mu   = cstn[2];
874    PhiG = cstn[4];
875    tp   = "Nix";
876    p    = 1;
877    phi  = jet(f,3);
878    fr   = f - phi;
879    RFlg = 1,2,3;
880    alpha = coeffs(fr, x(1))[1,1];
881    beta = (fr - alpha) / x(1);
882    debug_log(3, "f    = ", Show(f));
883    debug_log(3, "fr   = ", Show(fr));
884    debug_log(3, "alpha= ", Show(alpha));
885    debug_log(3, "beta = ", Show(beta));
886
887//---------------------------- begin of loop ----------------------------------
888    while(6*p<Mu) {
889      JetId = x(2)^(3*p+1);
890      JetId = phi + x(2)^(3*p+1);
891      //--------------------------- step 59(k) --------------------------------
892      w     = weight(JetId);
893      fk    = jet(fr, 3*w[1], w);
894      if(fk!=0) { return(printresult(60,f, "Q["+string(6*p+4)+"]", cstn, p)); }
895
896      JetId = phi + x(1)*x(2)^(2*p+1);
897      w     = weight(JetId);
898      fk    = jet(fr, 3*w[1], w);
899      if(fk!=0) { return(printresult(61,f, "Q["+string(6*p+5)+"]", cstn, p)); }
900
901      JetId = phi + x(2)^(3*p+2);
902      w     = weight(JetId);
903      fk    = jet(fr, 3*w[1], w);
904      if(fk!=0) { return(printresult(62,f, "Q["+string(6*p+6)+"]", cstn, p)); }
905
906      //--------------------------- step 63(k) --------------------------------
907      p++;
908      JetId = phi + x(2)^(3*p);
909      w     = weight(JetId);
910      fk    = jet(f, 3*w[1], w);
911      JetId = std(jacob(fk));
912      Dim   = dim(JetId);
913      Mult  = mult(JetId);
914      if(Dim==0) { return(printresult(64, f, "Q["+string(p)+",0]", cstn, p)); }
915      if(Dim==1) {
916        if(Mult == 1) {
917           return(printresult(65, f, "Q["+string(p)+","+string(Mu-(6*p+2))+"]",
918                  cstn, p));
919        }
920        if(Mult == 2) {
921          fk    = jet(fr, 3*w[1], w);
922          f_tmp = Coeffs(phi, x(1))[4,1] *x(1)^3+fk;
923          v    = Faktorisiere(f, f_tmp, 3 , p, RFlg);
924          f    = v[1];
925          Phi  = v[2];
926          PhiG = Phi(PhiG);
927          cstn[4] = PhiG;
928          fr = f - phi;
929        }
930      }
931    }
932    // Should never reach this line
933    return(printresult(59, f, "error!", cstn, -1));
934}
935
936///////////////////////////////////////////////////////////////////////////////
937static proc Funktion66 (poly f, list cstn)
938{
939//---------------------------- initialisation ---------------------------------
940    int kx = 1; // Koordinate x
941    int ky = 2; // Koordinate y
942    int kz = 3; // Koordinate z
943    poly f3 = jet(f, 3);
944    ideal JetId;
945
946    debug_log(1, "   Step 66");
947    debug_log(2, "F3=", Show(f3));
948    poly fx = diff(f3, x(kx));
949    JetId = jacob(fx);
950    JetId = std(JetId);
951    "nach x:",Show(fx), "  Id=", JetId, "  Dim=", dim(JetId);
952
953    poly fy = diff(f3, x(ky));
954    JetId = jacob(fx);
955    JetId = std(JetId);
956    "nach y:",Show(fy), "  Id=", JetId, "  Dim=", dim(JetId);
957
958    poly fz = diff(f3, x(kz));
959    JetId = jacob(fx);
960    JetId = std(JetId);
961    "nach z:",Show(fz), "  Id=", JetId, "  Dim=", dim(JetId);
962    return(printresult(1, 66, "error!", cstn, -1));
963}
964
965///////////////////////////////////////////////////////////////////////////////
966static proc Funktion82 (poly f, list cstn)
967{
968//---------------------------- initialisation ---------------------------------
969    poly   f3, b1, b2, b3;
970    int    i, kx, ky, kz, Fall;
971    ideal  Jfsyz, B;
972    intvec kv = 1,2,3;
973    matrix Mat;
974    map    PhiG, VERT;
975    list   v;
976    def ring_top=basering;
977    debug_log(1, "   Step 82");
978
979    f3 = jet(f,3);
980    kx = 1;     // Koordinate x
981    ky = 2;     // Koordinate y
982    kz = 3;     // Koordinate z
983    B  = maxideal(1);
984    Jfsyz = jacob(f3);
985    PhiG  = cstn[4];
986    Fall  = 2;
987
988//------------------ find coordinatechange that f3 ~ g(x,z) -------------------
989    if (diff(f3, x(1)) == 0) { kx, ky = swap(kx, ky); }
990    if (diff(f3, x(2)) == 0) {  }
991    if (diff(f3, x(3)) == 0) { kz, ky = swap(kz, ky); }
992    if ( (diff(f3, x(1)) != 0) && (diff(f3, x(2)) != 0) &&
993          (diff(f3, x(3)) != 0) ) {
994      Mat = matrix(syz(Jfsyz));
995      b1  = Mat[1,1];
996      b2  = Mat[2,1];
997      b3  = Mat[3,1];
998
999      if( b1 != 0) {
1000        VERT = ring_top,b1*x(kx), b2*x(kx)+x(ky), b3*x(kx) + x(kz);
1001        kx, ky = swap(kx, ky);
1002      }
1003      else {
1004        if(b2!=0) { VERT = ring_top,x(kx)+b1*x(ky),b2*x(ky),b3*x(ky)+x(kz); }
1005        else {
1006          if(b3!=0) { VERT = ring_top,x(kx)+b1*x(kz),x(ky)+b2*x(kz),b3*x(kz); }
1007          else { VERT = ring_top,B; }
1008        }
1009      }
1010      f    = VERT(f);
1011      PhiG = VERT(PhiG);
1012      cstn[4] = PhiG;
1013    }
1014    VERT = ring_top,x(kx),x(ky),x(kz);
1015    f    = VERT(f);
1016    PhiG = VERT(PhiG);
1017    cstn[4] = PhiG;
1018    f3   = jet(f,3);
1019
1020    if( (f3-subst(f3, x(kx), 0)) == 0) { kx, ky = swap(kx, ky); }
1021    if( (f3-subst(f3, x(kz), 0)) == 0) { kz, ky = swap(kz, ky); }
1022
1023//------------ find coordinatechange for f3 ~ x3+xz2, if possible  ------------
1024    matrix C = coeffs(f3, x(kx));
1025    if(size(C) == 3) { C = coeffs(f3, x(kz)); kx,kz=swap(kx, kz); }
1026    if(C[1,1] == 0 && C[3,1] == 0) { Fall = 1; }
1027    if(C[1,1] != 0 && C[3,1] != 0 ) { Fall = 3; }
1028    if(C[1,1] == 0 && C[3,1] != 0 ) { Fall = 2; }
1029    if(C[1,1] != 0 && C[3,1] == 0 ) { Fall = 2; kx,kz=swap(kx, kz); }
1030
1031    if(Fall == 1) {
1032      VERT=ring_top,x(kx),x(ky),x(kz);
1033    }
1034    if(Fall == 2) {
1035       v = tschirnhaus(f3/x(kz), x(kx));
1036       b1, VERT = v[1..2];
1037    }
1038    if(Fall == 3) {
1039      v = tschirnhaus(f3/x(kx), x(kx));
1040      b1, VERT = [1..2];
1041      debug_log(2, "B1=", Show(jet(VERT(f),3)));
1042      v = tschirnhaus(f3/x(kz), x(kz));
1043      b2, VERT = [1..2];
1044      debug_log(2, "B2=", Show(jet(VERT(f),3)));
1045    }
1046    f    = VERT(f);
1047    PhiG = VERT(PhiG);
1048    cstn[4] = PhiG;
1049    f3   = jet(f,3);
1050
1051//------------- if f3 ~ x3+xz2 then continue with classification  -------------
1052    C = coeffs(f3, x(1));
1053    if( C[1,1] == 0 && C[2,1] != 0 && C[3,1] == 0 && C[4,1] != 0 ) {
1054      return(Funktion83(f, cstn));
1055    }
1056    return(printresult(82, f, "error!", cstn, -1));
1057}
1058
1059///////////////////////////////////////////////////////////////////////////////
1060static proc Isomorphie_s82_z (poly f, poly fk, int p)
1061{
1062//---------------------------- initialisation ---------------------------------
1063    poly   Relation, a, b;
1064    ideal  Jfsyz, B;
1065    matrix Mat;
1066    map    VERT;
1067    list   v;
1068    def ring_top=basering;
1069
1070    debug_log(1, "      Isomorphie 82/90 z");
1071    debug_log(2, "tt=", Show(fk));
1072    Jfsyz = fk, diff(fk, x(3));
1073    Mat   = matrix(syz(Jfsyz));
1074    Relation = -2 * Mat[2,1] / Mat[1,1];
1075    a = Coeff(Relation, x(3), x(3));
1076    b = Coeff(Relation, x(2), x(2)^p);
1077    B = maxideal(1);
1078    B[rvar(x(3))] = x(3)-b*x(2)^p;
1079    VERT = ring_top,B;
1080    v    = VERT(f), VERT;
1081    debug_log(2, VERT);
1082    debug_log(2, "      z res=", Show(VERT(fk)));
1083    return(v);
1084}
1085
1086///////////////////////////////////////////////////////////////////////////////
1087static proc Isomorphie_s82_x (poly f, poly fk, int p)
1088{
1089//---------------------------- initialisation ---------------------------------
1090    matrix Mat;
1091    poly   Relation, a, b;
1092    ideal  Jfsyz, B;
1093    map    VERT;
1094    list   v;
1095    def ring_top=basering;
1096
1097    debug_log(1, "      Isomorphie 82/90 x");
1098    debug_log(2, "tt=", Show(fk));
1099    Jfsyz = fk, diff(fk, x(1));
1100    Mat   = matrix(syz(Jfsyz));
1101    Relation = -3 * Mat[2,1] / Mat[1,1];
1102    a = Coeff(Relation, x(1), x(1));
1103    b = Coeff(Relation, x(2), x(2)^p);
1104    B = maxideal(1);
1105    B[rvar(x(1))] = x(1)-b*x(2)^p;
1106    VERT = ring_top,B;
1107    v    = VERT(f), VERT;
1108    debug_log(2, VERT);
1109    debug_log(2, "      x res=", Show(VERT(fk)));
1110
1111    return(v);
1112}
1113
1114///////////////////////////////////////////////////////////////////////////////
1115static proc Funktion83 (poly f, list cstn)
1116{
1117//---------------------------- initialisation ---------------------------------
1118    poly   fk, phi, a, b;
1119    ideal  JetId, Jf, B;
1120    int    k, Dim, Mult;
1121    intvec w;
1122    map    PhiG, Phi;
1123    list   v;
1124    matrix Mat;
1125    def ring_top=basering;
1126    debug_log(1, "      Step 83");
1127
1128    k    = 1;
1129    PhiG = cstn[4];
1130//---------------------------- begin of loop ----------------------------------
1131    while(k<10) {
1132      phi   = jet(f, 3);
1133      //--------------------------- step 83(k) --------------------------------
1134      JetId = x(1)^3 + x(3)^3 + x(2)^(3*k+1);
1135      fk    = jet(f- phi, 3*w[1], weight(JetId)) ;
1136      if(fk!=0) { return(printresult(84,f,"U["+string(12*k)+"]",cstn,4*k-3)); }
1137
1138      JetId = x(1)^3 + x(3)^3 + x(1)*x(2)^(2*k+1);
1139      fk    = jet(f, 3*w[1], weight(JetId)) ;
1140      //--------------------------- step 85(k) --------------------------------
1141      if ( fk != phi ) {
1142        Jf   = std(jacob(fk));
1143        Dim  = dim(Jf);
1144        if(Dim==0) {return(printresult(86,f,"U["+string(k)+",0]",cstn,4*k-2));}
1145        if(Dim==1) {return(printresult(87,f,"U["+string(k)+",p]",cstn,4*k-2));}
1146      }
1147
1148      //--------------------------- step 88(k) --------------------------------
1149      JetId = x(1)^3 + x(3)^3 + x(2)^(3*k+2);
1150      fk    = jet(f- phi, 3*w[1], weight(JetId)) ;
1151      if(fk!=0) {return(printresult(89,f,"U["+string(12*k+4)+"]",cstn,4*k-2));}
1152
1153      //--------------------------- step 90(k) --------------------------------
1154      k++;
1155      JetId = x(1)^3 + x(3)^3 + x(2)^(3*k);
1156      fk    = jet(f, 3*w[1], weight(JetId)) ;
1157      Jf    = std(jacob(fk));
1158      Dim   = dim(Jf);
1159      Mult  = mult(Jf);
1160      if ( Dim == 0 ) { return(printresult(83, f, "NoClass", cstn, -1)); }
1161      if ( Dim == 1 ) {
1162        if ( Mult == 4 ) {
1163          if( fk - phi != 0) { // b!=0  und/oder b'!=0
1164            if( Coeff(fk,x(1)*x(2), x(1)^2*x(2)^k) == 0 ) { // b=0 und b'!=0
1165              a    = (fk - Coeff(fk, x(1), x(1)^3)*x(1)^3) / x(1);
1166              v    = Isomorphie_s82_z(f, a, k);
1167            }
1168            else {
1169              if( Coeff(fk,x(1)*x(2)*x(3), x(1)*x(2)^k*x(3)) == 0 ){
1170                        // b!=0 und b'=0
1171                a    = subst(fk, x(3), 0);
1172                v    = Isomorphie_s82_x(f, a, k);
1173              }
1174              else {
1175                a = Coeff(fk,x(1)*x(2)*x(3), x(1)*x(2)^k*x(3));
1176                b = Coeff(fk,x(2)*x(3), x(2)^(2*k)*x(3));
1177                B = maxideal(1);
1178                B[rvar(x(1))] = x(1)-b/a*x(2)^k;
1179                Phi  = ring_top,B;
1180                f    = Phi(f);
1181                PhiG = Phi(PhiG);
1182                cstn[4] = PhiG;
1183                fk   = jet(f, 3*w[1], w) ;
1184                a    = (fk - Coeff(fk, x(1), x(1)^3)*x(1)^3) / x(1);
1185                v    = Isomorphie_s82_z(f, a, k);
1186              } // ende else b!=0 und b'=0
1187            } // ende else b=0 und b'!=0
1188            f, Phi  = v[1..2];
1189            PhiG    = Phi(PhiG);
1190            cstn[4] = PhiG;
1191          } //ende fk-phi!=0
1192        } // ende mult=4
1193        else { return(printresult(83, f, "NoClass", cstn, -1)); }
1194      } // ende dim=1
1195      else { return(printresult(83, f, "NoClass", cstn, -1)); }
1196    } // ENDE While
1197    return(printresult(83, f, "error!", cstn, -1));
1198}
1199
1200///////////////////////////////////////////////////////////////////////////////
1201static proc Funktion97 (poly f, list cstn)
1202{
1203//---------------------------- initialisation ---------------------------------
1204    poly   f3, l1, l2, a, b, c, prod;
1205    ideal  Jfsyz, Jf, B;
1206    int    k, i, pt, kx, ky, kz, Dim, Mult, Mu;
1207    matrix Mat;
1208    map    PhiG, VERT;
1209    def ring_top=basering;
1210    debug_log(1, "   Step 97");
1211
1212    Mu   = cstn[2];
1213    PhiG = cstn[4];
1214    kx   = 1;   // Koordinate x
1215    ky   = 2;   // Koordinate y
1216    kz   = 3;   // Koordinate z
1217    B    = maxideal(1); // Abbildungs-ideal
1218    pt   = 2;
1219    f3   = jet(f, 3);
1220    k    = 1;
1221
1222//--------------------------- compute f3 ~ x2y --------------------------------
1223    // vertausche 2 Koordinaten sodass d2f/dx2 <>0 ist.
1224    for(i=1;i<4;i=i+1) {
1225      if(diff(diff(f3, x(i)), x(i)) != 0) { kx = i; i=4; }
1226    }
1227    if(kx == 2) { ky = 1; kz = 3; }
1228    if(kx == 3) { ky = 2; kz = 1; }
1229    //-------------------------- compute -l1*l2 -------------------------------
1230    f3    = jet(f, 3);
1231    Jfsyz = f3, diff(f3, x(kx));
1232    Mat   = matrix(syz(Jfsyz));
1233    if(deg(Mat[2,1])>1) {
1234      Jfsyz = f3, Mat[2,1];
1235      Mat   = matrix(syz(Jfsyz));
1236
1237      // berechen Abb. sodass f=x2*l2
1238      l1 = Mat[2,1];
1239      a  = Coeff(l1, x(kx), x(kx));
1240      l1 =  l1 / number(a);
1241      b  = Coeff(l1, x(ky), x(ky));
1242      c  = Coeff(l1, x(kz), x(kz));
1243      B[rvar(x(kx))] = x(kx) - b * x(ky) - c * x(kz);
1244      VERT = ring_top, B;
1245      f    = VERT(f);
1246      PhiG = VERT(PhiG);
1247      cstn[4] = PhiG;
1248
1249      f3 = jet(f, 3);
1250      l2 = f3 / x(kx)^2;
1251
1252      // sorge dafuer, dass b<>0 ist.
1253      b = Coeff(l2, x(ky), x(ky));
1254      if( b== 0) { ky, kz = swap(ky, kz); }
1255
1256      // Koordinaten-Transf. s.d. f=x2y
1257      b  = Coeff(l2, x(ky), x(ky));
1258      l2 =  l2 / number(b);
1259      a  = Coeff(l2, x(kx), x(kx));
1260      c  = Coeff(l2, x(kz), x(kz));
1261      B  = maxideal(1);
1262      B[rvar(x(ky))] = -a * x(kx) + x(ky) - c * x(kz);
1263      VERT = ring_top, B;
1264      f    = VERT(f);
1265      PhiG = VERT(PhiG);
1266      cstn[4] = PhiG;
1267    }
1268
1269//------------------------------- step 98 ---------------------------------
1270    Jfsyz = x(kx)^2*x(ky) + x(ky)^4 + x(kz)^4;
1271    a     = jet(f, 8, weight(Jfsyz));
1272    Jf    = std(jacob(a));
1273    Dim   = dim(Jf);
1274    Mult  = mult(Jf);
1275    if( Dim == 0) { return(printresult(99, f, "V[1,0]", cstn, 3)); }
1276    if( Dim == 1) {
1277      if(Mult==1) {return(printresult(100,f,"V[1,"+string(Mu-15)+"]",cstn,3));}
1278      if(Mult==2){return(printresult(101,f,"V#[1,"+string(Mu-15)+"]",cstn,3));}
1279    }
1280    " Dim=",Dim," Mult=",Mult;
1281    return(printresult(102, f, "NoClass", cstn, -1));
1282}
1283
1284///////////////////////////////////////////////////////////////////////////////
1285proc tschirnhaus (poly f, poly x)
1286"USAGE:    tschirnhaus();"
1287{
1288//---------------------------- initialisation ---------------------------------
1289    poly   b;
1290    ideal  B;
1291    int    n, j, hc;
1292    matrix cf;
1293    intvec z;
1294    string s;
1295    list   v;
1296    map    Phi, EH;
1297    def ring_top=basering;
1298
1299    n    = nvars(basering);
1300    cf   = coeffs(f, x);
1301    hc   = nrows(cf) - 1; // hoechster exponent von x_i
1302    b    = cf[hc+1,1];    // koeffizient von x_i^hc
1303    B    = maxideal(1);
1304    z[n] = 0;
1305    EH   = ring_top, z;
1306    Phi  = ring_top, B;
1307    v[1] = f;
1308    if ( EH(b) != 0)    // pruefe ob der Koeff von x_i^hc
1309    { B[rvar(x)] = x -1*(cf[hc,1]/(hc*b));
1310      v[1] = Phi(f);
1311    }
1312    v[2] = Phi;
1313    return(v);
1314}
1315
1316///////////////////////////////////////////////////////////////////////////////
1317static proc Isomorphie_s17 (poly f, poly fk, int k, int ct, list #)
1318{
1319//---------------------------- initialisation ---------------------------------
1320    ideal  Jfsyz, JetId, bb;
1321    poly   a, b, c, d, Relation, an, bn;
1322    int    B,C, alpha, beta, gamma, g;
1323    matrix Matx, Maty;
1324    map    PhiG, VERT;
1325    list   v;
1326    def ring_top=basering;
1327
1328    if(size(#)==1) { PhiG = #[1]; }
1329    else { PhiG = ring_top,maxideal(1); }
1330    bb = maxideal(1);
1331
1332    // Ziel: bestimme a,b,c,d sodass  fk = (ax+by^k)^3(cx+dy) gilt.
1333    debug_log(2, "Isomorphie_s17:");
1334    debug_log(2, "Faktor: f=",Show(f)," Jet=",Show(fk)," k=",k, "cnt=", ct);
1335
1336    if( k == 1) {
1337      Jfsyz = fk, diff(fk, x(1));
1338      Matx  = matrix(syz(Jfsyz));
1339      Jfsyz = fk, diff(fk, x(2));
1340      Maty  = matrix(syz(Jfsyz));
1341
1342      a = Coeff(fk, x(1), x(1)^4);
1343      b = Coeff(fk, x(2), x(2)^4);
1344      c = Coeff(fk, x(1)*x(2), x(1)^3*x(2));
1345      d = Coeff(fk, x(1)*x(2), x(1)*x(2)^3);
1346
1347      if( (a != 0) && (b != 0) ) {
1348        B = -int(Coeff(Matx[1,1], x(2), x(2)));
1349        C = -int(Coeff(Maty[1,1], x(1), x(1)));
1350        alpha = int(Coeff(Matx[2,1], x(1), x(1)^2));
1351        beta  = int(Coeff(Matx[2,1], x(1)*x(2), x(1)*x(2)));
1352        gamma = int(Coeff(Matx[2,1], x(2), x(2)^2));
1353
1354        bb[rvar(x(1))] = x(1) - (2*gamma / (B - beta))*x(2);
1355        bb[rvar(x(2))] = x(2) - ((C - beta) / (2*gamma))*x(1);
1356        VERT     = ring_top,bb;
1357        Relation = VERT(f);
1358        fk       = jet(Relation, 4);
1359
1360        an = Coeff(fk, x(1), x(1)^4);
1361        bn = Coeff(fk, x(2), x(2)^4);
1362        if( (an != 0) & (bn != 0) ) { VERT=ring_top,x(1),(x(2) + a*x(1))/ b; }
1363        f    = VERT(f);
1364        fk   = jet(f, 4);
1365        PhiG = VERT(PhiG);
1366
1367        a = Coeff(fk, x(1), x(1)^4);
1368        b = Coeff(fk, x(2), x(2)^4);
1369        c = Coeff(fk, x(1)*x(2), x(1)^3*x(2));
1370        d = Coeff(fk, x(1)*x(2), x(1)*x(2)^3);
1371        Jfsyz = fk, diff(fk, x(1));
1372        Matx  = matrix(syz(Jfsyz));
1373        Jfsyz = fk, diff(fk, x(2));
1374        Maty  = matrix(syz(Jfsyz));
1375      }
1376
1377      if( (a == 0) || (b == 0) ) {
1378        if( a == 0) {
1379          if( c == 0) { // y3(ax+by)
1380            Relation = - Matx[2,1] / Matx[1,1];
1381            a = Coeff(Relation, x(1), x(1));
1382            b = Coeff(Relation, x(2), x(2));
1383            VERT=ring_top,a*x(2)^k - b*x(1), x(1);
1384          }
1385          else { // (ax+by)^3y
1386            Relation = - 3*Matx[2,1] / Matx[1,1];
1387            a = Coeff(Relation, x(1), x(1));
1388            b = Coeff(Relation, x(2), x(2));
1389            VERT=ring_top,a*x(1) - b*x(2), x(2);
1390          }
1391        }
1392        else {
1393          if( d == 0) { // x3(ax+by)
1394            Relation = - Maty[2,1] / Maty[1,1];
1395            a = Coeff(Relation, x(1), x(1));
1396            b = Coeff(Relation, x(2), x(2));
1397            VERT=ring_top,x(1), b*x(2)^k - a*x(1);
1398          }
1399          else { // x(ax+by)^3
1400            Relation = - 3*Maty[2,1] / Maty[1,1];
1401            a = Coeff(Relation, x(1), x(1));
1402            b = Coeff(Relation, x(2), x(2));
1403            VERT=ring_top,x(2), b*x(1) - a*x(2);
1404          }
1405        }
1406        f    = VERT(f);
1407        PhiG = VERT(PhiG);
1408      }
1409      else {  //      "Weder b noch a sind 0";
1410        if(ct > 5) { v[1]=f; v[2]=PhiG; return(v); }
1411        fk = jet(f, 4);
1412        return(Isomorphie_s17(f, fk, k, ct+1, PhiG));
1413      }
1414    }
1415    else {  // k >1
1416      a     = fk/x(2);
1417      Jfsyz = a, diff(a, x(1));
1418      Matx  = matrix(syz(Jfsyz));
1419      Relation = -3 * Matx[2,1] / Matx[1,1];
1420      a    = Coeff(Relation, x(1), x(1));
1421      b    = Coeff(Relation, x(2), x(2)^k);
1422      VERT = basering,x(1)-b*x(2)^k,x(2);
1423      f    = VERT(f);
1424      PhiG = VERT(PhiG);
1425      JetId = x(1)^3*x(2) + x(2)^(3*k+1);
1426      fk = jet(f, 3*k+1, weight(JetId));
1427    }
1428    v = f, PhiG;
1429    debug_log(2, "Isomorphie_s17: done");
1430    debug_log(2, "Faktor: f=",Show(f)," Jet=",Show(fk)," k=",k);
1431
1432    return(v);
1433}
1434
1435///////////////////////////////////////////////////////////////////////////////
1436static proc printresult (int step, poly f, string typ, list cstn, int m)
1437{
1438//---------------------------- initialisation ---------------------------------
1439  int corank, Mu, K;
1440  list v;
1441
1442  corank, Mu, K = cstn[1..3];
1443  debug_log(0,"   Arnold step number "+string(step));
1444  if( @DeBug >= 0 )
1445  {
1446    "The singularity";
1447    "   "+Show(jet(f, K))+"";
1448    if( typ != "error!" && typ != "NoClass" )
1449    {
1450      "is R-equivalent to "+typ+".";
1451    }
1452    if ( typ == "NoClass" )
1453    {
1454      "is not in Arnolds list.";
1455    }
1456//    if(K>=0)  { "  det = "+string(K); }
1457    if(Mu>=0) { "   Milnor number = "+string(Mu); }
1458    if(m>=0)  { "   modality      = "+string(m); }
1459  }
1460  v = f, typ, corank, cstn[4];
1461  return(v);
1462}
1463
1464///////////////////////////////////////////////////////////////////////////////
1465static proc Funktion47 (poly f, list cstn)
1466{
1467  int corank = cstn[1];
1468  int Mu = cstn[2];
1469  int K  = cstn[3];
1470  string s = "The Singularity "+Show(jet(f, K));
1471  string tp="";
1472//  return(printresult(47, f, tp, cstn, -1));
1473
1474  s = s +" has 4-jet equal to zero. (F47), mu="+string(Mu);
1475
1476  s; // +"  ("+SG_Typ+")";
1477  s = "No further classification available.";
1478  s;
1479  return(Show(f), tp, corank);
1480}
1481
1482///////////////////////////////////////////////////////////////////////////////
1483static proc Funktion91 (poly f, list cstn, int k)
1484{
1485  string tp  = "U*[k,0]";
1486  return(printresult(91, f, tp, cstn, -1));
1487}
1488
1489///////////////////////////////////////////////////////////////////////////////
1490static proc Funktion92 (poly f, list cstn, int k)
1491{
1492  string tp  = "UP[k]";
1493  return(printresult(92, f, tp, cstn, -1));
1494}
1495
1496///////////////////////////////////////////////////////////////////////////////
1497static proc Funktion93 (poly f, list cstn, int k)
1498{
1499  string tp  = "UQ[k]";
1500  return(printresult(93, f, tp, cstn, -1));
1501}
1502
1503///////////////////////////////////////////////////////////////////////////////
1504static proc Funktion94 (poly f, list cstn, int k)
1505{
1506  string tp  = "UR[k]";
1507  return(printresult(94, f, tp, cstn, -1));
1508}
1509
1510///////////////////////////////////////////////////////////////////////////////
1511static proc Funktion95 (poly f, list cstn, int k)
1512{
1513  string tp  = "US[k]";
1514  return(printresult(95, f, tp, cstn, -1));
1515}
1516
1517///////////////////////////////////////////////////////////////////////////////
1518static proc Funktion96 (poly f, list cstn, int k)
1519{
1520  string tp  = "UT[k]";
1521  return(printresult(96, f, tp, cstn, -1));
1522}
1523
1524///////////////////////////////////////////////////////////////////////////////
1525proc morsesplit(poly f)
1526"
1527USAGE:    morsesplit(f);        f=poly
1528RETURN:   Normal form of f in M^3 after application of the splitting lemma
1529COMPUTE:  apply the splitting lemma (generalized Morse lemma) to f
1530EXAMPLE:  example morsesplit; shows an example"
1531{
1532//---------------------------- initialisation ---------------------------------
1533  poly f_out;
1534  int  n, K, Mu, corank;
1535  list v;
1536
1537  if(defined(@ringdisplay) != 0 ) { kill @ringdisplay; }
1538  string @ringdisplay = "setring "+nameof(basering);
1539  export @ringdisplay;
1540
1541  def ring_ext=basering;
1542
1543  n = nvars(basering);
1544
1545  // if trace/debug mode not set, do it!
1546  init_debug();
1547  K, Mu, corank = basicinvariants(f);
1548  ring ring_top=char(basering),(x(1..n)),(c,ds);
1549
1550  map Conv=ring_ext,maxideal(1);
1551  setring ring_top;
1552  v = Morse(jet(Conv(f),K), K, corank, 0);
1553  poly f_out = v[1];
1554  setring ring_ext;
1555  map ConvUp = ring_top, maxideal(1);
1556  return(ConvUp(f_out));
1557}
1558example
1559{ "EXAMPLE"; echo=2;
1560   ring r=0,(x,y,z),ds;
1561   export r;
1562   init_debug(1);
1563   poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
1564   poly g=morsesplit(f);
1565   g;
1566}
1567
1568///////////////////////////////////////////////////////////////////////////////
1569static proc Coeffs (list #)
1570{
1571  matrix m=matrix(coeffs(#[1],#[2]), deg(#[1])+1, 1);
1572  return(m);
1573}
1574
1575///////////////////////////////////////////////////////////////////////////////
1576static proc Morse(poly f, int K, int corank, int ShowPhi)
1577{
1578//---------------------------- initialisation ---------------------------------
1579  poly   fc, f2, a, P, Beta, fi;
1580  ideal  Jfx, B;
1581  int    n, i, j, k, Rang, Done;
1582  matrix Mat;
1583  map    Id, Psi, Phi, PhiG;
1584  intvec Abb, RFlg;
1585  list   v;
1586
1587  fi = f;
1588  n = nvars(basering);
1589  init_debug();
1590
1591  def ring_top=basering;
1592
1593  debug_log(3, "Split the polynomial below using determinacy: ", string(K));
1594  debug_log(3, Show(fi));
1595
1596  for( j=1; j<=n ; j++) { Abb[j] = 0; }
1597
1598  RFlg = GetRf(fi, n);
1599  debug_log(2, "Permutations:", RFlg );
1600  PhiG=ring_top,maxideal(1);
1601
1602//----------------- find quadratic term, if there is only one -----------------
1603  B = maxideal(1);
1604  if(corank == (n-1))
1605  {
1606    Done = 0;
1607    f2   = jet(fi, 2);
1608    j    = 1;
1609    Jfx  = f2, diff(f2, x(j));
1610    while(j<=n && (diff(f2, x(j))==0))
1611    {
1612      j++;
1613      Jfx = f2, diff(f2, x(j));
1614    }
1615    Mat  = matrix(syz(Jfx));
1616    Beta = 2*Mat[2,1]/Mat[1,1];
1617    for( j=1; j<=n ; j++)
1618    {
1619      f2 = Coeff(Beta, x(RFlg[j]), x(RFlg[j]));
1620      if(f2!=0)
1621      {
1622        k = RFlg[j];
1623        break;
1624      }
1625    }
1626    for( j=1; j<=n ; j=j+1)
1627    {
1628      f2 = Coeff(Beta, x(j), x(j));
1629      if(j == k) { B[rvar(x(j))] = (2*f2*x(j)-Beta) / number(f2); }
1630    }
1631    Phi  =ring_top,B;
1632    fi   = Phi(fi);
1633    PhiG = Phi(PhiG);
1634  }
1635  if( ShowPhi > 1) { PhiG; }
1636
1637//------------------------ compute spliting lemma -----------------------------
1638  fc = fi;
1639  i  = 1;              // Index fuer Variablen wird bearbeitet
1640  while( i <= n)
1641  {
1642    Phi=ring_top,maxideal(1);
1643    debug_log(6, "Pruefe Variable x(" +string(RFlg[i]) + ")");
1644    debug_log(6, "--------------------");
1645    j  = i + 1;        // setze j fuer evtle Verschiebung
1646    f2 = jet(fc,2);
1647    debug_log(6, "Rechne 2-Jet =" , string(f2));
1648    if( (f2 - subst(f2, x(RFlg[i]), 0)) == 0 ) { Abb[RFlg[i]] = 1; }
1649    if( (f2 - subst(f2, x(RFlg[i]), 0)) != 0 )
1650    {
1651      while( (j<=n) || (i==n) )
1652      {
1653        debug_log(6, "Pruefe 2-Jet mit Wert : " + string(jet(fc,2)));
1654        a = Coeff(jet(fc,2), x(RFlg[i]), x(RFlg[i])^2);
1655        debug_log(6,"Koeffizient von x(" + string(RFlg[i]) + ")^2 ist:", a);
1656        if( (a != 0) || (i==n) )
1657        {
1658          debug_log(6, "BREAK!!!!!!!!!!!!!!");
1659          break;
1660        }
1661        debug_log(6,"Verschiebe evtl Variable x(",string(RFlg[j]),") um x(",
1662                     string(RFlg[i]), ")");
1663        B = maxideal(1);
1664        for( k=1; k<=n ; k=k+1)
1665        {
1666          if(k==RFlg[j]) { B[rvar(x(k))] = x(k) + x(RFlg[i]); }
1667        }
1668        Phi = ring_top,B;
1669        fc  = Phi(fi);
1670        j++;
1671      }               // Ende while( (j<=n) || (i==n) )
1672
1673      debug_log(6, "Moegliche Verschiebung fertig!");
1674      PhiG = Phi(PhiG);
1675      if( ShowPhi > 1) { "NachVersch.:"; Phi; }
1676
1677      if( (j<=n) || (i==n))
1678      {
1679        P = Coeff(fc, x(RFlg[i]), x(RFlg[i]));
1680        debug_log(6, "Koeffizient von x("+string(RFlg[i])+") ist: ",
1681                      string(P));
1682        if(P != 0)
1683        {
1684          debug_log(6, "1 Koeffizient von x("+string(RFlg[i])+") ist: ",
1685                       string(P));
1686          debug_log(6, "a=" + string(a));
1687          P = P / number (2 * a);
1688          debug_log(6, "2 Koeffizient von x("+string(RFlg[i])+") ist: ",
1689                       string(P));
1690          B = maxideal(1);
1691          for( k=1; k<=n ; k=k+1) {if(k==RFlg[i]) {B[rvar(x(k))]=x(k)-P;}}
1692          Phi =ring_top,B;
1693          debug_log(6, "Quadratische-Ergaenzung durch:", Phi);
1694          fi   = Phi(fc);
1695          PhiG = Phi(PhiG);
1696          fc   = jet(fi,K);
1697          P    = Coeff(fc, x(RFlg[i]), x(RFlg[i]));
1698          if( P != 0)
1699          {
1700            fi = fc;
1701            continue;
1702          }
1703        }     // Ende if(P != 0)
1704              // Fertig mit Quadratischer-Ergaenzung
1705      }               // Ende if( (j<=n) || (i==n))
1706    }                 // Ende if( (f2 - subst(f2, x(RFlg[i]), 0)) != 0 )
1707
1708    fi = fc;
1709    i++;
1710    debug_log(6, "++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
1711  }
1712  debug_log(6, "Ende  ---------------------------------------------------");
1713
1714//--------------------------- collect results ---------------------------------
1715  if( ShowPhi > 0 )
1716  {
1717    "Abbildung innerhalb des Morse-Lemmas:";
1718    PhiG;
1719    "Vergleich:";
1720    "PhiG(f)= " + Show(jet(PhiG(f), K));
1721    "fi     = " + Show(fi);
1722  }
1723
1724  Rang = 0;
1725  B    = maxideal(1);
1726  for( i=1; i<=n ; i++) { if(Abb[i] != 1) { Rang ++; B[rvar(x(i))] = 0; } }
1727  Phi  = ring_top,B;
1728  PhiG = Phi(PhiG);
1729  fi   = Phi(fi);
1730  v    = fi, PhiG;
1731  debug_log(2, "rank determined with Morse rg=", Rang);
1732  debug_log(1, "Residual singularity f=",Show(fi));
1733  return(v);
1734}
1735
1736///////////////////////////////////////////////////////////////////////////////
1737static proc Coeff(poly f, list #)
1738{
1739//---------------------------- initialisation ---------------------------------
1740  poly   a, term;
1741  int    n, i;
1742  matrix K;
1743
1744  n     = nvars(basering);
1745  i     = 1;
1746  term  = #[2];
1747  K     = coef(f, #[1]);
1748
1749  while( (i<=ncols(K)) && (K[1,i] != term) )
1750  { i++;
1751    if(i>ncols(K)) { break; }
1752  }
1753  if(i<=ncols(K)) { a = K[2,i]; }
1754  if(i>ncols(K)) { a = 0; }
1755
1756  return(a);
1757}
1758
1759///////////////////////////////////////////////////////////////////////////////
1760static proc ReOrder(poly f)
1761{
1762//---------------------------- initialisation ---------------------------------
1763  poly  result;
1764  ideal B = maxideal(1);
1765  int   i, n, Ctn, Ctv;
1766  map   conv;
1767
1768  n = nvars(basering);
1769  Ctv = 1;              // Zahl der Vorhandenen Variablen
1770  Ctn = n;              // Zahl der Nicht-Vorhandenen Variablen
1771  def ring_top=basering;
1772
1773  for( i=1; i<=n; i=i+1)
1774  { result = subst(f,x(i), 0) - f;
1775    if( result != 0 ) { B[rvar(x(i))] = x(Ctv); Ctv++; }
1776    else { B[rvar(x(i))] = x(Ctn); Ctn--; }
1777  }
1778
1779  conv = ring_top,B;
1780  return(conv);
1781}
1782
1783///////////////////////////////////////////////////////////////////////////////
1784proc quickclass(poly f)
1785"
1786USAGE:    quickclass(f);         f=poly
1787RETURN:   Normal form of f in Arnold's list
1788REMARK:   try to determine the normal form of f by invariants, mainly by
1789          computing the Hilbert function of the Milnor algebra,
1790          no coordinate change is needed (see also proc 'milnorcode').
1791EXAMPLE:  example quickclass; shows an example"
1792{
1793//---------------------------- initialisation ---------------------------------
1794  string Typ;
1795  int    cnt, K, Mu, corank;
1796  list   v;
1797  def ring_top=basering;
1798  // check basic condition on the basering.
1799  if(checkring()) { return(f); }
1800  if( f==0 ) {
1801    "Normal form : 0";
1802    return(f);
1803  }
1804  if( jet(f,0)!=0 ) {
1805    "Normal form : 1";
1806    return(f);
1807  }
1808  K, Mu, corank = basicinvariants(f);
1809  if(Mu<0) {
1810    debug_log(0, "The Milnor number of the function is infinite.");
1811    return(f);
1812  }
1813
1814  // Do the classification of f
1815  // typ: list of types matching the milnorcode
1816  // cnt: number of matches found
1817  v = HKclass(milnorcode(f));
1818  Typ, cnt = v[1..2];
1819  "Singularity R-equivalent to :",Typ;
1820  if(cnt==0) {
1821    "Hilbert polynomial not recognised. Milnor code = ", milnorcode(f);
1822    return();
1823  }
1824  if(cnt==1) {
1825    debug_log(1,"Getting normal form from database.");
1826    "normal form :",A_L(Typ);
1827    return(A_L(Typ));
1828  }
1829  // Hier nun der Fall cnt>1
1830  "Hilbert-Code of Jf^2";
1831  "We have ", cnt, "cases to test";
1832  Cubic(f);
1833  return(v);
1834}
1835example
1836{ "EXAMPLE:"; echo=2;
1837   ring r=0,(x,y,z),ds;
1838   poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
1839   quickclass(f);
1840}
1841
1842///////////////////////////////////////////////////////////////////////////////
1843proc milnorcode (poly f, list #)
1844"USAGE:    milnorcode(f[,e]); f=poly, e=int
1845RETURN:   intvec, coding the Hilbert function of the e-th Milnor algebra
1846          of f, i.e. of basering/(jacob(f)^e) (default e=1), according
1847          to proc Hcode
1848EXAMPLE:  example milnorcode; shows an example"
1849{
1850  int  e=1;
1851  if(size(#)==1) { e=#[1]; }
1852  ideal jf=std(jacob(f)^e);
1853  intvec v=hilb(jf,2);
1854
1855  return(Hcode(v));
1856}
1857example
1858{ "EXAMPLE:"; echo=2;
1859  ring r=0,(x,y,z),ds;
1860  poly f=x2y+y3+z2;
1861  milnorcode(f);
1862  milnorcode(f,2);  // a big second argument may result in memory overflow
1863}
1864
1865///////////////////////////////////////////////////////////////////////////////
1866proc Hcode (intvec v)
1867"USAGE:    Hcode(v); v=intvec
1868RETURN:   intvec, coding v according to the number of successive
1869          repetitions of an entry
1870EXAMPLE:  example Hcode; shows an example."
1871{
1872  int    col, len, i, cur, cnt, maxcoef, nlen;
1873  intvec hil1, hil2;
1874
1875  col      = 1;
1876  len      = size(v);
1877  v[len+1] = 0;
1878
1879  init_debug();
1880  debug_log(1, "Hcode:", v );
1881
1882  for(i=1; i<=len; i++) { if( v[i] > maxcoef) { maxcoef = v[i]; } }
1883
1884  nlen       = 2*maxcoef - 1;
1885  hil1[nlen] = 0;
1886  hil2[nlen] = 0;
1887
1888  for(i=1; i<=nlen; i++)
1889  { if( i > maxcoef) { hil2[i] = 2*maxcoef-i; }
1890    else { hil2[i] = i; }
1891  }
1892
1893  for(i=1; i<=nlen; i++)
1894  { cnt=0;
1895    while( (col<=len) && (v[col] == hil2[i]) )
1896    { cnt++; col++; }
1897    hil1[i] = cnt;
1898  }
1899  return(hil1);
1900}
1901example
1902{ "EXAMPLE:"; echo=2;
1903  intvec v1 = 1,3,5,5,2;
1904  Hcode(v1);
1905  intvec v2 = 1,2,3,4,4,4,4,4,4,4,3,2,1;
1906  Hcode(v2);
1907}
1908
1909///////////////////////////////////////////////////////////////////////////////
1910static proc Cubic (poly f)
1911{
1912//---------------------------- initialisation ---------------------------------
1913  poly  f3;
1914  ideal Jf, Jf1, Jf2;
1915  int   Dim, Mult;
1916
1917  f3 = jet(f, 3);
1918  if( jet(f,2) != 0) { return("2-jet non zero"); }
1919  if( f3 == 0 ) { return("null form"); }
1920
1921  Jf1  = jacob(f3);
1922  Jf   = std(Jf1);
1923  Dim  = dim(Jf);
1924  Mult = mult(Jf);
1925
1926  if(Dim == 0) { return("P[8]:smooth cubic"); } // x3 + y3 + z3 + axyz
1927  if(Dim == 1) {
1928    if(Mult == 2) {
1929      Jf2  = wedge(jacob(Jf1),3-Dim), Jf1;
1930      Jf2  = std(Jf2);
1931      Dim  = dim(Jf2);
1932      if (Dim == 0) { return("R:conic + line"); }       // x3 + xyz
1933      if (Dim == 1) { return("Q:cuspidal cubic"); }  // x3 + yz2
1934    }
1935    if(Mult == 3) {
1936      Jf2 = wedge(jacob(Jf1),3-Dim), Jf1;
1937      Jf2 = std(Jf2);
1938      Dim = dim(Jf2);
1939      if(Dim == 0) { return("T:three lines"); } // xyz
1940      if(Dim == 1) { return("S:conic + tangent"); }     // x2z + yz2
1941    }
1942    if(Mult == 4) { return("U:three concurrent lines"); }       // x3 + xz2
1943  }
1944  if(Dim == 2) {
1945    if(Mult == 1) { return("V:doubleline + line"); }    // x2y
1946    if(Mult == 2) { return("V': tripple line"); }       // x3
1947  }
1948  if(Dim == 3) { return("P[9]:nodal cubic"); }  // x3 + y3 + xyz
1949
1950  return("");
1951}
1952
1953///////////////////////////////////////////////////////////////////////////////
1954static proc parity  (int e)
1955"USAGE:    parity()"
1956{
1957  int r = e/2;
1958  if( 2*r == e ) { return(0); }
1959  return(1);
1960}
1961
1962///////////////////////////////////////////////////////////////////////////////
1963static proc HKclass (intvec sg)
1964{
1965//---------------------------- initialisation ---------------------------------
1966  int    cnt = 0;
1967  string SG_Typ = "";
1968  list   v;
1969
1970  // if trace/debug mode not set, do it!
1971  init_debug();
1972  debug_log(1, "Milnor code : ", sg );
1973  if(size(sg) == 1) { v ="A["+string(sg[1])+"]", 1; return(v); }
1974  if(size(sg) == 3) { return(HKclass3(sg, SG_Typ, cnt)); }
1975  if(size(sg) == 5) { return(HKclass5(sg, SG_Typ, cnt)); }
1976  if(size(sg) == 7) { return(HKclass7(sg, SG_Typ, cnt)); }
1977  debug_log(1, "No solution found." );
1978  v = "", 0;
1979  return(v);
1980}
1981
1982///////////////////////////////////////////////////////////////////////////////
1983static proc HKclass3 (intvec sg, string SG_Typ, int cnt)
1984{
1985  list v;
1986
1987  if(sg[1] == 1) { v = HKclass3_teil_1(sg, SG_Typ, cnt); }
1988  debug_log(6, "HKclass3: ", v[1], " cnt=", v[2]);
1989  return(v);
1990}
1991
1992///////////////////////////////////////////////////////////////////////////////
1993static proc HKclass3_teil_1 (intvec sg, string SG_Typ, int cnt)
1994{
1995  int  k, r, s;
1996  list v;
1997
1998  debug_log(2, "entering HKclass3_teil_1", sg);
1999  if(sg[2]==1) { SG_Typ=SG_Typ+" D[k]=D["+string(sg[3]+3)+"]";cnt++;} // D[k]
2000  if(sg[2]>=1) {
2001    if( parity(sg[2])) { // sg[2] ist ungerade
2002      if(sg[2]<=sg[3]) {
2003        k = (sg[2]+1)/2;
2004        if(k>1) {
2005          cnt++;
2006          SG_Typ=SG_Typ+" J[k,r]=J["+string(k)+","+string(sg[3]+1-2*k)+"]";
2007        }// J[k,r]
2008      }
2009      if(sg[2]==sg[3]+2) {                              // E[6k+2]
2010        k = (sg[2]-1)/2;
2011        if(k>0) {cnt++; SG_Typ=SG_Typ+" E[6k+2]=E[" + string(6*k+2) + "]"; }
2012      }
2013    }
2014    else {              // sg[2] ist gerade
2015      if( sg[2] == sg[3]+1) {                           // E[6k]
2016        k = sg[2]/2; cnt++; SG_Typ=SG_Typ + " E[6k]=E[" + string(6*k) + "]"; }
2017      if( sg[2] == sg[3]) {                             // E[6k+1]
2018        k=sg[2]/2; cnt++; SG_Typ=SG_Typ+" E[6k+1]=E["+string(6*k+1)+"]"; }
2019    }
2020  }
2021
2022  debug_log(2, "finishing HKclass3_teil_1");
2023  debug_log(6, "HKclass3: ", SG_Typ, " cnt=", cnt);
2024  v = SG_Typ, cnt;
2025  return(v);
2026}
2027
2028///////////////////////////////////////////////////////////////////////////////
2029static proc HKclass5 (intvec sg, string SG_Typ, int cnt)
2030{
2031  list v;
2032
2033  if(sg[1] == 1 && sg[2] == 1) { v = HKclass5_teil_1(sg, SG_Typ,cnt); }
2034  if(sg[1] == 1 && sg[2] == 0) { v = HKclass5_teil_2(sg, SG_Typ,cnt); }
2035  debug_log(6, "HKclass5: ", v[1], " cnt=", v[2]);
2036  return(v);
2037}
2038
2039///////////////////////////////////////////////////////////////////////////////
2040static proc HKclass5_teil_1 (intvec sg, string SG_Typ, int cnt)
2041{
2042  int  k, r, s;
2043  list v;
2044
2045  debug_log(2, "entering HKclass5_teil_1", sg);
2046  if(parity(sg[3])) {  // Dritte Stelle soll ungerade sein
2047    k = (sg[3]+1)/2;
2048    if(sg[3] > sg[4]) {
2049      k--;
2050      if( (sg[4]==sg[5]) && (sg[3] == sg[4]+1) && k>0 ) { // W[12k+6]
2051        SG_Typ = SG_Typ + " W[12k+6]=W["+string(12*k+6)+"]"; cnt++;
2052      }
2053      if( (sg[3]==sg[5]) && (sg[3] == sg[4]+2) && k>0 ) { // W[12k+5]
2054        SG_Typ = SG_Typ + " W[12k+5]=W["+string(12*k+5)+"]"; cnt++;
2055      }
2056    }
2057    else {  // sg[3] <= sg[4]
2058      if( (sg[3]==sg[4]) && (sg[5] >= sg[3]) ) {
2059        r = sg[5] - sg[4];
2060        SG_Typ=SG_Typ +" X[k,r]=X["+string(k)+","+string(r)+"]"; cnt++;
2061      }
2062      if( (sg[3]==1) && (sg[4]==3) && (sg[5]>=sg[4])){    // Z[1,r]
2063        r = sg[5] - sg[4];
2064        SG_Typ = SG_Typ + " Z[1,r]=Z[1,"+string(r)+"]"; cnt++;
2065      }
2066
2067      if( sg[4] == sg[5]) {
2068        if(parity(sg[4])) {                                  // Z[k,r,0]
2069          r = (sg[4] - sg[3])/2;
2070          if( r>0 ) { cnt++;
2071            SG_Typ = SG_Typ + " Z[k,r,0]=Z["+string(k)+","+string(r)+",0]";
2072          }
2073        }
2074        else {                                                // Z[k,12k+6r]
2075          r = (sg[4] - 2*k)/2; cnt++;
2076          SG_Typ = SG_Typ+" Z[k,12k+6r]=Z["+string(k)+","+string(12*k+6*r)+"]";
2077        }
2078      }
2079
2080      if( parity(sg[4]) ) {  // 4. Stelle ist ungerade
2081        if(sg[4] == sg[5]+2) {                              // Z[k,12k+6r+1]
2082          r = (sg[4]-2*k-1)/2; cnt++;
2083          SG_Typ=SG_Typ+" Z[k,12k+6r+1]=Z["+string(k)+",";
2084          SG_Typ=SG_Typ+string(12*k+6*r+1)+"]";
2085       }
2086       if( (sg[5]>sg[4]) && (sg[4]>sg[3]) ) {           // Z[k,r,s]
2087          r = (sg[4] - sg[3])/2; cnt++;
2088          s = sg[5] - sg[4];
2089          SG_Typ = SG_Typ + " Z[k,r,s]=";
2090          SG_Typ = SG_Typ + "Z["+string(k)+","+string(r)+","+string(s)+"]";
2091        }
2092      }
2093      else {  // 4. Stelle ist gerade
2094        if( sg[4] == sg[5]+1) {                             // Z[k,12k+6r-1]
2095          r = (sg[4] - 2*k)/2; cnt++;
2096          SG_Typ=SG_Typ+" Z[k,12k+6r-1]=Z["+string(k)+",";
2097          SG_Typ=SG_Typ+string(12*k+6*r-1)+"]";
2098        }
2099      }
2100
2101      if(sg[4]>sg[3]) {                                     // Y[k,r,s]
2102        r = sg[4] - sg[3];
2103        s = sg[5] - sg[3] + r;
2104        if( s<0 ) { s = -s; }
2105        SG_Typ = SG_Typ + " Y[k,r,s]="; cnt++;
2106        SG_Typ = SG_Typ + "Y["+string(k)+","+string(r)+","+string(s)+"]";
2107      }
2108    }
2109  }
2110  else {  // Dritte Stelle soll gerade sein
2111    k = sg[3]/2;
2112    // sortiere verschiedene W's
2113    if(k>0) {
2114      if( (sg[4]==2*k-1) && (sg[4]==sg[5]) ) {  // W[12k]
2115        SG_Typ = SG_Typ + " W[12k]=W["+string(12*k)+"]"; cnt++;
2116      }
2117      if( (sg[4]==2*k-1) && (sg[3]==sg[5]) ) {  // W[12k+1]
2118        SG_Typ = SG_Typ + " W[12k+1]=W["+string(12*k+1)+"]"; cnt++;
2119      }
2120      if( (sg[4]==2*k) && (sg[5]>=sg[4]) ) {    // W[k,r]
2121        r = sg[5] - sg[4];
2122        SG_Typ=SG_Typ+" W[k,r]=W["+string(k)+","+string(r)+"]"; cnt++;
2123      }
2124      if( (sg[5]==2*k-1) && (sg[4]>sg[3]) ) {  // W#[k,2r-1]
2125        r = sg[4] - sg[3]; cnt++;
2126        SG_Typ = SG_Typ + " W#[k,2r-1]=W["+string(k)+","+string(2*r-1)+"]";
2127      }
2128      if( (sg[5]==2*k) && (sg[4]>sg[3]) ) {  // W#[k,2r]
2129        r = sg[4] - sg[3]; cnt++;
2130        SG_Typ = SG_Typ + " W#[k,2r]=W["+string(k)+","+string(2*r)+"]";
2131      }
2132    }   // ENDIF k>0
2133  }
2134  debug_log(2, "finishing HKclass5_teil_1");
2135  debug_log(6, "HKclass5_teil_1: ", SG_Typ, " cnt=", cnt);
2136  v = SG_Typ, cnt;
2137  return(v);
2138}
2139
2140///////////////////////////////////////////////////////////////////////////////
2141static proc HKclass5_teil_2 (intvec sg, string SG_Typ, int cnt)
2142{
2143  int k, r, s;
2144  list v;
2145
2146  debug_log(2, "entering HKclass5_teil_2", sg);
2147  // finde T[p,q,r]
2148  k = sg[3] + 1;
2149  r = sg[4] + k;
2150  s = sg[5] + r - 1;
2151  if(k>2 && r>2 && s>2) {                               // T[k,r,s]
2152    cnt++;
2153    SG_Typ = SG_Typ + " T[k,r,s]=T["+string(k)+","+string(r)+","+string(s)+"]";
2154  }
2155
2156  // finde weitere Moeglicjkeiten.
2157  if(sg[3]==2) {  // Q[...]
2158    if(parity(sg[4])) { // 4. Stelle ist ungerade.
2159      if(sg[4]==sg[5]) {                                // Q[6k+4]
2160        k=(sg[4]+1)/2; cnt++; SG_Typ=SG_Typ+" Q[6k+4]=Q["+string(6*k+4)+"]";
2161      }
2162      if(sg[4]+1==sg[5]) {                      // Q[6k+5]
2163        k=sg[5]/2; cnt++; SG_Typ=SG_Typ+" Q[6k+5]=Q["+string(6*k+5)+"]";
2164      }
2165    }
2166    else { // 4. Stelle ist gerade.
2167      if(sg[4]==sg[5]+1) {                      // Q[6k+6]
2168        k=sg[4]/2; cnt++; SG_Typ=SG_Typ+" Q[6k+6]=Q["+string(6*k+6)+"]";
2169      }
2170      if(sg[4]<sg[5]) {                 // Q[k,r]
2171        k = (sg[4]+2)/2;
2172        if(k>=2) {
2173          r=sg[5]+1-2*k; cnt++;
2174          SG_Typ=SG_Typ+" Q[k,r]=Q["+string(k)+","+string(r)+"]";
2175        }
2176      }
2177    }
2178  }
2179  else {           // S[...]
2180    if(parity(sg[3])) {  // 3. Stelle ist ungerade.
2181      k = (sg[3]-1)/2;
2182      if(sg[3]==sg[4]+3 && sg[3]==sg[5]+2) {    // S[12k-1]
2183        cnt++; SG_Typ = SG_Typ + " S[12k-1]=S["+string(12*k-1)+"]";
2184      }
2185      if(sg[3]==sg[4]+3 && sg[3]==sg[5]+1) {    // s[12k]
2186        cnt++; SG_Typ = SG_Typ + " S[12k]=S["+string(12*k)+"]";
2187      }
2188      if(sg[3]==sg[4]+2 && sg[5]>=sg[4]+1) {    // S[k,r]
2189        r = sg[5] - 2*k; cnt++;
2190        SG_Typ = SG_Typ + " S[k,r]=S["+string(k)+","+string(r)+"]";
2191      }
2192      if(sg[3]==sg[5]+2 && sg[4]>=sg[5]) {              // S#[k,2r-1]
2193        r = sg[4] - 2*k + 1; cnt++;
2194        SG_Typ = SG_Typ + " S#[k,2r-1]=S#["+string(k)+","+string(2*r-1)+"]";
2195      }
2196      if(sg[3]==sg[5]+1 && sg[4]>=sg[5]) {              // S#[k,2r]
2197        r = sg[4] - 2*k + 1; cnt++;
2198        SG_Typ = SG_Typ + " S#[k,2r]=S#["+string(k)+","+string(2*r)+"]";
2199      }
2200    }
2201    else { // 3. Stelle ist gerade.
2202      if(sg[3]==sg[5]+1 && sg[5]==sg[4]+3) {    // S[12k+4]
2203        k = (sg[3]-2)/2; cnt++;
2204        SG_Typ = SG_Typ + " S[12k+4]=S["+string(12*k+4)+"]";
2205      }
2206      if(sg[3]==sg[5]+2 && sg[5]==sg[4]+1) {    // S[12k+5]
2207        k = (sg[3]-2)/2; cnt++;
2208        SG_Typ = SG_Typ + " S[12k+5]=S["+string(12*k+5)+"]";
2209      }
2210    }
2211  }
2212  debug_log(2, "finishing HKclass5_teil_2");
2213  debug_log(6, "HKclass5_teil_2: ", SG_Typ, " cnt=", cnt);
2214  v = SG_Typ, cnt;
2215  return(v);
2216}
2217
2218///////////////////////////////////////////////////////////////////////////////
2219static proc HKclass7 (intvec sg, string SG_Typ, int cnt)
2220{
2221  list v;
2222
2223  if(sg[1]==1 && sg[2]==0 && sg[3]==1) {
2224      v=HKclass7_teil_1(sg, SG_Typ, cnt);
2225  }
2226  else {
2227      v[1]="not in list";
2228      v[2]=0;
2229  }
2230  debug_log(6, "HKclass7: ", v[1], " cnt=", v[2]);
2231  return(v);
2232}
2233
2234///////////////////////////////////////////////////////////////////////////////
2235static proc HKclass7_teil_1 (intvec sg, string SG_Typ, int cnt)
2236{
2237  int k, r, s;
2238  list v;
2239
2240  debug_log(2, "entering HKclass7_teil_1", sg);
2241  if(sg[4] == 2) {                                      // V[...]
2242    if(sg[5] == 0 && sg[6] == 1 && sg[7]>0) {   // V[1,r]
2243      r = sg[7] - 1; cnt++; SG_Typ = SG_Typ + " V[1,r]=V[1,"+string(r)+"]";
2244    }
2245    if(sg[5] == 1 && sg[7] == 1) {                      // V#[1,2r-1]
2246      r=sg[6]+1; cnt++; SG_Typ=SG_Typ+" V#[1,2r-1]=V#[1,"+string(2*r-1)+"]";
2247    }
2248    if(sg[5] == 1 && sg[7] == 2) {                      // V#[1,2r]
2249      r=sg[6]+1; cnt++; SG_Typ=SG_Typ+" V#[1,2r]=V#[1,"+string(2*r)+"]";
2250    }
2251  }
2252  //            Moegliche U[...]'s
2253  k = sg[4];
2254  if(sg[5]==2*k-1 && sg[6]==0 && sg[7]==sg[5]) {        // U[12k]
2255    cnt++;SG_Typ = SG_Typ + " U[12k]=U["+string(12*k)+"]";
2256  }
2257  if(sg[5]==2*k && sg[6]==0 && sg[7]==sg[5]) {  // U[12k+4]
2258    cnt++;SG_Typ = SG_Typ + " U[12k+4]=U["+string(12*k+4)+"]";
2259  }
2260  if(sg[5]==2*k-1 && sg[6]>0 && sg[7]==sg[5]) { // U[k,2r-1]
2261    r=sg[6]-1; cnt++;
2262    SG_Typ=SG_Typ+" U[k,2r-1]=U["+string(k)+","+string(2*r-1)+"]";
2263  }
2264  if(sg[5]==2*k-1 && sg[6]>0 && sg[7]==2*k) {   // U[k,2r]
2265    r = sg[6]; cnt++;
2266    SG_Typ = SG_Typ + " U[k,2r]=U["+string(k)+","+string(2*r)+"]";
2267  }
2268  debug_log(2, "finishing HKclass7_teil_1");
2269  debug_log(6, "HKclass7_teil_1: ", SG_Typ, " cnt=", cnt);
2270  v = SG_Typ, cnt;
2271  return(v);
2272}
2273
2274///////////////////////////////////////////////////////////////////////////////
2275proc singularity(string typ, list #)
2276"USAGE:    singularity(t, l); t=string (name of singularity),
2277          l=list of integers/polynomials (indices/parmeters of singularity)
2278COMPUTE:  get the singularity named by type t from the database.
2279          list l is as follows: @*
2280          l= k [,r [,s [,a [,b [,c [,d]..]: k,r,s=int   a,b,c,d=poly. @*
2281          The name of the dbm-databasefile is: NFlist.[dir,pag].
2282          The file is found in the current directory. If it does not
2283          exist, please run the script MakeDBM first.
2284RETURN:   Normal form and corank of the singularity named by type t and its
2285          index (indices) l.
2286EXAMPLE:  example singularity; shows an example"
2287{
2288  poly a1, a2, a3, a4, f;
2289  int k, r, s;
2290  int len = size(#);
2291  list v, ret;
2292
2293  classify_init();
2294  ret = 0, 0;
2295  k = #[1];
2296  if(len>=2) { r = #[2]; }
2297  else { r = 0; }
2298  if(len>=3) { s = #[3]; }
2299  else { s = 0; }
2300  if( k<0 || r<0 || s<0) {
2301    "Initial condition failed: k>=0; r>=0; s>=0";
2302    "k="+string(k)+" r="+string(r)+"   s="+string(s);
2303    return(ret);
2304  }
2305  int crk;
2306
2307  init_debug();
2308  def ring_top=basering;
2309
2310  if(len>=4) { a1 = #[4]; }
2311  else { a1=1; }
2312  if(len>=5) { a2 = #[5]; }
2313  else { a2=1; }
2314  if(len>=6) { a3 = #[6]; }
2315  else { a3=1; }
2316  if(len>=7) { a4 = #[7]; }
2317  else { a4=1; }
2318
2319  debug_log(4, "Values: len=", string(len), " k=", string(k), " r=",
2320        string(r));
2321  if(defined(RingNF) != 0 ) { kill RingNF; }
2322  ring RingNF=char(basering),(x,y,z),(c,ds);
2323  poly f;
2324  map Conv=ring_top,maxideal(1);
2325  v = Singularitaet(typ, k, r, s, Conv(a1), Conv(a2), Conv(a4), Conv(a4));
2326  f, crk = v[1..2];
2327  debug_log(2, "Info=", f );
2328  setring ring_top;
2329  if(defined(Phi) != 0 ) { kill Phi; }
2330  map Phi=RingNF,maxideal(1);
2331
2332  ret = Phi(f), crk;
2333  return(ret);
2334}
2335example
2336{ "EXAMPLE"; echo=2;
2337  ring r=0,(x,y,z),(c,ds);
2338  init_debug(0);
2339  singularity("E[6k]",6);
2340  singularity("T[k,r,s]", 3, 7, 5);
2341  poly f=y;
2342  singularity("J[k,r]", 4, 0, 0, f);
2343}
2344
2345///////////////////////////////////////////////////////////////////////////////
2346static proc Singularitaet (string typ,int k,int r,int s,poly a,poly b,poly c,poly d)
2347{
2348  list   v;
2349  string DBMPATH=system("getenv","DBMPATH");
2350  string DatabasePath, Database, S, Text, Tp;
2351  poly   f, f1;
2352  int    crk, Mu, ret;
2353  intvec MlnCd;
2354
2355  if( DBMPATH != "" ) { DatabasePath = DBMPATH+"/NFlist"; }
2356  else { DatabasePath = "NFlist"; }
2357  Database="DBM: ",DatabasePath;
2358
2359  link dbmLink=Database;
2360  debug_log(2, "Opening Singalarity-database: ", newline, Database);
2361  Tp = read(dbmLink, typ);
2362  debug_log(2,"DBMread(", typ, ")=", Tp, ".");
2363  if( Tp != "(null)" && Tp !="" ) {
2364    string Key = "I_", typ;
2365    S = "f = ", Tp, ";";
2366    debug_log(2,"S=", S, " Tp=", Tp, "Key=", Key);
2367    execute(S);
2368    execute(read(dbmLink, Key)+";");
2369    debug_log(1, "Polynom f=", f,  "  crk=", crk, "  Mu=", Mu,
2370                " MlnCd=", MlnCd);
2371    v = f, crk, Mu, MlnCd;
2372  }
2373  else {
2374    v = 0, 0, 0, 0;
2375  }
2376  close(dbmLink);
2377  return(v);
2378}
2379
2380///////////////////////////////////////////////////////////////////////////////
2381proc RandomPolyK (int M, string Typ)
2382"USAGE:    RandomPolyK(M, Typ)"
2383{
2384//---------------------------- initialisation ---------------------------------
2385  int    n, b, i, k, r, s, crk;
2386  ideal  B;
2387  map    Phi;
2388  string txt, Tp;
2389  list   v;
2390
2391  def ring_ext=basering;
2392  n=4;
2393  if(M<5) { M = 5; }
2394
2395  k = random(1, M);
2396  r = random(-5, 2*M);
2397  s = random(-5, 2*M);
2398  if(r<0) { r = 0; }
2399  if(s<0) { s = 0; }
2400
2401  ring RgAnf=char(basering),(x,y,z,t),(c,ds);
2402  poly f;
2403
2404  v = singularity(Typ, k, r, s);
2405  f, crk = v[1..2];
2406//  f = f +t2;
2407  if(crk==1) { f = f + y2 + z2; }
2408  if(crk==2) { f = f + z2; }
2409  txt="RandomPoly-Series: gewaehlt fall "+Typ+" mit";
2410  txt=txt+" f="+string(f);
2411  txt;
2412  setring ring_ext;
2413  B = maxideal(1);
2414
2415  r=1;
2416  for(i=n; i>0; i--,r++) {
2417//  for(i=1; i<=n; i=i+1)
2418    B[rvar(x(r))] = x(i);
2419    if(i>2 && random(1,10)<3) { B[rvar(x(r))] = B[rvar(x(r))] + x(i-1); }
2420//    if(i==1 && random(1,10)<4) { B[rvar(x(r))] = B[rvar(x(r))]- x(n); }
2421    if(i>0) {
2422      for(b=3; b<5; b=b+1) {
2423        // B[rvar(x(r))] = B[rvar(x(r))] + random(0,9) * x(i)^(b+2);
2424        if(random(1,20)<3) {
2425          B[rvar(x(r))] = B[rvar(x(r))] - random(-2,2)*x(b)^2;
2426        }
2427      }
2428    }
2429  }
2430  Phi=RgAnf, B;
2431  Phi;
2432  poly fr=Phi(f);
2433  fr = fr+(x(1)+x(2))^2;
2434//  return(Phi(f));
2435  return(fr);
2436}
2437
2438///////////////////////////////////////////////////////////////////////////////
2439proc debug_log (int level, list #)
2440"USAGE:    debug_log(level,li); level=int, li=comma separated \"message\" list
2441COMPUTE:  print \"messages\" if level>=@DeBug.
2442          useful for user-defined trace messages.
2443EXAMPLE:  example debug_log; shows an example
2444SEE ALSO: init_debug
2445"
2446{
2447   int len = size(#);
2448//   int printresult = printlevel - level +1;
2449//   if(level>1) {
2450//     dbprint(printresult, "Debug:("+ string(level)+ "): ", #[2..len]);
2451//   }
2452//   else { dbprint(printresult, #[1..len]); }
2453   if( defined(@DeBug) == 0 ) { init_debug(); }
2454   if(@DeBug>=level) {
2455      if(level>1) { "Debug:("+ string(level)+ "): ", #[1..len]; }
2456      else { #[1..len]; }
2457   }
2458}
2459example
2460{ "EXAMPLE:"; echo=2;
2461  example init_debug;
2462}
2463
2464///////////////////////////////////////////////////////////////////////////////
2465proc init_debug(list #)
2466"USAGE:    init_debug([level]);  level=int
2467COMPUTE:  Set the global variable @DeBug to level. The variable @DeBug is
2468          used by the function debug_log(level, list of strings) to know
2469          when to print the list of strings. init_debug() reports only
2470          changes of @DeBug.
2471NOTE:     The procedure init_debug(n); is usefull as trace-mode. n may
2472          range from 0 to 10, higher values of n give more information.
2473EXAMPLE:  example init_debug; shows an example"
2474{
2475  int newDebug=0;
2476  if( defined(@DeBug) != 0 ) { newDebug = @DeBug; }
2477
2478  if( size(#) > 0 ) {
2479    newDebug=#[1];
2480  }
2481  else {
2482    string s=system("getenv", "SG_DEBUG");
2483    if( s != "" && defined(@DeBug)==0) {
2484      s="newDebug="+s;
2485      execute(s);
2486    }
2487  }
2488  if( defined(@DeBug) == 0) {
2489    int @DeBug = newDebug;
2490    export @DeBug;
2491    if(@DeBug>0) { "Debugging level is set to ", @DeBug; }
2492  }
2493  else {
2494    if( (size(#) == 0) && (newDebug < @DeBug) ) { return(); }
2495    if( @DeBug != newDebug) {
2496      int oldDebug = @DeBug;
2497      @DeBug = newDebug;
2498      if(@DeBug>0) { "Debugging level change from ", oldDebug, " to ",@DeBug; }
2499      else {
2500        if( @DeBug==0 && oldDebug>0 ) { "Debugging switched off."; }
2501      }
2502    }
2503  }
2504  printlevel = @DeBug;
2505}
2506example
2507{ "EXAMPLE:"; echo=2;
2508  init_debug();
2509  debug_log(1,"no trace information printed");
2510  init_debug(1);
2511  debug_log(1,"some trace information");
2512  init_debug(2);
2513  debug_log(2,"nice for debugging scripts");
2514  init_debug(0);
2515}
2516
2517///////////////////////////////////////////////////////////////////////////////
2518proc basicinvariants(poly f)
2519"USAGE:    basicinvariants(f);   f = poly
2520COMPUTE:  Compute basic invariants of f: an upper bound d for the
2521          determinacy, the milnor number mu and the corank c of f
2522RETURN:   intvec: d, mu, c
2523EXAMPLE:  example basicinvariants; shows an example"
2524{
2525  intvec v;
2526  ideal Jfs = std(jacob(f));
2527  v[1] = system("HC")+1;
2528  v[2] = vdim(Jfs);
2529  v[3] = corank(f);
2530  if( v[2]<v[1] ) { v[1] = v[2]+1; }
2531  return(v);
2532}
2533example
2534{ "EXAMPLE:"; echo=2;
2535   ring r=0,(x,y,z),ds;
2536   basicinvariants((x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3);
2537}
2538
2539///////////////////////////////////////////////////////////////////////////////
2540proc corank(poly f)
2541"USAGE:    corank(f);   f=poly
2542RETURN:   the corank of the Hessian matrix of f, of type int
2543REMARK:   corank(f) is the number of variables occuring in the residual
2544          singularity after applying 'morsesplit' to f
2545EXAMPLE:  example corank; shows an example"
2546{
2547  matrix M = jacob(jacob(jet(f,2)));
2548  list lba = bareiss(M);
2549  int cr = nvars(basering) - size(module(lba[1]));
2550  return(cr);
2551}
2552example
2553{ "EXAMPLE:"; echo=2;
2554  ring r=0,(x,y,z),ds;
2555  poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
2556  corank(f);
2557}
2558///////////////////////////////////////////////////////////////////////////////
2559static proc Faktorisiere(poly f, poly fk, int pt, int k, intvec RFlg)
2560{
2561//---------------------------- initialisation ---------------------------------
2562  poly   a, b, Relation;
2563  ideal  B, Jfsyz;
2564  map    PhiG, VERT;
2565  matrix Mat;
2566  list   v;
2567  def    ring_top=basering;
2568
2569  // Ziel: bestimme a,b sodass  fk = (ax+by^k)^pt gilt.
2570  B    = maxideal(1);
2571  PhiG = ring_top,B;
2572  debug_log(2, "Faktor: f=",Show(f)," Jet=",Show(fk)," k=",k," exp=",pt);
2573
2574//----------------------- compute role of x and y -----------------------------
2575  Jfsyz = fk, diff(fk, x(1));
2576  Mat   = matrix(syz(Jfsyz));
2577  if( (fk-subst(fk,x(1),0)) != 0  &&  (fk-subst(fk,x(2),0)) != 0 ) {
2578    // Wenn k>0 ist die Wahl fuer x & y bereits getroffen
2579    // sonst bestimmen x und y
2580    Jfsyz    = fk, diff(fk, x(1));
2581    Mat      = matrix(syz(Jfsyz));
2582    Relation = -pt * Mat[2,1] / Mat[1,1];
2583    a = Coeff(Relation, x(1), x(1));
2584    b = Coeff(Relation, x(2), x(2)^k);
2585    B = maxideal(1);
2586    if( (RFlg[1]==1 && k==1) || k>1) { B[rvar(x(1))] = x(1)-b*x(2)^k; }
2587    else { B[rvar(x(2))] = x(2)-b*x(1)^k; }
2588    VERT = basering,B;
2589    f    = VERT(f);
2590    PhiG = VERT(PhiG);
2591  }
2592
2593//------------------- permutation of x and y, if needed -----------------------
2594  if( k==1 ) {
2595    debug_log(2, "Fak-7:",Show(f)," jet=",Show(fk));
2596    if(Coeff(jet(f, pt), x(1), x(1)^pt) == 0) {
2597      VERT = basering,x(2),x(1);
2598      f    = VERT(f);
2599      PhiG = VERT(PhiG);
2600    }
2601  }
2602  debug_log(2, "Fak-8:",Show(f)," jet=",Show(fk));
2603  debug_log(6, "Faktorisiere liefert: f=", Show(f));
2604  v[1] = f;
2605  v[2] = PhiG;
2606  return(v);
2607}
2608
2609///////////////////////////////////////////////////////////////////////////////
2610static proc Teile(poly f, poly fk)
2611{
2612  ideal  Jfsyz = f, fk;
2613  poly   Relation;
2614  matrix Mat = matrix(syz(Jfsyz));
2615  Relation   = -1 * Mat[2,1]/Mat[1,1];
2616  return(Relation);
2617}
2618
2619///////////////////////////////////////////////////////////////////////////////
2620static proc GetRf (poly fi, int n)
2621"USAGE:    GetRf();"
2622{
2623//---------------------------- initialisation ---------------------------------
2624  int    j, k, l1, l1w;
2625  matrix Koef;
2626  intvec RFlg;
2627
2628  RFlg[n] = 0;
2629  intvec Haeufigkeit = RFlg;
2630
2631  for( j=1; j<=n ; j=j+1) {
2632    Koef = coef(fi, x(j));
2633    Haeufigkeit[j] = ncols(Koef);
2634    if(Coeff(fi, x(j),0) == 0) { Haeufigkeit[j] = Haeufigkeit[j] + 1;}
2635  }
2636  for( j=n; j>0 ; j=j-1) {
2637    l1  = 0;
2638    l1w = 0;
2639    for(k=1;k<=n;k=k+1) { if(Haeufigkeit[k]>l1w) { l1=k; l1w=Haeufigkeit[k];}}
2640    RFlg[j]        = l1;
2641    Haeufigkeit[l1] = 0;
2642  }
2643  debug_log(2, "Permutations:", RFlg);
2644  return(RFlg);
2645}
2646
2647///////////////////////////////////////////////////////////////////////////////
2648static proc Show(poly g)
2649{
2650  string s;
2651  def ring_save=basering;
2652
2653  execute(@ringdisplay);
2654  map showpoly=ring_save,maxideal(1);
2655  s = string(showpoly(g));
2656  setring ring_save;
2657  return (s);
2658}
2659
2660///////////////////////////////////////////////////////////////////////////////
2661static proc checkring
2662{
2663  int CH = char(basering);
2664  if(CH >= 2 && CH<=13) {
2665    "Ring has characteristic ",CH;
2666    "Characteristic other than 0 or 0<char<13 is not yet implemented";
2667    return(1);
2668  }
2669  return(0);  // characteristic of ring is OK, return (0)
2670}
2671
2672///////////////////////////////////////////////////////////////////////////////
2673static proc DecodeNormalFormString (string S_in)
2674"USAGE:    DecodeNormalFormString"
2675{
2676//---------------------------- initialisation ---------------------------------
2677  int    C_eq, a, b, i, t, k, r, s;
2678  string s1, s2, s3, s4, s_in, Typ;
2679  list v = "Error", 0, 0, 0;
2680
2681  C_eq = find(S_in, "=")+1;
2682  s_in = S_in[C_eq,30];
2683  debug_log(2, "Decode:");
2684
2685  debug_log(2, "S_in=", S_in,"  s_in=",s_in );
2686  a = find(s_in, "[")+1;
2687  b = find(s_in, "]")-1;
2688  t = 1;
2689  k = 0;
2690  r = 0;
2691  s = 0;
2692
2693  if(a<0 || b<0) { return(v); }
2694  Typ = s_in[1..a-1];
2695  s1  = s_in[a..b];
2696  debug_log(6, "Suche Type:", Typ);
2697  //---------------------- decode between brackets ----------------------------
2698  if( find(s1, ",") == 0) {
2699    debug_log(8, "  Number of columns: 0");
2700    s2 = "k = "+s1+";";
2701    execute(s2);
2702    if( (Typ=="A[") || (Typ=="D[") ) { s3 = "k"; }
2703    if( Typ == "E[") { t = 6; }
2704    if( Typ == "W[") { t = 12; }
2705    if( Typ == "Q[") { t = 6; }
2706    if( Typ == "Z[") { t = 6; }
2707    if( Typ == "U[") { t = 12; }
2708    if( t > 1 ) {
2709      i = k;
2710      k = k/t;
2711      b = i - t*k;
2712      if( (s1 == "Q[") && (b==0) ) { k=k-1; b=6; }
2713      if(Typ == "Z[") {
2714        if(b==0) { k=k-1; b=6; }
2715        if(b==1) { k=k-1; b=7; }
2716      }
2717      if( b == 0 ) { s3 = string(t)+"k"; }
2718      else { s3 = string(t)+"k+"+string(b); }
2719    }
2720    if( Typ == "S[") {
2721      i = k+1;
2722      k = i/12;
2723      b = i - 12*k;
2724      if( b == 1 ) { s3 = "k"; }
2725      else {
2726        if(b==0) { s3 = "12k"+string(b-1); }
2727        else { s3 = "12k+"+string(b-1); }
2728      }
2729    }
2730    s2 = Typ + s3 +"]";
2731  }  // es kommt mindestens ein komma vor...
2732  //----------------------- more than 1 paramater -----------------------------
2733  else {
2734    b  = find(s1, ",");
2735    s2 = "k = ",s1[1..b-1],";";
2736    execute(s2);
2737    s1 = s1[b+1..size(s1)];
2738    if(find(s1, ",") == 0) {
2739      debug_log(8, "  Number of columns 1");
2740      s2 = "r = "+s1+";";
2741      execute(s2);
2742      s4 = "r";
2743      s3 = "k";
2744      if(r==0) { s4 = string(0); }
2745      if(k==0 && Typ=="Z[") { s3 = string(1); }
2746      if(Typ[2] == "#") {
2747        i = r+1;
2748        r = i/2;
2749        b = i - 2*r;
2750        if( b == 1 ) { s4 = "2r"; }
2751        else { s4 = "2r-1"; }
2752      }
2753      s2 = Typ + s3 + "," + s4 +"]";
2754    }  // es kommt mindestens zwei komma vor...
2755    //----------------------- get third parameter -----------------------------
2756    else {
2757      debug_log(8, "  Number of columns >=2");
2758      debug_log(2, "Y[k,r,s] / Z[k,r,s] / T[k,r,s]");
2759      b  = find(s1, ",");
2760      s2 = "r = ",s1[1..b-1],";";
2761      execute(s2);
2762      s2 = "s = ",s1[b+1..size(s1)],";";
2763      execute(s2);
2764      if(Typ=="Y[") { s2 = "Y[k,r,s]"; }
2765      if(Typ=="Z[") { s2 = "Z[k,r,s]"; }
2766      if(Typ=="T[") { s2 = "T[k,r,s]"; }
2767    }
2768  }
2769  debug_log(2, "Looking for Normalform of ",s2,"with (k,r,s) = (",
2770        k,",",r,",", s, ")" );
2771  v = s2, k, r, s;
2772  return(v);
2773}
2774
2775///////////////////////////////////////////////////////////////////////////////
2776proc A_L
2777"USAGE:    A_L(f);  f poly
2778          A_L(s);  s string, the name of the singularity
2779COMPUTE:  the normal form of f in Arnold's list of singularities in case 1,
2780          in case 2 nothing has to be computed.
2781RETURN:   A_L(f): compute via 'milnorcode' the class of f and return the normal
2782          form of f found in the database.
2783          A_L(\"name\"): get the normal form from the database for the
2784          singularity given by its name.
2785EXAMPLE:  example A_L; shows an example"
2786{
2787  // if trace/debug mode not set, do it!
2788  init_debug();
2789
2790  if( typeof(#[1]) == "string" ) {
2791    if(checkring()) { return(#[1]); }
2792    return(normalform(#[1]));
2793  }
2794  if( typeof(#[1]) == "poly" ) {
2795    if(checkring()) { return(#[1]); }
2796    return(quickclass(#[1]));
2797  }
2798
2799}
2800example
2801{ "EXAMPLE:"; echo=2;
2802  ring r=0,(a,b,c),ds;
2803  poly f=A_L("E[13]");
2804  f;
2805  A_L(f);
2806}
2807
2808///////////////////////////////////////////////////////////////////////////////
2809proc normalform(string s_in)
2810"USAGE:    normalform(s);  s=string
2811RETURN:   Arnold's normal form of singularity with name s
2812EXAMPLE:  example normalform; shows an example."
2813{
2814  string Typ;
2815  int    k, r, s, crk;
2816  int    n, i;
2817  poly   f;
2818  list   v;
2819  def ring_ext = basering;
2820  n = nvars(basering);
2821  ring ring_top=char(basering),(x(1..n)),(c,ds);
2822
2823  if(checkring()) { return(s_in); }
2824  if(nvars(basering)<=1) {
2825    "We need at least 2 variables in basering, you have",nvars(basering),".";
2826    return();
2827  }
2828  // if trace/debug mode not set, do it!
2829  init_debug();
2830
2831  v = DecodeNormalFormString(s_in);
2832  Typ, k, r, s = v[1..4];
2833  if(Typ=="Error") { return(0); }
2834  v = singularity(Typ, k, r, s);
2835  poly f_out;
2836  f_out, crk = v[1..2];
2837  if(crk>1) { for(i=crk+1;i<=n;i=i+1) { f_out = f_out + x(i)^2; } }
2838  setring ring_ext;
2839  map conv_top2ext=ring_top,maxideal(1);
2840  f = conv_top2ext(f_out);
2841//  f, crk = v[1..2];
2842  return(f);
2843}
2844example
2845{ "EXAMPLE:"; echo=2;
2846  ring r=0,(a,b,c),ds;
2847  normalform("E[13]");
2848}
2849
2850///////////////////////////////////////////////////////////////////////////////
2851proc swap
2852"USAGE:    swap(a,b);
2853RETURN:   b,a if a,b is the input (any type)"
2854{
2855  return(#[2],#[1]);
2856}
2857example
2858{ "EXAMPLE:"; echo=2;
2859  swap("variable1","variable2");
2860}
2861
2862///////////////////////////////////////////////////////////////////////////////
2863proc Setring(int c, string name)
2864"USAGE:    "
2865{
2866  string s="ring "+name+"=0,(x(1.."+ string(c) +")),(c,ds);";
2867  return(s);
2868}
2869
2870///////////////////////////////////////////////////////////////////////////////
2871proc internalfunctions()
2872"USAGE:   internalfunctions();
2873RETURN:  nothing, display names of internal procedures of classify.lib
2874EXAMPLE: no example"
2875{ "   Internal functions for the classification using Arnold's method,";
2876 "   the function numbers correspond to numbers in Arnold's classifier:";
2877 "Klassifiziere(poly f);      //determine the type of the singularity f
2878  Funktion1bis (poly f, list cstn)
2879  Funktion3 (poly f, list cstn)
2880  Funktion6 (poly f, list cstn)
2881  Funktion13 (poly f, list cstn)
2882  Funktion17 (poly f, list cstn)
2883  Funktion25 (poly f, list cstn)
2884  Funktion40 (poly f, list cstn, int k)
2885  Funktion47 (poly f, list cstn)
2886  Funktion50 (poly f, list cstn)
2887  Funktion58 (poly fin, list cstn)
2888  Funktion59 (poly f, list cstn)
2889  Funktion66 (poly f, list cstn)
2890  Funktion82 (poly f, list cstn)
2891  Funktion83 (poly f, list cstn)
2892  Funktion91 (poly f, list cstn, int k)
2893  Funktion92 (poly f, list cstn, int k)
2894  Funktion93 (poly f, list cstn, int k)
2895  Funktion94 (poly f, list cstn, int k)
2896  Funktion95 (poly f, list cstn, int k)
2897  Funktion96 (poly f, list cstn, int k)
2898  Funktion97 (poly f, list cstn)
2899  Isomorphie_s82_x (poly f, poly fk, int k)
2900  Isomorphie_s82_z (poly f, poly fk, int k)
2901  Isomorphie_s17 (poly f, poly fk, int k, int ct)
2902  printresult (string f,string typ,int Mu,int m,int corank,int K)
2903  ";
2904  "   Internal functions for the classifcation by invariants:
2905  Cubic (poly f)
2906  parity (int e)             //return the parity of e
2907  HKclass (intvec i)
2908  HKclass3( intvec i, string SG_Typ, int cnt)
2909  HKclass3_teil_1 (intvec i, string SG_Typ, int cnt)
2910  HKclass5 (intvec i, string SG_Typ, int cnt)
2911  HKclass5_teil_1 (intvec i, string SG_Typ, int cnt)
2912  HKclass5_teil_2 (intvec i, string SG_Typ, int cnt)
2913  HKclass7 (intvec i, string SG_Typ, int cnt)
2914  HKclass7_teil_1 (intvec i, string SG_Typ, int cnt)
2915  ";
2916  "   Internal functions for the Morse-splitting lemma:
2917  Morse(poly fi, int K, int corank)  //splitting lemma itself
2918  Coeffs (list #)
2919  Coeff
2920  ";
2921  "   Internal functions providing tools:
2922  ReOrder(poly f)
2923  Singularitaet(string typ,int k,int r,int s,poly a,poly b,poly c,poly d)
2924  RandomPolyK
2925  Faktorisiere(poly f, poly g, int p, int k)   compute g = (ax+by^k)^p
2926  Teile(poly f, poly g);             //divides f by g
2927  GetRf(poly f, int n);
2928  Show(poly f);
2929  checkring();
2930  DecodeNormalFormString(string s);
2931  Setring(int n, string ringname);
2932  ";
2933}
2934example
2935{
2936  "EXAMPLE"; echo=2;
2937  internalfunctions();
2938}
2939
2940///////////////////////////////////////////////////////////////////////////////
2941// E n d   O f   F i l e
2942
Note: See TracBrowser for help on using the repository browser.