# source:git/Singular/LIB/classify.lib@5c8c19

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