source: git/Singular/LIB/classify.lib @ a42239

// $Id: classify.lib,v 1.20 1998-04-03 23:18:15 krueger Exp $
///////////////////////////////////////////////////////////////////////////////

LIBRARY:  classify.lib  PROCEDURES FOR THE ARNOLD-CLASSIFIER OF SINGULARITIES  

   A library for classifying isolated hypersurface singularities w.r.t. right
   equivalence, based on the determinator of singularities by V.I. Arnold.
   Author: Kai Krueger, krueger@mathematik.uni-kl.de
   last modified: 08.11.1997

basicinvariants(f);  computes Milnor number, determinacy-bound and corank of f
classify(f);         normal form of poly f determined with Arnold's method
corank(f);           computes the corank of f (i.e. of the Hessian of f) 
Hcode(v);            coding of intvec v acoording to the number repetitions 
init_debug([n]);     print trace and debugging information depending on int n  
internalfunctions(); display names of internal procedures of this library 
milnorcode(f[,e]);   Hilbert poly of [e-th] Milnor algebra coded with Hcode
morsesplit(f);       residual part of f after applying the splitting lemma
quickclass(f)        normal form of f determined by invariants (milnorcode)
singularity(s,[]);   normal form of singularity given by its name s and index
swap (a,b);          returns b,a
tschirnhaus(f,v);    Tschirnhaus transformation of f w.r.t. variable v
A_L(s/f)              shortcut for quickclass(f) or normalform(s)
normalform(s);       normal form of singularity given by its name s
debug_log (lev,[])   print trace and debugging information w.r.t level>@DeBug
           (parameters in square brackets [] are optional)

LIB "inout.lib";
LIB "elim.lib";
LIB "sing.lib";

///////////////////////////////////////////////////////////////////////////////
proc classify_init
{
  string s;
  link l="DBM:r NFlist";
  s = read(l,"VERSION");
  if (s == "" ) {
    "(classify.lib): Need to create database...";
    LIB"makedbm.lib";
    create_sing_dbm();
  }
  close(l);
  l="DBM:r NFlist";
  s = read(l,"VERSION");
  //"(classify.lib): Creation done. Current version:", s;
}
///////////////////////////////////////////////////////////////////////////////

proc classify (poly f_in)
USAGE:    classify(f);  f=poly
COMPUTE:  normal form and singularity type of f with respect to right
          equivalence, as given in the book "Singularities of 
          differentiables maps, Volume I" by V.I. Arnold, S.M. Gusein-Zade,
          A.N. Varchenko
RETURN:   normal form of f, of type poly
REMARK:   This version of classify is only alpha. Please send bugs and
          comments to: "Kai Krueger" <krueger@mathematik.uni-kl.de>
          Be sure to have at least Singular version 0.9.3, better 1.0.1
          Updates can be found under:
          URL=http://www.mathematik.uni-kl.de/~krueger/Singular/
NOTE:     type init_debug(n); (0 <= n <= 10) in order to get intermediate 
          information, higher values of n give more information.
          The proc creates several global objects with names all starting
          with @, hence there should be no name conflicts 
EXAMPLE:  example classify; shows an example
{
//---------------------------- initialisation ---------------------------------
  poly   f_out;
  int    n, i, corank, show_nf;
  string s2;
  list   v;
  def ring_ext = basering; 

  init_debug();                    // initialize trace/debug mode
  if(checkring()) { return(f_in); }
  n = nvars(basering);
  show_nf  = 1;                // return normal form if set to '1'

  // define new ring
  ring ring_top=char(basering),(x(1..n)),(c,ds);

  map conv_ext2top=ring_ext,maxideal(1);

  if(defined(@ringdisplay)!=0) { kill @ringdisplay; }
  string @ringdisplay = "ring_ext";
  export @ringdisplay;

  v = Klassifiziere(conv_ext2top(f_in));
  poly f_out = v[1];
  s2 = v[2];          // s2: Typ des Polynoms f z.b: E[18]
  corank = v[3];

//---------------- collect results and create return-value --------------------
  if( s2=="error!" || s2=="NoClass") { 
      setring ring_ext;
      return(f_in);
  }

  if(show_nf==1) {
    poly f_nf = normalform(s2);
    for(i=corank+1;i<=n;i=i+1) { f_nf = f_nf + x(i)^2; }
    debug_log(2, "Normal form NF(f)=", f_nf);
  }
  if(corank>1) { for(i=corank+1;i<=n;i=i+1) { f_out = f_out + x(i)^2; } }
  setring ring_ext;
  map conv_top2ext=ring_top,maxideal(1);

  if(show_nf == 1) { return(conv_top2ext(f_nf)); }
  else { return(conv_top2ext(f_out)); }
}
example
{"EXAMPLE"; echo=2;
   ring r=0,(x,y,z),ds;
   poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
   classify(f);
   init_debug(3);
   classify(f);
}

