singular-data 1:4.1.0-p3+ds-2build1 / usr / share / singular / LIB / classify.lib

///////////////////////////////////////////////////////////////////////////////
version="version classify.lib 4.0.2.0 Apr_2015 "; // $Id: e6af0122579da218b042132258501ba6c0bfad8f $
category="Singularities";
info="
LIBRARY:  classify.lib  Arnold Classifier of Singularities
AUTHOR:   Kai Krueger, krueger@mathematik.uni-kl.de

OVERVIEW:
   A library for classifying isolated hypersurface singularities w.r.t. right
   equivalence, based on the determinator of singularities by V.I. Arnold.

SEE ALSO:  realclassify_lib

PROCEDURES:
 basicinvariants(f); computes Milnor number, determinacy-bound and corank of
 classify(f);        normal form of polynomial 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 polynomial 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
 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
 swap(a,b);          swaps the arguments
 modality(f);        modality of the singularity
 complexSingType(f); complex type of the singularity as a string
 prepRealclassify(f);
                     the modality and the complex type of the singularity at
                     once
           (parameters in square brackets [] are optional)
";

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

///////////////////////////////////////////////////////////////////////////////
proc classify_init
{
  string s;
  link l="DBM:r NFlist";
  s = read(l,"VERSION");
  if (s == "" )
  {
    if (printlevel > 0)
    {
      "(classify.lib): Need to create database...";
    }
    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 differentiable
          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 beta. Please send bugs and
          comments to: \"Kai Krueger\" <krueger@mathematik.uni-kl.de> @*
          Be sure to have at least Singular version 1.0.1.
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;

  // Namespaces:
  if(defined(A_Z)==0) { proc A_Z = General::A_Z; export A_Z; }
  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));
  if(typeof(v[1])=="poly")
  {
     poly f_out = v[1];
     s2 = v[2];          // s2: Typ des Polynoms f z.b: E[18]
     corank = v[3];
  }
  else
  {
     s2="NoClass";
  }

//---------------- 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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 ------------------------------
    if(typeof(v_class[1])!="poly") { return(v); }
    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);
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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-3)+"]", 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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 = v[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));
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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));
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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*number(gamma) / (B - beta))*x(2);
        bb[rvar(x(2))] = x(2) - ((C - number(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);
}

///////////////////////////////////////////////////////////////////////////////
static proc printresult (int step, poly f, string typ, list cstn, int m)
{
  if(defined(onlyreturninvariants))
  {
    if(onlyreturninvariants == 1)
    {
      Modality = m;
      Type = typ;
    }
  }
//---------------------------- 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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));
  string tp="";
//  return(printresult(47, f, tp, cstn, -1));

  s = s +" has 4-jet equal to zero. (F47), mu="+string(Mu);

  s; // +"  ("+SG_Typ+")";
  s = "No further classification available.";
  s;
  return(Show(f), tp, corank);
}

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

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

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

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

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

///////////////////////////////////////////////////////////////////////////////
static 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;
}

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

///////////////////////////////////////////////////////////////////////////////
static 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, "Split the polynomial below using determinacy: ", string(K));
  debug_log(3, Show(fi));

  for( j=1; j<=n ; j++) { Abb[j] = 0; }

  RFlg = GetRf(fi, n);
  debug_log(2, "Permutations:", 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++;
      Jfx = f2, diff(f2, x(j));
    }
    Mat  = matrix(syz(Jfx));
    Beta = 2*Mat[2,1]/Mat[1,1];
    for( j=1; j<=n ; j++)
    {
      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++;
      }               // 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++;
    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, "Residual singularity f=",Show(fi));
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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 algebra,
          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);
}

///////////////////////////////////////////////////////////////////////////////
static 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("");
}

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

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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) div 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) div 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] div 2; cnt++; SG_Typ=SG_Typ + " E[6k]=E[" + string(6*k) + "]";
      }
      if( sg[2] == sg[3])
      {                             // E[6k+1]
        k=sg[2] div 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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) div 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]) div 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) div 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) div 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]) div 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) div 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] div 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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) div 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] div 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] div 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) div 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) div 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) div 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) div 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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);
  }
  else
  {
      v[1]="not in list";
      v[2]=0;
  }
  debug_log(6, "HKclass7: ", v[1], " cnt=", v[2]);
  return(v);
}

