source: git/Singular/LIB/classify.lib @ 5011fd

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