///////////////////////////////////////////////////////////////////////////////
proc Klassifiziere (poly f)
{ 
//--------------------------- initialisation ----------------------------------
    string s1;
    int    n, cnt, corank_f, K, Mu;
    list   v, cstn;
    map    PhiG;
    def ring_top = basering; 

    n = nvars(basering);    // Zahl der Variablen des aktuellen Rings.
    PhiG = ring_top, maxideal(1);
    cstn[4] = PhiG;
    if( defined(@ringdisplay) == 0) { 
      string @ringdisplay;               // Define always 'ringdisplay' to be
      export @ringdisplay;               // able to run 'Show(f)'
    }
    @ringdisplay = "setring RingB;";
    if(defined(RingB)!=0) { kill RingB; }
    execute ("ring RingB="+charstr(basering)+",("+A_Z("x", n)+"),(c,ds);");
    export RingB;
    setring ring_top;

//---------------------- compute basciinvariants ------------------------------
    if(jet(f,0) != 0 ) { 
      cstn[1] = corank(f); cstn[2]=-1; cstn[3]=1;
      return(printresult(1, f, "a unit", cstn, -1));
    }

    debug_log(1, "Computing Basicinvariants of f ...");
    K, Mu, corank_f = basicinvariants(f);
    debug_log(0, "About the singularity :");
    debug_log(0, "          Milnor number(f)   = "+string(Mu));
    debug_log(0, "          Corank(f)          = "+string(corank_f));
    debug_log(0, "          Determinacy       <= "+string(K));
    cstn[1] = corank_f;
    cstn[2] = Mu;
    cstn[3] = K;

    if( Mu == 0) { 
      cstn[1]=1; cstn[3]=1;
      return(printresult(1, f, "A[0]", cstn, 0)); }

    if(Mu<0) {
      debug_log(0, "The Milnor number of the function is infinite.");
      debug_log(0, "The singularity is not in Arnolds list.");
      return(printresult(1, 1, "error!", cstn, -1));
    }

    f = jet(f, K);
    v = HKclass(milnorcode(f));
    if(v[2]>0) { debug_log(0, "Guessing type via Milnorcode: ", v[1]);}
    else { 
      debug_log(0, "Hilbert polynomial not recognised. Milnor code = ",
                milnorcode(f));
    }
    debug_log(0, "");
    debug_log(0, "Computing normal form ...");

//------------ step 1, classification according to corank(f) ------------------
    if(corank_f == 0) { 
       return(printresult(2, f, "A["+string(Mu)+"]", cstn, 0)); }
    if(corank_f == 1) { 
       return(printresult(2, f, "A["+string(Mu)+"]", cstn, 0)); }
    cstn[4] = 0;
    if(corank_f == 2) { return(Funktion1bis(f, cstn)); }
    if(corank_f == 3) { return(Funktion1bis(f, cstn)); }
    return(printresult(105, f, "NoClass", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion1bis (poly f, list cstn)
{
//---------------------------- initialisation ---------------------------------
   def ring_top=basering;
   poly   g;
   int    n, corank, K;
   map    conv, PhiG;
   string s1;
   list   v_res, v_class, v, iv;
  
   corank = cstn[1];
   K = cstn[3];
   n = nvars(basering);

//-------------------- apply morsesplit if n>corank ---------------------------
   if( n > corank) {
     debug_log(0, 
      "I have to apply the splitting lemma. This will take some time....:-)");
     v_res = Morse(f, K, corank, 0);
     g = v_res[1];
     PhiG = v_res[2];

     conv = ReOrder(g);
     g = conv(g);
     PhiG = conv(PhiG);

     if(defined(RingB) != 0 ) { kill RingB; }
     ring ring_rest=char(basering),(x(1..corank)),(c,ds);
  
     map MapReduce=ring_top,maxideal(1);
     poly G = MapReduce(g);
     map PhiG=ring_top,maxideal(1);// Konstruiere Id auf r
     PhiG = MapReduce(PhiG);
  
     execute("ring RingB="+charstr(basering)+",("+A_Z("x",corank)+"),(c,ds);");
     export RingB;
     setring ring_rest;
   }
   else { 
     ring ring_rest=char(basering),(x(1..corank)),(c,ds);
     map  PhiG=ring_top,maxideal(1);
     poly G = PhiG(f);
   }
  
//--------------------- step 1 of Arnold now finished -------------------------
    map phi = ring_rest,maxideal(1);
    cstn[4] = phi;
    if(corank == 2) { v_class = Funktion3(G, cstn); }
    else {
      if(corank == 3) { v_class = Funktion50(G, cstn); }
      else { v_class = printresult(1, f, "error!", cstn, -1); }
    }
//-------------------------- classification done ------------------------------
    poly f_result = v_class[1];
    v[2] = v_class[2];
    v[3] = v_class[3];
    map Phi = v_class[4];
    PhiG = Phi(PhiG);
    setring ring_top;
    if(defined(conv)!=0) { kill conv; }
    map conv=ring_rest,maxideal(1);
    v[1] = conv(f_result);
    return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion3 (poly f, list cstn)
{ 
//---------------------------- initialisation ---------------------------------
    poly f3 = jet(f, 3);
    ideal Jf;
    int Dim, Mult, Mu;
    debug_log(1, "Step 3");

    if( f3 == 0 ) { return(Funktion13(f, cstn)); }

    // f3 ~ x3 , x2y+y3 , x2y
    Jf = std(jacob(f3));
    Dim = dim(Jf);
    if(Dim == 0) { return(printresult(4, f, "D[4]", cstn, 0)); }

    Mult = mult(Jf);
    Mu = cstn[2];
    if(Dim == 1) {
      if( Mult == 1) { return(printresult(5,f,"D["+string(Mu)+"]", cstn, 0)); }
      if( Mult == 2) { return(Funktion6(f, cstn));}         // series E,J
      debug_log(0, "dimension 1 und deg != 1, 2 => error, ", 
                        "this should never occur");
      return(printresult(3, f, "error!", cstn, -1));
      // Should never reach this line
    }
    // Should never reach this line
    return(printresult(3, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion6 (poly f, list cstn)
{ // Arnold's steps 6-12
//---------------------------- initialisation ---------------------------------
    poly   f3, fk;
    ideal  JetId, Jf;
    int    k, Dim, Mult, n, Mu;
    map    PhiG, Phi;
    intvec RFlg;
    list   v;

    def ring_top=basering;
    f3   = jet(f, 3);            // 3-Jet von f
    n    = nvars(basering);
    Mu   = cstn[2];
    PhiG = cstn[4];
    k    = 1;
    debug_log(1, "   Step 6");

    RFlg = GetRf(f, n);
    v    = Faktorisiere(f, f3, 3, 1, RFlg);
    f    = v[1];
    Phi  = v[2];
    PhiG = Phi(PhiG);
    cstn[4] = PhiG;

//---------------------------- begin of loop ----------------------------------
    while( (6*k) <= Mu ) {
      JetId = x(1)^3+x(2)^(3*k);
      fk    = jet(f, 3*k, weight(JetId));
      //--------------------------- step 6(k) ---------------------------------
      if( fk == Coeff(fk,x(1), x(1)^3)*x(1)^3 ) {
        JetId = x(1)^3+x(2)^(3*k+1);           // check jet x3,y3k+1  : E[6k]
        fk    = jet(f, 3*(3*k+1), weight(JetId));
        if( Coeff(fk,x(2),x(2)^(3*k+1)) != 0 ) {
          return(printresult(7, f, "E["+string(6*k)+"]", cstn, k-1));
        }

        JetId = x(1)^3+x(1)*x(2)^(2*k+1);      // check jet x3,xy2k+1 : E[6k+1]
        fk    = jet(f, 3*(2*k+1), weight(JetId));
        if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(2*k+1)) != 0 ) {
          return(printresult(8, f,"E["+string(6*k+1)+"]", cstn, k-1));
        }

        JetId = x(1)^3+x(2)^(3*k+2);           // check jet x3,y3k+1  : E[6k+2]
        fk    = jet(f, 3*(3*k+2), weight(JetId));
        if( Coeff(fk,x(2),x(2)^(3*k+2)) != 0 ) {
          return(printresult(9, f,"E["+string(6*k+2)+"]", cstn, k-1));
        }

        //------------------------- step 10(k) --------------------------------
        k++;
        JetId = x(1)^3+x(2)^(3*k);
        fk    = jet(f, 3*k, weight(JetId));
        Jf    = std(jacob(fk));
        Dim   = dim(Jf);

        if(Dim==0) { return(printresult(11,f,"J["+string(k)+",0]",cstn,k-1)); }
        Mult = mult(Jf);
        if( Dim ==1  && Mult==1) { 
          return(printresult(12,f,"J["+string(k)+","+string(Mu - 6*k +2)+"]",
                 cstn, k-1));
        }
        if( Dim == 1  && Mult == 2) {
          if(Coeff(fk, x(2), x(2)^(3*k)) != 0) {
            v    = Faktorisiere(f, fk, 3, k, RFlg);
            f    = v[1];
            Phi  = v[2];
            PhiG = Phi(PhiG);
            cstn[4] = PhiG;
          }
        }
      }
    }
    // Should never reach this line
    return(printresult(6, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion13 (poly f, list cstn)
{ 
//---------------------------- initialisation ---------------------------------
    poly f4;
    ideal Jf;
    int Dim, Mult, Mu;

    debug_log(1, "   Step 13");
    Mu = cstn[2];
    f4 = jet(f, 4);
    if( f4 == 0 ) { return(Funktion47(f, cstn)); }

    // f4 ~ x4+ax2y2+y4, x4+x2y2, x2y2, x3y, x4
    Jf  = std(jacob(f4));
    Dim = dim(Jf);

    if(Dim==0) { return(printresult(14,f,"X[9] = X[1,0] = T[2,4,4]",cstn,1)); }
    Mult = mult(Jf);
    if( Dim == 1) {
      if( Mult == 1 ) { 
        return(printresult(15, f, 
              "X[1,"+string(Mu-9)+"] = T[2,4,"+string(Mu-5)+"]", cstn, 1));
      }
      if( Mult == 2 ) {
        Jf = Jf, jacob(Jf);
        Jf = std(Jf);
        Dim = dim(Jf);
        if(Dim==0){return(printresult(16,f,"Y[1,p,q] = T[2,4+p,4+q]",cstn,1));}
        if( Dim == 1 ) { return(Funktion17(f, cstn)); }
      }
      if( Mult == 3 ) { return(Funktion25(f, cstn)); }
    } 
    // Should never reach this line
    return(printresult(13, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion17 (poly f, list cstn)
{ // Analog zu Fumktion 6, Kombination 17-24
//---------------------------- initialisation ---------------------------------
    poly  fk, ft;
    ideal JetId, Jf;
    int   p, Dim, Mult, Mu;
    list  v;
    map   PhiG, Phi;

    def ring_top=basering;
    debug_log(1, "      Step 17");
    Mu   = cstn[2];
    PhiG = cstn[4];
    fk   = jet(f, 4);
    p    = 1;

//---------------------------- begin of loop ----------------------------------
    while( 3*p<= Mu) {
      debug_log(1, "         Step 18("+string(p)+")");
      v   = Isomorphie_s17(f, fk, p, 1);
      f, Phi  = v[1..2];
      PhiG    = Phi(PhiG);
      cstn[4] = PhiG;

      if ( p>1) {
        JetId = x(1)^3*x(2) + x(2)^(3*p);
        fk = jet(f, 3*p, weight(JetId));
      }
      //--------------------------- step 18(p) --------------------------------
      JetId = x(1)^3*x(2) + x(2)^(3*p+2);     // check jet x3y,y3k+2  : Z[6p+5]
      fk = jet(f, 3*(3*p+2), weight(JetId));
      if( Coeff(fk, x(2), x(2)^(3*p+2)) != 0) { 
        return(printresult(19,f, "Z["+string(6*p+5)+"]", cstn, p));
      }

      JetId = x(1)^3*x(2)+x(1)*x(2)^(2*p+2);  // check jet x3y,xy2k+2 : Z[6p+6]
      fk = jet(f, 2*(3*p+2)+1, weight(JetId));
      if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(2*p+2)) != 0) { 
        return(printresult(20, f, "Z["+string(6*p+6)+"]", cstn, p));
      }

      JetId = x(1)^3*x(2) + x(2)^(3*p+3);     // check jet x3y,y3k+3  : Z[6p+7]
      fk = jet(f, 3*(3*p+3), weight(JetId));
      if( Coeff(fk, x(2), x(2)^(3*p+3)) != 0) {
        return(printresult(21, f, "Z["+string(6*p+7)+"]", cstn, p));
      }

      //---------------------------- step 22 ----------------------------------
      p = p+1;
      JetId = x(1)^3*x(2) + x(2)^(3*p+1);
      fk   = jet(f, 3*p+1, weight(JetId));
      ft   = Teile(fk, x(2));
      Jf   = std(jacob(ft));
      Dim  = dim(Jf);
      Mult = mult(Jf);
      if(Dim==0) { return(printresult(23,f,"Z["+string(p-1)+",0]", cstn, p)); }
      if( Mult == 1 ) { 
         return(printresult(24, f, "Z["+string(p-1)+","+string(Mu-3-6*p)+"]", 
                cstn, p));
      }
    }
    // Should never reach this line
    return(printresult(17, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion25 (poly f, list cstn)
{ // Analog zu Fumktion 6, Kombination 25-46
//---------------------------- initialisation ---------------------------------
    poly   fk, ft;
    ideal  JetId, Jf;
    int    k, Dim, Mult, Mu;
    map    PhiG, Phi;
    intvec RFlg;
    list   v;
    def ring_top=basering;

    debug_log(1, "      Step 25");
    Mu   = cstn[2];
    PhiG = cstn[4];
    fk   = jet(f, 4);
    k    = 1;
    RFlg = GetRf(f, 2);

//---------------------------- begin of loop ----------------------------------
    while (k<Mu) {
      v    =  Faktorisiere(f, fk, 4 , k, RFlg);
      f, Phi  = v[1..2];
      PhiG    = Phi(PhiG);
      cstn[4] = PhiG;

      //--------------------------- step 26(k) --------------------------------
      JetId = x(1)^4 + x(2)^(4*k+1);          // check jet x4,y4k+1  : W[12k]
      fk    = jet(f, 4*(4*k+1), weight(JetId));
      if( Coeff(fk, x(2), x(2)^(4*k+1)) != 0) {
        return(printresult(27, f, "W["+string(12*k)+"]", cstn, 3*k-2));
      }

      JetId = x(1)^4 + x(1)*x(2)^(3*k+1);     // check jet x4,xy3k+1 : W[12k+1]
      fk    = jet(f, 4*(3*k+1), weight(JetId));
      if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(3*k+1)) != 0) {
        return(printresult(28, f, "W["+string(12*k+1)+"]", cstn, 3*k-2));
      }

      //--------------------------- step 29(k) --------------------------------
      JetId = x(1)^4 + x(2)^(4*k+2);
      fk    = jet(f, 2*(4*k+2), weight(JetId));
      if( Coeff(fk, x(2), x(2)^(4*k+2)) != 0) { 
        Jf  = std(jacob(fk));
        Dim = dim(Jf);
        if(Dim==0) {return(printresult(30,f,"W["+string(k)+",0]",cstn,3*k-1));}
        if(Dim==1) { 
           return(printresult(32, f, 
                  "W#["+string(k)+","+string(Mu-12*k-2)+"]", cstn, 3*k-1));
        }
        return(printresult(29, f, "error!", cstn, -1));
      }
      else {
        // x^4 oder x^2(x^2+x(2)^2k+1)
        ft  = Teile(fk, x(1)^2);
        Jf  = std(jacob(ft));
        Dim = dim(Jf);
        if( Dim == 0 ) { 
           return(printresult(31, f, "W["+string(k)+","+string(Mu-12*k-3)+"]",
                  cstn, 3*k-1));
        }
        if( Dim != 1 ) { return(printresult(29, f, "error!", cstn, -1)); }

        //-------------------------- step 33(k) -------------------------------
        JetId = x(1)^4 + x(1)*x(2)^(3*k+2);   // check jet x4,xy3k+2 : W[12k+5]
        fk    = jet(f, 4*(3*k+2), weight(JetId));
        if( Coeff(fk, x(1)*x(2), x(1)*x(2)^(3*k+2)) != 0) {
          return(printresult(34, f,"W["+string(12*k+5)+"]", cstn, 3*k-1));
        }

        JetId = x(1)^4 + x(2)^(4*k+3);        // check jet x4,y4k+3  : W[12k+6]
        fk    = jet(f, 4*(4*k+3), weight(JetId));
        if( Coeff(fk, x(2), x(2)^(4*k+3)) != 0){
          return(printresult(35, f,"W["+string(12*k+6)+"]", cstn, 3*k-1));
        }

        //-------------------------- step 36(k) -------------------------------
        k = k+1;
        JetId = x(1)^4 + x(2)^(4*k);
        fk    = jet(f, (4*k), weight(JetId));
        Jf    = std(jacob(fk));
        Dim   = dim(Jf);
        Mult  = mult(Jf);
        if(Dim==0) {return(printresult(37,f,"X["+string(k)+",0]",cstn,3*k-1));}
        if(Dim==1) { 
          if(Mult==1) { 
             return(printresult(38, f,"X["+string(k)+","+string(Mu-12*k+3)+"]",
                    cstn, 3*k-1));
          }
          if(Mult==2) { 
            ft  = Teile(fk, x(1)^2);
            Jf  = std(jacob(ft));
            Dim = dim(Jf);
            if( Dim == 0) { return(Funktion40(f, cstn, k)); }
            if( Dim == 1) {
               return(printresult(39, f, "Y["+string(k)+",r,s]", cstn,3*k-1));
            }
          }
          if(Mult!=3) { 
            return(printresult(36, f, "error!", cstn, -1)); }
        }
        else { return(printresult(36, f, "error!", cstn, -1)); }
      }
    }  // Ende der While-Schleife
    // Should never reach this line
    return(printresult(25, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion40 (poly f, list cstn, int k)
{
//---------------------------- initialisation ---------------------------------
    int    r, kr, rr, sr, oldDebug;
    poly   fk, f2, a, b, c;
    ideal  JetId, Jfsyz;
    string Typ, RestRing, s1;
    list   v1, v2;
    def ring_top=basering;

    debug_log(1, "         Step 40" );

//------------------------------ compute f2 -----------------------------------
    JetId = x(1)^4 + x(2)^(4*k);
    fk  = jet(f, (4*k), weight(JetId));
    f2  = -fk / x(1)^3;
    Jfsyz = f - fk, x(1)^3, f2;
    matrix Mat = matrix(syz(Jfsyz));
    a = Mat[2,1] / Mat[1,1] - Mat[2,2];
    b = - Mat[3,1] / Mat[1,1] + Mat[3,2];
    ring tmp_ring=char(basering), (x(1),x(2)),(c,ds);
    map map_top2tmp=ring_top,maxideal(1);
    oldDebug = @DeBug;
    init_debug(-1);
//------------------------------ classify f2 ----------------------------------
    v1=Klassifiziere(map_top2tmp(b));
    init_debug(oldDebug);
    Typ = v1[2];
    v2 = DecodeNormalFormString(Typ);
    Typ, kr, rr, sr = v2[1..4];
    r   = kr-k;
    setring ring_top;
    if( Typ == "E[6k]" ) { 
       return(printresult(42, f, "Z["+string(k)+","+string(12*k+6*r-1)+"]", 
              cstn, 3*k+r-2));
    }
    if( Typ == "E[6k+1]" ) { 
       return(printresult(43, f, "Z["+string(k)+","+string(12*k+6*r)+"]",
              cstn, 3*k+r-2));
    }
    if( Typ == "E[6k+2]" ) {
       return(printresult(44, f, "Z["+string(k)+","+string(12*k+6*r+1)+"]",
              cstn, 3*k+r-2));
    }
    if( Typ == "J[k,0]" ) {
       return(printresult(45, f, "Z["+string(k)+","+string(r)+",0]",
              cstn, 3*k+r-2));
    }
    if( Typ == "J[k,r]" ) {
       return(printresult(46,f,"Z["+string(k)+","+string(r)+","+string(rr)+"]",
              cstn, 3*k+r-2));
    }
    // Should never reach this line
    return(printresult(40, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion50 (poly f, list cstn)
{ 
//---------------------------- initialisation ---------------------------------
    poly  f3;
    ideal Jf, Jf1, Jf2;
    int   Dim, Mult, Mu;
    debug_log(1, "Step 50");

    f3 = jet(f, 3);
    if( f3 == 0 ) { return(printresult(104, f, "NoClass", cstn, -1)); }

    // f3 ~ 
    Jf1  = jacob(f3);
    Jf   = std(Jf1);
    Dim  = dim(Jf);
    Mult = mult(Jf);
    Mu   = cstn[2];

    if(Dim == 0) { 
       return(printresult(51, f, "P[8] = T[3,3,3]", cstn, 1));
    } // x3 + y3 + z3 + axyz
    if(Dim == 1) { 
      if (Mult == 1) {
        return(printresult(52, f,"P["+string(Mu)+"] = T[3,3,"+string(Mu-5)+"]",
               cstn, 1));
      } // x3 + y3 + xyz
      if(Mult == 2) {
        Jf2 = wedge(jacob(Jf1),3-Dim), Jf1;
        Jf2 = std(Jf2);
        Dim = dim(Jf2);
        if (Dim==0) { return(printresult(54,f,"R[p,q] = T[3,p,q]", cstn, 1)); }
        if (Dim==1) { return(Funktion58(f, cstn)); }  // x3 + yz2
      }
      if(Mult == 3) {
        Jf2 = wedge(jacob(Jf1),3-Dim), Jf1;
        Jf2 = std(Jf2);
        Dim = dim(Jf2);
        if(Dim == 0) { return(printresult(56, f, "T[p,q,r]", cstn, 1)); }
        if(Dim == 1) { return(Funktion66(f, cstn)); }   // x2z + yz2
      }
      if(Mult == 4) { return(Funktion82(f, cstn)); }    // x3 + xz2
    }
    if(Dim == 2) {
      if(Mult == 1) { return(Funktion97(f, cstn)); }    // x2y
      if(Mult == 2) { return(printresult(103,f,"NoClass", cstn, -1));}
    }

    // Should never reach this line
    return(printresult(50, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion58 (poly fin, list cstn)
{
//---------------------------- initialisation ---------------------------------
    poly  f, f3, a, b, phi, b1, b2, b3, fa, fb, fc;
    ideal B, Jf3, J1, J2, S, Jfsyz;
    int   kx, ky, kz;
    string tp;
    matrix M[2][3];
    matrix C[2][3], D;
    list v;
    map  PhiG, VERT;
    def ring_top=basering;
    debug_log(1, "   Step 58");

    f    = fin;
    f3   = jet(f, 3);
    PhiG = cstn[4];
    tp   = "Nix";
    kx   = 1;     // Koordinate x
    ky   = 2;     // Koordinate y
    kz   = 3;     // Koordinate z
    B    = maxideal(1);     // ideal fuer Abbildungen
    Jf3  = jacob(f3);
    S    = sat(Jf3, maxideal(1))[1];
    J1   = diff(S[1], x(kx)), diff(S[1], x(ky)), diff(S[1], x(kz)),
           diff(S[2], x(kx)), diff(S[2], x(ky)), diff(S[2], x(kz));
    M    = J1;
    J2   = minor(M, 2), S;
 
    //------------------ determine coordinate named 'x' -----------------------
    S  = sat(J2, maxideal(1))[1];
    J1 = Coeff(S[1], x(1), x(1)), Coeff(S[1], x(2), x(2)),
          Coeff(S[1], x(3), x(3)), Coeff(S[2], x(1), x(1)),
          Coeff(S[2], x(2), x(2)), Coeff(S[2], x(3), x(3));
    C  = J1;
    D  = syz(C);
    kill C;

    b1 = D[1,1];
    b2 = D[2,1];
    b3 = D[3,1];

    debug_log(6, "f3,s1=", Show(f3));
    if( b1 != 0) {
      VERT=ring_top,-1*b1*x(1), -1*b2*x(1)+x(2), -1*b3*x(1) + x(3); 
      kx=1; ky=2; kz=3;
    }
    else {
      if( b2 != 0) {
        VERT=ring_top, x(1) + -1*b1*x(2), -1*b2*x(2), -1*b3*x(2) + x(3);
        kx=2; ky=1; kz=3;
      }
      else {
        if( b3 != 0) {
          VERT=ring_top,x(1) + -1*b1*x(3), x(2) + -1*b2*x(3), -1*b3*x(3);
          kx=3; ky=1; kz=2;
        }
      }
    }
    f       = VERT(f);
    PhiG    = VERT(PhiG);
    cstn[4] = PhiG;
    debug_log(6, VERT);
    f3 = jet(f,3);
    debug_log(6, "f3,s2=", Show(f3));

    //---------------- compute f_2 such that j3f = xf_2+f_3 -------------------
    debug_log(6, "1) x=", kx, "  y=", ky, "  z=", kz );
    matrix C = Coeffs(f3, x(kx)); 
    fb=C[2,1];  // Coeff von x^1
    fc=C[1,1];  // Coeff von x^0
    if(diff(fb, x(ky)) != 0) {
      Jfsyz = fb, diff(fb, x(ky));
      matrix Mat = matrix(syz(Jfsyz));
      B = maxideal(1);     // setze Abbildungs-ideal
      if( nrows(Mat) == 2) {
        poly Relation = -2 * Mat[2,1] / Mat[1,1];
        b = Coeff(Relation, x(kz), x(kz));
        B[rvar(x(ky))] = x(ky)-b*x(kz);
      } 
      else {
        Jfsyz = fb, diff(fb, x(kz));
        Mat = matrix(syz(Jfsyz));
        poly Relation = -2 * Mat[2,1];
        a = Coeff(Relation, x(ky), x(ky));
        B[rvar(x(kz))] = x(kz)-a*x(kz);
        ky, kz = swap(ky, kz);
      }
      VERT = ring_top, B;
      f    = VERT(f);
      PhiG = VERT(PhiG);
      cstn[4] = PhiG;
      f3 = jet(f,3);
      kill Mat;
    } 
    else { ky,kz = swap(ky,kz); }

    //------- compute tschirnhaus for 'z' and get f3=f_1(x,y,z)y^2+z^3 --------
    C    = Coeffs(f3, x(kx)); 
    fb   = C[2,1];  // Coeff von x^1
    fc   = C[1,1];  // Coeff von x^0
    v    = tschirnhaus(fc, x(kz));
    fc, VERT = v[1..2];
    f    = VERT(f);
    PhiG = VERT(PhiG);
    cstn[4] = PhiG;
    f3 = jet(f,3);
  
    //------------------- compute f_1 and get f3=xy^2+z^3 ---------------------
    fb = (f3 - 1*(Coeffs(f3, x(kz))[4,1])*x(kz)^3)/(x(ky)^2);
    fc=(x(kx)-1*(Coeffs(fb,x(ky))[2,1])*x(ky)-1*(Coeffs(fb,x(kz))[2,1])*x(kz));
    fa = Coeffs(fb, x(kx))[2,1];
    if ( fa != 0 ) {
      B = maxideal(1);
      B[rvar(x(kx))] = fc / fa;
      VERT = ring_top, B;
      f    = VERT(f);
      PhiG = VERT(PhiG);
      cstn[4] = PhiG;
      f3   = jet(f,3);
    }

    //--------------------- permutation of x,y,z  -----------------------------
    if(Coeffs(f3, x(1))[4,1]!=0) { 
      kx=1;
      if(Coeffs(f3, x(2))[3,1]==0) { ky=2; kz=3; }
      else { ky=3; kz=2; }
    }
    else {
      if(Coeffs(f3, x(2))[4,1]!=0) { 
        kx=2;
        if(Coeffs(f3, x(3))[3,1]==0) { ky=3; kz=1; }
        else { ky=1; kz=3; }
      }
      else { 
        kx=3;
        if(Coeffs(f3, x(1))[3,1]==0) { ky=1; kz=2; }
        else { ky=2; kz=1; }
      }
    }
    VERT = ring_top, x(kx), x(ky), x(kz);
    f    = VERT(f);
    PhiG = VERT(PhiG);
    cstn[4] = PhiG;
    f3   = jet(f,3);
    return(Funktion59(f, cstn));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion59 (poly f, list cstn)
{
//---------------------------- initialisation ---------------------------------
    poly   phi, fr, fk, alpha, beta, f_tmp;
    ideal JetId;
    int    p, Dim, Mult, Mu;
    string tp;
    list   v;
    map    PhiG, Phi;
    intvec w, RFlg;
    def ring_top=basering;
    debug_log(1, "      Step 59");

    Mu   = cstn[2];
    PhiG = cstn[4];
    tp   = "Nix";
    p    = 1;
    phi  = jet(f,3);
    fr   = f - phi;
    RFlg = 1,2,3;
    alpha = coeffs(fr, x(1))[1,1];
    beta = (fr - alpha) / x(1);
    debug_log(3, "f    = ", Show(f));
    debug_log(3, "fr   = ", Show(fr));
    debug_log(3, "alpha= ", Show(alpha));
    debug_log(3, "beta = ", Show(beta));

//---------------------------- begin of loop ----------------------------------
    while(6*p<Mu) {
      JetId = x(2)^(3*p+1);
      JetId = phi + x(2)^(3*p+1);
      //--------------------------- step 59(k) --------------------------------
      w     = weight(JetId);
      fk    = jet(fr, 3*w[1], w);
      if(fk!=0) { return(printresult(60,f, "Q["+string(6*p+4)+"]", cstn, p)); }

      JetId = phi + x(1)*x(2)^(2*p+1);
      w     = weight(JetId);
      fk    = jet(fr, 3*w[1], w);
      if(fk!=0) { return(printresult(61,f, "Q["+string(6*p+5)+"]", cstn, p)); }

      JetId = phi + x(2)^(3*p+2);
      w     = weight(JetId);
      fk    = jet(fr, 3*w[1], w);
      if(fk!=0) { return(printresult(62,f, "Q["+string(6*p+6)+"]", cstn, p)); }

      //--------------------------- step 63(k) --------------------------------
      p++;
      JetId = phi + x(2)^(3*p);
      w     = weight(JetId);
      fk    = jet(f, 3*w[1], w);
      JetId = std(jacob(fk));
      Dim   = dim(JetId);
      Mult  = mult(JetId);
      if(Dim==0) { return(printresult(64, f, "Q["+string(p)+",0]", cstn, p)); }
      if(Dim==1) { 
        if(Mult == 1) {
           return(printresult(65, f, "Q["+string(p)+","+string(Mu-(6*p+2))+"]",
                  cstn, p));
        }
        if(Mult == 2) { 
          fk    = jet(fr, 3*w[1], w);
          f_tmp = Coeffs(phi, x(1))[4,1] *x(1)^3+fk;
          v    = Faktorisiere(f, f_tmp, 3 , p, RFlg);
          f    = v[1];
          Phi  = v[2];
          PhiG = Phi(PhiG);
          cstn[4] = PhiG;
          fr = f - phi;
        }
      }
    }
    // Should never reach this line
    return(printresult(59, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion66 (poly f, list cstn)
{
//---------------------------- initialisation ---------------------------------
    int kx = 1; // Koordinate x
    int ky = 2; // Koordinate y
    int kz = 3; // Koordinate z
    poly f3 = jet(f, 3);
    ideal JetId;

    debug_log(1, "   Step 66");
    debug_log(2, "F3=", Show(f3));
    poly fx = diff(f3, x(kx));
    JetId = jacob(fx);
    JetId = std(JetId);
    "nach x:",Show(fx), "  Id=", JetId, "  Dim=", dim(JetId);

    poly fy = diff(f3, x(ky));
    JetId = jacob(fx);
    JetId = std(JetId);
    "nach y:",Show(fy), "  Id=", JetId, "  Dim=", dim(JetId);

    poly fz = diff(f3, x(kz));
    JetId = jacob(fx);
    JetId = std(JetId);
    "nach z:",Show(fz), "  Id=", JetId, "  Dim=", dim(JetId);
    return(printresult(1, 66, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion82 (poly f, list cstn)
{
//---------------------------- initialisation ---------------------------------
    poly   f3, b1, b2, b3;
    int    i, kx, ky, kz, Fall;
    ideal  Jfsyz, B;
    intvec kv = 1,2,3;
    matrix Mat;
    map    PhiG, VERT;
    list   v;
    def ring_top=basering;
    debug_log(1, "   Step 82");

    f3 = jet(f,3);
    kx = 1;     // Koordinate x
    ky = 2;     // Koordinate y
    kz = 3;     // Koordinate z
    B  = maxideal(1);
    Jfsyz = jacob(f3);
    PhiG  = cstn[4];
    Fall  = 2;

//------------------ find coordinatechange that f3 ~ g(x,z) -------------------
    if (diff(f3, x(1)) == 0) { kx, ky = swap(kx, ky); }
    if (diff(f3, x(2)) == 0) {  }
    if (diff(f3, x(3)) == 0) { kz, ky = swap(kz, ky); }
    if ( (diff(f3, x(1)) != 0) && (diff(f3, x(2)) != 0) &&
          (diff(f3, x(3)) != 0) ) {
      Mat = matrix(syz(Jfsyz));
      b1  = Mat[1,1];
      b2  = Mat[2,1];
      b3  = Mat[3,1];

      if( b1 != 0) {
        VERT = ring_top,b1*x(kx), b2*x(kx)+x(ky), b3*x(kx) + x(kz);
        kx, ky = swap(kx, ky);
      }
      else {
        if(b2!=0) { VERT = ring_top,x(kx)+b1*x(ky),b2*x(ky),b3*x(ky)+x(kz); }
        else {
          if(b3!=0) { VERT = ring_top,x(kx)+b1*x(kz),x(ky)+b2*x(kz),b3*x(kz); }
          else { VERT = ring_top,B; }
        }
      }
      f    = VERT(f);
      PhiG = VERT(PhiG);
      cstn[4] = PhiG;
    }
    VERT = ring_top,x(kx),x(ky),x(kz);
    f    = VERT(f);
    PhiG = VERT(PhiG);
    cstn[4] = PhiG;
    f3   = jet(f,3);

    if( (f3-subst(f3, x(kx), 0)) == 0) { kx, ky = swap(kx, ky); }
    if( (f3-subst(f3, x(kz), 0)) == 0) { kz, ky = swap(kz, ky); }

//------------ find coordinatechange for f3 ~ x3+xz2, if possible  ------------
    matrix C = coeffs(f3, x(kx)); 
    if(size(C) == 3) { C = coeffs(f3, x(kz)); kx,kz=swap(kx, kz); }
    if(C[1,1] == 0 && C[3,1] == 0) { Fall = 1; }
    if(C[1,1] != 0 && C[3,1] != 0 ) { Fall = 3; }
    if(C[1,1] == 0 && C[3,1] != 0 ) { Fall = 2; }
    if(C[1,1] != 0 && C[3,1] == 0 ) { Fall = 2; kx,kz=swap(kx, kz); }

    if(Fall == 1) {
      VERT=ring_top,x(kx),x(ky),x(kz);
    }
    if(Fall == 2) { 
       v = tschirnhaus(f3/x(kz), x(kx));
       b1, VERT = [1..2];
    }
    if(Fall == 3) { 
      v = tschirnhaus(f3/x(kx), x(kx));
      b1, VERT = [1..2];
      debug_log(2, "B1=", Show(jet(VERT(f),3)));
      v = tschirnhaus(f3/x(kz), x(kz));
      b2, VERT = [1..2];
      debug_log(2, "B2=", Show(jet(VERT(f),3)));
    }
    f    = VERT(f);
    PhiG = VERT(PhiG);
    cstn[4] = PhiG;
    f3   = jet(f,3);

//------------- if f3 ~ x3+xz2 then continue with classification  -------------
    C = coeffs(f3, x(1)); 
    if( C[1,1] == 0 && C[2,1] != 0 && C[3,1] == 0 && C[4,1] != 0 ) { 
      return(Funktion83(f, cstn));
    }
    return(printresult(82, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Isomorphie_s82_z (poly f, poly fk, int p)
{
//---------------------------- initialisation ---------------------------------
    poly   Relation, a, b;
    ideal  Jfsyz, B;
    matrix Mat;
    map    VERT;
    list   v;
    def ring_top=basering;

    debug_log(1, "      Isomorphie 82/90 z");
    debug_log(2, "tt=", Show(fk));
    Jfsyz = fk, diff(fk, x(3));
    Mat   = matrix(syz(Jfsyz));
    Relation = -2 * Mat[2,1] / Mat[1,1];
    a = Coeff(Relation, x(3), x(3));
    b = Coeff(Relation, x(2), x(2)^p);
    B = maxideal(1);
    B[rvar(x(3))] = x(3)-b*x(2)^p;
    VERT = ring_top,B;
    v    = VERT(f), VERT;
    debug_log(2, VERT);
    debug_log(2, "      z res=", Show(VERT(fk)));
    return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc Isomorphie_s82_x (poly f, poly fk, int p)
{
//---------------------------- initialisation ---------------------------------
    matrix Mat;
    poly   Relation, a, b;
    ideal  Jfsyz, B;
    map    VERT;
    list   v;
    def ring_top=basering;

    debug_log(1, "      Isomorphie 82/90 x");
    debug_log(2, "tt=", Show(fk));
    Jfsyz = fk, diff(fk, x(1));
    Mat   = matrix(syz(Jfsyz));
    Relation = -3 * Mat[2,1] / Mat[1,1];
    a = Coeff(Relation, x(1), x(1));
    b = Coeff(Relation, x(2), x(2)^p);
    B = maxideal(1);
    B[rvar(x(1))] = x(1)-b*x(2)^p;
    VERT = ring_top,B;
    v    = VERT(f), VERT;
    debug_log(2, VERT);
    debug_log(2, "      x res=", Show(VERT(fk)));

    return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion83 (poly f, list cstn)
{
//---------------------------- initialisation ---------------------------------
    poly   fk, phi, a, b;
    ideal  JetId, Jf, B;
    int    k, Dim, Mult;
    intvec w;
    map    PhiG, Phi;
    list   v;
    matrix Mat;
    def ring_top=basering;
    debug_log(1, "      Step 83");

    k    = 1;
    PhiG = cstn[4];
//---------------------------- begin of loop ----------------------------------
    while(k<10) {
      phi   = jet(f, 3);
      //--------------------------- step 83(k) --------------------------------
      JetId = x(1)^3 + x(3)^3 + x(2)^(3*k+1);
      fk    = jet(f- phi, 3*w[1], weight(JetId)) ;
      if(fk!=0) { return(printresult(84,f,"U["+string(12*k)+"]",cstn,4*k-3)); }

      JetId = x(1)^3 + x(3)^3 + x(1)*x(2)^(2*k+1);
      fk    = jet(f, 3*w[1], weight(JetId)) ;
      //--------------------------- step 85(k) --------------------------------
      if ( fk != phi ) {
        Jf   = std(jacob(fk));
        Dim  = dim(Jf);
        if(Dim==0) {return(printresult(86,f,"U["+string(k)+",0]",cstn,4*k-2));}
        if(Dim==1) {return(printresult(87,f,"U["+string(k)+",p]",cstn,4*k-2));}
      }

      //--------------------------- step 88(k) --------------------------------
      JetId = x(1)^3 + x(3)^3 + x(2)^(3*k+2);
      fk    = jet(f- phi, 3*w[1], weight(JetId)) ;
      if(fk!=0) {return(printresult(89,f,"U["+string(12*k+4)+"]",cstn,4*k-2));}

      //--------------------------- step 90(k) --------------------------------
      k++;
      JetId = x(1)^3 + x(3)^3 + x(2)^(3*k);
      fk    = jet(f, 3*w[1], weight(JetId)) ;
      Jf    = std(jacob(fk));
      Dim   = dim(Jf);
      Mult  = mult(Jf);
      if ( Dim == 0 ) { return(printresult(83, f, "NoClass", cstn, -1)); }
      if ( Dim == 1 ) {
        if ( Mult == 4 ) {
          if( fk - phi != 0) { // b!=0  und/oder b'!=0 
            if( Coeff(fk,x(1)*x(2), x(1)^2*x(2)^k) == 0 ) { // b=0 und b'!=0
              a    = (fk - Coeff(fk, x(1), x(1)^3)*x(1)^3) / x(1);
              v    = Isomorphie_s82_z(f, a, k);
            } 
            else {
              if( Coeff(fk,x(1)*x(2)*x(3), x(1)*x(2)^k*x(3)) == 0 ){ 
                        // b!=0 und b'=0
                a    = subst(fk, x(3), 0);
                v    = Isomorphie_s82_x(f, a, k);
              }
              else {
                a = Coeff(fk,x(1)*x(2)*x(3), x(1)*x(2)^k*x(3));
                b = Coeff(fk,x(2)*x(3), x(2)^(2*k)*x(3));
                B = maxideal(1);
                B[rvar(x(1))] = x(1)-b/a*x(2)^k;
                Phi  = ring_top,B;
                f    = Phi(f);
                PhiG = Phi(PhiG);
                cstn[4] = PhiG;
                fk   = jet(f, 3*w[1], w) ;
                a    = (fk - Coeff(fk, x(1), x(1)^3)*x(1)^3) / x(1);
                v    = Isomorphie_s82_z(f, a, k);
              } // ende else b!=0 und b'=0
            } // ende else b=0 und b'!=0
            f, Phi  = v[1..2];
            PhiG    = Phi(PhiG);
            cstn[4] = PhiG;
          } //ende fk-phi!=0
        } // ende mult=4
        else { return(printresult(83, f, "NoClass", cstn, -1)); }
      } // ende dim=1
      else { return(printresult(83, f, "NoClass", cstn, -1)); }
    } // ENDE While
    return(printresult(83, f, "error!", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion97 (poly f, list cstn)
{
//---------------------------- initialisation ---------------------------------
    poly   f3, l1, l2, a, b, c, prod;
    ideal  Jfsyz, Jf, B;
    int    k, i, pt, kx, ky, kz, Dim, Mult, Mu;
    matrix Mat;
    map    PhiG, VERT;
    def ring_top=basering;
    debug_log(1, "   Step 97");

    Mu   = cstn[2];
    PhiG = cstn[4];
    kx   = 1;   // Koordinate x
    ky   = 2;   // Koordinate y
    kz   = 3;   // Koordinate z
    B    = maxideal(1); // Abbildungs-ideal
    pt   = 2;
    f3   = jet(f, 3);
    k    = 1;

//--------------------------- compute f3 ~ x2y --------------------------------
    // vertausche 2 Koordinaten sodass d2f/dx2 <>0 ist.
    for(i=1;i<4;i=i+1) {
      if(diff(diff(f3, x(i)), x(i)) != 0) { kx = i; i=4; }
    }
    if(kx == 2) { ky = 1; kz = 3; }
    if(kx == 3) { ky = 2; kz = 1; }
    //-------------------------- compute -l1*l2 -------------------------------
    f3    = jet(f, 3);
    Jfsyz = f3, diff(f3, x(kx));
    Mat   = matrix(syz(Jfsyz));
    if(deg(Mat[2,1])>1) {
      Jfsyz = f3, Mat[2,1];
      Mat   = matrix(syz(Jfsyz));

      // berechen Abb. sodass f=x2*l2
      l1 = Mat[2,1];
      a  = Coeff(l1, x(kx), x(kx));
      l1 =  l1 / number(a);
      b  = Coeff(l1, x(ky), x(ky));
      c  = Coeff(l1, x(kz), x(kz));
      B[rvar(x(kx))] = x(kx) - b * x(ky) - c * x(kz);
      VERT = ring_top, B;
      f    = VERT(f);
      PhiG = VERT(PhiG);
      cstn[4] = PhiG;

      f3 = jet(f, 3);
      l2 = f3 / x(kx)^2;

      // sorge dafuer, dass b<>0 ist.
      b = Coeff(l2, x(ky), x(ky));
      if( b== 0) { ky, kz = swap(ky, kz); }
  
      // Koordinaten-Transf. s.d. f=x2y
      b  = Coeff(l2, x(ky), x(ky));
      l2 =  l2 / number(b);
      a  = Coeff(l2, x(kx), x(kx));
      c  = Coeff(l2, x(kz), x(kz));
      B  = maxideal(1);
      B[rvar(x(ky))] = -a * x(kx) + x(ky) - c * x(kz);
      VERT = ring_top, B;
      f    = VERT(f);
      PhiG = VERT(PhiG);
      cstn[4] = PhiG;
    } 

//------------------------------- step 98 ---------------------------------
    Jfsyz = x(kx)^2*x(ky) + x(ky)^4 + x(kz)^4;
    a     = jet(f, 8, weight(Jfsyz));
    Jf    = std(jacob(a));
    Dim   = dim(Jf);
    Mult  = mult(Jf);
    if( Dim == 0) { return(printresult(99, f, "V[1,0]", cstn, 3)); }
    if( Dim == 1) {
      if(Mult==1) {return(printresult(100,f,"V[1,"+string(Mu-15)+"]",cstn,3));}
      if(Mult==2){return(printresult(101,f,"V#[1,"+string(Mu-15)+"]",cstn,3));}
    }
    " Dim=",Dim," Mult=",Mult;
    return(printresult(102, f, "NoClass", cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc tschirnhaus (poly f, poly x)
USAGE:    tschirnhaus();
{
//---------------------------- initialisation ---------------------------------
    poly   b;
    ideal  B;
    int    n, j, hc;
    matrix cf;
    intvec z;
    string s;
    list   v;
    map    Phi, EH;
    def ring_top=basering;

    n    = nvars(basering);
    cf   = coeffs(f, x);
    hc   = nrows(cf) - 1; // hoechster exponent von x_i
    b    = cf[hc+1,1];    // koeffizient von x_i^hc
    B    = maxideal(1);
    z[n] = 0;
    EH   = ring_top, z;
    Phi  = ring_top, B;
    v[1] = f;
    if ( EH(b) != 0)    // pruefe ob der Koeff von x_i^hc 
    { B[rvar(x)] = x -1*(cf[hc,1]/(hc*b));
      v[1] = Phi(f);
    }
    v[2] = Phi;
    return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc Isomorphie_s17 (poly f, poly fk, int k, int ct, list #)
{
//---------------------------- initialisation ---------------------------------
    ideal  Jfsyz, JetId, bb;
    poly   a, b, c, d, Relation, an, bn;
    int    B,C, alpha, beta, gamma, g;
    matrix Matx, Maty;
    map    PhiG, VERT;
    list   v;
    def ring_top=basering;

    if(size(#)==1) { PhiG = #[1]; }
    else { PhiG = ring_top,maxideal(1); }
    bb = maxideal(1);

    // Ziel: bestimme a,b,c,d sodass  fk = (ax+by^k)^3(cx+dy) gilt.
    debug_log(2, "Isomorphie_s17:");
    debug_log(2, "Faktor: f=",Show(f)," Jet=",Show(fk)," k=",k, "cnt=", ct);

    if( k == 1) {
      Jfsyz = fk, diff(fk, x(1));
      Matx  = matrix(syz(Jfsyz));
      Jfsyz = fk, diff(fk, x(2));
      Maty  = matrix(syz(Jfsyz));

      a = Coeff(fk, x(1), x(1)^4);
      b = Coeff(fk, x(2), x(2)^4);
      c = Coeff(fk, x(1)*x(2), x(1)^3*x(2));
      d = Coeff(fk, x(1)*x(2), x(1)*x(2)^3);

      if( (a != 0) && (b != 0) ) {  
        B = -int(Coeff(Matx[1,1], x(2), x(2)));
        C = -int(Coeff(Maty[1,1], x(1), x(1)));
        alpha = int(Coeff(Matx[2,1], x(1), x(1)^2));
        beta  = int(Coeff(Matx[2,1], x(1)*x(2), x(1)*x(2)));
        gamma = int(Coeff(Matx[2,1], x(2), x(2)^2));

        bb[rvar(x(1))] = x(1) - (2*gamma / (B - beta))*x(2);
        bb[rvar(x(2))] = x(2) - ((C - beta) / (2*gamma))*x(1);
        VERT     = ring_top,bb;
        Relation = VERT(f);
        fk       = jet(Relation, 4);

        an = Coeff(fk, x(1), x(1)^4);
        bn = Coeff(fk, x(2), x(2)^4);
        if( (an != 0) & (bn != 0) ) { VERT=ring_top,x(1),(x(2) + a*x(1))/ b; }
        f    = VERT(f);
        fk   = jet(f, 4);
        PhiG = VERT(PhiG);

        a = Coeff(fk, x(1), x(1)^4);
        b = Coeff(fk, x(2), x(2)^4);
        c = Coeff(fk, x(1)*x(2), x(1)^3*x(2));
        d = Coeff(fk, x(1)*x(2), x(1)*x(2)^3);
        Jfsyz = fk, diff(fk, x(1));
        Matx  = matrix(syz(Jfsyz));
        Jfsyz = fk, diff(fk, x(2));
        Maty  = matrix(syz(Jfsyz));
      }

      if( (a == 0) || (b == 0) ) {
        if( a == 0) {
          if( c == 0) { // y3(ax+by)
            Relation = - Matx[2,1] / Matx[1,1];
            a = Coeff(Relation, x(1), x(1));
            b = Coeff(Relation, x(2), x(2));
            VERT=ring_top,a*x(2)^k - b*x(1), x(1);
          }
          else { // (ax+by)^3y
            Relation = - 3*Matx[2,1] / Matx[1,1];
            a = Coeff(Relation, x(1), x(1));
            b = Coeff(Relation, x(2), x(2));
            VERT=ring_top,a*x(1) - b*x(2), x(2);
          }
        }
        else {
          if( d == 0) { // x3(ax+by)
            Relation = - Maty[2,1] / Maty[1,1];
            a = Coeff(Relation, x(1), x(1));
            b = Coeff(Relation, x(2), x(2));
            VERT=ring_top,x(1), b*x(2)^k - a*x(1);
          }
          else { // x(ax+by)^3
            Relation = - 3*Maty[2,1] / Maty[1,1];
            a = Coeff(Relation, x(1), x(1));
            b = Coeff(Relation, x(2), x(2));
            VERT=ring_top,x(2), b*x(1) - a*x(2);
          }
        }
        f    = VERT(f);
        PhiG = VERT(PhiG);
      } 
      else {  //      "Weder b noch a sind 0";
        if(ct > 5) { v[1]=f; v[2]=PhiG; return(v); }
        fk = jet(f, 4);
        return(Isomorphie_s17(f, fk, k, ct+1, PhiG));
      }
    } 
    else {  // k >1
      a     = fk/x(2);
      Jfsyz = a, diff(a, x(1));
      Matx  = matrix(syz(Jfsyz));
      Relation = -3 * Matx[2,1] / Matx[1,1];
      a    = Coeff(Relation, x(1), x(1));
      b    = Coeff(Relation, x(2), x(2)^k);
      VERT = basering,x(1)-b*x(2)^k,x(2);
      f    = VERT(f);
      PhiG = VERT(PhiG);
      JetId = x(1)^3*x(2) + x(2)^(3*k+1);
      fk = jet(f, 3*k+1, weight(JetId));
    }
    v = f, PhiG;
    debug_log(2, "Isomorphie_s17: done");
    debug_log(2, "Faktor: f=",Show(f)," Jet=",Show(fk)," k=",k);

    return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc printresult (int step, poly f, string typ, list cstn, int m)
{
//---------------------------- initialisation ---------------------------------
  int corank, Mu, K;
  list v;

  corank, Mu, K = cstn[1..3];
  debug_log(0,"   Arnold step number "+string(step));
  if( @DeBug >= 0 ) {
    "The singularity";
    "   `"+Show(jet(f, K))+"'";
    if( typ != "error!" && typ != "NoClass" ) {
      "is R-equivalent to "+typ+".";
    }
    if ( typ == "NoClass" ) {
      "is not in Arnolds list.";
    }
//    if(K>=0)  { "  det = "+string(K); }
    if(Mu>=0) { "   Milnor number = "+string(Mu); }
    if(m>=0)  { "   modality      = "+string(m); }
  }
  v = f, typ, corank, cstn[4];
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion47 (poly f, list cstn)
{
    int corank = cstn[1];
    int Mu = cstn[2];
    int K  = cstn[3];
    string s = "The Singularity '";+Show(jet(f, K), corank, K));
    string tp="";
//    return(printresult(47, f, tp, cstn, -1));

    s = s +"' has 4-jet equal to zero. (F47), mu="+string(Mu);
  
    s; // +"  ("+SG_Typ+")";
    return(Show(f), tp, corank);
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion91 (poly f, list cstn, int k)
{
    string tp  = "U*[k,0]";
    return(printresult(91, f, tp, cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion92 (poly f, list cstn, int k)
{
    string tp  = "UP[k]";
    return(printresult(92, f, tp, cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion93 (poly f, list cstn, int k)
{
    string tp  = "UQ[k]";
    return(printresult(93, f, tp, cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion94 (poly f, list cstn, int k)
{
    string tp  = "UR[k]";
    return(printresult(94, f, tp, cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion95 (poly f, list cstn, int k)
{
    string tp  = "US[k]";
    return(printresult(95, f, tp, cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc Funktion96 (poly f, list cstn, int k)
{
    string tp  = "UT[k]";
    return(printresult(96, f, tp, cstn, -1));
}

///////////////////////////////////////////////////////////////////////////////
proc morsesplit(poly f)
USAGE:    morsesplit(f);        f=poly
RETURN:   Normal form of f in M^3 after application of the splitting lemma
COMPUTE:  apply the splitting lemma (generalized Morse lemma) to f
EXAMPLE:  example morsesplit; shows an example
{
//---------------------------- initialisation ---------------------------------
  poly f_out;
  int  n, K, Mu, corank;
  list v;

  if(defined(@ringdisplay) != 0 ) { kill @ringdisplay; }
  string @ringdisplay = "setring "+nameof(basering);
  export @ringdisplay;

  def ring_ext=basering;

  n = nvars(basering);

  // if trace/debug mode not set, do it!
  init_debug();
  K, Mu, corank = basicinvariants(f);
  ring ring_top=char(basering),(x(1..n)),(c,ds);

  map Conv=ring_ext,maxideal(1);
  setring ring_top;
  v = Morse(jet(Conv(f),K), K, corank, 0);
  poly f_out = v[1];
  setring ring_ext;
  map ConvUp = ring_top, maxideal(1);
  return(ConvUp(f_out));
}
example
{ "EXAMPLE"; echo=2;
   ring r=0,(x,y,z),ds;
   export r;
   init_debug(1);
   poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
   poly g=morsesplit(f);
   g;
}

///////////////////////////////////////////////////////////////////////////////
proc Coeffs (list #)
{
  matrix m=matrix(coeffs(#[1],#[2]), deg(#[1])+1, 1);
  return(m);
}

///////////////////////////////////////////////////////////////////////////////
proc Morse(poly f, int K, int corank, int ShowPhi)
{ 
//---------------------------- initialisation ---------------------------------
    poly   fc, f2, a, P, Beta, fi;
    ideal  Jfx, B;
    int    n, i, j, k, Rang, Done;
    matrix Mat;
    map    Id, Psi, Phi, PhiG;
    intvec Abb, RFlg;
    list   v;

    fi = f;
    n = nvars(basering);
    init_debug();

    def ring_top=basering;
  
    debug_log(3, "Spalte folgendes Polynom mit Bestimmtheit: ", string(K));
    debug_log(3, Show(fi));
  
    for( j=1; j<=n ; j++) { Abb[j] = 0; }
  
    RFlg = GetRf(fi, n);
    debug_log(2, "Reihenfolge fuer Vertauschungen:", RFlg );
    PhiG=ring_top,maxideal(1);
  
//----------------- find quadratic term, if there is only one -----------------
    B = maxideal(1);
    if(corank == (n-1)) {
      Done = 0;
      f2   = jet(fi, 2);
      j    = 1;
      Jfx  = f2, diff(f2, x(j));
      while(j<=n && (diff(f2, x(j))==0)) {
        j   = j+1;
        Jfx = f2, diff(f2, x(j));
      }
      Mat  = matrix(syz(Jfx));
      Beta = 2*Mat[2,1]/Mat[1,1];
      for( j=1; j<=n ; j=j+1) {
        f2 = Coeff(Beta, x(RFlg[j]), x(RFlg[j]));
        if(f2!=0) {
          k = RFlg[j];
          break;
        }
      }
      for( j=1; j<=n ; j=j+1) {
        f2 = Coeff(Beta, x(j), x(j));
        if(j == k) { B[rvar(x(j))] = (2*f2*x(j)-Beta) / number(f2); }
      }
      Phi  =ring_top,B;
      fi   = Phi(fi);   
      PhiG = Phi(PhiG);
    }
    if( ShowPhi > 1) { PhiG; }
  
//------------------------ compute spliting lemma -----------------------------
    fc = fi;
    i  = 1;              // Index fuer Variablen wird bearbeitet
    while( i <= n) {
      Phi=ring_top,maxideal(1);
      debug_log(6, "Pruefe Variable x(" +string(RFlg[i]) + ")");
      debug_log(6, "--------------------");
      j  = i + 1;        // setze j fuer evtle Verschiebung
      f2 = jet(fc,2);
      debug_log(6, "Rechne 2-Jet =" , string(f2));
      if( (f2 - subst(f2, x(RFlg[i]), 0)) == 0 ) { Abb[RFlg[i]] = 1; }
      if( (f2 - subst(f2, x(RFlg[i]), 0)) != 0 ) {
        while( (j<=n) || (i==n) ) {
          debug_log(6, "Pruefe 2-Jet mit Wert : " + string(jet(fc,2)));
          a = Coeff(jet(fc,2), x(RFlg[i]), x(RFlg[i])^2);
          debug_log(6,"Koeffizient von x(" + string(RFlg[i]) + ")^2 ist:", a);
          if( (a != 0) || (i==n) ) { 
            debug_log(6, "BREAK!!!!!!!!!!!!!!");
            break;
          }
          debug_log(6,"Verschiebe evtl Variable x(",string(RFlg[j]),") um x(",
                       string(RFlg[i]), ")");
          B = maxideal(1);
          for( k=1; k<=n ; k=k+1) {
            if(k==RFlg[j]) { B[rvar(x(k))] = x(k) + x(RFlg[i]); }
          }
          Phi = ring_top,B;
          fc  = Phi(fi);
          j   = j + 1;
        }               // Ende while( (j<=n) || (i==n) )
  
        debug_log(6, "Moegliche Verschiebung fertig!");
        PhiG = Phi(PhiG);
        if( ShowPhi > 1) { "NachVersch.:"; Phi; }
  
        if( (j<=n) || (i==n)) {
          P = Coeff(fc, x(RFlg[i]), x(RFlg[i]));
          debug_log(6, "Koeffizient von x("+string(RFlg[i])+") ist: ",
                        string(P));
          if(P != 0) {
            debug_log(6, "1 Koeffizient von x("+string(RFlg[i])+") ist: ",
                         string(P));
            debug_log(6, "a=" + string(a));
            P = P / number (2 * a);
            debug_log(6, "2 Koeffizient von x("+string(RFlg[i])+") ist: ",
                         string(P));
            B = maxideal(1);
            for( k=1; k<=n ; k=k+1) {if(k==RFlg[i]) {B[rvar(x(k))]=x(k)-P;}}
            Phi =ring_top,B;
            debug_log(6, "Quadratische-Ergaenzung durch:", Phi);
            fi   = Phi(fc);
            PhiG = Phi(PhiG);
            fc   = jet(fi,K);
            P    = Coeff(fc, x(RFlg[i]), x(RFlg[i]));
            if( P != 0) {
              fi = fc;
              continue;
            }
          }     // Ende if(P != 0)
                // Fertig mit Quadratischer-Ergaenzung
        }               // Ende if( (j<=n) || (i==n))
      }                 // Ende if( (f2 - subst(f2, x(RFlg[i]), 0)) != 0 )
  
      fi = fc;
      i  = i + 1;
      debug_log(6, "++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
    }
    debug_log(6, "Ende  ---------------------------------------------------");
  
//--------------------------- collect results ---------------------------------
    if( ShowPhi > 0 ) {
      "Abbildung innerhalb des Morse-Lemmas:";
      PhiG;
      "Vergleich:";
      "PhiG(f)= " + Show(jet(PhiG(f), K));
      "fi     = " + Show(fi);
    }
  
    Rang = 0;
    B    = maxideal(1);
    for( i=1; i<=n ; i++) { if(Abb[i] != 1) { Rang ++; B[rvar(x(i))] = 0; } }
    Phi  = ring_top,B;
    PhiG = Phi(PhiG);
    fi   = Phi(fi);
    v    = fi, PhiG;
    debug_log(2, "rank determined with Morse rg=", Rang);
    debug_log(1, "Rest singularity f=",Show(fi));
    return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc Coeff(poly f, list #)
{ 
//---------------------------- initialisation ---------------------------------
  poly   a, term;
  int    n, i;
  matrix K;

  n     = nvars(basering);
  i     = 1;
  term  = #[2];
  K     = coef(f, #[1]);

  while( (i<=ncols(K)) && (K[1,i] != term) )
  { i++;
    if(i>ncols(K)) { break; }
  }
  if(i<=ncols(K)) { a = K[2,i]; }
  if(i>ncols(K)) { a = 0; }

  return(a);
}

///////////////////////////////////////////////////////////////////////////////
proc ReOrder(poly f)
{
//---------------------------- initialisation ---------------------------------
  poly  result;
  ideal B = maxideal(1);
  int   i, n, Ctn, Ctv;
  map   conv;

  n = nvars(basering);
  Ctv = 1;              // Zahl der Vorhandenen Variablen
  Ctn = n;              // Zahl der Nicht-Vorhandenen Variablen
  def ring_top=basering;

  for( i=1; i<=n; i=i+1)
  { result = subst(f,x(i), 0) - f;
    if( result != 0 ) { B[rvar(x(i))] = x(Ctv); Ctv++; } 
    else { B[rvar(x(i))] = x(Ctn); Ctn--; }
  }

  conv = ring_top,B;
  return(conv);
}

///////////////////////////////////////////////////////////////////////////////
proc quickclass(poly f);
USAGE:    quickclass(f);         f=poly
RETURN:   Normal form of f in Arnold's list
REMARK:   try to determine the normal form of f by invariants, mainly by
          computing the Hilbert function of the Milnor albegra,
          no coordinate change is needed (see also proc 'milnorcode').
EXAMPLE:  example quickclass; shows an example
{
//---------------------------- initialisation ---------------------------------
  string Typ;
  int    cnt, K, Mu, corank;
  list   v;
  def ring_top=basering;
  // check basic condition on the basering.
  if(checkring()) { return(f); }
  if( f==0 ) { 
    "Normal form : 0";
    return(f);
  }
  if( jet(f,0)!=0 ) { 
    "Normal form : 1";
    return(f);
  }
  K, Mu, corank = basicinvariants(f);
  if(Mu<0) {
    debug_log(0, "The Milnor number of the function is infinite.");
    return(f);
  }

  // Do the classification of f
  // typ: list of types matching the milnorcode
  // cnt: number of matches found
  v = HKclass(milnorcode(f));
  Typ, cnt = v[1..2];
  "Singularity R-equivalent to :",Typ;
  if(cnt==0) {
    "Hilbert polynomial not recognised. Milnor code = ", milnorcode(f);
    return();
  }
  if(cnt==1) { 
    debug_log(1,"Getting normal form from database.");
    "normal form :",A_L(Typ);
    return(A_L(Typ));
  }
  // Hier nun der Fall cnt>1
  "Hilbert-Code of Jf^2";
  "We have ", cnt, "cases to test";
  Cubic(f);
  return(v);
}
example
{ "EXAMPLE:"; echo=2;
   ring r=0,(x,y,z),ds;
   poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
   quickclass(f);
}

///////////////////////////////////////////////////////////////////////////////
proc milnorcode (poly f, list #)
USAGE:    milnorcode(f[,e]); f=poly, e=int
RETURN:   intvec, coding the Hilbert function of the e-th Milnor algebra
          of f, i.e. of basering/(jacob(f)^e) (default e=1), according
          to proc Hcode
EXAMPLE:  example milnorcode; shows an example
{
  int  e=1;
  if(size(#)==1) { e=#[1]; }
  ideal jf=std(jacob(f)^e);
  intvec v=hilb(jf,2);
  
  return(Hcode(v));
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),ds;
  poly f=x2y+y3+z2;
  milnorcode(f);
  milnorcode(f,2);  // a big second argument may result in memory overflow
}

///////////////////////////////////////////////////////////////////////////////
proc Hcode (intvec v)
USAGE:    Hcode(v); v=intvec
RETURN:   intvec, coding v according to the number of successive
          repetitions of an entry
EXAMPLE:  example Hcode; shows an example.
{
  int    col, len, i, cur, cnt, maxcoef, nlen;
  intvec hil1, hil2;

  col      = 1;
  len      = size(v);
  v[len+1] = 0;

  init_debug();
  debug_log(1, "Hcode:", v );

  for(i=1; i<=len; i++) { if( v[i] > maxcoef) { maxcoef = v[i]; } }

  nlen       = 2*maxcoef - 1;
  hil1[nlen] = 0;
  hil2[nlen] = 0;

  for(i=1; i<=nlen; i++)
  { if( i > maxcoef) { hil2[i] = 2*maxcoef-i; }
    else { hil2[i] = i; }
  }

  for(i=1; i<=nlen; i++)
  { cnt=0;
    while( (col<=len) && (v[col] == hil2[i]) )
    { cnt++; col++; }
    hil1[i] = cnt;
  }
  return(hil1);
}
example
{ "EXAMPLE:"; echo=2;
  intvec v1 = 1,3,5,5,2;
  Hcode(v1);
  intvec v2 = 1,2,3,4,4,4,4,4,4,4,3,2,1;
  Hcode(v2);
}

///////////////////////////////////////////////////////////////////////////////
proc Cubic (poly f)
{ 
//---------------------------- initialisation ---------------------------------
  poly  f3;
  ideal Jf, Jf1, Jf2;
  int   Dim, Mult;

  f3 = jet(f, 3);
  if( jet(f,2) != 0) { return("2-jet non zero"); }
  if( f3 == 0 ) { return("null form"); }

  Jf1  = jacob(f3);
  Jf   = std(Jf1);
  Dim  = dim(Jf);
  Mult = mult(Jf);

  if(Dim == 0) { return("P[8]:smooth cubic"); } // x3 + y3 + z3 + axyz
  if(Dim == 1) { 
    if(Mult == 2) {
      Jf2  = wedge(jacob(Jf1),3-Dim), Jf1;
      Jf2  = std(Jf2);
      Dim  = dim(Jf2);
      if (Dim == 0) { return("R:conic + line"); }       // x3 + xyz
      if (Dim == 1) { return("Q:cuspidal cubic"); }  // x3 + yz2
    }
    if(Mult == 3) {
      Jf2 = wedge(jacob(Jf1),3-Dim), Jf1;
      Jf2 = std(Jf2);
      Dim = dim(Jf2);
      if(Dim == 0) { return("T:three lines"); } // xyz
      if(Dim == 1) { return("S:conic + tangent"); }     // x2z + yz2
    }
    if(Mult == 4) { return("U:three concurrent lines"); }       // x3 + xz2
  }
  if(Dim == 2) {
    if(Mult == 1) { return("V:doubleline + line"); }    // x2y
    if(Mult == 2) { return("V': tripple line"); }       // x3
  }
  if(Dim == 3) { return("P[9]:nodal cubic"); }  // x3 + y3 + xyz

  return("");
}

///////////////////////////////////////////////////////////////////////////////
proc parity  (int e)
USAGE:    parity()
{ 
  int r = e/2;
  if( 2*r == e ) { return(0); }
  return(1);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass (intvec sg)
{
//---------------------------- initialisation ---------------------------------
  int    cnt = 0;
  string SG_Typ = "";
  list   v;
  
  // if trace/debug mode not set, do it!
  init_debug();
  debug_log(1, "Milnor code : ", sg );
  if(size(sg) == 1) { v ="A["+string(sg[1])+"]", 1; return(v); }
  if(size(sg) == 3) { return(HKclass3(sg, SG_Typ, cnt)); }
  if(size(sg) == 5) { return(HKclass5(sg, SG_Typ, cnt)); }
  if(size(sg) == 7) { return(HKclass7(sg, SG_Typ, cnt)); }
  debug_log(1, "No solution found." );
  v = "", 0;
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass3 (intvec sg, string SG_Typ, int cnt)
{ 
  list v;

  if(sg[1] == 1) { v = HKclass3_teil_1(sg, SG_Typ, cnt); }
  debug_log(6, "HKclass3: ", v[1], " cnt=", v[2]);
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass3_teil_1 (intvec sg, string SG_Typ, int cnt)
{ 
  int  k, r, s;
  list v;

  debug_log(2, "entering HKclass3_teil_1", sg);
  if(sg[2]==1) { SG_Typ=SG_Typ+" D[k]=D["+string(sg[3]+3)+"]";cnt++;} // D[k]
  if(sg[2]>=1) {
    if( parity(sg[2])) { // sg[2] ist ungerade
      if(sg[2]<=sg[3]) {
        k = (sg[2]+1)/2; 
        if(k>1) { 
          cnt++;
          SG_Typ=SG_Typ+" J[k,r]=J["+string(k)+","+string(sg[3]+1-2*k)+"]";
        }// J[k,r]
      }
      if(sg[2]==sg[3]+2) {                              // E[6k+2]
        k = (sg[2]-1)/2;
        if(k>0) {cnt++; SG_Typ=SG_Typ+" E[6k+2]=E[" + string(6*k+2) + "]"; }
      }
    }
    else {              // sg[2] ist gerade
      if( sg[2] == sg[3]+1) {                           // E[6k]
        k = sg[2]/2; cnt++; SG_Typ=SG_Typ + " E[6k]=E[" + string(6*k) + "]"; }
      if( sg[2] == sg[3]) {                             // E[6k+1]
        k=sg[2]/2; cnt++; SG_Typ=SG_Typ+" E[6k+1]=E["+string(6*k+1)+"]"; }
    }
  }

  debug_log(2, "finishing HKclass3_teil_1");
  debug_log(6, "HKclass3: ", SG_Typ, " cnt=", cnt);
  v = SG_Typ, cnt;
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass5 (intvec sg, string SG_Typ, int cnt)
{ 
  list v;

  if(sg[1] == 1 && sg[2] == 1) { v = HKclass5_teil_1(sg, SG_Typ,cnt); }
  if(sg[1] == 1 && sg[2] == 0) { v = HKclass5_teil_2(sg, SG_Typ,cnt); }
  debug_log(6, "HKclass5: ", v[1], " cnt=", v[2]);
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass5_teil_1 (intvec sg, string SG_Typ, int cnt)
{
  int  k, r, s;
  list v;

  debug_log(2, "entering HKclass5_teil_1", sg);
  if(parity(sg[3])) {  // Dritte Stelle soll ungerade sein
    k = (sg[3]+1)/2;
    if(sg[3] > sg[4]) {
      k--;
      if( (sg[4]==sg[5]) && (sg[3] == sg[4]+1) && k>0 ) { // W[12k+6]
        SG_Typ = SG_Typ + " W[12k+6]=W["+string(12*k+6)+"]"; cnt++;
      }
      if( (sg[3]==sg[5]) && (sg[3] == sg[4]+2) && k>0 ) { // W[12k+5]
        SG_Typ = SG_Typ + " W[12k+5]=W["+string(12*k+5)+"]"; cnt++;
      }
    }
    else {  // sg[3] <= sg[4]
      if( (sg[3]==sg[4]) && (sg[5] >= sg[3]) ) {
        r = sg[5] - sg[4];
        SG_Typ=SG_Typ +" X[k,r]=X["+string(k)+","+string(r)+"]"; cnt++;
      } 
      if( (sg[3]==1) && (sg[4]==3) && (sg[5]>=sg[4])){    // Z[1,r]
        r = sg[5] - sg[4];
        SG_Typ = SG_Typ + " Z[1,r]=Z[1,"+string(r)+"]"; cnt++;
      }

      if( sg[4] == sg[5]) {
        if(parity(sg[4])) {                                  // Z[k,r,0]
          r = (sg[4] - sg[3])/2;
          if( r>0 ) { cnt++;
            SG_Typ = SG_Typ + " Z[k,r,0]=Z["+string(k)+","+string(r)+",0]";
          }
        } 
        else {                                                // Z[k,12k+6r]
          r = (sg[4] - 2*k)/2; cnt++;
          SG_Typ = SG_Typ+" Z[k,12k+6r]=Z["+string(k)+","+string(12*k+6*r)+"]";
        }
      }

      if( parity(sg[4]) ) {  // 4. Stelle ist ungerade
        if(sg[4] == sg[5]+2) {                              // Z[k,12k+6r+1]
          r = (sg[4]-2*k-1)/2; cnt++;
          SG_Typ=SG_Typ+" Z[k,12k+6r+1]=Z["+string(k)+",";
          SG_Typ=SG_Typ+string(12*k+6*r+1)+"]";
       }
       if( (sg[5]>sg[4]) && (sg[4]>sg[3]) ) {           // Z[k,r,s]
          r = (sg[4] - sg[3])/2; cnt++;
          s = sg[5] - sg[4];
          SG_Typ = SG_Typ + " Z[k,r,s]=";
          SG_Typ = SG_Typ + "Z["+string(k)+","+string(r)+","+string(s)+"]";
        }
      }
      else {  // 4. Stelle ist gerade
        if( sg[4] == sg[5]+1) {                             // Z[k,12k+6r-1]
          r = (sg[4] - 2*k)/2; cnt++;
          SG_Typ=SG_Typ+" Z[k,12k+6r-1]=Z["+string(k)+",";
          SG_Typ=SG_Typ+string(12*k+6*r-1)+"]";
        }
      }

      if(sg[4]>sg[3]) {                                     // Y[k,r,s]
        r = sg[4] - sg[3];
        s = sg[5] - sg[3] + r;
        if( s<0 ) { s = -s; }
        SG_Typ = SG_Typ + " Y[k,r,s]="; cnt++;
        SG_Typ = SG_Typ + "Y["+string(k)+","+string(r)+","+string(s)+"]";
      }
    }
  }
  else {  // Dritte Stelle soll gerade sein
    k = sg[3]/2;
    // sortiere verschiedene W's
    if(k>0) {
      if( (sg[4]==2*k-1) && (sg[4]==sg[5]) ) {  // W[12k]
        SG_Typ = SG_Typ + " W[12k]=W["+string(12*k)+"]"; cnt++;
      }
      if( (sg[4]==2*k-1) && (sg[3]==sg[5]) ) {  // W[12k+1]
        SG_Typ = SG_Typ + " W[12k+1]=W["+string(12*k+1)+"]"; cnt++;
      }
      if( (sg[4]==2*k) && (sg[5]>=sg[4]) ) {    // W[k,r]
        r = sg[5] - sg[4];
        SG_Typ=SG_Typ+" W[k,r]=W["+string(k)+","+string(r)+"]"; cnt++;
      }
      if( (sg[5]==2*k-1) && (sg[4]>sg[3]) ) {  // W#[k,2r-1]
        r = sg[4] - sg[3]; cnt++;
        SG_Typ = SG_Typ + " W#[k,2r-1]=W["+string(k)+","+string(2*r-1)+"]";
      }
      if( (sg[5]==2*k) && (sg[4]>sg[3]) ) {  // W#[k,2r]
        r = sg[4] - sg[3]; cnt++;
        SG_Typ = SG_Typ + " W#[k,2r]=W["+string(k)+","+string(2*r)+"]";
      }
    }   // ENDIF k>0
  }
  debug_log(2, "finishing HKclass5_teil_1");
  debug_log(6, "HKclass5_teil_1: ", SG_Typ, " cnt=", cnt);
  v = SG_Typ, cnt;
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass5_teil_2 (intvec sg, string SG_Typ, int cnt)
{
  int k, r, s;
  list v;

  debug_log(2, "entering HKclass5_teil_2", sg);
  // finde T[p,q,r]
  k = sg[3] + 1;
  r = sg[4] + k;
  s = sg[5] + r - 1;
  if(k>2 && r>2 && s>2) {                               // T[k,r,s]
    cnt++;
    SG_Typ = SG_Typ + " T[k,r,s]=T["+string(k)+","+string(r)+","+string(s)+"]";
  }

  // finde weitere Moeglicjkeiten.
  if(sg[3]==2) {  // Q[...]
    if(parity(sg[4])) { // 4. Stelle ist ungerade.
      if(sg[4]==sg[5]) {                                // Q[6k+4]
        k=(sg[4]+1)/2; cnt++; SG_Typ=SG_Typ+" Q[6k+4]=Q["+string(6*k+4)+"]";
      }
      if(sg[4]+1==sg[5]) {                      // Q[6k+5]
        k=sg[5]/2; cnt++; SG_Typ=SG_Typ+" Q[6k+5]=Q["+string(6*k+5)+"]";
      }
    }
    else { // 4. Stelle ist gerade.
      if(sg[4]==sg[5]+1) {                      // Q[6k+6]
        k=sg[4]/2; cnt++; SG_Typ=SG_Typ+" Q[6k+6]=Q["+string(6*k+6)+"]";
      }
      if(sg[4]<sg[5]) {                 // Q[k,r]
        k = (sg[4]+2)/2; 
        if(k>=2) {
          r=sg[5]+1-2*k; cnt++;
          SG_Typ=SG_Typ+" Q[k,r]=Q["+string(k)+","+string(r)+"]";
        }
      }
    }
  }
  else {           // S[...]
    if(parity(sg[3])) {  // 3. Stelle ist ungerade.
      k = (sg[3]-1)/2;
      if(sg[3]==sg[4]+3 && sg[3]==sg[5]+2) {    // S[12k-1]
        cnt++; SG_Typ = SG_Typ + " S[12k-1]=S["+string(12*k-1)+"]";
      }
      if(sg[3]==sg[4]+3 && sg[3]==sg[5]+1) {    // s[12k]
        cnt++; SG_Typ = SG_Typ + " S[12k]=S["+string(12*k)+"]";
      }
      if(sg[3]==sg[4]+2 && sg[5]>=sg[4]+1) {    // S[k,r]
        r = sg[5] - 2*k; cnt++;
        SG_Typ = SG_Typ + " S[k,r]=S["+string(k)+","+string(r)+"]";
      }
      if(sg[3]==sg[5]+2 && sg[4]>=sg[5]) {              // S#[k,2r-1]
        r = sg[4] - 2*k + 1; cnt++;
        SG_Typ = SG_Typ + " S#[k,2r-1]=S#["+string(k)+","+string(2*r-1)+"]";
      }
      if(sg[3]==sg[5]+1 && sg[4]>=sg[5]) {              // S#[k,2r]
        r = sg[4] - 2*k + 1; cnt++;
        SG_Typ = SG_Typ + " S#[k,2r]=S#["+string(k)+","+string(2*r)+"]";
      }
    }
    else { // 3. Stelle ist gerade.
      if(sg[3]==sg[5]+1 && sg[5]==sg[4]+3) {    // S[12k+4]
        k = (sg[3]-2)/2; cnt++;
        SG_Typ = SG_Typ + " S[12k+4]=S["+string(12*k+4)+"]";
      }
      if(sg[3]==sg[5]+2 && sg[5]==sg[4]+1) {    // S[12k+5]
        k = (sg[3]-2)/2; cnt++;
        SG_Typ = SG_Typ + " S[12k+5]=S["+string(12*k+5)+"]";
      }
    }
  }
  debug_log(2, "finishing HKclass5_teil_2");
  debug_log(6, "HKclass5_teil_2: ", SG_Typ, " cnt=", cnt);
  v = SG_Typ, cnt;
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass7 (intvec sg, string SG_Typ, int cnt)
{ 
  list v;

  if(sg[1]==1 && sg[2]==0 && sg[3]==1) { v=HKclass7_teil_1(sg, SG_Typ, cnt); }
  debug_log(6, "HKclass7: ", v[1], " cnt=", v[2]);
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc HKclass7_teil_1 (intvec sg, string SG_Typ, int cnt)
{
  int k, r, s;
  list v;

  debug_log(2, "entering HKclass7_teil_1", sg);
  if(sg[4] == 2) {                                      // V[...]
    if(sg[5] == 0 && sg[6] == 1 && sg[7]>0) {   // V[1,r]
      r = sg[7] - 1; cnt++; SG_Typ = SG_Typ + " V[1,r]=V[1,"+string(r)+"]";
    }
    if(sg[5] == 1 && sg[7] == 1) {                      // V#[1,2r-1]
      r=sg[6]+1; cnt++; SG_Typ=SG_Typ+" V#[1,2r-1]=V#[1,"+string(2*r-1)+"]";
    }
    if(sg[5] == 1 && sg[7] == 2) {                      // V#[1,2r]
      r=sg[6]+1; cnt++; SG_Typ=SG_Typ+" V#[1,2r]=V#[1,"+string(2*r)+"]";
    }
  }
  //            Moegliche U[...]'s
  k = sg[4];
  if(sg[5]==2*k-1 && sg[6]==0 && sg[7]==sg[5]) {        // U[12k]
    cnt++;SG_Typ = SG_Typ + " U[12k]=U["+string(12*k)+"]";
  }
  if(sg[5]==2*k && sg[6]==0 && sg[7]==sg[5]) {  // U[12k+4]
    cnt++;SG_Typ = SG_Typ + " U[12k+4]=U["+string(12*k+4)+"]";
  }
  if(sg[5]==2*k-1 && sg[6]>0 && sg[7]==sg[5]) { // U[k,2r-1]
    r=sg[6]-1; cnt++;
    SG_Typ=SG_Typ+" U[k,2r-1]=U["+string(k)+","+string(2*r-1)+"]";
  }
  if(sg[5]==2*k-1 && sg[6]>0 && sg[7]==2*k) {   // U[k,2r]
    r = sg[6]; cnt++;
    SG_Typ = SG_Typ + " U[k,2r]=U["+string(k)+","+string(2*r)+"]";
  }
  debug_log(2, "finishing HKclass7_teil_1");
  debug_log(6, "HKclass7_teil_1: ", SG_Typ, " cnt=", cnt);
  v = SG_Typ, cnt;
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc singularity(string typ, list #)
USAGE:    singularity(t, l); t=string (name of singularity), 
          l=list of integers (index/indices of singularity)
COMPUTE:  get the Singularity named by type t from the database.
          list l is as follows:
          l= k [,r [,s [,a [,b [,c [,d]]]]]] k,r,s=int   a,b,c,d=poly
          The name of the dbm-databasefile is: NFlist.[dir,pag]
          The file is found in the current directory. If it does not
          exists, please run the script MakeDBM first.
RETURN:   Normal form and corank of the singularity named by type t and its
          index (indices) l
EXAMPLE:  example singularity; shows an example
{
  poly a1, a2, a3, a4, f;
  int k, r, s;
  int len = size(#);
  list v, ret;

  classify_init();
  ret = 0, 0;
  k = #[1];
  if(len>=2) { r = #[2]; }
  else { r = 0; }
  if(len>=3) { s = #[3]; }
  else { s = 0; }
  if( k<0 || r<0 || s<0) { 
    "Initial condition failed: k>=0; r>=0; s>=0";
    "k="+string(k)+" r="+string(r)+"   s="+string(s);
    return(ret);
  }
  int crk;

  init_debug();
  def ring_top=basering;

  if(len>=4) { a1 = #[4]; } 
  else { a1=1; }
  if(len>=5) { a2 = #[5]; } 
  else { a2=1; }
  if(len>=6) { a3 = #[6]; } 
  else { a3=1; }
  if(len>=7) { a4 = #[7]; } 
  else { a4=1; }

  debug_log(4, "Values: len=", string(len), " k=", string(k), " r=", 
        string(r));
  if(defined(RingNF) != 0 ) { kill RingNF; }
  ring RingNF=char(basering),(x,y,z),(c,ds);
  poly f;
  map Conv=ring_top,maxideal(1);
  v = Singularitaet(typ, k, r, s, Conv(a1), Conv(a2), Conv(a4), Conv(a4));
  f, crk = v[1..2];
  debug_log(2, "Info=", f );
  setring ring_top;
  if(defined(Phi) != 0 ) { kill Phi; }
  map Phi=RingNF,maxideal(1);

  ret = Phi(f), crk;
  return(ret);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,y,z),(c,ds);
  init_debug(0);
  singularity("E[6k]",6);
  singularity("T[k,r,s]", 3, 7, 5);
  poly f=y;
  singularity("J[k,r]", 4, 0, 0, f);
}

///////////////////////////////////////////////////////////////////////////////
proc Singularitaet (string typ,int k,int r,int s,poly a,poly b,poly c,poly d)
{
  list   v;  
  string DBMPATH=system("getenv","DBMPATH");
  string DatabasePath, Database, S, Text, Tp;
  poly   f, f1;
  int    crk, Mu, ret;
  intvec MlnCd;

  if( DBMPATH != "" ) { DatabasePath = DBMPATH+"/NFlist"; }
  else { DatabasePath = "NFlist"; }
  Database="DBM: ",DatabasePath;

  link dbmLink=Database;
  debug_log(2, "Opening Singalarity-database: ", newline, Database);
  Tp = read(dbmLink, typ); 
  debug_log(2,"DBMread(", typ, ")=", Tp, ".");
  if( Tp != "(null)" && Tp !="" ) {
    string Key = "I_", typ;
    S = "f = ", Tp, ";";
    debug_log(2,"S=", S, " Tp=", Tp, "Key=", Key);
    execute S;
    execute read(dbmLink, Key)+";";
    debug_log(1, "Polynom f=", f,  "  crk=", crk, "  Mu=", Mu, 
                " MlnCd=", MlnCd);
    v = f, crk, Mu, MlnCd;
  }
  else {
    v = 0, 0, 0, 0;
  }
  close(dbmLink);
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc RandomPolyK (int M, string Typ)
USAGE:    RandomPolyK(M, Typ)
{
//---------------------------- initialisation ---------------------------------
  int    n, b, i, k, r, s, crk;
  ideal  B;
  map    Phi;
  string txt, Tp;
  list   v;

  def ring_ext=basering;
  n=4;
  if(M<5) { M = 5; }

  k = random(1, M);
  r = random(-5, 2*M);
  s = random(-5, 2*M);
  if(r<0) { r = 0; }
  if(s<0) { s = 0; }

  ring RgAnf=char(basering),(x,y,z,t),(c,ds);
  poly f;

  v = singularity(Typ, k, r, s);
  f, crk = v[1..2];
//  f = f +t2;
  if(crk==1) { f = f + y2 + z2; }
  if(crk==2) { f = f + z2; }
  txt="RandomPoly-Series: gewaehlt fall `"+Typ+"' mit";
  txt=txt+" f="+string(f);
  txt;
  setring ring_ext;
  B = maxideal(1);

  r=1;
  for(i=n; i>0; i--,r++) {
//  for(i=1; i<=n; i=i+1)
    B[rvar(x(r))] = x(i);
    if(i>2 && random(1,10)<3) { B[rvar(x(r))] = B[rvar(x(r))] + x(i-1); }
//    if(i==1 && random(1,10)<4) { B[rvar(x(r))] = B[rvar(x(r))]- x(n); }
    if(i>0) {
      for(b=3; b<5; b=b+1) {
        // B[rvar(x(r))] = B[rvar(x(r))] + random(0,9) * x(i)^(b+2);
        if(random(1,20)<3) {
          B[rvar(x(r))] = B[rvar(x(r))] - random(-2,2)*x(b)^2;
        }
      }
    }
  }
  Phi=RgAnf, B;
  Phi;
  poly fr=Phi(f);
  fr = fr+(x(1)+x(2))^2;
//  return(Phi(f));
  return(fr);
}

///////////////////////////////////////////////////////////////////////////////
proc debug_log (int level, list #)
USAGE:    debug_log(level,li); level=int, li=comma separated "message" list
COMPUTE:  print "messages" if level>=@DeBug.
          useful for user-defined trace messages.
EXAMPLE:  example debug_log; shows an example
SEE ALSO: init_debug();
{ 
   int len = size(#);
   if( defined(@DeBug) == 0 ) { init_debug(); }
   if(@DeBug>=level) { 
      if(level>1) { "Debug:("+ string(level)+ "): ", #[1..len]; }
      else { #[1..len]; }
   }
}
example
{ "EXAMPLE:"; echo=2;
  example init_debug;
}

///////////////////////////////////////////////////////////////////////////////
proc init_debug(list #)
USAGE:    init_debug([level]);  level=int
COMPUTE:  Set the global variable @DeBug to level. The variable @DeBug is 
          used by the function debug_log(level, list of strings) to know 
          when to print the list of strings. init_debug() reports only 
          changes of @DeBug.
NOTE:     The procedure init_debug(n); is usefull as trace-mode. n may
          range from 0 to 10, higher values of n give more information.
EXAMPLE:  example init_debug; shows an example
{ 
  int newDebug=0;
  if( defined(@DeBug) != 0 ) { newDebug = @DeBug; }

  if( size(#) > 0 ) { 
    newDebug=#[1];
  }
  else {
    string s=system("getenv", "SG_DEBUG");
    if( s != "" && defined(@DeBug)==0) {
      s="newDebug="+s;
      execute s;
    }
  }
  if( defined(@DeBug) == 0) { 
    int @DeBug = newDebug;
    export @DeBug;
    if(@DeBug>0) { "Debugging level is set to ", @DeBug; }
  }
  else {
    if( (size(#) == 0) && (newDebug < @DeBug) ) { return(); }
    if( @DeBug != newDebug) {
      int oldDebug = @DeBug;
      @DeBug = newDebug;
      if(@DeBug>0) { "Debugging level change from ", oldDebug, " to ",@DeBug; }
      else {
        if( @DeBug==0 && oldDebug>0 ) { "Debugging switched off."; }
      }
    }
  }
}
example
{ "EXAMPLE:"; echo=2;
  init_debug();
  debug_log(1,"no trace information printed");
  init_debug(1);
  debug_log(1,"some trace information");
  init_debug(2);
  debug_log(2,"nice for debugging scripts");
  init_debug(0);
}

///////////////////////////////////////////////////////////////////////////////
proc basicinvariants(poly f)
USAGE:    basicinvariants(f);   f = poly
COMPUTE:  Compute basic invariants of f: an upper bound d for the
          determinacy, the milnor number mu and the corank c of f
RETURN:   intvec: d, mu, c
EXAMPLE:  example basicinvariants; shows an example
{ 
  intvec v;
  ideal Jfs = std(jacob(f));
  v[1] = system("HC")+1;
  v[2] = vdim(Jfs);
  v[3] = corank(f);
  if( v[2]<v[1] ) { v[1] = v[2]+1; }
  return(v);
}
example
{ "EXAMPLE:"; echo=2;
   ring r=0,(x,y,z),ds;
   basicinvariants((x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3);
}

///////////////////////////////////////////////////////////////////////////////
proc corank(poly f)
USAGE:    corank(f);   f=poly
RETURN:   the corank of the Hessian matrix of f, of type int
REMARK:   corank(f) is the number of variables occuring in the residual
          singulartity after applying 'morsesplit' to f
EXAMPLE:  example corank; shows an example
{
  matrix M = jacob(jacob(jet(f,2)));
  int cr = nvars(basering) - size(module(transpose(bareiss(M))));
  return(cr);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),ds;
  poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
  corank(f);
}
///////////////////////////////////////////////////////////////////////////////
proc Faktorisiere(poly f, poly fk, int pt, int k, intvec RFlg)
{
//---------------------------- initialisation ---------------------------------
  poly   a, b, Relation;
  ideal  B, Jfsyz;
  map    PhiG, VERT;
  matrix Mat;
  list   v;
  def    ring_top=basering;

  // Ziel: bestimme a,b sodass  fk = (ax+by^k)^pt gilt.
  B    = maxideal(1);
  PhiG = ring_top,B;
  debug_log(2, "Faktor: f=",Show(f)," Jet=",Show(fk)," k=",k," exp=",pt);

//----------------------- compute role of x and y -----------------------------
  Jfsyz = fk, diff(fk, x(1));
  Mat   = matrix(syz(Jfsyz));
  if( (fk-subst(fk,x(1),0)) != 0  &&  (fk-subst(fk,x(2),0)) != 0 ) {
    // Wenn k>0 ist die Wahl fuer x & y bereits getroffen
    // sonst bestimmen x und y
    Jfsyz    = fk, diff(fk, x(1));
    Mat      = matrix(syz(Jfsyz));
    Relation = -pt * Mat[2,1] / Mat[1,1];
    a = Coeff(Relation, x(1), x(1));
    b = Coeff(Relation, x(2), x(2)^k);
    B = maxideal(1);
    if( (RFlg[1]==1 && k==1) || k>1) { B[rvar(x(1))] = x(1)-b*x(2)^k; }
    else { B[rvar(x(2))] = x(2)-b*x(1)^k; }
    VERT = basering,B;
    f    = VERT(f);
    PhiG = VERT(PhiG);
  }

//------------------- permutation of x and y, if needed -----------------------
  if( k==1 ) {
    debug_log(2, "Fak-7:",Show(f)," jet=",Show(fk));
    if(Coeff(jet(f, pt), x(1), x(1)^pt) == 0) {
      VERT = basering,x(2),x(1);
      f    = VERT(f);
      PhiG = VERT(PhiG);
    }
  }
  debug_log(2, "Fak-8:",Show(f)," jet=",Show(fk));
  debug_log(6, "Faktorisiere liefert: f=", Show(f));
  v[1] = f;
  v[2] = PhiG;
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc Teile(poly f, poly fk)
{
  ideal  Jfsyz = f, fk;
  poly   Relation;
  matrix Mat = matrix(syz(Jfsyz));
  Relation   = -1 * Mat[2,1]/Mat[1,1];
  return(Relation);
}

///////////////////////////////////////////////////////////////////////////////
proc GetRf (poly fi, int n)
USAGE:    GetRf();
{ 
//---------------------------- initialisation ---------------------------------
  int    j, k, l1, l1w;
  matrix Koef;
  intvec RFlg;

  RFlg[n] = 0;
  intvec Haeufigkeit = RFlg;

  for( j=1; j<=n ; j=j+1) {
    Koef = coef(fi, x(j));
    Haeufigkeit[j] = ncols(Koef);
    if(Coeff(fi, x(j),0) == 0) { Haeufigkeit[j] = Haeufigkeit[j] + 1;}
  }
  for( j=n; j>0 ; j=j-1) {
    l1  = 0;
    l1w = 0;
    for(k=1;k<=n;k=k+1) { if(Haeufigkeit[k]>l1w) { l1=k; l1w=Haeufigkeit[k];}}
    RFlg[j]        = l1;
    Haeufigkeit[l1] = 0;
  }
  debug_log(2, "Reihenfolge fuer Vertauschungen:", RFlg);
  return(RFlg);
}

///////////////////////////////////////////////////////////////////////////////
proc Show(poly g)
{ 
  string s;
  def ring_save=basering;

  execute @ringdisplay;
  map showpoly=ring_save,maxideal(1);
  s = string(showpoly(g));
  setring ring_save;
  return (s);
}

///////////////////////////////////////////////////////////////////////////////
proc checkring
{
  int CH = char(basering);
  if(CH >= 2 && CH<=13) {
    "Ring has characteristic ",CH;
    "Characteristic other than 0 or 0<char<13 is not yet implemented";
    return(1);
  }
  return(0);  // characteristic of ring is OK, return (0)
}

///////////////////////////////////////////////////////////////////////////////
proc DecodeNormalFormString (string S_in)
USAGE:    DecodeNormalFormString
{
//---------------------------- initialisation ---------------------------------
  int    C_eq, a, b, i, t, k, r, s;
  string s1, s2, s3, s4, s_in, Typ;
  list v = "Error", 0, 0, 0;

  C_eq = find(S_in, "=")+1;
  s_in = S_in[C_eq,30];
  debug_log(2, "Decode:");

  debug_log(2, "S_in=", S_in,"  s_in=",s_in );
  a = find(s_in, "[")+1;
  b = find(s_in, "]")-1;
  t = 1;
  k = 0;
  r = 0;
  s = 0;

  if(a<0 || b<0) { return(v); }
  Typ = s_in[1..a-1];
  s1  = s_in[a..b];
  debug_log(6, "Suche Type:", Typ);
  //---------------------- decode between brackets ----------------------------
  if( find(s1, ",") == 0) {
    debug_log(8, "  Number of columns: 0");
    s2 = "k = "+s1+";";
    execute s2;
    if( (Typ=="A[") || (Typ=="D[") ) { s3 = "k"; }
    if( Typ == "E[") { t = 6; }
    if( Typ == "W[") { t = 12; }
    if( Typ == "Q[") { t = 6; }
    if( Typ == "Z[") { t = 6; }
    if( Typ == "U[") { t = 12; }
    if( t > 1 ) {
      i = k;
      k = k/t;
      b = i - t*k;
      if( (s1 == "Q[") && (b==0) ) { k=k-1; b=6; }
      if(Typ == "Z[") {
        if(b==0) { k=k-1; b=6; }
        if(b==1) { k=k-1; b=7; }
      }
      if( b == 0 ) { s3 = string(t)+"k"; }
      else { s3 = string(t)+"k+"+string(b); }
    }
    if( Typ == "S[") {
      i = k+1;
      k = i/12;
      b = i - 12*k;
      if( b == 1 ) { s3 = "k"; }
      else { 
        if(b==0) { s3 = "12k"+string(b-1); }
        else { s3 = "12k+"+string(b-1); }
      }
    }
    s2 = Typ + s3 +"]";
  }  // es kommt mindestens ein komma vor...
  //----------------------- more than 1 paramater -----------------------------
  else {
    b  = find(s1, ",");
    s2 = "k = ",s1[1..b-1],";";
    execute s2;
    s1 = s1[b+1..size(s1)];
    if(find(s1, ",") == 0) {
      debug_log(8, "  Number of columns 1");
      s2 = "r = "+s1+";";
      execute s2;
      s4 = "r";
      s3 = "k";
      if(r==0) { s4 = string(0); }
      if(k==0 && Typ=="Z[") { s3 = string(1); }
      if(Typ[2] == "#") { 
        i = r+1;
        r = i/2;
        b = i - 2*r;
        if( b == 1 ) { s4 = "2r"; } 
        else { s4 = "2r-1"; }
      }
      s2 = Typ + s3 + "," + s4 +"]";
    }  // es kommt mindestens zwei komma vor...
    //----------------------- get third parameter -----------------------------
    else { 
      debug_log(8, "  Number of columns >=2");
      debug_log(2, "Y[k,r,s] / Z[k,r,s] / T[k,r,s]");
      b  = find(s1, ",");
      s2 = "r = ",s1[1..b-1],";";
      execute s2;
      s2 = "s = ",s1[b+1..size(s1)],";";
      execute s2;
      if(Typ=="Y[") { s2 = "Y[k,r,s]"; }
      if(Typ=="Z[") { s2 = "Z[k,r,s]"; }
      if(Typ=="T[") { s2 = "T[k,r,s]"; }
    }
  }
  debug_log(2, "Looking for Normalform of ",s2,"with (k,r,s) = (",
        k,",",r,",", s, ")" );
  v = s2, k, r, s;
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
proc A_L
USAGE:    A_L(f);         f=poly
          A_L("name");    type=string
COMPUTE:  the normal form in Arnold's list of the singularity given either
          by a polynomial f or by its name.
RETURN:   A_L(f): compute via 'milnorcode' the class of f and
          return the normal form of f found in the database.
          A_L("name"): Get the normal form from the database for
          the singularity given by its name.
EXAMPLE:  example A_L; shows an example
{ 
  // if trace/debug mode not set, do it!
  init_debug();

  if( typeof(#[1]) == "string" ) {
    if(checkring()) { return(#[1]); }
    return(normalform(#[1]));
  }
  if( typeof(#[1]) == "poly" ) {
    if(checkring()) { return(#[1]); }
    return(quickclass(#[1]));
  }
  
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(a,b,c),ds;  
  poly f=A_L("E[13]");
  f;
  A_L(f);
}

///////////////////////////////////////////////////////////////////////////////
proc normalform(string s_in)
USAGE:    normalform(s);  s=string
RETURN:   Arnold's normal form of singularity with name s
EXAMPLE:  example normalform; shows an example.
{ 
  string Typ;
  int    k, r, s, crk;
  poly   f;
  list   v;

  if(checkring()) { return(s_in); }
  if(nvars(basering)<=1) {
    "We need at least 2 variables in basering, you have",nvars(basering),".";
    return();
  }
  // if trace/debug mode not set, do it!
  init_debug();

  v = DecodeNormalFormString(s_in);
  Typ, k, r, s = v[1..4];
  if(Typ=="Error") { return(0); }
  v = singularity(Typ, k, r, s);
  f, crk = v[1..2];
  return(f);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(a,b,c),ds;  
  normalform("E[13]");
}

///////////////////////////////////////////////////////////////////////////////
proc swap 
USAGE:    swap(a,b);
RETURN:   b,a if b,a is the input (any type)
{
  return(#[2],#[1]); 
}
example
{ "EXAMPLE:"; echo=2;
  swap("variable1","variable2");
}

///////////////////////////////////////////////////////////////////////////////
proc Setring(int c, string name)
USAGE:    
{
  string s="ring "+name+"=0,(x(1.."+ string(c) +")),(c,ds);";
  return(s);
}

///////////////////////////////////////////////////////////////////////////////
proc internalfunctions
USAGE:   internalfunctions();
RETURN:  nothing, display names of internal procedures of classify.lib
EXAMPLE: no example
{ "   Internal functions for the classification using Arnold's method,";
  "   the function numbers correspond to numbers in Arnold's classifier:";
 "Klassifiziere(poly f);             determine the type of the singularity f
  Funktion1bis (poly f, list cstn)
  Funktion3 (poly f, list cstn)
  Funktion6 (poly f, list cstn)
  Funktion13 (poly f, list cstn)
  Funktion17 (poly f, list cstn)
  Funktion25 (poly f, list cstn)
  Funktion40 (poly f, list cstn, int k)
  Funktion47 (poly f, list cstn)
  Funktion50 (poly f, list cstn)
  Funktion58 (poly fin, list cstn)
  Funktion59 (poly f, list cstn)
  Funktion66 (poly f, list cstn)
  Funktion82 (poly f, list cstn)
  Funktion83 (poly f, list cstn)
  Funktion91 (poly f, list cstn, int k)
  Funktion92 (poly f, list cstn, int k)
  Funktion93 (poly f, list cstn, int k)
  Funktion94 (poly f, list cstn, int k)
  Funktion95 (poly f, list cstn, int k)
  Funktion96 (poly f, list cstn, int k)
  Funktion97 (poly f, list cstn)
  Isomorphie_s82_x (poly f, poly fk, int k)
  Isomorphie_s82_z (poly f, poly fk, int k)
  Isomorphie_s17 (poly f, poly fk, int k, int ct)
  printresult (string f, string typ, int Mu, int m, int corank, int K)
  ";
  "   Internal functions for the classifcation by invariants:
  Cubic (poly f)
  parity (int e)               return the parity of e
  HKclass (intvec i)
  HKclass3( intvec i, string SG_Typ, int cnt)
  HKclass3_teil_1 (intvec i, string SG_Typ, int cnt)
  HKclass5 (intvec i, string SG_Typ, int cnt)
  HKclass5_teil_1 (intvec i, string SG_Typ, int cnt)
  HKclass5_teil_2 (intvec i, string SG_Typ, int cnt)
  HKclass7 (intvec i, string SG_Typ, int cnt)
  HKclass7_teil_1 (intvec i, string SG_Typ, int cnt)
  ";
  "   Internal functions for the Morse-splitting lemma:
  Morse(poly fi, int K, int corank)          Splittinglemma itself
  Coeffs (list #)
  Coeff
  ";
  "   Internal functions providing tools:
  ReOrder(poly f)
  Singularitaet (string typ,int k,int r,int s,poly a,poly b,poly c,poly d)
  RandomPolyK
  Faktorisiere(poly f, poly g, int p, int k)     compute g = (ax+by^k)^p
  Teile(poly f, poly g);                  Teilt f durch g.
  GetRf(poly f, int n);
  Show(poly f);
  checkring();
  DecodeNormalFormString(string s);
  Setring(int n, string ringname);
  ";
}
///////////////////////////////////////////////////////////////////////////////
// E n d   O f   F i l e