///////////////////////////////////////////////////////////////////////////////
static 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/polynomials (indices/parmeters 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
          exist, 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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(#);
//   int printresult = printlevel - level +1;
//   if(level>1)
//   {
//     dbprint(printresult, "Debug:("+ string(level)+ "): ", #[2..len]);
//   }
//   else { dbprint(printresult, #[1..len]); }
   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."; }
      }
    }
  }
  printlevel = @DeBug;
}
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
          singularity after applying 'morsesplit' to f
EXAMPLE:  example corank; shows an example"
{
  matrix M = jacob(jacob(jet(f,2)));
  list lba = bareiss(M);
  int cr = nvars(basering) - size(module(lba[1]));
  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);
}
///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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, "Permutations:", RFlg);
  return(RFlg);
}

///////////////////////////////////////////////////////////////////////////////
static 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);
}

///////////////////////////////////////////////////////////////////////////////
static 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)
}

///////////////////////////////////////////////////////////////////////////////
static 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 div 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 parameter -----------------------------
  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 div 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],";";
      if ((s2[5]>"0") && (s2[5]<="9")) { execute(s2); }
      s2 = "s = ",s1[b+1..size(s1)],";";
      if ((s2[5]>"0") && (s2[5]<="9")) { 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(s);  s string, the name of the singularity
COMPUTE:  the normal form of f in Arnold's list of singularities in case 1,
          in case 2 nothing has to be computed.
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;
  int    n, i;
  poly   f;
  list   v;
  def ring_ext = basering;
  n = nvars(basering);
  ring ring_top=char(basering),(x(1..n)),(c,ds);

  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);
  poly f_out;
  f_out, crk = v[1..2];
  if(crk>1) { for(i=crk+1;i<=n;i=i+1) { f_out = f_out + x(i)^2; } }
  setring ring_ext;
  map conv_top2ext=ring_top,maxideal(1);
  f = conv_top2ext(f_out);
//  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 a,b 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)  //splitting lemma 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);             //divides f by g
  GetRf(poly f, int n);
  Show(poly f);
  checkring();
  DecodeNormalFormString(string s);
  Setring(int n, string ringname);
  ";
}
example
{
  "EXAMPLE"; echo=2;
  internalfunctions();
}

///////////////////////////////////////////////////////////////////////////////
proc prepRealclassify(poly f)
"
USAGE:    prepRealclassify(f); f poly
RETURN:   a list, containing the modality of the singularity and the type of
          the singularity as a string
          @* This procedure is needed in realclassify.lib in order to avoid
          classify() being called more than once.
EXAMPLE:  example prepRealclassify; shows an example"
{
  exportinvariants(f);
  return(list(Modality, Type));
}
example
{
  "EXAMPLE"; echo = 2;
  ring r = 0, (x,y,z), ds;
  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
  prepRealclassify(f);
}

///////////////////////////////////////////////////////////////////////////////
proc modality(poly f)
"
USAGE:    modality(f); f poly
RETURN:   the modality of the singularity
EXAMPLE:  example modality; shows an example"
{
  exportinvariants(f);
  return(Modality);
}
example
{
  "EXAMPLE"; echo = 2;
  ring r = 0, (x,y,z), ds;
  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
  modality(f);
}

///////////////////////////////////////////////////////////////////////////////
proc complexSingType(poly f)
"
USAGE:    complexSingType(f); f poly
RETURN:   the type of the singularity as a string
EXAMPLE:  example complexSingType; shows an example"
{
  exportinvariants(f);
  return(Type);
}
example
{
  "EXAMPLE"; echo = 2;
  ring r = 0, (x,y,z), ds;
  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
  complexSingType(f);
}

///////////////////////////////////////////////////////////////////////////////
/* some of the invariants will be exported during the computation of
 * classify(f) such that they are accessible at for all procs in the library
 */
static proc exportinvariants(poly f)
{
  init_debug(-1);
  if(!defined(onlyreturninvariants))
  {
    int onlyreturninvariants;
    export(onlyreturninvariants);
  }
  onlyreturninvariants = 1;
  if(!defined(Modality))
  {
    int Modality;
    export(Modality);
  }
  if(!defined(Type))
  {
    string Type;
    export(Type);
  }
  f = classify(f);
}

///////////////////////////////////////////////////////////////////////////////
// E n d   O f   F i l e