version="$Id: hnoether.lib,v 1.29.2.14 2002/10/21 14:41:56 lossen Exp $";
category="Singularities";
info="
LIBRARY:  hnoether.lib   Hamburger-Noether (Puiseux) Development
AUTHOR:   Martin Lamm,   lamm@mathematik.uni-kl.de

OVERVIEW:
 A library for computing the Hamburger-Noether, resp. Puiseux, development
 of a plane curve singularity following [Campillo, A.: Algebroid curves
 in positive characteristic, Springer LNM 813 (1980)]. @*
 The library contains also procedures for computing the (topological)
 numerical invariants of plane curve singularities.

MAIN PROCEDURES:
 hnexpansion(f);            Hamburger-Noether (H-N) development of f
 sethnering(hn);            changes to the hnering created by hnexpansion
 develop(f [,n]);           H-N development of irreducible curves
 extdevelop(hne,n);         extension of the H-N development hne of f
 parametrisation(f [,x]);   a parametrization of f 
 displayHNE(hne);           display H-N development as an ideal
 invariants(f);             invariants of f, e.g. the characteristic exponents
 displayInvariants(f);      display invariants of f
 multsequence(f);           sequence of multiplicities
 displayMultsequence(f);    display sequence of multiplicities
 intersection(hne1,hne2);   intersection multiplicity of two curves
 stripHNE(hne);             reduce amount of memory consumed by hne
 is_irred(f);               test if f is irreducible
 delta(f);                  delta invariant of f
 newtonpoly(f);             (local) Newton polygon of f
 is_NND(f);                 test if f is Newton non-degenerate


AUXILIARY PROCEDURES:
 puiseux2generators(m,n);   convert Puiseux pairs to generators of semigroup
 separateHNE(hne1,hne2);    number of quadratic transf. needed for separation
 squarefree(f);             a squarefree divisor of the poly f
 allsquarefree(f,l);        the maximal squarefree divisor of the poly f
 further_hn_proc();         show further procedures useful for interactive use

KEYWORDS: Hamburger-Noether expansion; Puiseux expansion; curve singularities
";

// HNdevelop(f);              Hamburger-Noether (H-N) development of f
// reddevelop(f);             H-N development of reducible curves
// essdevelop(f);             H-N development of essential branches
// multiplicities(hne);       multiplicities of blowed up curves


///////////////////////////////////////////////////////////////////////////////
LIB "primitiv.lib";
LIB "inout.lib";
LIB "sing.lib";

///////////////////////////////////////////////////////////////////////////////

proc further_hn_proc()
"USAGE: further_hn_proc();
NOTE:  The library @code{hnoether.lib} contains some more procedures which
       are not shown when typing @code{help hnoether.lib;}. They may be useful
       for interactive use (e.g. if you want to do the calculation of an HN
       development \"by hand\" to see the intermediate results), and they
       can be enumerated by calling @code{further_hn_proc()}. @*
       Use @code{help <procedure>;} for detailed information about each of
       them.
"
{
 "
 The following procedures are also part of `hnoether.lib':

 getnm(f);           intersection pts. of Newton polygon with axes
 T_Transform(f,Q,N); returns f(y,xy^Q)/y^NQ (f: poly, Q,N: int)
 T1_Transform(f,d,M); returns f(x,y+d*x^M)  (f: poly,d:number,M:int)
 T2_Transform(f,d,M,N,ref);   a composition of T1 & T
 koeff(f,I,J);       gets coefficient of indicated monomial of poly f
 redleit(f,S,E);     restriction of monomials of f to line (S-E)
 leit(f,n,m);        special case of redleit (for irred. polynomials)
 testreducible(f,n,m); tests whether f is reducible
 charPoly(f,M,N);    characteristic polynomial of f
 find_in_list(L,p);  find int p in list L
 get_last_divisor(M,N); last divisor in Euclid's algorithm
 factorfirst(f,M,N); try to factor f without `factorize'
 factorlist(L);      factorize a list L of polynomials
 referencepoly(D);   a polynomial f s.t. D is the Newton diagram of f";

//       static procedures not useful for interactive use:
// polytest(f);        tests coefficients and exponents of poly f
// extractHNEs(H,t);   extracts output H of HN to output of reddevelop
// HN(f,grenze);       recursive subroutine for reddevelop
// constructHNEs(...); subroutine for HN
}
example
{ echo=2;
  further_hn_proc();
}
///////////////////////////////////////////////////////////////////////////////

proc getnm (poly f)
"USAGE:   getnm(f); f bivariate polynomial
RETURN:  intvec(n,m) : (0,n) is the intersection point of the Newton
         polygon of f with the y-axis, n=-1 if it doesn't exist
         (m,0) is its intersection point with the x-axis,
         m=-1 if this point doesn't exist
ASSUME:  ring has ordering `ls' or `ds'
EXAMPLE: example getnm; shows an example
"
{
 // assume being called by develop ==> ring ordering is ls (ds would also work)
 return(ord(subst(f,var(1),0)),ord(subst(f,var(2),0)));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y),ds;
   poly f = x5+x4y3-y2+y4;
   getnm(f);
}
///////////////////////////////////////////////////////////////////////////////

proc leit (poly f, int n, int m)
"USAGE:   leit(f,n,m);  poly f, int n,m
RETURN:  all monomials on the line from (0,n) to (m,0) in the Newton diagram
EXAMPLE: example leit;  shows an example
"
{
 return(jet(f,m*n,intvec(n,m))-jet(f,m*n-1,intvec(n,m)))
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y),ds;
   poly f = x5+x4y3-y2+y4;
   leit(f,2,5);
}
///////////////////////////////////////////////////////////////////////////////
proc testreducible (poly f, int n, int m)
"USAGE:   testreducible(f,n,m);  f poly, n,m int
RETURN:  1 if there are points in the Newton diagram below the line (0,n)-(m,0)
         0 else
EXAMPLE: example testreducible;  shows an example
"
{
 return(size(jet(f,m*n-1,intvec(n,m))) != 0)
}
example
{ "EXAMPLE:"; echo = 2;
  ring rg=0,(x,y),ls;
  testreducible(x2+y3-xy4,3,2);
}
///////////////////////////////////////////////////////////////////////////////
proc T_Transform (poly f, int Q, int N)
"USAGE:   T_Transform(f,Q,N);  f poly, Q,N int
RETURN:  f(y,xy^Q)/y^NQ   if x,y are the ring variables
NOTE:    this is intended for irreducible power series f
EXAMPLE: example T_Transform;  shows an example
"
{
 map T = basering,var(2),var(1)*var(2)^Q;
 return(T(f)/var(2)^(N*Q));
}
example
{ "EXAMPLE:"; echo = 2;
  ring exrg=0,(x,y),ls;
  export exrg;
  T_Transform(x3+y2-xy3,1,2);
  kill exrg;
}
///////////////////////////////////////////////////////////////////////////////
proc T1_Transform (poly f, number d, int Q)
"USAGE:   T1_Transform(f,d,Q);  f poly, d number, Q int
RETURN:  f(x,y+d*x^Q)   if x,y are the ring variables
EXAMPLE: example T1_Transform;  shows an example
"
{
 map T1 = basering,var(1),var(2)+d*var(1)^Q;
 return(T1(f));
}
example
{ "EXAMPLE:"; echo = 2;
  ring exrg=0,(x,y),ls;
  export exrg;
  T1_Transform(y2-2xy+x2+x2y,1,1);
  kill exrg;
}
///////////////////////////////////////////////////////////////////////////////

proc T2_Transform (poly f, number d, int M, int N, poly refpoly)
"USAGE:   T2_Transform(f,d,M,N,ref); f poly, d number; M,N int; ref poly
RETURN:  list: poly T2(f,d',M,N), number d' in \{ d, 1/d \}
ASSUME:  ref has the same Newton polygon as f (but can be simpler)
         for this you can e.g. use the proc `referencepoly' or simply f again
COMMENT: T2 is a composition of T_Transform and T1_Transform; the exact
         definition can be found in  Rybowicz: `Sur le calcul des places ...'
         or in  Lamm: `Hamburger-Noether-Entwicklung von Kurvensingularitaeten'
SEE ALSO: T_Transform, T1_Transform, referencepoly
EXAMPLE: example T2_Transform;  shows an example
"
{
 //---------------------- compute gcd and extgcd of N,M -----------------------
  int ggt=gcd(M,N);
  M=M/ggt; N=N/ggt;
  list ts=extgcd(M,N);
  int tau,sigma=ts[2],-ts[3];
  int s,t;
  poly xp=var(1);
  poly yp=var(2);
  poly hilf;
  if (sigma<0) { tau=-tau; sigma=-sigma;}
 // es gilt: 0<=tau<=N, 0<=sigma<=M, |N*sigma-M*tau| = 1 = ggT(M,N)
  if (N*sigma < M*tau) { d = 1/d; }
 //--------------------------- euklid. Algorithmus ----------------------------
  int R;
  int M1,N1=M,N;
  for ( R=M1%N1; R!=0; ) { M1=N1; N1=R; R=M1%N1;}
  int Q=M1 / N1;
  map T1 = basering,xp,yp+d*xp^Q;
  map Tstar=basering,xp^(N-Q*tau)*yp^tau,xp^(M-sigma*Q)*yp^sigma;
  if (defined(HNDebugOn)) {
   "Trafo. T2: x->x^"+string(N-Q*tau)+"*y^"+string(tau)+", y->x^"
    +string(M-sigma*Q)+"*y^"+string(sigma);
   "delt =",d,"Q =",Q,"tau,sigma =",tau,sigma;
  }
 //------------------- Durchfuehrung der Transformation T2 --------------------
  f=Tstar(f);
  refpoly=Tstar(refpoly);
 //--- dividiere f so lange durch x & y, wie die Div. aufgeht, benutze ein  ---
 //--- Referenzpolynom mit gleichem Newtonpolynom wie f zur Beschleunigung: ---
  for (hilf=refpoly/xp; hilf*xp==refpoly; hilf=refpoly/xp) {refpoly=hilf; s++;}
  for (hilf=refpoly/yp; hilf*yp==refpoly; hilf=refpoly/yp) {refpoly=hilf; t++;}
  f=f/(xp^s*yp^t);
  return(list(T1(f),d));
}
example
{ "EXAMPLE:"; echo = 2;
  ring exrg=0,(x,y),ds;
  export exrg;
  poly f=y2-2x2y+x6-x5y+x4y2;
  T2_Transform(f,1/2,4,1,f);
  T2_Transform(f,1/2,4,1,referencepoly(newtonpoly(f,1)));
  // if  size(referencepoly) << size(f)  the 2nd example would be faster
  referencepoly(newtonpoly(f,1));
  kill exrg;
}
///////////////////////////////////////////////////////////////////////////////

proc koeff (poly f, int I, int J)
"USAGE:   koeff(f,I,J); f bivariate polynomial, I,J integers
RETURN:  if f = sum(a(i,j)*x^i*y^j), then koeff(f,I,J)= a(I,J) (of type number)
NOTE:    J must be in the range of the exponents of the 2nd ring variable
EXAMPLE: example koeff;  shows an example
"
{
  matrix mat = coeffs(coeffs(f,var(2))[J+1,1],var(1));
  if (size(mat) <= I) { return(0);}
  else { return(leadcoef(mat[I+1,1]));}
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  koeff(x2+2xy+3xy2-x2y-2y3,1,2);
}
///////////////////////////////////////////////////////////////////////////////

proc squarefree (poly f)
"USAGE:  squarefree(f);  f poly
ASSUME:  f is a bivariate polynomial (in the first 2 ring variables).
RETURN:  poly, a squarefree divisor of f.
NOTE:    Usually, the return value is the greatest squarefree divisor, but
         there is one exception: factors with a p-th root, p the
         characteristic of the basering, are lost.
SEE ALSO: allsquarefree
EXAMPLE: example squarefree; shows some examples.
"
{
 //----------------- Wechsel in geeigneten Ring & Variablendefinition ---------
  def altring = basering;
  int e;
  int gcd_ok=1;
  string mipl="0";
  if (size(parstr(altring))==1) { mipl=string(minpoly); }
 //---- test: char = (p^k,a) (-> gcd not implemented) or (p,a) (gcd works) ----
  if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
    string tststr=charstr(basering);
    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
    gcd_ok=(tststr==string(char(basering)));
  }
  execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;");
  if ((gcd_ok!=0) && (mipl!="0")) { execute("minpoly="+mipl+";"); }
  poly f=fetch(altring,f);
  poly dif,g,l;
  if ((char(basering)==0) and (charstr(basering)!=string(char(basering)))
      and (mipl!="0")) {
    gcd_ok=0;                   // since Singular 1.2 gcd no longer implemented
  }
  if (gcd_ok!=0) {
 //--------------------- Berechne f/ggT(f,df/dx,df/dy) ------------------------
    dif=diff(f,x);
    if (dif==0) { g=f; }        // zur Beschleunigung
    else { g=gcd(f,dif); }
    if (g!=1) {                 // sonst schon sicher, dass f quadratfrei
     dif=diff(f,y);
     if (dif!=0) { g=gcd(g,dif); }
    }
    if (g!=1) {
     e=0;
     if (g==f) { l=1; }         // zur Beschleunigung
     else {
       module m=syz(ideal(g,f));
       if (deg(m[2,1])>0) {
         "!! The Singular command 'syz' has returned a wrong result !!";
         l=1;                   // Division f/g muss aufgehen
       }
       else { l=m[1,1]; }
     }
    }
    else { e=1; }
  }
  else {
 //------------------- Berechne syz(f,df/dx) oder syz(f,df/dy) ----------------
 //-- Achtung: Ist f reduzibel, koennen Faktoren mit Ableitung Null verloren --
 //-- gehen! Ist aber nicht weiter schlimm, weil char (p^k,a) nur im irred.  --
 //-- Fall vorkommen kann. Wenn f nicht g^p ist, wird auf jeden Fall         --
 //------------------------ ein Faktor gefunden. ------------------------------
    dif=diff(f,x);
    if (dif == 0) {
     dif=diff(f,y);
     if (dif==0) { e=2; l=1; } // f is of power divisible by char of basefield
     else { l=syz(ideal(dif,f))[1,1];  // x^p+y^(p-1) abgedeckt
            if (subst(f,x,0)==0) { l=l*x; }
            if (deg(l)==deg(f))  { e=1;}
            else {e=0;}
     }
    }
    else { l=syz(ideal(dif,f))[1,1];
           if (subst(f,y,0)==0) { l=l*y; }
           if (deg(l)==deg(f))  { e=1;}
           else {e=0;}
    }
  }
 //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses --------
  setring altring;
  if (e==1) { return(f); }    // zur Beschleunigung
  else {
   poly l=fetch(rsqrf,l);
   return(l);
  }
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=3,(x,y),dp;
 squarefree((x3+y)^2);
 squarefree((x+y)^3*(x-y)^2); // Warning: (x+y)^3 is lost
 squarefree((x+y)^4*(x-y)^2); // result is (x+y)*(x-y)
}
///////////////////////////////////////////////////////////////////////////////

proc allsquarefree (poly f, poly l)
"USAGE : allsquarefree(f,g);  f,g poly
ASSUME: g is the output of @code{squarefree(f)}.
RETURN: the greatest squarefree divisor of f.
NOTE  : This proc uses factorize to get the missing factors of f not in g and,
        therefore, may be slow.
SEE ALSO: squarefree
EXAMPLE: example allsquarefree;  shows an example
"
{
 //------------------------ Wechsel in geeigneten Ring ------------------------
 def altring = basering;
 string mipl="0";
 if (size(parstr(altring))==1) { mipl=string(minpoly); }
 if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
   string tststr=charstr(basering);
   tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
   if (tststr!=string(char(basering))) {
     " Sorry -- not implemented for this ring (gcd doesn't work)";
     return(l);
   }
 }
 execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;");
 if (mipl!="0") { execute("minpoly="+mipl+";"); }
 poly f=fetch(altring,f);
 poly l=fetch(altring,l);
 //---------- eliminiere bereits mit squarefree gefundene Faktoren ------------
 poly g=l;
 while (deg(g)!=0) {
   f=syz(ideal(g,f))[1,1];                         // f=f/g;
   g=gcd(f,l);
 }                                                 // jetzt f=h^p
 //--------------- Berechne uebrige Faktoren mit factorize --------------------
 if (deg(f)>0) {
  g=1;
//*CL old:  ideal factf=factorize(f,1);
//*         for (int i=1; i<=size(factf); i++) { g=g*factf[i]; }
  ideal factf=factorize(f)[1];
  for (int i=2; i<=size(factf); i++) { g=g*factf[i]; }
  poly testp=squarefree(g);
  if (deg(testp)<deg(g)) {
    "!! factorize has not worked correctly !!";
    if (testp==1) {" We cannot proceed ..."; g=1;}
    else {" But we could recover some factors..."; g=testp;}
  }
  l=l*g;
 }
 //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses --------
 setring altring;
 l=fetch(rsqrf,l);
 return(l);
}
example
{ "EXAMPLE:"; echo = 2;
  ring exring=7,(x,y),dp;
  poly f=(x+y)^7*(x-y)^8;
  poly g=squarefree(f);
  g;                      // factor x+y lost, since characteristic=7
  allsquarefree(f,g);     // all factors (x+y)*(x-y) found
}
///////////////////////////////////////////////////////////////////////////////

proc is_irred (poly f)
"USAGE:   is_irred(f); f poly
ASSUME:  f is a squarefree bivariate polynomial (in the first 2 ring
         variables).
RETURN:  int (0 or 1): @*
         - @code{is_irred(f)=1} if f is irreducible as a formal power
         series over the algebraic closure of its coefficient field (f
         defines an analytically irreducible curve at zero), @*
         - @code{is_irred(f)=0} otherwise.
NOTE:    0 and units in the ring of formal power series are considered to be
         not irreducible.
KEYWORDS: irreducible power series
EXAMPLE: example is_irred;  shows an example
"
{
  int pl=printlevel;
  printlevel=-1;
  list hnl=develop(f,-1);
  printlevel=pl;
  return(hnl[5]);
}
example
{ "EXAMPLE:"; echo = 2;
  ring exring=0,(x,y),ls;
  is_irred(x2+y3);
  is_irred(x2+y2);
  is_irred(x2+y3+1);
}
///////////////////////////////////////////////////////////////////////////////

static proc polytest(poly f)
"USAGE : polytest(f); f poly in x and y
RETURN: a monomial of f with |coefficient| > 16001
          or exponent divisible by 32003, if there is one
        0 else (in this case computing a squarefree divisor
                in characteristic 32003 could make sense)
NOTE:   this procedure is only useful in characteristic zero, because otherwise
        there is no appropriate ordering of the leading coefficients
"
{
 poly verbrecher=0;
 intvec leitexp;
 for (; (f<>0) and (verbrecher==0); f=f-lead(f)) {
  if ((leadcoef(f)<-16001) or (leadcoef(f)>16001)) {verbrecher=lead(f);}
  leitexp=leadexp(f);
  if (( ((leitexp[1] % 32003) == 0)   and (leitexp[1]<>0))
     or ( ((leitexp[2] % 32003) == 0) and (leitexp[2]<>0)) )
       {verbrecher=lead(f);}
 }
 return(verbrecher);
}

//////////////////////////////////////////////////////////////////////////////


proc develop
"USAGE:   develop(f [,n]); f poly, n int
ASSUME:  f is a bivariate polynomial (in the first 2 ring variables) and
         irreducible as power series (for reducible f use @code{hnexpansion}).
RETURN:  list @code{L} with:
@texinfo
@table @asis
@item @code{L[1]}; matrix:
         Each row contains the coefficients of the corresponding line of the
         Hamburger-Noether expansion (HNE). The end of the line is marked in
         the matrix by the first ring variable (usually x).
@item @code{L[2]}; intvec:
         indicating the length of lines of the HNE
@item @code{L[3]}; int:
         0  if the 1st ring variable was transversal (with respect to f), @*
         1  if the variables were changed at the beginning of the
            computation, @*
        -1  if an error has occurred.
@item @code{L[4]}; poly:
         the transformed polynomial of f to make it possible to extend the
         Hamburger-Noether development a posteriori without having to do
         all the previous calculation once again (0 if not needed)
@item @code{L[5]}; int:
         1  if the curve has exactly one branch (i.e., is irreducible), @*
         0  else (i.e., the curve has more than one HNE, or f is not valid).
@end table
@end texinfo
DISPLAY: The (non zero) elements of the HNE (if not called by another proc).
NOTE:    The optional parameter @code{n} affects only the computation of
         the LAST line of the HNE. If it is given, the HN-matrix @code{L[1]}
         will have at least @code{n} columns. @*
         Otherwise, the number of columns will be chosen minimal such that the
         matrix contains all necessary information (i.e., all lines of the HNE
         but the last (which is in general infinite) have place). @*
         If @code{n} is negative, the algorithm is stopped as soon as the
         computed information is sufficient for @code{invariants(L)}, but the
         HN-matrix @code{L[1]} may still contain undetermined elements, which
         are marked with the 2nd variable (of the basering). @*
         For time critical computations it is recommended to use
         @code{ring ...,(x,y),ls} as basering - it increases the algorithm's
         speed. @*
         If @code{printlevel>=0} comments are displayed (default is
         @code{printlevel=0}).
SEE ALSO: hnexpansion, extdevelop, displayHNE
EXAMPLES: example develop;         shows an example
          example paramametrize;   shows an example for using the 2nd parameter
"
{
 //--------- Abfangen unzulaessiger Ringe: 1) nur eine Unbestimmte ------------
 poly f=#[1];
 if (size(#) > 1) {int maxspalte=#[2];}
 else             {int maxspalte= 1 ; }
 if (nvars(basering) < 2) {
   " Sorry. I need two variables in the ring.";
   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0),0));}
 if (nvars(basering) > 2) {
   dbprint(printlevel-voice+2,
   " Warning! You have defined too many variables!
 All variables except the first two will be ignored!"
           );
 }
 
 string namex=varstr(1); string namey=varstr(2);
 list return_error=matrix(maxideal(1)[2]),intvec(0),int(-1),poly(0),int(0);

 //------------- 2) mehrere Unbestimmte, weitere unzulaessige Ringe -----------
 // Wir koennen einheitlichen Rueckgabewert benutzen, aus dem ersichtlich ist,
 // dass ein Fehler aufgetreten ist: return_error.
 //----------------------------------------------------------------------------

 if (charstr(basering)=="real") {
  " The algorithm doesn't work with 'real' as coefficient field.";
                     // denn : map from characteristic -1 to -1 not implemented
  return(return_error);
 }
 if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
 //-- teste, ob char = (p^k,a) (-> a primitiv; erlaubt) oder (p,a[,b,...]) ----
    string tststr=charstr(basering);
    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
    int primit=(tststr==string(char(basering)));
    if (primit!=0) {
      " Such extensions of Z/p are not implemented.";
      " Please try (p^k,a) as ground field or use `hnexpansion'.";
      return(return_error);
    }
 }
 //---- Ende der unzulaessigen Ringe; Ringwechsel in einen guenstigen Ring: ---

 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");

 def altring = basering;
 if (ringwechsel) {
   string mipl=string(minpoly);
   execute("ring guenstig = ("+charstr(altring)+"),(x,y),ls;");
   if ((char(basering)==0) && (mipl!="0")) {
     execute("minpoly="+mipl+";");
   }}
 else { def guenstig=basering; }
 export guenstig;

 //-------------------------- Initialisierungen -------------------------------
 map m=altring,x,y;
 if (ringwechsel) { poly f=m(f); }
 if (defined(HNDebugOn))
 {"received polynomial: ",f,", where x =",namex,", y =",namey;}

 int M,N,Q,R,l,e,hilf,eps,getauscht,Abbruch,zeile,exponent,Ausgabe;

 // Werte von Ausgabe: 0 : normale HNE-Matrix,
 // 1 : Fehler aufgetreten - Matrix (namey) zurueck
 // 2 : Die HNE ist eine Nullzeile - Matrix (0) zurueck
 // int maxspalte=1; geaendert: wird jetzt am Anfang gesetzt

 int minimalHNE=0;          // Flag fuer minimale HNE-Berechnung
 int einzweig=1;            // Flag fuer Irreduzibilit"at
 intvec hqs;                // erhaelt die Werte von h(zeile)=Q;

 if (maxspalte<0) {
   minimalHNE=1;
   maxspalte=1;
 }

 number c,delt;
 int p = char(basering);
 string ringchar=charstr(basering);
 map xytausch = basering,y,x;
 if ((p!=0) and (ringchar != string(p))) {
                            // coefficient field is extension of Z/pZ
   execute("int n_elements="+
           ringchar[1,size(ringchar)-size(parstr(basering))-1]+";");
                            // number of elements of actual ring
   number generat=par(1);   // generator of the coefficient field of the ring
 }


 //========= Abfangen von unzulaessigen oder trivialen Eingaben ===============
 //------------ Nullpolynom oder Einheit im Potenzreihenring: -----------------
 if (f == 0) {
   dbprint(printlevel+1,"The given polynomial is the zero-polynomial !");
   Abbruch=1; Ausgabe=1;
 }
 else {
   intvec nm = getnm(f);
   N = nm[1]; M = nm[2]; // Berechne Schnittpunkte Newtonpolygon mit Achsen
   if (N == 0) {
     dbprint(printlevel+1,"The given polynomial is a unit as power series !");
     Abbruch=1; Ausgabe=1;
   }
   else {
    if (N == -1) {
      if ((voice==2) && (printlevel > -1)) { "The HNE is x = 0"; }
      Abbruch=1; Ausgabe=2; getauscht=1;
      if (M <> 1) { einzweig=0; }
    }
    else {
     if (M == -1) {
       if ((voice==2) && (printlevel > -1)) { "The HNE is y = 0"; }
       Abbruch=1; Ausgabe=2;
       if (N <> 1) { einzweig=0; }
   }}}
 }
 //--------------------- Test auf Quadratfreiheit -----------------------------
 if (Abbruch==0) {

 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring guenstig
 //----------------------------------------------------------------------------

  if ((p==0) and (size(charstr(basering))==1)) {
   int testerg=(polytest(f)==0);
   ring zweitring = 32003,(x,y),dp;
   map polyhinueber=guenstig,x,y;     // fetch geht nicht
   poly f=polyhinueber(f);
   poly test_sqr=squarefree(f);
   if (test_sqr != f) {
    if (printlevel>0) {
      "Most probably the given polynomial is not squarefree. But the test was";
      "made in characteristic 32003 and not 0 to improve speed. You can";
      "(r) redo the test in char 0 (but this may take some time)";
      "(c) continue the development, if you're sure that the polynomial",
      "IS squarefree";
      if (testerg==1) {
        "(s) continue the development with a squarefree factor (*)";}
      "(q) or just quit the algorithm (default action)";
      "";"Please enter the letter of your choice:";
      string str=read("")[1];
    }
    else { string str="r"; }      // printlevel <= 0: non-interactive behaviour
    setring guenstig;
    map polyhinueber=zweitring,x,y;
    if (str=="r") {
      poly test_sqr=squarefree(f);
      if (test_sqr != f) {
       if (printlevel>0) { "The given polynomial is in fact not squarefree."; }
       else              { "The given polynomial is not squarefree!"; }
       "I'll continue with the radical.";
       if (printlevel>0) { pause("Hit RETURN to continue:"); }
       f=test_sqr;
      }
      else {
       dbprint(printlevel,
        "everything is ok -- the polynomial is squarefree in char(k)=0");
      }
    }
    else {
      if ((str=="s") and (testerg==1)) {
       "(*) attention: it could be that the factor is only one in char 32003!";
        f=polyhinueber(test_sqr);
      }
      else {
        if (str<>"c") {
          setring altring;kill guenstig;kill zweitring;
          return(return_error);}
        else { "if the algorithm doesn't terminate, you were wrong...";}
    }}
    kill zweitring;
    nm = getnm(f);             // N,M haben sich evtl. veraendert
    N = nm[1]; M = nm[2];      // Berechne Schnittpunkte Newtonpoly mit Achsen
    if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
   }
   else {
     setring guenstig;
     kill zweitring;
   }
  }
 // ------------------- Fall Charakteristik > 0 -------------------------------
  else {
   poly test_sqr=squarefree(f);
   if (test_sqr == 1) {
    "The given polynomial is of the form g^"+string(p)+", therefore",
    "reducible.";"Please try again.";
    setring altring;
    kill guenstig;
    return(return_error);}
   if (test_sqr != f) {
    "The given polynomial is not squarefree. I'll continue with the radical.";
    if (p != 0)
     {"But if the polynomial contains a factor of the form g^"+string(p)+",";
      "this factor will be lost.";}
    if (printlevel>0) { pause("Hit RETURN to continue:"); }
    f=test_sqr;
    nm = getnm(f);              // N,M haben sich veraendert
    N = nm[1]; M = nm[2];       // Berechne Schnittpunkte Newtonpoly mit Achsen
    if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
   }

  }                             // endelse(p==0)

  if (N==0) {
    " Sorry. The remaining polynomial is a unit in the power series ring...";
    setring altring;kill guenstig;return(return_error);
  }
 //---------------------- gewaehrleiste, dass x transvers ist -----------------
  if (M < N)
  { f = xytausch(f);            // Variablentausch : x jetzt transvers
    getauscht = 1;              // den Tausch merken
    M = M+N; N = M-N; M = M-N;  // M, N auch vertauschen
  }
  if (defined(HNDebugOn)) {
   if (getauscht) {"x<->y were exchanged; poly is now ",f;}
   else           {"x , y were not exchanged";}
   "M resp. N are now",M,N;
  }
 }                              // end(if Abbruch==0)

 ideal a(0);
 while (Abbruch==0) {

 //================= Beginn der Schleife (eigentliche Entwicklung) ============

 //------------------- ist das Newtonpolygon eine gerade Linie? ---------------
  if (testreducible(f,N,M)) {
    dbprint(printlevel+1," The given polynomial is not irreducible");
    kill guenstig;
    setring altring;
    return(return_error);       // Abbruch der Prozedur!
  }
  R = M%N;
  Q = M / N;

 //-------------------- Fall Rest der Division R = 0 : ------------------------
  if (R == 0) {
    c = koeff(f,0,N);
    if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
    e = gcd(M,N);
 //----------------- Test, ob leitf = c*(y^N - delta*x^(m/e))^e ist -----------
    if (p==0) {
      delt = koeff(f,M/ e,N - N/ e) / (-1*e*c);
      if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
         leit(f,N,M)," = ",c,"* (y -",delt,"* x^"+string(M/ e)+")^",e," ?";}
      if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) {
        dbprint(printlevel+1," The given polynomial is reducible !");
        Abbruch=1; Ausgabe=1; }
    }
    else {                     // p!=0
      if (e%p != 0) {
        delt = koeff(f,M/ e,N - N/ e) / (-1*e*c);
        if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
           leit(f,N,M)," = ",c,"* (y -",delt,"* x^"+string(M/ e)+")^",e," ?";}
        if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) {
           dbprint(printlevel+1," The given polynomial is reducible !");
           Abbruch=1; Ausgabe=1; }
      }

      else {                   // e%p == 0
        eps = e;
        for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
        if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;}
        delt = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);

        if ((ringchar != string(p)) and (delt != 0)) {
 //- coeff. field is not Z/pZ => we`ve to correct delta by taking (p^l)th root-
          if (delt == generat) {exponent=1;}
          else {
           if (delt == 1) {exponent=0;}
           else {
            exponent=pardeg(delt);

 //-- an dieser Stelle kann ein Fehler auftreten, wenn wir eine transzendente -
 //-- Erweiterung von Z/pZ haben: dann ist das hinzuadjungierte Element kein  -
 //-- Erzeuger der mult. Gruppe, d.h. in Z/pZ (a) gibt es i.allg. keinen      -
 //-- Exponenten mit z.B. a2+a = a^exp                                        -
 //----------------------------------------------------------------------------
          }}
          delt = generat^(extgcd(n_elements-1,p^l)[3]*exponent);
        }

        if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
          leit(f,N,M)," = ",c,"* (y^"+string(N/ e),"-",delt,"* x^"
          +string(M/ e)+")^",e,"  ?";}
        if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) {
          dbprint(printlevel+1," The given polynomial is reducible !");
          Abbruch=1; Ausgabe=1; }
      }
    }
    if (Abbruch == 0) {
      f = T1_Transform(f,delt,M/ e);
      dbprint(printlevel-voice+2,"a("+string(zeile)+","+string(Q)+") = "+string(delt)); 
      a(zeile)[Q]=delt;
      if (defined(HNDebugOn)) {"transformed polynomial: ",f;}}

      nm=getnm(f); N=nm[1]; M=nm[2];        // Neuberechnung des Newtonpolygons
  }
 //--------------------------- Fall R > 0 : -----------------------------------
  else {
    dbprint(printlevel-voice+2, "h("+string(zeile)+ ") ="+string(Q)); 
    hqs[zeile+1]=Q;                  // denn zeile beginnt mit dem Wert 0
    a(zeile)[Q+1]=x;                 // Markierung des Zeilenendes der HNE
    maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
                                     // Anpassung der Sp.zahl der HNE-Matrix
    f = T_Transform(f,Q,N);
    if (defined(HNDebugOn)) {"transformed polynomial: ",f;}
    zeile=zeile+1;
 //------------ Bereitstellung von Speicherplatz fuer eine neue Zeile: --------
    ideal a(zeile);
    M=N;N=R;
  }

 //--------------- schneidet das Newtonpolygon beide Achsen? ------------------
  if (M==-1) {
     dbprint(printlevel-voice+2,"The HNE is finite!"); 
     a(zeile)[Q+1]=x;   // Markiere das Ende der Zeile
     hqs[zeile+1]=Q;
     maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
     if (N <> 1) { einzweig=0; }
     f=0;               // transformiertes Polynom wird nicht mehr gebraucht
     Abbruch=1;
  }
  else {if (M<N) {"Something has gone wrong: M<N";}}
  if(defined(HNDebugOn)) {"new M,N:",M,N;}

 //----------------- Abbruchbedingungen fuer die Schleife: --------------------
  if ((N==1) and (Abbruch!=1) and ((M > maxspalte) or (minimalHNE==1))
      and (size(a(zeile))>0))
 //----------------------------------------------------------------------------
 // Abbruch, wenn die Matrix so voll ist, dass eine neue Spalte angefangen
 // werden muesste und die letzte Zeile nicht nur Nullen enthaelt
 // oder wenn die Matrix nicht voll gemacht werden soll (minimale Information)
 //----------------------------------------------------------------------------
   { Abbruch=1; hqs[zeile+1]=-1;
     if (maxspalte < ncols(a(zeile))) { maxspalte=ncols(a(zeile));}
     if ((minimalHNE==1) and (M <= maxspalte)) {
 // teile param mit, dass Eintraege der letzten Zeile nur teilw. richtig sind:-
       hqs[zeile+1]=-M;
 //------------- markiere den Rest der Zeile als unbekannt: -------------------
       for (R=M; R <= maxspalte; R++) { a(zeile)[R]=y;}
     }                  // R wird nicht mehr gebraucht
   }
 //========================= Ende der Schleife ================================

 }
 setring altring;
 if (Ausgabe == 0) {
 //-------------------- Ergebnis in den alten Ring transferieren: -------------
   map zurueck=guenstig,maxideal(1)[1],maxideal(1)[2];
   matrix amat[zeile+1][maxspalte];
   ideal uebergabe;
   for (e=0; e<=zeile; e=e+1) {
     uebergabe=zurueck(a(e));
     if (ncols(uebergabe) > 1) {
      amat[e+1,1..ncols(uebergabe)]=uebergabe;}
     else {amat[e+1,1]=uebergabe[1];}
   }
   if (ringwechsel) {
     if (nvars(altring)==2) { f=fetch(guenstig,f); }
     else                   { f=zurueck(f); }
   }
 }

 kill guenstig;
 if ((einzweig==0) && (voice==2) && (printlevel > -1)) {
    "// Note: The curve is reducible, but we were able to compute a HNE.";
    "// This means the result is only one of several existing HNE's.";
 }
 if (Ausgabe == 0) { return(list(amat,hqs,getauscht,f,einzweig));}
 if (Ausgabe == 1) { return(return_error);}             // error has occurred
 if (Ausgabe == 2) { return(list(matrix(ideal(0,x)),intvec(1),getauscht,
                                 poly(0),einzweig));}   // HNE is x=0 or y=0
}
example
{ "EXAMPLE:"; echo = 2;
  ring exring = 7,(x,y),ds;
  list hne=develop(4x98+2x49y7+x11y14+2y14);
  print(hne[1]);
  // therefore the HNE is:
  // z(-1)= 3*z(0)^7 + z(0)^7*z(1),
  // z(0) = z(1)*z(2),       (there is 1 zero in the 2nd row before x)
  // z(1) = z(2)^3*z(3),     (there are 3 zeroes in the 3rd row)
  // z(2) = z(3)*z(4),
  // z(3) = -z(4)^2 + 0*z(4)^3 +...+ 0*z(4)^8 + ?*z(4)^9 + ...
  // (the missing x in the last line indicates that it is not complete.)
  hne[2];
  parametrisation(hne);
  // parametrization:   x(t)= -t^14+O(t^21),  y(t)= -3t^98+O(t^105)
  // (the term -t^109 in y may have a wrong coefficient)
  displayHNE(hne);
}

///////////////////////////////////////////////////////////////////////////////
//               procedures to extract information out of HNE                //
///////////////////////////////////////////////////////////////////////////////

proc parametrisation
"USAGE: parametrisation(INPUT [,x]); INPUT list or poly, x int (optional)
ASSUME: INPUT is either a bivariate polynomial f defining a plane curve 
        singularity, or it is the output of @code{hnexpansion(f[,\"ess\"])}, 
        or of @code{develop(f)}, or of @code{extdevelop(develop(f),n)}, 
        or  the list @{hne} in the ring created by @code{hnexpansion(f)} 
        respectively one entry thereof.
RETURN: a list L containing a parametrization L[i] for each branch f[i] of f 
        in the following format: @*
        - if only the list INPUT is given, L[i] is an ideal of two polynomials 
        p[1],p[2]: if the HNE of was finite  then f[i](p[1],p[2])=0; if not, 
        the \"real\" parametrization will be two power series and p[1],p[2] are
        truncations of these series.@*
        - if the optional parameter x is given, L[i] is itself a list: 
        L[i][1] is the parametrization ideal as above and L[i][2] is an intvec 
        with two entries indicating the highest degree up to which the 
        coefficients of the monomials in L[i][1] are exact (entry -1 means that
        the corresponding parametrization is exact).
NOTE:   If the basering has only 2 variables, the first variable is chosen
        as indefinite. Otherwise, the 3rd variable is chosen. @*
        In case the Hamburger-Noether expansion of the curve f is needed
        for other purposes as well it is better to calculate this first
        with the aid of @code{hnexpansion} and use it as input instead of
        the polynomial itself.
SEE ALSO: develop, extdevelop
KEYWORDS: parametrization
EXAMPLE: example parametrisation;     shows an example
         example develop;         shows another example
"
{
  ////////////////////////////////////////////////////////////////////////
  // Decide on the type of input
  ////////////////////////////////////////////////////////////////////////
  // Do the case where the input is a polynomial
  if (typeof(#[1])=="poly")
  {
    list HNEXPANSION=hnexpansion(#[1]);
    if (size(#)==1)
    {
      return(parametrisation(HNEXPANSION));
    }
    else
    {
      return(parametrisation(HNEXPANSION,1));
    }
  }
  // Do the case where the input is not a polynomial.
  int zusatz;  
  if (typeof(#[1])=="list")
  {
    if (typeof(#[1][1])=="ring")
    { // Input is a HNEring and extra input x exists.
      zusatz=1;
      def HNE_RING=#[1][1];
    }
    else
    {
      if (typeof(#[1][1])=="list")
      { // Input is a reducible HN-expansion and extra input x exists.
        zusatz=1;
        list hne=#[1];
      }
      else
      {
        if (typeof(#[size(#)])=="list")  
        { // Input is a reducible HN-expansion and no extra input exists.
          list hne=#;
        }
        else  
        { // Input is an irreducible HN-expansion and extra input x exists 
          list hne;
          hne[1]=#[1];
          zusatz=1;
        }
      }
    }
  }
  if (typeof(#[1])=="matrix")
  {
    list hne;
    hne[1]=#;
  }
  if (typeof(#[1])=="ring")
  { // Input is a HNEring and no extra input exists.
    def HNE_RING=#[1];
  }
  ////////////////////////////////////////////////////////////////////////////
  // Calculate the parametrization.
  if (defined(HNE_RING))
  {
    def rettering=basering;
    setring HNE_RING;
  }
  int r=size(hne);
  list ErGeBnIs;
  // Calculate the parametrization for each branch with the aid of param.
  for (int lauf=1;lauf<=r;lauf++)
  {
    if (zusatz==1)
    {
      ErGeBnIs[lauf]=param(hne[lauf],1);
    }
    else
    {
      ErGeBnIs[lauf]=param(hne[lauf]);
    }
  }
  // Map the parametrization to the basering, if necessary, and return it.
  if (defined(HNE_RING))
  {
    setring rettering;
    list erg=fetch(HNE_RING,ErGeBnIs);
    kill HNE_RING;
    // If the basering has 3 variables, choose the third variable for the parametrization.
    if (nvars(rettering)>=3)
    {
      for (lauf=1;lauf<=r;lauf++)
      {
        if (zusatz==1)
        {
          erg[lauf][1]=subst(erg[lauf][1],var(1),var(3));
        }
        else
        {
          erg[lauf]=subst(erg[lauf],var(1),var(3));
        }
      }
    }
    return(erg);
  }
  else
  {
    return(ErGeBnIs);
  }
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y,t),ds;
 // 1st Example: input is a polynomial
 poly g=(x2-y3)*(x3-y5);
 parametrisation(g);
 // 2nd Example: input is the ring of a Hamburger-Noether expansion
 poly h=x2-y2-y3;
 list hn=hnexpansion(h);
 parametrisation(h,1);
 // 3rd Example: input is a Hamburger-Noether expansion
 poly f=x3+2xy2+y2;
 list hne=develop(f);
 list hne_extended=extdevelop(hne,10);
            //   compare the matrices ...
 print(hne[1]);
 print(hne_extended[1]);
            // ... and the resulting parametrizations:
 parametrisation(hne);
 parametrisation(hne_extended);
 parametrisation(hne_extended,0);
}


proc param
"USAGE:  param(L [,x]); L list, x any type (optional)
ASSUME:  L is the output of @code{develop(f)}, or of
        @code{extdevelop(develop(f),n)}, or one entry in the list @code{hne}
        in the ring created by @code{hnexpansion(f[,\"ess\"])}.
RETURN: a parametrization for f in the following format: @*
        - if only the list L is given, the result is an ideal of two
        polynomials p[1],p[2]: if the HNE was finite then f(p[1],p[2])=0};
        if not, the \"real\" parametrization will be two power series and
        p[1],p[2] are truncations of these series.@*
        - if the optional parameter x is given, the result is a list l:
        l[1]=param(L) (ideal) and l[2]=intvec with two entries indicating
        the highest degree up to which the coefficients of the monomials in
        l[1] are exact (entry -1 means that the corresponding parametrization
        is exact).
NOTE:   If the basering has only 2 variables, the first variable is chosen
        as indefinite. Otherwise, the 3rd variable is chosen.
SEE ALSO: parametrization, develop, extdevelop
KEYWORDS: parametrization
EXAMPLE: example param;     shows an example
         example develop;   shows another example
"
{
 //-------------------------- Initialisierungen -------------------------------
 if (typeof(#[1])=="list") {
   list Li=#[1];
   matrix m=Li[1];
   intvec v=Li[2];
   int switch=Li[3];
   int return_list=1;
 }
 else {
   matrix m=#[1];
   intvec v=#[2];
   int switch=#[3];
   int return_list=0;
 }
 if (switch==-1) {
   "An error has occurred in develop, so there is no HNE.";
   return(ideal(0,0));
 }
 int fehler,fehlervor,untergrad,untervor,beginn,i,zeile,hilf;

 if (nvars(basering) > 2) { poly z(size(v)+1)=var(3); }
 else                     { poly z(size(v)+1)=var(1); }
 poly z(size(v));
 zeile=size(v);
 //------------- Parametrisierung der untersten Zeile der HNE -----------------
 if (v[zeile] > 0) {
   fehler=0;           // die Parametrisierung wird exakt werden
   for (i=1; i<=v[zeile]; i++) {
     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
 }}
 else {
   untervor=1;         // = Untergrad der vorhergehenden Zeile
   if (v[zeile]==-1) {
     fehler=ncols(m)+1;
     for (i=1; i<=ncols(m); i++) {
       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
                       // = Untergrad der aktuellen Zeile
   }}
   else {
     fehler= -v[zeile];
     for (i=1; i<-v[zeile]; i++) {
       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
   }}
 }
 //------------- Parametrisierung der restlichen Zeilen der HNE ---------------
 for (zeile=size(v)-1; zeile>0; zeile--) {
   poly z(zeile);
   beginn=0;             // Beginn der aktuellen Zeile
   for (i=1; i<=v[zeile]; i++) {
     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
     if ((beginn==0) and (m[zeile,i]!=0)) { beginn=i;}
   }
   z(zeile)=z(zeile) + z(zeile+1)^v[zeile] * z(zeile+2);
   if (beginn==0) {
     if (fehler>0) {     // damit fehler=0 bleibt bei exakter Param.
     fehlervor=fehler;   // Fehler der letzten Zeile
     fehler=fehler+untergrad*(v[zeile]-1)+untervor;   // Fehler dieser Zeile
     hilf=untergrad;
     untergrad=untergrad*v[zeile]+untervor;
     untervor=hilf;}}    // untervor = altes untergrad
   else {
     fehlervor=fehler;
     fehler=fehler+untergrad*(beginn-1);
     untervor=untergrad;
     untergrad=untergrad*beginn;
 }}
 //--------------------- Ausgabe der Fehlerabschaetzung -----------------------
 if (switch==0) {
   if (fehler>0) {
     if (fehlervor>0) {
       if ((voice==2) && (printlevel > -1)) {
         "// ** Warning: result is exact up to order",fehlervor-1,"in",
         maxideal(1)[1],"and",fehler-1,"in",maxideal(1)[2],"!";
     }}
     else {
       if ((voice==2) && (printlevel > -1)) {
         "// ** Warning: result is exact up to order",fehler-1,"in",
         maxideal(1)[2],"!";
     }}
   }
   if (return_list==0) { return(ideal(z(2),z(1))); }
   else   { return(list(ideal(z(2),z(1)),intvec(fehlervor-1,fehler-1))); }
 }
 else {
   if (fehler>0) {
     if (fehlervor>0) {
      if ((voice==2) && (printlevel > -1)) {
        "// ** Warning: result is exact up to order",fehler-1,"in",
        maxideal(1)[1],"and",fehlervor-1,"in",maxideal(1)[2],"!";
     }}
     else {
      if ((voice==2) && (printlevel > -1)) {
        "// ** Warning: result is exact up to order",fehler-1,"in",
        maxideal(1)[1],"!";
     }}
   }
   if (return_list==0) { return(ideal(z(1),z(2))); }
   else   { return(list(ideal(z(1),z(2)),intvec(fehler-1,fehlervor-1))); }
 }
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y,t),ds;
 poly f=x3+2xy2+y2;
 list hne=develop(f);
 list hne_extended=extdevelop(hne,10);
            //   compare the matrices ...
 print(hne[1]);
 print(hne_extended[1]);
            // ... and the resulting parametrizations:
 param(hne);
 param(hne_extended);
 param(hne_extended,0);
}
///////////////////////////////////////////////////////////////////////////////

proc invariants
"USAGE:   invariants(INPUT); INPUT list or poly
ASSUME:  INPUT is the output of @code{develop(f)}, or of
         @code{extdevelop(develop(f),n)}, or one entry in the list @code{hne} 
         of the HNEring created by @code{hnexpansion}.
RETURN:  list, if INPUT contains a valid HNE:
@format
    invariants(INPUT)[1]:  intvec    (characteristic exponents)
    invariants(INPUT)[2]:  intvec    (generators of the semigroup)
    invariants(INPUT)[3]:  intvec    (Puiseux pairs, 1st components)
    invariants(INPUT)[4]:  intvec    (Puiseux pairs, 2nd components)
    invariants(INPUT)[5]:  int       (degree of the conductor)
    invariants(INPUT)[6]:  intvec    (sequence of multiplicities)
@end format
         an empty list, if INPUT contains no valid HNE.
ASSUME:  INPUT is bivariate polynomial f or the output of @code{hnexpansion(f[,\"ess\"])}, 
         or the list @code{hne} in the HNEring created by @code{hnexpansion}.
RETURN:  list INV, such that INV[i] is the output of @code{invariants(develop(f[i]))}
         as above, where f[i] is the ith branch of the curve f, and the last 
         entry contains further invariants of f in the format:
@format
    INV[i][1]    : intvec    (characteristic exponents)
    INV[i][2]    : intvec    (generators of the semigroup)
    INV[i][3]    : intvec    (Puiseux pairs, 1st components)
    INV[i][4]    : intvec    (Puiseux pairs, 2nd components)
    INV[i][5]    : int       (degree of the conductor)
    INV[i][6]    : intvec    (sequence of multiplicities)
    INV[last][1] : intmat    (contact matrix of the branches)
    INV[last][2] : intmat    (intersection multiplicities of the branches)
    INV[last][3] : int       (delta invariant of f)
@end format
NOTE:    In case the Hamburger-Noether expansion of the curve f is needed
         for other purposes as well it is better to calculate this first
         with the aid of @code{hnexpansion} and use it as input instead of
         the polynomial itself.
SEE ALSO: develop, displayInvariants, multsequence, intersection
KEYWORDS: characteristic exponents; semigroup of values; Puiseux pairs;
          conductor, degree; multiplicities, sequence of
EXAMPLE:  example invariants; shows an example
"
{
 //---- INPUT = poly, or HNEring, or hne of reducible curve  --------------------
  if (typeof(#[1])!="matrix")
  {
    if (typeof(#[1])=="poly")
    {
      list HNEXPANSION=hnexpansion(#[1]);
      return(invariants(HNEXPANSION));
    }
    if (typeof(#[1])=="ring")
    {
      def H_N_E_R_I_N_G=#[1];
      def rette_ring=basering;
      setring H_N_E_R_I_N_G;
    }
    if (typeof(#[1])=="list")
    {
      list hne=#;
    }
    list ErGeBnIs;
    for (int lauf=1;lauf<=size(hne);lauf++)
    {
      ErGeBnIs[lauf]=invariants(hne[lauf]);
    }
    // Calculate the intersection matrix and the intersection multiplicities.
    intmat contact[size(hne)][size(hne)];
    intmat intersectionmatrix[size(hne)][size(hne)];
    int Lauf;
    for (lauf=1;lauf<=size(hne);lauf++)
    {
      for (Lauf=lauf+1;Lauf<=size(hne);Lauf++)
      {
        contact[lauf,Lauf]=separateHNE(hne[lauf],hne[Lauf]);
        contact[Lauf,lauf]=contact[lauf,Lauf];
        intersectionmatrix[lauf,Lauf]=intersection(hne[lauf],hne[Lauf]);
        intersectionmatrix[Lauf,lauf]=intersectionmatrix[lauf,Lauf];
      }
    }
    // Calculate the delta invariant.
    int inters;
    int del=ErGeBnIs[size(hne)][5]/2;
    for(lauf=1;lauf<=size(hne)-1;lauf++)
    {
      del=del+ErGeBnIs[lauf][5]/2;
      for(Lauf=lauf+1;Lauf<=size(hne);Lauf++)
      {
         inters=inters+intersectionmatrix[lauf,Lauf];
      }
    }
    del=del+inters;
    list LAST=contact,intersectionmatrix,del;
    ErGeBnIs[size(hne)+1]=LAST;
    if (defined(H_N_E_R_I_N_G))
    {
      setring rette_ring;
      kill H_N_E_R_I_N_G;
    }
    return(ErGeBnIs);
  }
 //-------------------------- Initialisierungen -------------------------------
 matrix m=#[1];
 intvec v=#[2];
 int switch=#[3];
 list ergebnis;
 if (switch==-1) {
   "An error has occurred in develop, so there is no HNE.";
   return(ergebnis);
 }
 intvec beta,s,svorl,ordnung,multseq,mpuiseux,npuiseux,halbgr;
 int genus,zeile,i,j,k,summe,conductor,ggT;
 string Ausgabe;
 int nc=ncols(m); int nr=nrows(m);
 ordnung[nr]=1;
         // alle Indizes muessen (gegenueber [Ca]) um 1 erhoeht werden,
         // weil 0..r nicht als Wertebereich erlaubt ist (aber nrows(m)==r+1)

 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
 for (zeile=nr; zeile>1; zeile--) {
   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
      k=1;
      while (m[zeile,k]==0) {k++;}
      ordnung[zeile-1]=k*ordnung[zeile]; // vgl. auch Def. von untergrad in
      genus++;                           // proc param
      svorl[genus]=zeile;} // werden gerade in umgekehrter Reihenfolge abgelegt
   else {
      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
 }}
 //----------------- charakteristische Exponenten (beta) ----------------------
 s[1]=1;
 for (k=1; k <= genus; k++) { s[k+1]=svorl[genus-k+1];} // s[2]==s(1), u.s.w.
 beta[1]=ordnung[1]; //charakt. Exponenten: Index wieder verschoben
 for (k=1; k <= genus; k++) {
   summe=0;
   for (i=1; i <= s[k]; i++) {summe=summe+v[i]*ordnung[i];}
   beta[k+1]=summe+ordnung[s[k]]+ordnung[s[k]+1]-ordnung[1];
 }
 //--------------------------- Puiseuxpaare -----------------------------------
 int produkt=1;
 for (i=1; i<=genus; i++) {
   ggT=gcd(beta[1],beta[i+1]*produkt);
   mpuiseux[i]=beta[i+1]*produkt / ggT;
   npuiseux[i]=beta[1] / ggT;
   produkt=produkt*npuiseux[i];
 }
 //---------------------- Grad des Konduktors ---------------------------------
 summe=1-ordnung[1];
 if (genus > 0) {
   for (i=2; i <= genus+1; i++) {
     summe=summe + beta[i] * (ordnung[s[i-1]] - ordnung[s[i]]);
   }                              // n.b.: Indizierung wieder um 1 verschoben
 }
 conductor=summe;
 //------------------- Erzeuger der Halbgruppe: -------------------------------
 halbgr=puiseux2generators(mpuiseux,npuiseux);

 //------------------- Multiplizitaetensequenz: -------------------------------
 k=1;
 for (i=1; i<size(v); i++) {
   for (j=1; j<=v[i]; j++) {
     multseq[k]=ordnung[i];
     k++;
 }}
 multseq[k]=1;
 //--- fuelle die Multipl.seq. mit den notwendigen Einsen auf -- T.Keilen ----
 int tester=k;
 while((multseq[tester]==1) and (tester>1))
 {
   tester=tester-1;
 } 
 if ((multseq[tester]!=1) and (multseq[tester]!=k-tester))
 {
   for (i=k+1; i<=tester+multseq[tester]; i++)
   {
     multseq[i]=1;
   }   
 }
 //--- Ende T.Keilen --- 06.05.02
 //------------------------- Rueckgabe ----------------------------------------
 ergebnis=beta,halbgr,mpuiseux,npuiseux,conductor,multseq;
 return(ergebnis);
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y),dp;
 list hne=develop(y4+2x3y2+x6+x5y);
 list INV=invariants(hne);
 INV[1];                   // the characteristic exponents
 INV[2];                   // the generators of the semigroup of values
 INV[3],INV[4];            // the Puiseux pairs in packed form
 INV[5] / 2;               // the delta-invariant
 INV[6];                   // the sequence of multiplicities
                           // To display the invariants more 'nicely':
 displayInvariants(hne);
 /////////////////////////////
 INV=invariants((x2-y3)*(x3-y5));
 INV[1][1];                // the characteristic exponents of the first branch
 INV[2][6];                // the sequence of multiplicities of the second branch
 print(INV[size(INV)][1]);         // the contact matrix of the branches
 print(INV[size(INV)][2]);         // the intersection numbers of the branches
 INV[size(INV)][3];                // the delta invariant of the curve
}

///////////////////////////////////////////////////////////////////////////////

proc displayInvariants
"USAGE:  displayInvariants(INPUT); INPUT list or poly
ASSUME:  INPUT is a bivariate polynomial, or the output of @code{develop(f)}, or of
         @code{extdevelop(develop(f),n)}, or (one entry of) the list @code{hne}
         in the ring created by @code{hnexpansion(f[,\"ess\"])}.
RETURN:  none
DISPLAY: invariants of the corresponding branch, resp. of all branches,
         in a better readable form.
NOTE:    In case the Hamburger-Noether expansion of the curve f is needed
         for other purposes as well it is better to calculate this first
         with the aid of @code{hnexpansion} and use it as input instead of
         the polynomial itself.
SEE ALSO: invariants, intersection, develop, hnexpansion
EXAMPLE: example displayInvariants;  shows an example
"
{
  // INPUT = poly or ring
  if (typeof(#[1])=="poly") 
  {
    list HNEXPANSION=hnexpansion(#[1]);
    displayInvariants(HNEXPANSION);
    return();
  }
  if (typeof(#[1])=="ring")
  {
    def H_N_E_RING=#[1];
    def rettering=basering;
    setring H_N_E_RING;
    displayInvariants(hne);
    setring rettering;
    kill H_N_E_RING;
    return();
  }
  // INPUT = hne of a plane curve
 int i,j,k,mul;
 string Ausgabe;
 list ergebnis;
 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
 #=stripHNE(#);
 //-------------------- Ausgabe eines Zweiges ---------------------------------
 if (typeof(#[1])=="matrix") {
   ergebnis=invariants(#);
   if (size(ergebnis)!=0) {
    " characteristic exponents  :",ergebnis[1];
    " generators of semigroup   :",ergebnis[2];
    if (size(ergebnis[1])>1) {
     for (i=1; i<=size(ergebnis[3]); i++) {
       Ausgabe=Ausgabe+"("+string(ergebnis[3][i])+","
       +string(ergebnis[4][i])+")";
    }}
    " Puiseux pairs             :",Ausgabe;
    " degree of the conductor   :",ergebnis[5];
    " delta invariant           :",ergebnis[5]/2;
    " sequence of multiplicities:",ergebnis[6];
 }}
 //-------------------- Ausgabe aller Zweige ----------------------------------
 else {
  ergebnis=invariants(#);
  intmat contact=ergebnis[size(#)+1][1];
  intmat intersectionmatrix=ergebnis[size(#)+1][2];
  for (j=1; j<=size(#); j++) {
    " --- invariants of branch number",j,": ---";
    " characteristic exponents  :",ergebnis[j][1];
    " generators of semigroup   :",ergebnis[j][2];
    Ausgabe="";
    if (size(ergebnis[j][1])>1) {
     for (i=1; i<=size(ergebnis[j][3]); i++) {
       Ausgabe=Ausgabe+"("+string(ergebnis[j][3][i])+","
       +string(ergebnis[j][4][i])+")";
    }}
    " Puiseux pairs             :",Ausgabe;
    " degree of the conductor   :",ergebnis[j][5];
    " delta invariant           :",ergebnis[j][5]/2;
    " sequence of multiplicities:",ergebnis[j][6];
    "";
  }
  if (size(#)>1) 
  {
    " -------------- contact numbers : -------------- ";"";
    Ausgabe="branch |   ";
    for (j=size(#); j>1; j--) 
    {
      if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+"    "; }
      else                    { Ausgabe=Ausgabe+string(j)+"    "; }
    }   
    Ausgabe;
    Ausgabe="-------+";
    for (j=2; j<size(#); j++) { Ausgabe=Ausgabe+"------"; }
    Ausgabe=Ausgabe+"-----";
    Ausgabe;
  } 
  for (j=1; j<size(#); j++) 
  {
    if (size(string(j))==1) { Ausgabe="    "+string(j)+"  |"; }
    else                    { Ausgabe="   " +string(j)+"  |"; }
    for (k=size(#); k>j; k--) 
    {
      mul=contact[j,k];//separateHNE(#[j],#[k]);
      for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; }
      Ausgabe=Ausgabe+string(mul);
      if (k>j+1) { Ausgabe=Ausgabe+","; }
    }
    Ausgabe;
  }
  "";
  if (size(#)>1) 
  {
    " -------------- intersection multiplicities : -------------- ";"";
    Ausgabe="branch |   ";
    for (j=size(#); j>1; j--) 
    {
      if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+"    "; }
      else                    { Ausgabe=Ausgabe+string(j)+"    "; }
    }   
    Ausgabe;
    Ausgabe="-------+";
    for (j=2; j<size(#); j++) { Ausgabe=Ausgabe+"------"; }
    Ausgabe=Ausgabe+"-----";
    Ausgabe;
  } 
  for (j=1; j<size(#); j++) 
  {
    if (size(string(j))==1) { Ausgabe="    "+string(j)+"  |"; }
    else                    { Ausgabe="   " +string(j)+"  |"; }
    for (k=size(#); k>j; k--) 
    {
      mul=intersectionmatrix[j,k];//intersection(#[j],#[k]);
      for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; }
      Ausgabe=Ausgabe+string(mul);
      if (k>j+1) { Ausgabe=Ausgabe+","; }
    }
    Ausgabe;
  }
  "";
  " -------------- delta invariant of the curve : ",ergebnis[size(#)+1][3];
  
 }
 return();
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y),dp;
 list hne=develop(y4+2x3y2+x6+x5y);
 displayInvariants(hne);
}
///////////////////////////////////////////////////////////////////////////////

proc multiplicities
"USAGE:   multiplicities(L); L list
ASSUME:  L is the output of @code{develop(f)}, or of
         @code{extdevelop(develop(f),n)}, or one entry in the list @code{hne}
         in the ring created by @code{hnexpansion(f[,\"ess\"])}.
RETURN:  intvec of the different multiplicities that occur when successively
         blowing-up the curve singularity corresponding to f.
SEE ALSO: multsequence, develop
EXAMPLE: example multiplicities;  shows an example
"
{
 matrix m=#[1];
 intvec v=#[2];
 int switch=#[3];
 list ergebnis;
 if (switch==-1) {
   "An error has occurred in develop, so there is no HNE.";
   return(intvec(0));
 }
 intvec ordnung;
 int zeile,k;
 int nc=ncols(m); int nr=nrows(m);
 ordnung[nr]=1;
 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
 for (zeile=nr; zeile>1; zeile--) {
   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
      k=1;
      while (m[zeile,k]==0) {k++;}
      ordnung[zeile-1]=k*ordnung[zeile];
   }
   else {
      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
 }}
 return(ordnung);
}
example
{ "EXAMPLE:"; echo = 2;
  int p=printlevel; printlevel=-1;
  ring r=0,(x,y),dp;
  multiplicities(develop(x5+y7));
  // The first value is the multiplicity of the curve itself, here it's 5
  printlevel=p;
}
///////////////////////////////////////////////////////////////////////////////

proc puiseux2generators (intvec m, intvec n)
"USAGE:   puiseux2generators(m,n); m,n intvec
ASSUME:  m, resp. n, represent the 1st, resp. 2nd, components of Puiseux pairs
         (e.g., @code{m=invariants(L)[3]}, @code{n=invariants(L)[4]}).
RETURN:  intvec of the generators of the semigroup of values.
SEE ALSO: invariants
EXAMPLE: example puiseux2generators;  shows an example
"
{
 intvec beta;
 int q=1;
 //------------ glatte Kurve (eigentl. waeren m,n leer): ----------------------
 if (m==0) { return(intvec(1)); }
 //------------------- singulaere Kurve: --------------------------------------
 for (int i=1; i<=size(n); i++) { q=q*n[i]; }
 beta[1]=q; // == q_0
 m=1,m; n=1,n; // m[1] ist damit m_0 usw., genau wie beta[1]==beta_0
 for (i=2; i<=size(n); i++) {
  beta[i]=m[i]*q/n[i] - m[i-1]*q + n[i-1]*beta[i-1];
  q=q/n[i]; // == q_i
 }
 return(beta);
}
example
{ "EXAMPLE:"; echo = 2;
  // take (3,2),(7,2),(15,2),(31,2),(63,2),(127,2) as Puiseux pairs:
  puiseux2generators(intvec(3,7,15,31,63,127),intvec(2,2,2,2,2,2));
}
///////////////////////////////////////////////////////////////////////////////

proc intersection (list hn1, list hn2)
"USAGE:   intersection(hne1,hne2); hne1, hne2 lists
ASSUME:  hne1, hne2 represent a HNE (i.e., are the output of
        @code{develop(f)}, or of @code{extdevelop(develop(f),n)}, or
        one entry of the list @code{hne} in the ring created by
        @code{hnexpansion(f[,\"ess\"])}).
RETURN:  int, the intersection multiplicity of the branches corresponding to
         hne1 and hne2.
SEE ALSO: hnexpansion, displayInvariants
KEYWORDS: intersection multiplicity
EXAMPLE: example intersection;  shows an example
"
{
 //------------------ `intersect' ist schon reserviert ... --------------------
 int i,j,s,sum,schnitt,unterschied;
 matrix a1=hn1[1];
 matrix a2=hn2[1];
 intvec h1=hn1[2];
 intvec h2=hn2[2];
 intvec n1=multiplicities(hn1);
 intvec n2=multiplicities(hn2);
 if (hn1[3]!=hn2[3]) {
 //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern -
 //---------------- d.h. beide Kurven schneiden sich transversal --------------
   schnitt=n1[1]*n2[1];        // = mult(hn1)*mult(hn2)
 }
 else {
 //--------- die jeweils erste Zeile gehoert zum gleichen Parameter -----------
   unterschied=0;
   for (s=1; (h1[s]==h2[s]) && (s<size(h1)) && (s<size(h2))
              && (unterschied==0); s++) {
     for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]); i++) {;}
     if (i<=h1[s]) {
       unterschied=1;
       s--;                    // um s++ am Schleifenende wieder auszugleichen
     }
   }
   if (unterschied==0) {
     if ((s<size(h1)) && (s<size(h2))) {
       for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]) && (i<=h2[s]); i++) {;}
     }
     else {
 //-------------- Sonderfall: Unterschied in letzter Zeile suchen -------------
 // Beachte: Es koennen undefinierte Stellen auftreten, bei abbrechender HNE
 // muss die Ende-Markierung weg, h_[r] ist unendlich, die Matrix muss mit
 // Nullen fortgesetzt gedacht werden
 //----------------------------------------------------------------------------
       if (ncols(a1)>ncols(a2)) { j=ncols(a1); }
       else                     { j=ncols(a2); }
       unterschied=0;
       if ((h1[s]>0) && (s==size(h1))) {
         a1[s,h1[s]+1]=0;
         if (ncols(a1)<=ncols(a2)) { unterschied=1; }
       }
       if ((h2[s]>0) && (s==size(h2))) {
         a2[s,h2[s]+1]=0;
         if (ncols(a2)<=ncols(a1)) { unterschied=1; }
       }
       if (unterschied==1) {                   // mind. eine HNE war endlich
         matrix ma1[1][j]=a1[s,1..ncols(a1)];  // und bedarf der Fortsetzung
         matrix ma2[1][j]=a2[s,1..ncols(a2)];  // mit Nullen
       }
       else {
         if (ncols(a1)>ncols(a2)) { j=ncols(a2); }
         else                     { j=ncols(a1); }
         matrix ma1[1][j]=a1[s,1..j];          // Beschr. auf vergleichbaren
         matrix ma2[1][j]=a2[s,1..j];          // Teil (der evtl. y's enth.)
       }
       for (i=1; (ma1[1,i]==ma2[1,i]) && (i<j) && (ma1[1,i]!=var(2)); i++) {;}
       if (ma1[1,i]==ma2[1,i]) {
         "//** The two HNE's are identical!";
         "//** You have either tried to intersect a branch with itself,";
         "//** or the two branches have been developed separately.";
         "//   In the latter case use `extdevelop' to extend the HNE's until",
         "they differ.";
         return(-1);
       }
       if ((ma1[1,i]==var(2)) || (ma2[1,i]==var(2))) {
         "//** The two HNE's are (so far) identical. This is because they",
         "have been";
         "//** computed separately. I need more data; use `extdevelop' to",
         "extend them,";
         if (ma1[1,i]==var(2)) {"//** at least the first one.";}
         else                  {"//** at least the second one.";}
         return(-1);
       }
     }
   }
   sum=0;
   h1[size(h1)]=ncols(a1)+42;        // Ersatz fuer h1[r]=infinity
   h2[size(h2)]=ncols(a2)+42;
   for (j=1; j<s; j++) {sum=sum+h1[j]*n1[j]*n2[j];}
   if ((i<=h1[s]) && (i<=h2[s]))    { schnitt=sum+i*n1[s]*n2[s]; }
   if (i==h2[s]+1) { schnitt=sum+h2[s]*n1[s]*n2[s]+n2[s+1]*n1[s]; }
   if (i==h1[s]+1) { schnitt=sum+h1[s]*n2[s]*n1[s]+n1[s+1]*n2[s]; }
 // "s:",s-1,"i:",i,"S:",sum;
 }
 return(schnitt);
}
example
{
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  list hn=hnexpansion((x2-y3)*(x2+y3));
  def HNEring=hn[1];
  setring HNEring;
  intersection(hne[1],hne[2]);
}
///////////////////////////////////////////////////////////////////////////////

proc multsequence
"USAGE:   multsequence(INPUT); INPUT list or poly 
ASSUME:  INPUT is the output of @code{develop(f)}, or of @code{extdevelop(develop(f),n)}, 
         or one entry in the list @code{hne} of the ring created by @code{hnexpansion(f)}.
RETURN:  intvec corresponding to the multiplicity sequence of (a branch)
         of the curve (the same as @code{invariants(INPUT)[6]}).

ASSUME:  INPUT is a bivariate polynomial, or the output of @code{hnexpansion(f)},
         or the list @code{hne} in the ring created by @code{hnexpansion(f)}.
RETURN:  list of two integer matrices:
@texinfo
@table @asis
@item  @code{multsequence(INPUT)[1][i,*]}
   contains the multiplicities of the branches at their infinitely near point
   of 0 in its (i-1) order neighbourhood (i.e., i=1: multiplicity of the
   branches themselves, i=2: multiplicity of their 1st quadratic transformed,
   etc., @*
   Hence, @code{multsequence(INPUT)[1][*,j]} is the multiplicity sequence
   of branch j.
@item  @code{multsequence(INPUT)[2][i,*]}:
   contains the information which of these infinitely near points coincide.
@end table
@end texinfo
NOTE:  The order of elements of the list @code{hne} obtained from @code{hnexpansion(f[,\"ess\")}
       must not be changed (because then the coincident infinitely near points
       couldn't be grouped together, cf. meaning of 2nd intmat in example).
       Hence, it is not wise to compute the HNE of several polynomials
       separately, put them into a list INPUT and call @code{multsequence(INPUT)}. @*
       Use @code{displayMultsequence} to produce a better readable output for
       reducible curves on the screen. @*
       In case the Hamburger-Noether expansion of the curve f is needed
       for other purposes as well it is better to calculate this first
       with the aid of @code{hnexpansion} and use it as input instead of
       the polynomial itself.
SEE ALSO: displayMultsequence, develop, hnexpansion, separateHNE
KEYWORDS: multiplicity sequence
EXAMPLE: example multsequence;  shows an example
"
{
 //---- INPUT = poly, or HNEring --------------------
  if (typeof(#[1])=="poly")
  {
    list HNEXPANSION=hnexpansion(#[1]);
    return(multsequence(HNEXPANSION));
  }
  if (typeof(#[1])=="ring")
  {
    def H_N_ER_I_N_G=#[1];
    def ret_te_ring=basering;
    setring H_N_ER_I_N_G;
    list ErGeBnIs=multsequence(hne);
    setring ret_te_ring;
    kill H_N_ER_I_N_G;
    return(ErGeBnIs);
  }
 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
 #=stripHNE(#);
 int k,i,j;
 //----------------- Multiplizitaetensequenz eines Zweiges --------------------
 if (typeof(#[1])=="matrix") {
  intvec v=#[2];
  list ergebnis;
  if (#[3]==-1) {
    "An error has occurred in develop, so there is no HNE.";
   return(intvec(0));
  }
  intvec multips,multseq;
  multips=multiplicities(#);
  k=1;
  for (i=1; i<size(v); i++) {
    for (j=1; j<=v[i]; j++) {
      multseq[k]=multips[i];
      k++;
  }}
  multseq[k]=1;
  //--- fuelle die Multipl.seq. mit den notwendigen Einsen auf -- T.Keilen ----
  int tester=k;
  while((multseq[tester]==1) and (tester>1))
  {
    tester=tester-1;
  }
  if((multseq[tester]!=1) and (multseq[tester]!=k-tester))
  {
    for (i=k+1; i<=tester+multseq[tester]; i++)
    {
      multseq[i]=1;
    }
  }
  //--- Ende T.Keilen --- 06.05.02  
  return(multseq);
 }
 //---------------------------- mehrere Zweige --------------------------------
 else {
   list HNEs=#;
   int anzahl=size(HNEs);
   int maxlength=0;
   int bisher;
   intvec schnitt,ones;
   ones[anzahl]=0;
   ones=ones+1;                  // = 1,1,...,1
   for (i=1; i<anzahl; i++) {
     schnitt[i]=separateHNE(HNEs[i],HNEs[i+1]);
     j=size(multsequence(HNEs[i]));
     maxlength=maxlength*(j < maxlength) + j*(j >= maxlength);
     maxlength=maxlength*(schnitt[i]+1 < maxlength)
               + (schnitt[i]+1)*(schnitt[i]+1 >= maxlength);
   }
   j=size(multsequence(HNEs[anzahl]));
   maxlength=maxlength*(j < maxlength) + j*(j >= maxlength);

//-------------- Konstruktion der ersten zu berechnenden Matrix ---------------
   intmat allmults[maxlength][anzahl];
   for (i=1; i<=maxlength; i++)  { allmults[i,1..anzahl]=ones[1..anzahl]; }
   for (i=1; i<=anzahl; i++) {
     ones=multsequence(HNEs[i]);
     allmults[1..size(ones),i]=ones[1..size(ones)];
   }
//---------------------- Konstruktion der zweiten Matrix ----------------------
   intmat separate[maxlength][anzahl];
   for (i=1; i<=maxlength; i++) {
     k=1;
     bisher=0;
     if (anzahl==1) { separate[i,1]=1; }
     for (j=1; j<anzahl; j++)   {
       if (schnitt[j]<i) {
         separate[i,k]=j-bisher;
         bisher=j;
         k++;
       }
       separate[i,k]=anzahl-bisher;
     }
   }
  return(list(allmults,separate));
 }
}
example
{
  // -------- prepare for example ---------
  if (nameof(basering)=="HNEring") {
   def rettering=HNEring;
   kill HNEring;
  }
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  list hn=hnexpansion((x6-y10)*(x+y2-y3)*(x+y2+y3));   // 4 branches
  def HNEring=hn[1];
  setring HNEring;
  multsequence(hne[1]),"  |  ",multsequence(hne[2]),"  |  ",
  multsequence(hne[3]),"  |  ",multsequence(hne[4]);
  multsequence(hne);
  // The meaning of the entries of the 2nd matrix is as follows:
  displayMultsequence(hne);
  echo = 0;
  // --- restore HNEring if previously defined ---
  kill HNEring,r;
  if (defined(rettering)) {
   setring rettering;
   def HNEring=rettering;
   export HNEring;
  }
}
///////////////////////////////////////////////////////////////////////////////

proc displayMultsequence
"USAGE:   displayMultsequence(INPUT); INPUT list or poly
ASSUME:  INPUT is a bivariate polynomial, or the output of @code{develop(f)}, 
         or of @code{extdevelop(develop(f),n)}, or of of @code{hnexpansion(f[,\"ess\"])}, 
         or (one entry in) the list @code{hne} of the ring created by @code{hnexpansion(f[,\"ess \"])}.
RETURN:  nothing
DISPLAY: the sequence of multiplicities:
@format
 - if @code{INPUT=develop(f)} or @code{INPUT=extdevelop(develop(f),n)} or @code{INPUT=hne[i]}:
                      @code{a , b , c , ....... , 1}
 - if @code{INPUT=f} or @code{INPUT=hnexpansion(f[,\"ess\"])} or @code{INPUT=hne}:
                      @code{[(a_1, .... , b_1 , .... , c_1)],}
                      @code{[(a_2, ... ), ... , (... , c_2)],}
                      @code{ ........................................ ,}
                      @code{[(a_n),(b_n), ....., (c_n)]}
     with:
       @code{a_1 , ... , a_n} the sequence of multiplicities of the 1st branch,
       @code{[...]} the multiplicities of the j-th transformed of all branches,
       @code{(...)} indicating branches meeting in an infinitely near point.
@end format
NOTE:  The same restrictions for INPUT as in @code{multsequence} apply.@*
       In case the Hamburger-Noether expansion of the curve f is needed
       for other purposes as well it is better to calculate this first
       with the aid of @code{hnexpansion} and use it as input instead of
       the polynomial itself.
SEE ALSO: multsequence, develop, hnexpansion, separateHNE
EXAMPLE: example displayMultsequence;  shows an example
"
{
 //---- INPUT = poly, or HNEring --------------------
  if (typeof(#[1])=="poly")
  {
    list HNEXPANSION=hnexpansion(#[1]);
    displayMultsequence(HNEXPANSION);
    return();
  }
  if (typeof(#[1])=="ring")
  {
    def H_N_ER_I_N_G=#[1];
    def ret_te_ring=basering;
    setring H_N_ER_I_N_G;
    displayMultsequence(hne);
    setring ret_te_ring;
    kill H_N_ER_I_N_G;
    return();
  }
 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
 #=stripHNE(#);
 //----------------- Multiplizitaetensequenz eines Zweiges --------------------
 if (typeof(#[1])=="matrix") {
   if (#[3]==-1) {
     "An error has occurred in develop, so there is no HNE.";
   }
   else {
     "The sequence of multiplicities is  ",multsequence(#);
 }}
 //---------------------------- mehrere Zweige --------------------------------
 else {
   list multips=multsequence(#);
   int i,j,k,l;
   string output;
   for (i=1; i<=nrows(multips[1]); i++) {
     output="[";
     k=1;
     for (l=1; k<=ncols(multips[1]); l++) {
       output=output+"(";
       for (j=1; j<=multips[2][i,l]; j++) {
         output=output+string(multips[1][i,k]);
         k++;
         if (j<multips[2][i,l]) { output=output+","; }
       }
       output=output+")";
       if ((k-1) < ncols(multips[1])) { output=output+","; }
      }
     output=output+"]";
     if (i<nrows(multips[1])) { output=output+","; }
     output;
   }
 }
} // example multsequence; geht wegen echo nicht (muesste auf 3 gesetzt werden)
example
{
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  //// Example 1: Input = output of develop
  displayMultsequence(develop(x3-y5));
  //// Example 2: Input = bivariate polynomial
  displayMultsequence((x6-y10)*(x+y2-y3)*(x+y2+y3));
}
///////////////////////////////////////////////////////////////////////////////

proc separateHNE (list hn1,list hn2)
"USAGE:    separateHNE(hne1,hne2);  hne1, hne2 lists
ASSUME:   hne1, hne2 are HNEs (=output of
          @code{develop(f)}, @code{extdevelop(develop(f),n)}, or
          one entry in the list @code{hne} in the ring created by
          @code{hnexpansion(f[,\"ess\"])}.
RETURN:   number of quadratic transformations needed to separate both curves
          (branches).
SEE ALSO: develop, hnexpansion, multsequence, displayMultsequence
EXAMPLE:  example separateHNE;  shows an example
"
{
 int i,j,s,unterschied,separated;
 matrix a1=hn1[1];
 matrix a2=hn2[1];
 intvec h1=hn1[2];
 intvec h2=hn2[2];
 if (hn1[3]!=hn2[3]) {
 //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern -
 //---------------- d.h. beide Kurven schneiden sich transversal --------------
   separated=1;
 }
 else {
 //--------- die jeweils erste Zeile gehoert zum gleichen Parameter -----------
   unterschied=0;
   for (s=1; (h1[s]==h2[s]) && (s<size(h1)) && (s<size(h2))
              && (unterschied==0); s++) {
     for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]); i++) {;}
     if (i<=h1[s]) {
       unterschied=1;
       s--;                    // um s++ am Schleifenende wieder auszugleichen
     }
   }
   if (unterschied==0) {
     if ((s<size(h1)) && (s<size(h2))) {
       for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]) && (i<=h2[s]); i++) {;}
     }
     else {
 //-------------- Sonderfall: Unterschied in letzter Zeile suchen -------------
 // Beachte: Es koennen undefinierte Stellen auftreten, bei abbrechender HNE
 // muss die Ende-Markierung weg, h_[r] ist unendlich, die Matrix muss mit
 // Nullen fortgesetzt gedacht werden
 //----------------------------------------------------------------------------
       if (ncols(a1)>ncols(a2)) { j=ncols(a1); }
       else                     { j=ncols(a2); }
       unterschied=0;
       if ((h1[s]>0) && (s==size(h1))) {
         a1[s,h1[s]+1]=0;
         if (ncols(a1)<=ncols(a2)) { unterschied=1; }
       }
       if ((h2[s]>0) && (s==size(h2))) {
         a2[s,h2[s]+1]=0;
         if (ncols(a2)<=ncols(a1)) { unterschied=1; }
       }
       if (unterschied==1) {                   // mind. eine HNE war endlich
         matrix ma1[1][j]=a1[s,1..ncols(a1)];  // und bedarf der Fortsetzung
         matrix ma2[1][j]=a2[s,1..ncols(a2)];  // mit Nullen
       }
       else {
         if (ncols(a1)>ncols(a2)) { j=ncols(a2); }
         else                     { j=ncols(a1); }
         matrix ma1[1][j]=a1[s,1..j];          // Beschr. auf vergleichbaren
         matrix ma2[1][j]=a2[s,1..j];          // Teil (der evtl. y's enth.)
       }
       for (i=1; (ma1[1,i]==ma2[1,i]) && (i<j) && (ma1[1,i]!=var(2)); i++) {;}
       if (ma1[1,i]==ma2[1,i]) {
         "//** The two HNE's are identical!";
         "//** You have either tried to compare a branch with itself,";
         "//** or the two branches have been developed separately.";
         "//   In the latter case use `extdevelop' to extend the HNE's until",
         "they differ.";
         return(-1);
       }
       if ((ma1[1,i]==var(2)) || (ma2[1,i]==var(2))) {
         "//** The two HNE's are (so far) identical. This is because they",
         "have been";
         "//** computed separately. I need more data; use `extdevelop' to",
         "extend them,";
         if (ma1[1,i]==var(2)) {"//** at least the first one.";}
         else                  {"//** at least the second one.";}
         return(-1);
       }
     }
   }
   separated=i;
   for (j=1; j<s; j++) { separated=separated+h1[j]; }
 }
 return(separated);
}
example
{ "EXAMPLE:"; echo = 2;
  int p=printlevel; printlevel=-1;
  ring r=0,(x,y),dp;
  list hne1=develop(x);
  list hne2=develop(x+y);
  list hne3=develop(x+y2);
  separateHNE(hne1,hne2);  // two transversal lines
  separateHNE(hne1,hne3);  // one quadratic transform. gives 1st example
  printlevel=p;
}
///////////////////////////////////////////////////////////////////////////////

proc displayHNE(list ldev,list #)
"USAGE:   displayHNE(L[,n]); L list, n int
ASSUME:  L is the output of @code{develop(f)}, or of @code{exdevelop(f,n)},
         or of @code{hnexpansion(f[,\"ess\"])}, or (one entry in) the list 
         @code{hne} in the ring created by @code{hnexpansion(f[,\"ess\"])}.
RETURN:  - if only one argument is given, no return value, but
         display an ideal HNE of the following form:
     @example
     HNE[1]=-y+[]*z(0)^1+[]*z(0)^2+...+z(0)^<>*z(1)
     HNE[2]=-x+          []*z(1)^2+...+z(1)^<>*z(2)
     HNE[3]=             []*z(2)^2+...+z(2)^<>*z(3)
     .......             ..........................
     HNE[r+1]=           []*z(r)^2+[]*z(r)^3+......
     @end example
     where @code{x},@code{y} are the first 2 variables of the basering.
     The values of @code{[]} are the coefficients of the Hamburger-Noether
     matrix, the values of @code{<>} are represented by @code{x} in the
     HN-matrix.@*
     - if a second argument is given, create and export a new ring with
     name @code{displayring} containing an ideal @code{HNE} as described
     above.@*
     - if L corresponds to the output of @code{hnexpansion(f[,\"ess\"])}
     or to the list @code{hne} in the ring created by @code{hnexpansion(f[,\"ess\"])},
     @code{displayHNE(L[,n])} shows the HNE's of all branches of f in the form
     described above. The optional parameter is then ignored.
NOTE:  The 1st line of the above ideal (i.e., @code{HNE[1]}) means that
     @code{y=[]*z(0)^1+...}, the 2nd line (@code{HNE[2]}) means that
     @code{x=[]*z(1)^2+...}, so you can see which indeterminate
     corresponds to which line (it's also possible that @code{x} corresponds
     to the 1st line and @code{y} to the 2nd).

SEE ALSO: develop, hnexpansion
EXAMPLE: example displayHNE; shows an example
"
{
 if ((typeof(ldev[1])=="list") || (typeof(ldev[1])=="none")) {
   for (int i=1; i<=size(ldev); i++) {
     "// Hamburger-Noether development of branch nr."+string(i)+":";
     displayHNE(ldev[i]);"";
   }
   return();
 }
 //--------------------- Initialisierungen und Ringwechsel --------------------
 matrix m=ldev[1];
 intvec v=ldev[2];
 int switch=ldev[3];
 if (switch==-1) {
   "An error has occurred throughout the expansion, so there is no HNE.";
   return(ideal(0));
 }
 def altring=basering;
 /////////////////////////////////////////////////////////
 //  Change by T. Keilen 08.06.2002
 //  ring + ring does not work if one ring is an algebraic extension
/*
 if (parstr(basering)!="") {
   if (charstr(basering)!=string(char(basering))+","+parstr(basering)) {
     execute
      ("ring dazu=("+charstr(basering)+"),z(0.."+string(size(v)-1)+"),ls;");
   }
   else { ring dazu=char(altring),z(0..size(v)-1),ls; }
 }
 else   { ring dazu=char(altring),z(0..size(v)-1),ls; }
 def displayring=dazu+altring;
*/
 execute("ring displayring=("+charstr(basering)+"),(z(0.."+string(size(v)-1)+"),"+varstr(basering)+"),(ls("+string(size(v))+"),"+ordstr(basering)+");");
 // End change by T. Keilen
 //////////////////////////////////////////////////////////////
 setring displayring;
 if (size(#) != 0) {
    export displayring;
  }
 map holematrix=altring,0;        // mappt nur die Monome vom Grad Null
 matrix m=holematrix(m);
 //--------------------- Erzeuge Matrix n mit n[i,j]=z(j-1)^i -----------------
 int i;
 matrix n[ncols(m)][nrows(m)];
 for (int j=1; j<=nrows(m); j++) {
    for (i=1; i<=ncols(m); i++) { n[i,j]=z(j-1)^i; }
 }
 matrix displaymat=m*n;
 ideal HNE;
 for (i=1; i<nrows(m); i++) { HNE[i]=displaymat[i,i]+z(i)*z(i-1)^v[i]; }
 HNE[nrows(m)]=displaymat[nrows(m),nrows(m)];
 if (nrows(m)<2) { HNE[2]=z(0); }
 if (switch==0) {
    HNE[1] = HNE[1]-var(size(v)+2);
    HNE[2] = HNE[2]-var(size(v)+1);
 }
 else {
    HNE[1] = HNE[1]-var(size(v)+1);
    HNE[2] = HNE[2]-var(size(v)+2);
 }
if (size(#) == 0) {
   HNE;
   return();
 }
if (size(#) != 0) {
   "// basering is now 'displayring' containing ideal 'HNE'";
   keepring(displayring);
   export(HNE);
   return(HNE);
 }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  poly f=x3+2xy2+y2;
  list hn=develop(f);
  displayHNE(hn);
}
///////////////////////////////////////////////////////////////////////////////
//                      procedures for reducible curves                      //
///////////////////////////////////////////////////////////////////////////////

// proc newtonhoehne (poly f)
// USAGE:   newtonhoehne(f);   f poly
// ASSUME:  basering = ...,(x,y),ds  or ls
// RETURN:  list of intvec(x,y) of coordinates of the newtonpolygon of f
// NOTE:    This proc is only available in versions of Singular that know the
//         command  system("newton",f);  f poly
// {
// intvec nm = getnm(f);
//  if ((nm[1]>0) && (nm[2]>0)) { f=jet(f,nm[1]*nm[2],nm); }
//  list erg=system("newton",f);
//  int i; list Ausgabe;
//  for (i=1; i<=size(erg); i++) { Ausgabe[i]=leadexp(erg[i]); }
// return(Ausgabe);
// }
///////////////////////////////////////////////////////////////////////////////

proc newtonpoly (poly f, int #)
"USAGE:   newtonpoly(f);   f poly
ASSUME:  basering has exactly two variables; @*
         f is convenient, that is, f(x,0) != 0 != f(0,y).
RETURN:  list of intvecs (= coordinates x,y of the Newton polygon of f).
NOTE:    Procedure uses @code{execute}; this can be avoided by calling
         @code{newtonpoly(f,1)} if the ordering of the basering is @code{ls}. 
KEYWORDS: Newton polygon
EXAMPLE: example newtonpoly;  shows an example
"
{
  if (size(#)>=1)
  {
    if (typeof(#[1])=="int")
    {
      // this is done to avoid the "execute" command for procedures in
      //  hnoether.lib
      def is_ls=#[1];
    }
  }
  if (defined(is_ls)<=0)
  {
    def @Rold=basering;
    execute("ring @RR=("+charstr(basering)+"),("+varstr(basering)+"),ls;");
    poly f=imap(@Rold,f);
  }
  intvec A=(0,ord(subst(f,var(1),0)));
  intvec B=(ord(subst(f,var(2),0)),0);
  intvec C,H; list L;
  int abbruch,i;
  poly hilf;
  L[1]=A;
  f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
  if (defined(is_ls))
  {
    map xytausch=basering,var(2),var(1);
  }
  else
  {
    map xytausch=@RR,var(2),var(1); 
  }
  for (i=2; f!=0; i++) 
  {
     abbruch=0;
     while (abbruch==0) 
     {
        C=leadexp(f);          
        if(jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
        { 
           abbruch=1; 
        }       
        else 
        { 
           f=jet(f,-C[1]-1,intvec(-1,0)); 
        }
    }
    hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));
    H=leadexp(xytausch(hilf));
    A=H[2],H[1]; 
    L[i]=A;
    f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]-A[1]));
  }
  L[i]=B;
  if (defined(is_ls))
  {
    return(L);
  }
  else
  {
    setring @Rold;
    return(L);
  }
}
example
{
 "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),ls;
  poly f=x5+2x3y-x2y2+3xy5+y6-y7;
  newtonpoly(f);
}
///////////////////////////////////////////////////////////////////////////////

proc is_NND (poly f, list #)
"USAGE:   is_NND(f[,mu,NP]);   f poly, mu int, NP list of intvecs
ASSUME:  f is convenient, that is, f(x,0) != 0 != f(0,y);@*
         mu (optional) is Milnor number of f.@*
         NP (optional) is output of @code{newtonpoly(f)}. 
RETURN:  int: 1 if f in Newton non-degenerate, 0 otherwise.
SEE ALSO: newtonpoly
KEYWORDS: Newton non-degenerate; Newton polygon
EXAMPLE: example is_NND;  shows examples
"
{
  int i;
  int i_print=printlevel-voice+2;

  if (size(#)==0)
  {
    int mu=milnor(f);
    list NP=newtonpoly(f);
  }
  else
  {
    if (typeof(#[1])=="int")
    {
      def mu=#[1];
      def NP=#[2];
      for (i=1;i<=size(NP);i++)
      {
        if (typeof(NP[i])!="intvec")
        {
          print("third input cannot be Newton polygon ==> ignored ")
          NP=newtonpoly(f);
          i=size(NP)+1;
        }  
      }
    }
    else
    {
      print("second input cannot be Milnor number ==> ignored ")
      int mu=milnor(f);
      NP=newtonpoly(f);
    }
  }

  // computation of the Newton number:
  int s=size(NP);
  int nN=-NP[1][2]-NP[s][1]+1;
  intmat m[2][2];
  for(i=1;i<=s-1;i++)
  {
    m=NP[i+1],NP[i];
    nN=nN+det(m);
  }

  if(mu==nN)                   
  { // the Newton-polygon is non-degenerate
    return(1);
  }
  else
  {
    return(0);
  }
}
example
{
 "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),ls;
  poly f=x5+y3;
  is_NND(f);
  poly g=(x-y)^5+3xy5+y6-y7;
  is_NND(g);

  // if already computed, one should give the Minor number and Newton polygon
  // as second and third input: 
  int mu=milnor(g);
  list NP=newtonpoly(g);
  is_NND(g,mu,NP);
}


///////////////////////////////////////////////////////////////////////////////

proc charPoly(poly f, int M, int N)
"USAGE:  charPoly(f,M,N);  f bivariate poly,  M,N int: length and height
                          of Newton polygon of f, which has to be only one line
RETURN:  the characteristic polynomial of f
EXAMPLE: example charPoly;  shows an example
"
{
 poly charp;
 int Np=N/ gcd(M,N);
 f=subst(f,var(1),1);
 for(charp=0; f<>0; f=f-lead(f))
  { charp=charp+leadcoef(f)*var(2)^(leadexp(f)[2]/ Np);}
 return(charp);
}
example
{ "EXAMPLE:"; echo = 2;
  ring exring=0,(x,y),dp;
  charPoly(y4+2y3x2-yx6+x8,8,4);
  charPoly(y6+3y3x2-x4,4,6);
}
///////////////////////////////////////////////////////////////////////////////

proc find_in_list(list L,int p)
"USAGE:   find_in_list(L,p); L: list of intvec(x,y)
         (sorted in y: L[1][2]>=L[2][2]), int p >= 0
RETURN:  int i: L[i][2]=p if existent; otherwise i with L[i][2]<p if existent;
         otherwise i = size(L)+1;
EXAMPLE: example find_in_list;  shows an example
"
{
 int i;
 L[size(L)+1]=intvec(0,-1);          // falls p nicht in L[.][2] vorkommt
 for (i=1; L[i][2]>p; i++) {;}
 return(i);
}
example
{ "EXAMPLE:"; echo = 2;
  list L = intvec(0,4), intvec(1,2), intvec(2,1), intvec(4,0);
  find_in_list(L,1);
  L[find_in_list(L,2)];
}
///////////////////////////////////////////////////////////////////////////////

proc get_last_divisor(int M, int N)
"USAGE:   get_last_divisor(M,N); int M,N
RETURN:  int Q: M=q1*N+r1, N=q2*r1+r2, ..., ri=Q*r(i+1) (Euclidean alg.)
EXAMPLE: example get_last_divisor; shows an example
"
{
 int R=M%N; int Q=M / N;
 while (R!=0) {M=N; N=R; R=M%N; Q=M / N;}
 return(Q)
}
example
{ "EXAMPLE"; echo = 2;
  ring r=0,(x,y),dp;
  get_last_divisor(12,10);
}
///////////////////////////////////////////////////////////////////////////////
proc redleit (poly f,intvec S, intvec E)
"USAGE:   redleit(f,S,E);  f poly, S,E intvec(x,y)
         S,E are two different points on a line in the Newton diagram of f
RETURN:  poly g: all monomials of f which lie on or below that line
NOTE:    The main purpose is that if the line defined by S and E is part of the
         Newton polygon, the result is the quasihomogeneous leading form of f
         wrt. that line.
SEE ALSO: newtonpoly
EXAMPLE: example redleit;  shows an example
"
{
 if (E[1]<S[1]) { intvec H=E; E=S; S=H; } // S,E verkehrt herum eingegeben
 return(jet(f,E[1]*S[2]-E[2]*S[1],intvec(S[2]-E[2],E[1]-S[1])));
}
example
{ "EXAMPLE"; echo = 2;
  ring exring=0,(x,y),dp;
  redleit(y6+xy4-2x3y2+x4y+x6,intvec(3,2),intvec(4,1));
}
///////////////////////////////////////////////////////////////////////////////


proc extdevelop (list l, int Exaktheit)
"USAGE:   extdevelop(L,N); list L, int N
ASSUME:  L is the output of @code{develop(f)}, or of @code{extdevelop(l,n)},
         or one entry in the list @code{hne} in the ring created by
          @code{hnexpansion(f[,\"ess\"])}.
RETURN:  an extension of the Hamburger-Noether development of f as a list
         in the same format as L has (up to the last entry in the output
         of @code{develop(f)}).@*
         Type @code{help develop;}, resp. @code{help hnexpansion;} for more
         details.
NOTE:    The new HN-matrix will have at least N columns (if the HNE is not
         finite). In particular, if f is irreducible then (in most cases)
         @code{extdevelop(develop(f),N)} will produce the same result as
         @code{develop(f,N)}.@*
         If the matrix M of L has n columns then, compared with
         @code{parametrisation(L)}, @code{paramametrize(extdevelop(L,N))} will increase the
         exactness by at least (N-n) more significant monomials.
SEE ALSO: develop, hnexpansion, parametrisation
EXAMPLE: example extdevelop;  shows an example
"
{
 //------------ Initialisierungen und Abfangen unzulaessiger Aufrufe ----------
 matrix m=l[1];
 intvec v=l[2];
 int switch=l[3];
 if (nvars(basering) < 2) {
   " Sorry. I need two variables in the ring.";
   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0)));}
 if (switch==-1) {
   "An error has occurred in develop, so there is no HNE and no extension.";
   return(l);
 }
 poly f=l[4];
 if (f==0) {
   " No extension is possible";
   return(l);
 }
 int Q=v[size(v)];
 if (Q>0) {
   " The HNE was already exact";
   return(l);
 }
 else {
   if (Q==-1) { Q=ncols(m); }
   else { Q=-Q-1; }
 }
 int zeile=nrows(m);
 int spalten,i,M;
 ideal lastrow=m[zeile,1..Q];
 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");

 //------------------------- Ringwechsel, falls noetig ------------------------
 if (ringwechsel) {
  def altring = basering;
  int p = char(basering);
  if (charstr(basering)!=string(p)) {
     string tststr=charstr(basering);
     tststr=tststr[1..find(tststr,",")-1];     //-> "p^k" bzw. "p"
     if (tststr==string(p)) {
       if (size(parstr(basering))>1) {         // ring (p,a,..),...
        execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;");
       }
       else {                                  // ring (p,a),...
        string mipl=string(minpoly);
        ring extdguenstig=(p,`parstr(basering)`),(x,y),ls;
        if (mipl!="0") { execute("minpoly="+mipl+";"); }
       }
     }
     else {
       execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;");
     }
  }
  else {                               // charstr(basering)== p : no parameter
     ring extdguenstig=p,(x,y),ls;
  }
  export extdguenstig;
  map hole=altring,x,y;
 //----- map kann sehr zeitaufwendig sein, daher Vermeidung, wo moeglich: -----
  if (nvars(altring)==2) { poly f=fetch(altring,f); }
  else                   { poly f=hole(f);          }
  ideal a=hole(lastrow);
 }
 else { ideal a=lastrow; }
 list Newton=newtonpoly(f,1);
 int M1=Newton[size(Newton)-1][1];     // konstant
 number delt;
 if (Newton[size(Newton)-1][2]!=1) {
    " *** The transformed polynomial was not valid!!";}
 else {
 //--------------------- Fortsetzung der HNE ----------------------------------
  while (Q<Exaktheit) {
    M=ord(subst(f,y,0));
    Q=M-M1;
 //------ quasihomogene Leitform ist c*x^M1*y+d*x^(M1+Q) => delta=-d/c: -------
    delt=-koeff(f,M,0)/koeff(f,M1,1);
    a[Q]=delt;
    dbprint(printlevel-voice+2,"a("+string(zeile-1)+","+string(Q)+") = "+string(delt)); 
    if (Q<Exaktheit) {
     f=T1_Transform(f,delt,Q);
     if (defined(HNDebugOn)) { "transformed polynomial:",f; }
     if (subst(f,y,0)==0) {
       dbprint(printlevel-voice+2,"The HNE is finite!"); 
       a[Q+1]=x; Exaktheit=Q;
       f=0;                        // Speicherersparnis: f nicht mehr gebraucht
     }
    }
  }
 }
 //------- Wechsel in alten Ring, Zusammensetzung alte HNE + Erweiterung ------
 if (ringwechsel) {
  setring altring;
  map zurueck=extdguenstig,var(1),var(2);
  if (nvars(altring)==2) { f=fetch(extdguenstig,f); }
  else                   { f=zurueck(f);            }
  lastrow=zurueck(a);
 }
 else { lastrow=a; }
 if (ncols(lastrow)>ncols(m)) { spalten=ncols(lastrow); }
 else { spalten=ncols(m); }
 matrix mneu[zeile][spalten];
 for (i=1; i<nrows(m); i++) {
  mneu[i,1..ncols(m)]=m[i,1..ncols(m)];
 }
 mneu[zeile,1..ncols(lastrow)]=lastrow;
 if (lastrow[ncols(lastrow)]!=var(1)) {
  if (ncols(lastrow)==spalten) { v[zeile]=-1; }  // keine undefinierten Stellen
  else {
   v[zeile]=-Q-1;
   for (i=ncols(lastrow)+1; i<=spalten; i++) {
    mneu[zeile,i]=var(2);           // fuelle nicht def. Stellen der Matrix auf
 }}}
 else { v[zeile]=Q; }               // HNE war exakt
 if (ringwechsel)
 {
   if(system("with","Namespaces")) { kill Top::extdguenstig; }
   kill extdguenstig;
 }

 return(list(mneu,v,switch,f));
}
example
{
  if (defined(HNEring))
  {
    def save_r_i_n_g=HNEring;
    kill HNEring;
  }
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring exring=0,(x,y),dp;
  list hn=hnexpansion(x14-3y2x11-y3x10-y2x9+3y4x8+y5x7+3y4x6+x5*(-y6+y5)
                      -3y6x3-y7x2+y8);
  def HNEring=hn[1];
  setring HNEring;  echo=0;
  export(HNEring);  echo=2;
  print(hne[1][1]);    // HNE of 1st branch is finite
  print(extdevelop(hne[1],5)[1]);
  print(hne[2][1]);    // HNE of 2nd branch can be extended
  list ehne=extdevelop(hne[2],5);
  print(ehne[1]);      // new HN-matrix has 5 columns
  parametrisation(hne[2]);
  parametrisation(ehne);
  echo=0;
  if (defined(save_r_i_n_g))
  {
    kill HNEring;
    def HNEring=save_r_i_n_g;
  }
}
///////////////////////////////////////////////////////////////////////////////

proc stripHNE (list l)
"USAGE:   stripHNE(L);  L list
ASSUME:  L is the output of @code{develop(f)}, or of
         @code{extdevelop(develop(f),n)}, or (one entry of) the list
         @code{hne} in the ring created by @code{hnexpansion(f[,\"ess\"])}.
RETURN:  list in the same format as L, but all polynomials L[4], resp.
         L[i][4], are set to zero.
NOTE:    The purpose of this procedure is to remove huge amounts of data
         no longer needed. It is useful, if one or more of the polynomials
         in L consume much memory. It is still possible to compute invariants,
         parametrizations etc. with the stripped HNE(s), but it is not possible
         to use @code{extdevelop} with them.
SEE ALSO: develop, hnexpansion, extdevelop
EXAMPLE: example stripHNE;  shows an example
"
{
 list h;
 if (typeof(l[1])=="matrix") { l[4]=poly(0); }
 else {
  for (int i=1; i<=size(l); i++) {
    h=l[i];
    h[4]=poly(0);
    l[i]=h;
  }
 }
 return(l);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  list hne=develop(x2+y3+y4);
  hne;
  stripHNE(hne);
}
///////////////////////////////////////////////////////////////////////////////
static proc extractHNEs(list HNEs, int transvers)
"USAGE:  extractHNEs(HNEs,transvers);  list HNEs (output from HN),
        int transvers: 1 if x,y were exchanged, 0 else
RETURN: list of Hamburger-Noether-Extensions in the form of hne in hnexpansion
NOTE:   This procedure is only for internal purpose; examples don't make sense
"
{
 int i,maxspalte,hspalte,hnezaehler;
 list HNEaktu,Ergebnis;
 for (hnezaehler=1; hnezaehler<=size(HNEs); hnezaehler++) {
  maxspalte=0;
  HNEaktu=HNEs[hnezaehler];
  if (defined(HNDebugOn)) {"To process:";HNEaktu;}
  if (size(HNEaktu)!=size(HNEaktu[1])+2) {
     "The ideals and the hqs in HNEs[",hnezaehler,"] don't match!!";
     HNEs[hnezaehler];
  }
 //------------ ermittle  maximale Anzahl benoetigter Spalten: ----------------
  for (i=2; i<size(HNEaktu); i++) {
    hspalte=ncols(HNEaktu[i]);
    maxspalte=maxspalte*(hspalte < maxspalte)+hspalte*(hspalte >= maxspalte);
  }
 //------------- schreibe Ausgabe fuer hnezaehler-ten Zweig: ------------------
  matrix ma[size(HNEaktu)-2][maxspalte];
  for (i=1; i<=(size(HNEaktu)-2); i++) {
    if (ncols(HNEaktu[i+1]) > 1) {
      ma[i,1..ncols(HNEaktu[i+1])]=HNEaktu[i+1]; }
    else { ma[i,1]=HNEaktu[i+1][1];}
  }
  Ergebnis[hnezaehler]=list(ma,HNEaktu[1],transvers,HNEaktu[size(HNEaktu)]);
  kill ma;
 }
 return(Ergebnis);
}
///////////////////////////////////////////////////////////////////////////////

proc factorfirst(poly f, int M, int N)
"USAGE : factorfirst(f,M,N); f poly, M,N int
RETURN: number d: f=c*(y^(N/e) - d*x^(M/e))^e with e=gcd(M,N), number c fitting
        0 if d does not exist
EXAMPLE: example factorfirst;  shows an example
"
{
 number c = koeff(f,0,N);
 number delt;
 int eps,l;
 int p=char(basering);
 string ringchar=charstr(basering);

 if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
 int e = gcd(M,N);

 if (p==0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
 else {
   if (e%p != 0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
   else {
     eps = e;
     for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
     if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;}
     delt = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);

     if ((charstr(basering) != string(p)) and (delt != 0)) {
 //------ coefficient field is not Z/pZ => (p^l)th root is not identity -------
       delt=0;
       if (defined(HNDebugOn)) {
         "trivial factorization not implemented for",
         "parameters---I've to use 'factorize'";
       }
     }
   }
 }
 if (defined(HNDebugOn)) {"quasihomogeneous leading form:",f," = ",c,
        "* (y^"+string(N/ e),"-",delt,"* x^"+string(M/ e)+")^",e," ?";}
 if (f != c*(var(2)^(N/ e) - delt*var(1)^(M/ e))^e) {return(0);}
 else {return(delt);}
}
example
{ "EXAMPLE:"; echo = 2;
  ring exring=7,(x,y),dp;
  factorfirst(2*(y3-3x4)^5,20,15);
  factorfirst(x14+y7,14,7);
  factorfirst(x14+x8y3+y7,14,7);
}

///////////////////////////////////////////////////////////////////////////////
//         
//  the command HNdevelop is obsolete  --> here is the former help string:
//
///////////////////////////////////////////////////////////////////////////////
//
//ASSUME:  f is a bivariate polynomial (in the first 2 ring variables)
//CREATE:  ring with name @code{HNEring}, variables @code{x,y} and ordering
//         @code{ls} over a field extension of the current basering's ground
//         field. @*
//         Since the Hamburger-Noether development usually does not exist
//         in the originally given basering, @code{HNdevelop} always defines
//         @code{HNEring} and CHANGES to it. The field extension is chosen
//         minimally.
//RETURN:  list @code{L} of lists @code{L[i]} (corresponding to the output of
//         @code{develop(f[i])}, f[i] a branch of f, but the last entry being
//         omitted).
//@texinfo
//@table @asis
//@item @code{L[i][1]}; matrix:
//         Each row contains the coefficients of the corresponding line of the
//         Hamburger-Noether expansion (HNE) for f[i]. The end of the line is
//         marked in the matrix by the first ring variable (usually x).
//@item @code{L[i][2]}; intvec:
//         indicating the length of lines of the HNE
//@item @code{L[i][3]}; int:
//         0  if the 1st ring variable was transversal (with respect to f[i]), @*
//         1  if the variables were changed at the beginning of the
//            computation, @*
//        -1  if an error has occurred.
//@item @code{L[i][4]}; poly:
//         the transformed polynomial of f[i] to make it possible to extend the
//         Hamburger-Noether development a posteriori without having to do
//         all the previous calculation once again (0 if not needed)
//@end table
//@end texinfo
//NOTE:    @code{HNdevelop} decides which procedure (@code{develop} or
//         @code{reddevelop}) applies best to the given problem and calls it. @*
//         If f is known to be irreducible as a power series, @code{develop(f)}
//         should be chosen instead to avoid the change of basering. @*
//         If @code{printlevel>=2} comments are displayed (default is
//         @code{printlevel=0}).
//
//EXAMPLE: example HNdevelop;  shows an example
//
proc HNdevelop (poly f)
"USAGE:   HNdevelop(f); f poly
NOTE:     command is obsolete, use hnexpansion(f) instead.
SEE ALSO: hnexpansion, develop, extdevelop, param, displayHNE
"
{
 int irred=0;
 //--------- Falls Ring (p^k,a),...: Wechsel in (p,a),... + minpoly -----------
 if ((find(charstr(basering),string(char(basering)))!=1) &&
     (charstr(basering)<>"real")) {
   string strmip=string(minpoly);
   string strf=string(f);
   execute("ring tempr=("+string(char(basering))+","+parstr(basering)+"),("
           +varstr(basering)+"),dp;");
   execute("minpoly="+strmip+";");
   execute("poly f="+strf+";");
   list hne=reddevelop(f);
   if ((voice==2) && (printlevel > -1)) {
     "// Attention: The parameter",par(1),"has changed its meaning!";
     "// It need no longer be a generator of the cyclic group of unities!";
   }
 }
 else {
 //--- Falls Ring (0,a),... + minpoly : solange factorize nicht in Singular ---
 //------- implementiert ist, develop aufrufen (kann spaeter entfallen) -------
   if ((char(basering)==0) && (npars(basering)==1)) {
     if (string(minpoly)<>"0") { irred=1; }
   }
 //------------------ Aufruf der geeigneten Prozedur --------------------------
   if (irred==0) {
     list hne=pre_HN(f,0);       // = reddevelop(f);
     dbprint(printlevel-voice+2,
        "// result: "+string(size(hne))+" branch(es) successfully computed,",
        "//         basering has changed to HNEring");
   }
   else {
     def altring=basering;
     string strmip=string(minpoly);
     ring HNEring=(char(altring),`parstr(altring)`),(x,y),ls;
     execute("minpoly="+strmip+";");
     export HNEring;
     poly f=fetch(altring,f);
     list hn=develop(f,-1);
     list hne;
     if (hn[3] <> -1) {
       hne[1]=list(hn[1],hn[2],hn[3],hn[4]);
       if (hn[5] <> 1) {
   " ** WARNING : The curve is reducible, but only one branch could be found!";
       }
     }
     else { " ** Sorry -- could not find a HNE."; }
     dbprint(printlevel-voice+2,"// note: basering has changed to HNEring");
   }
 }
 keepring basering;
 return(hne);
}
example
{
  // -------- prepare for example ---------
  if (nameof(basering)=="HNEring") {
   def rettering=HNEring;
   kill HNEring;
  }
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  list hne=HNdevelop(x4-y6);
  nameof(basering);
  size(hne);           // number of branches
  print(hne[1][1]);    // HN-matrix of 1st branch
  param(hne[1]);       // parametrization of 1st branch
  param(hne[2]);       // parametrization of 2nd branch
  kill HNEring,r;
  echo = 0;
  // --- restore HNEring if previously defined ---
  if (defined(rettering)) {
   setring rettering;
   def HNEring=rettering;
   export HNEring;
  }
}

///////////////////////////////////////////////////////////////////////////////
//         
//  the command reddevelop is obsolete  --> here is the former help string:
//
///////////////////////////////////////////////////////////////////////////////
//ASSUME:  f is a bivariate polynomial (in the first 2 ring variables)
//CREATE:  ring with name @code{HNEring}, variables @code{x,y} and ordering
//         @code{ls} over a field extension of the current basering's ground
//         field. @*
//         Since the Hamburger-Noether development of a reducible curve
//         singularity usually does not exist in the originally given basering,
//         @code{reddevelop} always defines @code{HNEring} and CHANGES to it.
//         The field extension is chosen minimally.
//RETURN:  list @code{L} of lists @code{L[i]} (corresponding to the output of
//         @code{develop(f[i])}, f[i] a branch of f, but the last entry being
//         omitted).
//@texinfo
//@table @asis
//@item @code{L[i][1]}; matrix:
//         Each row contains the coefficients of the corresponding line of the
//         Hamburger-Noether expansion (HNE) for f[i]. The end of the line is
//         marked in the matrix by the first ring variable (usually x).
//@item @code{L[i][2]}; intvec:
//         indicating the length of lines of the HNE
//@item @code{L[i][3]}; int:
//         0  if the 1st ring variable was transversal (with respect to f[i]), @*
//         1  if the variables were changed at the beginning of the
//            computation, @*
//        -1  if an error has occurred.
//@item @code{L[i][4]}; poly:
//         the transformed polynomial of f[i] to make it possible to extend the
//         Hamburger-Noether development a posteriori without having to do
//         all the previous calculation once again (0 if not needed)
//@end table
//@end texinfo
//NOTE:    If @code{printlevel>=0} comments are displayed (default is
//         @code{printlevel=0}).
//
//EXAMPLE: example reddevelop;  shows an example
//
proc reddevelop (poly f)
"USAGE:   reddevelop(f); f poly
NOTE:     command is obsolete, use hnexpansion(f) instead.   
SEE ALSO: hnexpansion, develop, extdevelop, param, displayHNE
"
{
 list Ergebnis=pre_HN(f,0);
 if (size(Ergebnis)>0) {     // otherwise an error may have occurred
  dbprint(printlevel-voice+2,
   "// result: "+string(size(Ergebnis))+" branch(es) successfully computed,",
   "//         basering has changed to HNEring");
 }

 // ----- Lossen 10/02 : the branches have to be resorted to be able to
 // -----                display the multsequence in a nice way
 if (size(Ergebnis)>2)
 {  
   int i,j,k,m;
   list dummy;
   int nbsave;
   int no_br = size(Ergebnis);
   intmat nbhd[no_br][no_br];
   for (i=1;i<no_br;i++)
   {
     for (j=i+1;j<=no_br;j++)  
     { 
       nbhd[i,j]=separateHNE(Ergebnis[i],Ergebnis[j]);
       k=i+1;
       while ( (nbhd[i,k] >= nbhd[i,j]) and (k<j) )
       {
         k++;
       }
       if (k<j)  // branches have to be resorted 
       {
         dummy=Ergebnis[j];
         nbsave=nbhd[i,j];
         for (m=k; m<j; m++)
         {
           Ergebnis[m+1]=Ergebnis[m];
           nbhd[i,m+1]=nbhd[i,m];
         }
         Ergebnis[k]=dummy;
         nbhd[i,k]=nbsave;
       }
     }
   }
 }
 // -----
 
 keepring basering;
 return(Ergebnis);
}
example
{
  // -------- prepare for example ---------
  if (nameof(basering)=="HNEring")
  {
   def rettering=HNEring;
   kill HNEring;
  }
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring r = 32003,(x,y),dp;
  poly f = x25+x24-4x23-1x22y+4x22+8x21y-2x21-12x20y-4x19y2+4x20+10x19y
          +12x18y2-24x18y-20x17y2-4x16y3+x18+60x16y2+20x15y3-9x16y
          -80x14y3-10x13y4+36x14y2+60x12y4+2x11y5-84x12y3-24x10y5
          +126x10y4+4x8y6-126x8y5+84x6y6-36x4y7+9x2y8-1y9;
  list hne=reddevelop(f);
  size(hne);            // number of branches
  print(hne[1][1]);     // HN-matrix of 1st branch
  print(hne[4][1]);     // HN-matrix of 4th branch
  // a ring change was necessary, a is a parameter
  HNEring;
  kill HNEring,r;
  echo = 0;
  // --- restore HNEring if previously defined ---
  if (defined(rettering)) {
   setring rettering;
   def HNEring=rettering;
   export HNEring;
  }
}
///////////////////////////////////////////////////////////////////////////////

static proc pre_HN (poly f, int essential)
"NOTE: This procedure is only for internal use, it is called via
      reddevelop or essdevelop"
{
 def altring = basering;
 int p = char(basering);                 // Ringcharakteristik

 //-------------------- Tests auf Zulaessigkeit von basering ------------------
 if (charstr(basering)=="real") {
   " Singular cannot factorize over 'real' as ground field";
   return(list());
 }
 if (size(maxideal(1))<2) {
   " A univariate polynomial ring makes no sense !";
   return(list());
 }
 if (size(maxideal(1))>2) {
   dbprint(printlevel-voice+2,
   " Warning: all variables except the first two will be ignored!");
 }
 if (find(charstr(basering),string(char(basering)))!=1) {
   " ring of type (p^k,a) not implemented";
 //----------------------------------------------------------------------------
 // weder primitives Element noch factorize noch map "char p^k" -> "char -p"
 // [(p^k,a)->(p,a) ground field] noch fetch
 //----------------------------------------------------------------------------
   return(list());
 }
 //----------------- Definition eines neuen Ringes: HNEring -------------------
 string namex=varstr(1); string namey=varstr(2);
 if (string(char(altring))==charstr(altring)) {       // kein Parameter
   ring HNEring = char(altring),(x,y),ls;
   map m=altring,x,y;
   poly f=m(f);
   kill m;
 }
 else {
   string mipl=string(minpoly);
   if (mipl=="0") {
     " ** WARNING: No algebraic extension of this ground field is possible!";
     " ** We try to develop this polynomial, but if the need for an extension";
     " ** occurs during the calculation, we cannot proceed with the";
     " ** corresponding branches ...";
     execute("ring HNEring=("+charstr(basering)+"),(x,y),ls;");
 //--- ring ...=(char(.),`parstr()`),... geht nicht, wenn mehr als 1 Param. ---
   }
   else {
    string pa=parstr(altring);
    ring HNhelpring=p,`pa`,dp;
    execute("poly mipo="+mipl+";");  // Minimalpolynom in Polynom umgewandelt
    ring HNEring=(p,a),(x,y),ls;
    map getminpol=HNhelpring,a;
    mipl=string(getminpol(mipo));    // String umgewandelt mit 'a' als Param.
    execute("minpoly="+mipl+";");     // "minpoly=poly is not supported"
    kill HNhelpring, getminpol;
   }
   if (nvars(altring)==2) { poly f=fetch(altring,f); }
   else {
     map m=altring,x,y;
     poly f=m(f);
     kill m;
   }
 }
 export HNEring;

 if (defined(HNDebugOn))
 {"received polynomial: ",f,", with x =",namex,", y =",namey;}

 //----------------------- Variablendefinitionen ------------------------------
 int Abbruch,i,NullHNEx,NullHNEy;
 string str;
 list Newton,Ergebnis,hilflist;

 //====================== Tests auf Zulaessigkeit des Polynoms ================

 //-------------------------- Test, ob Einheit oder Null ----------------------
 if (subst(subst(f,x,0),y,0)!=0) {
   dbprint(printlevel+1,
           "The given polynomial is a unit in the power series ring!");
   keepring HNEring;
   return(list());                   // there are no HNEs
 }
 if (f==0) {
   dbprint(printlevel+1,"The given polynomial is zero!");
   keepring HNEring;
   return(list());                   // there are no HNEs
 }

 //-----------------------  Test auf Quadratfreiheit --------------------------

 if ((p==0) and (size(charstr(basering))==1)) {

 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring HNEring
 //----------------------------------------------------------------------------
  int testerg=(polytest(f)==0);
  ring zweitring = 32003,(x,y),dp;

  map polyhinueber=HNEring,x,y;         // fetch geht nicht
  poly f=polyhinueber(f);
  poly test_sqr=squarefree(f);
  if (test_sqr != f) {
   if (printlevel>0) {
     "Most probably the given polynomial is not squarefree. But the test was";
     "made in characteristic 32003 and not 0 to improve speed. You can";
     "(r) redo the test in char 0 (but this may take some time)";
     "(c) continue the development, if you're sure that the polynomial IS",
     "squarefree";
     if (testerg==1) {
       "(s) continue the development with a squarefree factor (*)";}
     "(q) or just quit the algorithm (default action)";
     "";"Please enter the letter of your choice:";
     str=read("")[1];             // reads one character
   }
   else { str="r"; }              // printlevel <= 0: non-interactive behaviour
   setring HNEring;
   map polyhinueber=zweitring,x,y;
   if (str=="r") {
     poly test_sqr=squarefree(f);
     if (test_sqr != f) {
      if (printlevel>0) { "The given polynomial is in fact not squarefree."; }
      else              { "The given polynomial is not squarefree!"; }
      "I'll continue with the radical.";
      f=test_sqr;
     }
     else {
      dbprint(printlevel,
       "everything is ok -- the polynomial is squarefree in characteristic 0");
     }
   }
   else {
     if ((str=="s") and (testerg==1)) {
       "(*)attention: it could be that the factor is only one in char 32003!";
       f=polyhinueber(test_sqr);
     }
     else {
       if (str<>"c") {
         setring altring;
         if(system("with","Namespaces")) { kill Top::HNEring; }
         kill HNEring;kill zweitring;
         return(list());}
       else { "if the algorithm doesn't terminate, you were wrong...";}
   }}
   kill zweitring;
   if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
  }
  else {
    setring HNEring;
    kill zweitring;
  }
 }
 //------------------ Fall Char > 0 oder Ring hat Parameter -------------------
 else {
  poly test_sqr=squarefree(f);
  if (test_sqr != f) {
   if (printlevel>0) {
    if (test_sqr == 1) {
     "The given polynomial is of the form g^"+string(p)+",";
     "therefore not squarefree.  You can:";
     " (q) quit the algorithm (recommended) or";
     " (f) continue with the full radical (using a factorization of the";
     "     pure power part; this could take much time)";
     "";"Please enter the letter of your choice:";
     str=read("")[1];
     if (str<>"f") { str="q"; }
    }
    else {
     "The given polynomial is not squarefree.";
     if (p != 0)
      {
       " You can:";
       " (c) continue with a squarefree divisor (but factors of the form g^"
       +string(p);
       "     are lost; this is recommended, takes no more time)";
       " (f) continue with the full radical (using a factorization of the";
       "     pure power part; this could take much time)";
       " (q) quit the algorithm";
       "";"Please enter the letter of your choice:";
       str=read("")[1];
       if ((str<>"f") && (str<>"q")) { str="c"; }
      }
     else { "I'll continue with the radical."; str="c"; }
    }                                // endelse (test_sqr!=1)
   }
   else {
     "//** Error: The given polynomial is not squarefree!";
     "//** Since the global variable `printlevel' has the value",printlevel,
       "we stop here.";
     "//   Either call me again with a squarefree polynomial f or assign";
     "            printlevel=1;";
     "//   before calling me with a non-squarefree f.";
     "//   If printlevel > 0, I will present to you some possibilities how to",
       "proceed.";
     str="q";
   }
   if (str=="q") {
    if(system("with","Namespaces")) { kill Top::HNEring; }
    setring altring;kill HNEring;
    return(list());
   }
   if (str=="c") { f=test_sqr; }
   if (str=="f") { f=allsquarefree(f,test_sqr); }
  }
  if (defined(HNDebugOn)) {"I continue with the polynomial",f; }

 }
 //====================== Ende Test auf Quadratfreiheit =======================
 if (subst(subst(f,x,0),y,0)!=0) {
   "Sorry. The remaining polynomial is a unit in the power series ring...";
   keepring HNEring;
   return(list());
 }
 //---------------------- Test, ob f teilbar durch x oder y -------------------
 if (subst(f,y,0)==0) {
   f=f/y; NullHNEy=1; }             // y=0 is a solution
 if (subst(f,x,0)==0) {
   f=f/x; NullHNEx=1; }             // x=0 is a solution

 Newton=newtonpoly(f,1);
 i=1; Abbruch=0;
 //----------------------------------------------------------------------------
 // finde Eckpkt. des Newtonpolys, der den Teil abgrenzt, fuer den x transvers:
 // Annahme: Newton ist sortiert, s.d. Newton[1]=Punkt auf der y-Achse,
 // Newton[letzt]=Punkt auf der x-Achse
 //----------------------------------------------------------------------------
 while ((i<size(Newton)) and (Abbruch==0)) {
  if ((Newton[i+1][1]-Newton[i][1])>=(Newton[i][2]-Newton[i+1][2]))
   {Abbruch=1;}
  else {i=i+1;}
 }
 int grenze1=Newton[i][2];
 int grenze2=Newton[i][1];
 //----------------------------------------------------------------------------
 // Stelle Ring bereit zur Uebertragung der Daten im Fall einer Koerperer-
 // weiterung. Definiere Objekte, die spaeter uebertragen werden.
 // Binde die Listen (azeilen,...) an den Ring (um sie nicht zu ueberschreiben
 // bei Def. in einem anderen Ring).
 // Exportiere Objekte, damit sie auch in der proc HN noch da sind
 //----------------------------------------------------------------------------
 ring HNE_noparam = char(altring),(a,x,y),ls;
 export HNE_noparam;
 poly f;
 list azeilen=ideal(0);
 list HNEs=ideal(0);
 list aneu=ideal(0);
 list faktoren=ideal(0);
 ideal deltais;
 poly delt;                   // nicht number, weil delta von a abhaengen kann
 export f,azeilen,HNEs,aneu,faktoren,deltais,delt;
 //----- hier steht die Anzahl bisher benoetigter Ringerweiterungen drin: -----
 int EXTHNEnumber=0; export EXTHNEnumber;
 setring HNEring;

 // ================= Die eigentliche Berechnung der HNE: =====================

 // ------- Berechne HNE von allen Zweigen, fuer die x transversal ist: -------
 if (defined(HNDebugOn))
   {"1st step: Treat Newton polygon until height",grenze1;}
 if (grenze1>0) {
  hilflist=HN(f,grenze1,1,essential);
  if (typeof(hilflist[1][1])=="ideal") { hilflist[1]=list(); }
 //- fuer den Fall, dass keine Zweige in transz. Erw. berechnet werden konnten-
  Ergebnis=extractHNEs(hilflist[1],0);
  if (hilflist[2]!=-1) {
    if (defined(HNDebugOn)) {" ring change in HN(",1,") detected";}
    poly transfproc=hilflist[2];
    map hole=HNE_noparam,transfproc,x,y;
    setring HNE_noparam;
    f=imap(HNEring,f);
    setring EXTHNEring(EXTHNEnumber);
    poly f=hole(f);
  }
 }
 if (NullHNEy==1) {
  Ergebnis=Ergebnis+list(list(matrix(ideal(0,x)),intvec(1),int(0),poly(0)));
 }
 // --------------- Berechne HNE von allen verbliebenen Zweigen: --------------
 if (defined(HNDebugOn))
    {"2nd step: Treat Newton polygon until height",grenze2;}
 if (grenze2>0) {
  map xytausch=basering,y,x;
  kill hilflist;
  def letztring=basering;
  if (EXTHNEnumber==0) { setring HNEring; }
  else                 { setring EXTHNEring(EXTHNEnumber); }
  list hilflist=HN(xytausch(f),grenze2,1,essential);
  if (typeof(hilflist[1][1])=="ideal") { hilflist[1]=list(); }
  if (not defined(Ergebnis)) {
 //-- HN wurde schon mal ausgefuehrt; Ringwechsel beim zweiten Aufruf von HN --
    if (defined(HNDebugOn)) {" ring change in HN(",1,") detected";}
    poly transfproc=hilflist[2];
    map hole=HNE_noparam,transfproc,x,y;
    setring HNE_noparam;
    list Ergebnis=imap(letztring,Ergebnis);
    setring EXTHNEring(EXTHNEnumber);
    list Ergebnis=hole(Ergebnis);
  }
  Ergebnis=Ergebnis+extractHNEs(hilflist[1],1);
 }
 if (NullHNEx==1) {
  Ergebnis=Ergebnis+list(list(matrix(ideal(0,x)),intvec(1),int(1),poly(0)));
 }
 //------------------- Loesche globale, nicht mehr benoetigte Objekte: --------
 if (EXTHNEnumber>0) {
  if(system("with","Namespaces")) { kill Top::HNEring; }
  kill HNEring;
  def HNEring=EXTHNEring(EXTHNEnumber);
  setring HNEring;
  export HNEring;
  kill EXTHNEring(1..EXTHNEnumber);
 }
 kill HNE_noparam;
 kill EXTHNEnumber;
 keepring basering;

 return(Ergebnis);
}

///////////////////////////////////////////////////////////////////////////////
//         
//  the command essdevelop is obsolete  --> here is the former help string:
//
///////////////////////////////////////////////////////////////////////////////
//ASSUME:  f is a bivariate polynomial (in the first 2 ring variables)
//CREATE:  ring with name @code{HNEring}, variables @code{x,y} and ordering
//         @code{ls} over a field extension of the current basering's ground 
//         field. @*
//         Since the Hamburger-Noether development of a reducible curve
//         singularity usually does not exist in the originally given basering,
//         @code{essdevelop} always defines @code{HNEring} and CHANGES to it.
//         The field extension is chosen minimally.
//RETURN:  list @code{L} of lists @code{L[i]} (corresponding to the output of
//         @code{develop(f[i])}, f[i] an \"essential\" branch of f, but the
//         last entry being omitted).@*
//         For more details type @code{help reddevelop;}.
//NOTE:    If the HNE needs a field extension, some of the branches will be
//         conjugate. In this case @code{essdevelop} reduces the computation to
//         one representative for each group of conjugate branches.@*
//         Note that the degree of each branch is in general less than the
//         degree of the field extension in which all HNEs can be put.@*
//         Use @code{reddevelop} or @code{HNdevelop} to compute a complete HNE,
//         i.e., a HNE for all branches.@*
//         If @code{printlevel>=0} comments are displayed (default is
//         @code{printlevel=0}).
//SEE ALSO: hnexpansion, develop, reddevelop, HNdevelop, extdevelop
//EXAMPLE: example essdevelop;  shows an example
proc essdevelop (poly f)
"USAGE:   essdevelop(f); f poly
NOTE:     command is obsolete, use hnexpansion(f,\"ess\") instead.
SEE ALSO: hnexpansion, develop, extdevelop, param
"
{
 list Ergebnis=pre_HN(f,1);
 dbprint(printlevel-voice+2,
    "// result: "+string(size(Ergebnis))+" branch(es) successfully computed;");
 if (string(minpoly) <> "0") {
   dbprint(printlevel-voice+2,
    "// note that conjugate branches are omitted and that the number",
    "// of branches represented by each remaining one may vary!");
 }
 dbprint(printlevel-voice+2,
    "// basering has changed to HNEring");
 keepring basering;
 return(Ergebnis);
}
example
{
  // -------- prepare for example ---------
  if (nameof(basering)=="HNEring") {
   def rettering=HNEring;
   kill HNEring;
  }
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring r=2,(x,y),dp;
  poly f=(x4+x2y+y2)*(x3+xy2+y3);
  // --------- compute all branches: ---------
  list hne=reddevelop(f);
  displayHNE(hne[1]);   // HN-matrix of 1st branch
  displayHNE(hne[4]);   // HN-matrix of 4th branch
  setring r;
  kill HNEring;
  // --- compute only one of conjugate branches: ---
  list hne=essdevelop(f);
  displayHNE(hne);
  // no. 1 of essdevelop represents no. 1 - 3 of reddevelop and
  // no. 2 of essdevelop represents no. 4 + 5 of reddevelop
  kill HNEring,r;
  echo = 0;
  // --- restore HNEring if previously defined ---
  if (defined(rettering)) {
   setring rettering;
   def HNEring=rettering;
   export HNEring;
  }
}

///////////////////////////////////////////////////////////////////////////////
static proc HN (poly f,int grenze, int Aufruf_Ebene, int essential)
"NOTE: This procedure is only for internal use, it is called via pre_HN"
{
 //---------- Variablendefinitionen fuer den unverzweigten Teil: --------------
 if (defined(HNDebugOn)) {"procedure HN",Aufruf_Ebene;}
 int Abbruch,ende,i,j,e,M,N,Q,R,zeiger,zeile,zeilevorher;
 intvec hqs;
 poly fvorher;
 list erg=ideal(0); list HNEs=ideal(0); // um die Listen an den Ring zu binden

 //-------------------- Bedeutung von Abbruch: --------------------------------
 //------- 0:keine Verzweigung | 1:Verzweigung,nicht fertig | 2:fertig --------
 //
 // Struktur von HNEs : Liste von Listen L (fuer jeden Zweig) der Form
 // L[1]=intvec (hqs), L[2],L[3],... ideal (die Zeilen (0,1,...) der HNE)
 // L[letztes]=poly (transformiertes f)
 //----------------------------------------------------------------------------
 list Newton;
 number delt;
 int p = char(basering);                // Ringcharakteristik
 list azeilen=ideal(0);
 ideal hilfid; list hilflist=ideal(0); intvec hilfvec;

 // ======================= der unverzweigte Teil: ============================
 while (Abbruch==0) {
  Newton=newtonpoly(f,1);
  zeiger=find_in_list(Newton,grenze);
  if (Newton[zeiger][2] != grenze)
    {"Didn't find an edge in the Newton polygon!";}
  if (zeiger==size(Newton)-1) {
    if (defined(HNDebugOn)) {"only one relevant side in Newton polygon";}
    M=Newton[zeiger+1][1]-Newton[zeiger][1];
    N=Newton[zeiger][2]-Newton[zeiger+1][2];
    R = M%N;
    Q = M / N;

 //-------- 1. Versuch: ist der quasihomogene Leitterm reine Potenz ? ---------
 //              (dann geht alles wie im irreduziblen Fall)
 //----------------------------------------------------------------------------
    e = gcd(M,N);
    delt=factorfirst(redleit(f,Newton[zeiger],Newton[zeiger+1])
                      /x^Newton[zeiger][1],M,N);
    if (delt==0) {
      if (defined(HNDebugOn)) {" The given polynomial is reducible !";}
      Abbruch=1;
    }
    if (Abbruch==0) {
 //-------------- f,zeile retten fuer den Spezialfall (###): ------------------
      fvorher=f;zeilevorher=zeile;
      if (R==0) {
 //------------- transformiere f mit T1, wenn kein Abbruch nachher: -----------
        if (N>1) { f = T1_Transform(f,delt,M/ e); }
        else     { ende=1; }
        if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delt;}
        azeilen[zeile+1][Q]=delt;
      }
      else {
 //------------- R > 0 : transformiere f mit T2 -------------------------------
        erg=T2_Transform(f,delt,M,N,referencepoly(Newton));
        f=erg[1];delt=erg[2];
 //------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ---------
        while (R!=0) {
         if (defined(HNDebugOn)) { "h("+string(zeile)+") =",Q; }
         hqs[zeile+1]=Q;         // denn zeile beginnt mit dem Wert 0
 //------------------ markiere das Zeilenende der HNE: ------------------------
         azeilen[zeile+1][Q+1]=x;
         zeile=zeile+1;
 //----------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ---------
         azeilen[zeile+1]=ideal(0);
         M=N; N=R; R=M%N; Q=M / N;
        }
        if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delt;}
        azeilen[zeile+1][Q]=delt;
      }
      if (defined(HNDebugOn)) {"transformed polynomial: ",f;}
      grenze=e;
 //----------------------- teste Abbruchbedingungen: --------------------------
      if (subst(f,y,0)==0) {              // <==> y|f
        dbprint(printlevel-voice+3,"finite HNE of one branch found");
           // voice abzufragen macht bei rekursiven procs keinen Sinn
        azeilen[zeile+1][Q+1]=x;
 //- Q wird nur in hqs eingetragen, wenn der Spezialfall nicht eintritt (s.u.)-
        Abbruch=2;
        if (grenze>1) {
         if (jet(f,1,intvec(0,1))==0) {
 //------ jet(...)=alle Monome von f, die nicht durch y2 teilbar sind ---------
 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
         else {
 //-------------------------- Spezialfall (###): ------------------------------
 // Wir haben das Problem, dass die HNE eines Zweiges hier abbricht, aber ein
 // anderer Zweig bis hierher genau die gleiche HNE hat, die noch weiter geht
 // Loesung: mache Transform. rueckgaengig und behandle f im Verzweigungsteil
 //----------------------------------------------------------------------------
          Abbruch=1;
          f=fvorher;zeile=zeilevorher;grenze=Newton[zeiger][2];
        }}
        else {f=0;}     // f nicht mehr gebraucht - spare Speicher
        if (Abbruch==2) { hqs[zeile+1]=Q; }
      }                 // Spezialfall nicht eingetreten
      else {
        if (ende==1) {
          dbprint(printlevel-voice+2,"HNE of one branch found");
          Abbruch=2; hqs[zeile+1]=-Q-1;}
      }
    }                   // end(if Abbruch==0)
  }                     // end(if zeiger...)
  else { Abbruch=1;}
 }                      // end(while Abbruch==0)

 // ===================== der Teil bei Verzweigung: ===========================

 if (Abbruch==1) {
 //---------- Variablendefinitionen fuer den verzweigten Teil: ----------------
  poly leitf,teiler,transformiert;
  list aneu=ideal(0);
  list faktoren;
  list HNEakut=ideal(0);
  ideal deltais;
  intvec eis;
  int zaehler,hnezaehler,zl,zl1,M1,N1,R1,Q1,needext;
  int numberofRingchanges,lastRingnumber,ringischanged,flag;
  string letztringname;

  zeiger=find_in_list(Newton,grenze);
  if (defined(HNDebugOn)) {
    "Branching part reached---Newton polygon :",Newton;
    "relevant part until height",grenze,", from",Newton[zeiger],"on";
  }
  azeilen=list(hqs)+azeilen; // hat jetzt Struktur von HNEs: hqs in der 1.Zeile

 //======= Schleife fuer jede zu betrachtende Seite des Newtonpolygons: =======
  for(i=zeiger; i<size(Newton); i++) {
   if (defined(HNDebugOn)) { "we consider side",Newton[i],Newton[i+1]; }
   M=Newton[i+1][1]-Newton[i][1];
   N=Newton[i][2]-Newton[i+1][2];
   R = M%N;
   Q = M / N;
   e=gcd(M,N);
   needext=1;
   letztringname=nameof(basering);
   lastRingnumber=EXTHNEnumber;
   faktoren=list(ideal(charPoly(redleit(f,Newton[i],Newton[i+1])
                       /(x^Newton[i][1]*y^Newton[i+1][2]),M,N)  ),
                 intvec(1));                  // = (zu faktoriserendes Poly, 1)

 //-- wechsle so lange in Ringerw., bis Leitform vollst. in Linearfakt. zerf.:-
   for (numberofRingchanges=1; needext==1; numberofRingchanges++) {
    leitf=redleit(f,Newton[i],Newton[i+1])/(x^Newton[i][1]*y^Newton[i+1][2]);
    delt=factorfirst(leitf,M,N);
    needext=0;
    if (delt==0) {

 //---------- Sonderbehandlung: faktorisiere einige Polynome ueber Q(a): -------
     if (charstr(basering)=="0,a") {
//*CL old: faktoren=factorize(charPoly(leitf,M,N),2);  // damit funktion. Bsp. Baladi 5
         faktoren=factorize(charPoly(leitf,M,N));  
     }
     else {
 //------------------ faktorisiere das charakt. Polynom: ----------------------
       if ((numberofRingchanges==1) or (essential==0)) {
         faktoren=factorlist(faktoren);
       }
       else {     // eliminiere alle konjugierten Nullstellen bis auf eine:
         ideal hilf_id;
         for (zaehler=1; zaehler<=size(faktoren[1]); zaehler++) {
           hilf_id=factorize(faktoren[1][zaehler])[1];
           if (size(hilf_id)>1) { faktoren[1][zaehler]=hilf_id[2]; }
           else                 { faktoren[1][zaehler]=hilf_id[1]; }
         }
       }
     }

     zaehler=1; eis=0;
     for (j=1; j<=size(faktoren[2]); j++) {
      teiler=faktoren[1][j];
      if (teiler/y != 0) {         // sonst war's eine Einheit --> wegwerfen!
        if (defined(HNDebugOn)) {"factor of leading form found:",teiler;}
        if (teiler/y2 == 0) {      // --> Faktor hat die Form cy+d
          deltais[zaehler]=-subst(teiler,y,0)/koeff(teiler,0,1); //=-d/c
          eis[zaehler]=faktoren[2][j];
          zaehler++;
        }
        else {
          dbprint(printlevel-voice+2,
             " Change of basering (field extension) necessary!");
          if (defined(HNDebugOn)) { teiler,"is not properly factored!"; }
          if (needext==0) { poly zerlege=teiler; }
          needext=1;
        }
      }
     }                             // end(for j)
    }
    else { deltais=ideal(delt); eis=e;}
    if (defined(HNDebugOn)) {"roots of char. poly:";deltais;
                             "with multiplicities:",eis;}
    if (needext==1) {
 //--------------------- fuehre den Ringwechsel aus: --------------------------
      ringischanged=1;
      if ((size(parstr(basering))>0) && string(minpoly)=="0") {
        " ** We've had bad luck! The HNE cannot completely be calculated!";
                                   // HNE in transzendenter Erw. fehlgeschlagen
        kill zerlege;
        ringischanged=0; break;    // weiter mit gefundenen Faktoren
      }
      if (parstr(basering)=="") {
        EXTHNEnumber++;
        splitring(zerlege,"EXTHNEring("+string(EXTHNEnumber)+")");
        poly transf=0;
        poly transfproc=0;
      }
      else {
        if (defined(translist)) { kill translist; } // Vermeidung einer Warnung
        if (numberofRingchanges>1) {  // ein Ringwechsel hat nicht gereicht
         list translist=splitring(zerlege,"",list(transf,transfproc,faktoren));
         poly transf=translist[1];
         poly transfproc=translist[2];
         list faktoren=translist[3];
        }
        else {
         if (defined(transfproc)) { // in dieser proc geschah schon Ringwechsel
          EXTHNEnumber++;
          list translist=splitring(zerlege,"EXTHNEring("
               +string(EXTHNEnumber)+")",list(a,transfproc));
          poly transf=translist[1];
          poly transfproc=translist[2];
         }
         else {
          EXTHNEnumber++;
          list translist=splitring(zerlege,"EXTHNEring("
               +string(EXTHNEnumber)+")",a);
          poly transf=translist[1];
          poly transfproc=transf;
        }}
      }
 //----------------------------------------------------------------------------
 // transf enthaelt jetzt den alten Parameter des Ringes, der aktiv war vor
 // Beginn der Schleife (evtl. also ueber mehrere Ringwechsel weitergereicht),
 // transfproc enthaelt den alten Parm. des R., der aktiv war zu Beginn der
 // Prozedur, und der an die aufrufende Prozedur zurueckgegeben werden muss
 // transf ist Null, falls der alte Ring keinen Parameter hatte,
 // das gleiche gilt fuer transfproc
 //----------------------------------------------------------------------------

 //------ Neudef. von Variablen, Uebertragung bisher errechneter Daten: -------
      poly leitf,teiler,transformiert;
      list aneu=ideal(0);
      ideal deltais;
      number delt;
      setring HNE_noparam;
      if (defined(letztring)) { kill letztring; }
      if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
      else                  { def letztring=HNEring; }
      f=imap(letztring,f);
      faktoren=imap(letztring,faktoren);
      setring EXTHNEring(EXTHNEnumber);
      map hole=HNE_noparam,transf,x,y;
      poly f=hole(f);
      if (not defined(faktoren)) {
        list faktoren=hole(faktoren);
      }
    }
   }    // end (while needext==1) bzw. for (numberofRingchanges)

   if (eis==0) { i++; continue; }
   if (ringischanged==1) {
    list erg,hilflist,HNEakut;            // dienen nur zum Sp. von Zwi.erg.
    ideal hilfid;
    erg=ideal(0); hilflist=erg; HNEakut=erg;

    hole=HNE_noparam,transf,x,y;
    setring HNE_noparam;
    azeilen=imap(letztring,azeilen);
    HNEs=imap(letztring,HNEs);

    setring EXTHNEring(EXTHNEnumber);
    list azeilen=hole(azeilen);
    list HNEs=hole(HNEs);
    kill letztring;
    ringischanged=0;
   }

 //============ Schleife fuer jeden gefundenen Faktor der Leitform: ===========
   for (j=1; j<=size(eis); j++) {
 //-- Mache Transf. T1 oder T2, trage Daten in HNEs ein, falls HNE abbricht: --

 //------------------------ Fall R==0: ----------------------------------------
    if (R==0) {
      transformiert = T1_Transform(f,number(deltais[j]),M/ e);
      if (defined(HNDebugOn)) {
        "a("+string(zeile)+","+string(Q)+") =",deltais[j];
        "transformed polynomial: ",transformiert;
      }
      if (subst(transformiert,y,0)==0) {
       dbprint(printlevel-voice+3,"finite HNE found");
       hnezaehler++;
 //------------ trage deltais[j],x ein in letzte Zeile, fertig: ---------------
       HNEakut=azeilen+list(poly(0));        // =HNEs[hnezaehler];
       hilfid=HNEakut[zeile+2]; hilfid[Q]=deltais[j]; hilfid[Q+1]=x;
       HNEakut[zeile+2]=hilfid;
       HNEakut[1][zeile+1]=Q;                // aktualisiere Vektor mit den hqs
       HNEs[hnezaehler]=HNEakut;
       if (eis[j]>1) {
        transformiert=transformiert/y;
        if (subst(transformiert,y,0)==0) {
 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
        else {
 //------ Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden --------
          eis[j]=eis[j]-1;
        }
       }
      }
    }
    else {
 //------------------------ Fall R <> 0: --------------------------------------
      erg=T2_Transform(f,number(deltais[j]),M,N,referencepoly(Newton));
      transformiert=erg[1];delt=erg[2];
      if (defined(HNDebugOn)) {"transformed polynomial: ",transformiert;}
      if (subst(transformiert,y,0)==0) {
       dbprint(printlevel-voice+3,"finite HNE found");
       hnezaehler++;
 //---------------- trage endliche HNE in HNEs ein: ---------------------------
       HNEakut=azeilen;           // dupliziere den gemeins. Anfang der HNE's
       zl=2;                      // (kommt schliesslich nach HNEs[hnezaehler])
 //----------------------------------------------------------------------------
 // Werte von:  zeile: aktuelle Zeilennummer der HNE (gemeinsamer Teil)
 //             zl   : die HNE spaltet auf; zeile+zl ist der Index fuer die
 // Zeile des aktuellen Zweigs; (zeile+zl-2) ist die tatsaechl. Zeilennr.
 // (bei 0 angefangen) der HNE  ([1] <- intvec(hqs), [2] <- 0. Zeile usw.)
 //----------------------------------------------------------------------------

 //---------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ------
       M1=M;N1=N;R1=R;Q1=M1/ N1;
       while (R1!=0) {
        if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; }
        HNEakut[1][zeile+zl-1]=Q1;
        HNEakut[zeile+zl][Q1+1]=x;
                                  // markiere das Zeilenende der HNE
        zl=zl+1;
 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
        HNEakut[zeile+zl]=ideal(0);

        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
       }
       if (defined(HNDebugOn)) {
         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delt;
       }
       HNEakut[zeile+zl][Q1]  =delt;
       HNEakut[zeile+zl][Q1+1]=x;
       HNEakut[1][zeile+zl-1] =Q1;     // aktualisiere Vektor mit hqs
       HNEakut[zeile+zl+1]=poly(0);
       HNEs[hnezaehler]=HNEakut;
 //-------------------- Ende der Eintragungen in HNEs -------------------------

       if (eis[j]>1) {
        transformiert=transformiert/y;
        if (subst(transformiert,y,0)==0) {
 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
         else {
 //--------- Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden -----
          eis[j]=eis[j]-1;
       }}
      }                           // endif (subst()==0)
    }                             // endelse (R<>0)

 //========== Falls HNE nicht abbricht: Rekursiver Aufruf von HN: =============
 //------------------- Berechne HNE von transformiert -------------------------
    if (subst(transformiert,y,0)!=0) {
     lastRingnumber=EXTHNEnumber;
     list HNerg=HN(transformiert,eis[j],Aufruf_Ebene+1,essential);
     if (HNerg[2]==-1) {          // kein Ringwechsel in HN aufgetreten
       aneu=HNerg[1];  }
     else {
       if (defined(HNDebugOn))
          {" ring change in HN(",Aufruf_Ebene+1,") detected";}
       list aneu=HNerg[1];
       poly transfproc=HNerg[2];

 //- stelle lokale Var. im neuen Ring wieder her und rette ggf. ihren Inhalt: -
       list erg,hilflist,faktoren,HNEakut;
       ideal hilfid;
       erg=ideal(0); hilflist=erg; faktoren=erg; HNEakut=erg;
       poly leitf,teiler,transformiert;

       map hole=HNE_noparam,transfproc,x,y;
       setring HNE_noparam;
       if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
       else                  { def letztring=HNEring; }
       HNEs=imap(letztring,HNEs);
       azeilen=imap(letztring,azeilen);
       deltais=imap(letztring,deltais);
       delt=imap(letztring,delt);
       f=imap(letztring,f);

       setring EXTHNEring(EXTHNEnumber);
       list HNEs=hole(HNEs);
       list azeilen=hole(azeilen);
       ideal deltais=hole(deltais);
       number delt=number(hole(delt));
       poly f=hole(f);
     }
     kill HNerg;
 //----------------------------------------------------------------------------
 // HNerg muss jedesmal mit "list" neu definiert werden, weil vorher noch nicht
 // ------- klar ist, ob der Ring nach Aufruf von HN noch derselbe ist --------

 //============= Verknuepfe bisherige HNE mit von HN gelieferten HNEs: ========
     if (R==0) {
       HNEs,hnezaehler=constructHNEs(HNEs,hnezaehler,aneu,azeilen,zeile,
                       deltais,Q,j);
     }
     else {
      for (zaehler=1; zaehler<=size(aneu); zaehler++) {
       hnezaehler++;
       HNEakut=azeilen;          // dupliziere den gemeinsamen Anfang der HNE's
       zl=2;                     // (kommt schliesslich nach HNEs[hnezaehler])
 //---------------- Trage Beitrag dieser Transformation T2 ein: ---------------
 //--------- Zur Bedeutung von zeile, zl: siehe Kommentar weiter oben ---------

 //--------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: -------
       M1=M;N1=N;R1=R;Q1=M1/ N1;
       while (R1!=0) {
        if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; }
        HNEakut[1][zeile+zl-1]=Q1;
        HNEakut[zeile+zl][Q1+1]=x;    // Markierung des Zeilenendes der HNE
        zl=zl+1;
 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
        HNEakut[zeile+zl]=ideal(0);
        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
       }
       if (defined(HNDebugOn)) {
         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delt;
       }
       HNEakut[zeile+zl][Q1]=delt;

 //--- Daten aus T2_Transform sind eingetragen; haenge Daten von HN an: -------
       hilfid=HNEakut[zeile+zl];
       for (zl1=Q1+1; zl1<=ncols(aneu[zaehler][2]); zl1++) {
        hilfid[zl1]=aneu[zaehler][2][zl1];
       }
       HNEakut[zeile+zl]=hilfid;
 //--- vorher HNEs[.][zeile+zl]<-aneu[.][2], jetzt [zeile+zl+1] <- [3] usw.: --
       for (zl1=3; zl1<=size(aneu[zaehler]); zl1++) {
         HNEakut[zeile+zl+zl1-2]=aneu[zaehler][zl1];
       }
 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
       hilfvec=HNEakut[1],aneu[zaehler][1];
       HNEakut[1]=hilfvec;
 //----------- weil nicht geht: liste[1]=liste[1],aneu[zaehler][1] ------------
       HNEs[hnezaehler]=HNEakut;
      }                     // end(for zaehler)
     }                      // endelse (R<>0)
    }                       // endif (subst()!=0)  (weiteres Aufblasen mit HN)

   }                        // end(for j) (Behandlung der einzelnen delta_i)

  }
  keepring basering;
  if (defined(transfproc)) { return(list(HNEs,transfproc)); }
  else                     { return(list(HNEs,poly(-1))); }
 // -1 als 2. Rueckgabewert zeigt an, dass kein Ringwechsel stattgefunden hat -
 }
 else {
  HNEs[1]=list(hqs)+azeilen+list(f); // f ist das transform. Poly oder Null
  keepring basering;
  return(list(HNEs,poly(-1)));
 //-- in dieser proc trat keine Verzweigung auf, also auch kein Ringwechsel ---
 }
}
///////////////////////////////////////////////////////////////////////////////

static proc constructHNEs (list HNEs,int hnezaehler,list aneu,list azeilen,
                    int zeile,ideal deltais,int Q,int j)
"NOTE: This procedure is only for internal use, it is called via HN"
{
  int zaehler,zl;
  ideal hilfid;
  list hilflist;
  intvec hilfvec;
  for (zaehler=1; zaehler<=size(aneu); zaehler++) {
     hnezaehler++;
     // HNEs[hnezaehler]=azeilen;            // dupliziere gemeins. Anfang
 //----------------------- trage neu berechnete Daten ein ---------------------
     hilfid=azeilen[zeile+2];
     hilfid[Q]=deltais[j];
     for (zl=Q+1; zl<=ncols(aneu[zaehler][2]); zl++) {
      hilfid[zl]=aneu[zaehler][2][zl];
     }
     hilflist=azeilen; hilflist[zeile+2]=hilfid;
 //------------------ haenge uebrige Zeilen von aneu[] an: --------------------
     for (zl=3; zl<=size(aneu[zaehler]); zl++) {
      hilflist[zeile+zl]=aneu[zaehler][zl];
     }
 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
     if (hilflist[1]==0) { hilflist[1]=aneu[zaehler][1]; }
     else { hilfvec=hilflist[1],aneu[zaehler][1]; hilflist[1]=hilfvec; }
     HNEs[hnezaehler]=hilflist;
  }
  return(HNEs,hnezaehler);
}
///////////////////////////////////////////////////////////////////////////////

proc referencepoly (list newton)
"USAGE:   referencepoly(newton);
         newton is list of intvec(x,y) which represents points in the Newton
         diagram (e.g. output of the proc newtonpoly)
RETURN:  a polynomial which has newton as Newton diagram
SEE ALSO: newtonpoly
EXAMPLE: example referencepoly;  shows an example
"
{
 poly f;
 for (int i=1; i<=size(newton); i++) {
   f=f+var(1)^newton[i][1]*var(2)^newton[i][2];
 }
 return(f);
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y),ds;
 referencepoly(list(intvec(0,4),intvec(2,3),intvec(5,1),intvec(7,0)));
}
///////////////////////////////////////////////////////////////////////////////

proc factorlist (list L)
"USAGE:   factorlist(L);   L a list in the format of `factorize'
RETURN:  the nonconstant irreducible factors of
         L[1][1]^L[2][1] * L[1][2]^L[2][2] *...* L[1][size(L[1])]^...
         with multiplicities (same format as factorize)
SEE ALSO: factorize
EXAMPLE: example factorlist;  shows an example
"
{
 // eine Sortierung der Faktoren eruebrigt sich, weil keine zwei versch.
 // red.Fakt. einen gleichen irred. Fakt. haben koennen (I.3.27 Diplarb.)
 int i,gross;
 list faktoren,hilf;
 ideal hil1,hil2;
 intvec v,w;
 for (i=1; (L[1][i] == jet(L[1][i],0)) && (i<size(L[1])); i++) {;}
 if (L[1][i] != jet(L[1][i],0)) {
   hilf=factorize(L[1][i]);
 // Achtung!!! factorize(..,2) wirft entgegen der Beschreibung nicht nur
 // konstante Faktoren raus, sondern alle Einheiten in der LOKALISIERUNG nach
 // der Monomordnung!!! Im Beispiel unten verschwindet der Faktor x+y+1, wenn
 // man ds statt dp als Ordnung nimmt!
   hilf[2]=hilf[2]*L[2][i];
   hil1=hilf[1];
   gross=size(hil1);
   if (gross>1) {
       // faktoren=list(hilf[1][2..gross],hilf[2][2..gross]);
       // -->  `? indexed object must have a name'
     v=hilf[2];
     faktoren=list(ideal(hil1[2..gross]),intvec(v[2..gross]));
   }
   else         { faktoren=hilf; }
 }
 else {
   faktoren=L;
 }

 for (i++; i<=size(L[2]); i++) {
 //------------------------- linearer Term -- irreduzibel ---------------------
   if (L[1][i] == jet(L[1][i],1)) {
     if (L[1][i] != jet(L[1][i],0)) {           // konst. Faktoren eliminieren
       hil1=faktoren[1];
       hil1[size(hil1)+1]=L[1][i];
       faktoren[1]=hil1;
       v=faktoren[2],L[2][i];
       faktoren[2]=v;
     }
   }
 //----------------------- nichtlinearer Term -- faktorisiere -----------------
   else {
     hilf=factorize(L[1][i]);
     hilf[2]=hilf[2]*L[2][i];
     hil1=faktoren[1];
     hil2=hilf[1];
     gross=size(hil2);
       // hil2[1] ist konstant, wird weggelassen:
     hil1[(size(hil1)+1)..(size(hil1)+gross-1)]=hil2[2..gross];
       // ideal+ideal does not work due to simplification;
       // ideal,ideal not allowed
     faktoren[1]=hil1;
     w=hilf[2];
     v=faktoren[2],w[2..gross];
     faktoren[2]=v;
   }
 }
 return(faktoren);
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y),ds;
 list L=ideal(x,(x-y)^2*(x+y+1),x+y),intvec(2,2,1);
 L;
 factorlist(L);
}
///////////////////////////////////////////////////////////////////////////////

proc deltaHNE(list hne)
"USAGE:   deltaHNE(L);  L list
NOTE:     command is obsolete, use hnexpansion(f,\"ess\") instead.   
SEE ALSO: delta, deltaLoc
"
{
   int i,j,inters;
   int n=size(hne);
   list INV;
   for (i=1;i<=n;i++)
   {
     INV[i]=invariants(hne[i]);
   }
   int del=INV[n][5]/2;
   for(i=1;i<=n-1;i++)
   {
      del=del+INV[i][5]/2;
      for(j=i+1;j<=n;j++)
      {
         inters=inters+intersection(hne[i],hne[j]);
      }
   }    
   return(del+inters);
}

///////////////////////////////////////////////////////////////////////////////

proc delta
"USAGE:  delta(INPUT);  INPUT a polynomial defining an isolated  plane curve
         singularity at 0, or the Hamburger-Noether expansion thereof, i.e.
         the output of @code{develop(f)}, or the output of @code{hnexpansion(f[,\"ess\"])},
         or (one of the entries of) the list @code{hne} in the ring created
         by @code{hnexpansion(f[,\"ess\"])}.
RETURN:  the delta invariant of the singularity at 0, the vector space
         dimension of R~/R, where R~ is the normalization of the
         singularity R=basering/f
NOTE:    In case the Hamburger-Noether expansion of the curve f is needed
         for other purposes as well it is better to calculate this first
         with the aid of @code{hnexpansion} and use it as input instead of
         the polynomial itself.
SEE ALSO: deltaLoc, invariants
KEYWORDS: delta invariant
EXAMPLE: example delta;  shows an example
"
{
  if (typeof(#[1])=="poly")
  { // INPUT = polynomial defining the curve
    list HNEXPANSION=hnexpansion(#[1]);
    return(delta(HNEXPANSION));
  }
  if (typeof(#[1])=="ring")
  { // INPUT = HNEring of curve
    def r_e_t_t_e_r_i_n_g=basering;
    def H_N_E_RING=#[1];
    setring H_N_E_RING;
    int del=delta(hne);
    setring r_e_t_t_e_r_i_n_g;
    kill H_N_E_RING;    
    return(del);
  }
  if (typeof(#[1])=="matrix")  
  { // INPUT = hne of an irreducible curve  
    return(invariants(#)[5]/2);
  }
  else 
  { // INPUT = hne of a reducible curve
    list INV=invariants(#);
    return(INV[size(INV)][3]);
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = 32003,(x,y),ds;
  poly f = x25+x24-4x23-1x22y+4x22+8x21y-2x21-12x20y-4x19y2+4x20+10x19y
           +12x18y2-24x18y-20x17y2-4x16y3+x18+60x16y2+20x15y3-9x16y
           -80x14y3-10x13y4+36x14y2+60x12y4+2x11y5-84x12y3-24x10y5
           +126x10y4+4x8y6-126x8y5+84x6y6-36x4y7+9x2y8-1y9;
  delta(f);
}

///////////////////////////////////////////////////////////////////////////////


proc hnexpansion(poly f,list #)
"USAGE:   hnexpansion(f); or hnexpansion(f,\"ess\");  f poly

USAGE:   hnexpansion(f); f poly
ASSUME:  f is a bivariate polynomial (in the first 2 ring variables)
CREATE:  ring with variables @code{x,y} and ordering @code{ls} over a 
         field extension of the current basering's ground field,         
         since the Hamburger-Noether development usually does not exist
         in the originally given basering. The field extension is chosen
         minimally.@*
         Moreover, in the ring a list @code{hne} of lists @code{hne[i]} is 
         created (corresponding to the output of @code{develop(f[i])}, 
         f[i] a branch of f, but the last entry being omitted).
@texinfo
@table @asis
@item @code{hne[i][1]}; matrix:
         Each row contains the coefficients of the corresponding line of the
         Hamburger-Noether expansion (HNE) for f[i]. The end of the line is
         marked in the matrix by the first ring variable (usually x).
@item @code{hne[i][2]}; intvec:
         indicating the length of lines of the HNE
@item @code{hne[i][3]}; int:
         0  if the 1st ring variable was transversal (with respect to f[i]), @*
         1  if the variables were changed at the beginning of the
            computation, @*
        -1  if an error has occurred.
@item @code{hne[i][4]}; poly:
         the transformed polynomial of f[i] to make it possible to extend the
         Hamburger-Noether development a posteriori without having to do
         all the previous calculation once again (0 if not needed)
@end table
@end texinfo
RETURN:  a list, say @code{hn}, containing the created ring 
NOTE:    to use the ring type: @code{def HNEring=hn[i]; setring HNEring;}.
         @*
         If f is known to be irreducible as a power series, @code{develop(f)}
         could be chosen instead to avoid the change of basering. @*
         Increasing  @code{printlevel} leads to more and more comments.

USAGE:   hnexpansion(f,\"ess\"); f poly
ASSUME:  f is a bivariate polynomial (in the first 2 ring variables)
CREATE:  ring with variables @code{x,y} and ordering @code{ls} over a 
         field extension of the current basering's ground field,         
         since the Hamburger-Noether development usually does not exist
         in the originally given basering. The field extension is chosen
         minimally.
         @*
         Moreover, in the ring a list @code{hne} of lists @code{hne[i]} is 
         created (corresponding to the output of @code{develop(f[i])}, f[i] an
         \"essential\" branch of f, but the last entry being omitted). See 
         @code{hnexpansion} above for more details.
RETURN:  a list, say @code{hn}, containing the created ring 
NOTE:    to use the ring type: @code{def hnering=hn[i]; setring hnering;}.
         @*
         Alternatively you may use the procedure sethnering and type: 
         @code{sethnering(hn);}
         @*
         If the HNE needs a field extension, some of the branches will be
         conjugate. In this case @code{hnexpansion(f,\"ess\")} reduces the 
         computation to one representative for each group of conjugate 
         branches.@*
         Note that the degree of each branch is in general less than the degree
         of the field extension in which all HNEs can be put.@*
         Use @code{hnexpansion(f)} to compute a complete HNE, i.e., a HNE for 
         all branches.@*
         Increasing  @code{printlevel} leads to more and more comments.        
SEE ALSO: develop, extdevelop, parametrisation, displayHNE
EXAMPLE: example hnexpansion;  shows an example
"
{
  def rettering=basering;
  if (defined(HNEring))
  {
    def @HNEring = HNEring;  
    kill HNEring;
  }
  if (size(#)==1)
  {
    list hne=essdevelop(f);
  }
  else
  {
    list hne=HNdevelop(f);
  }
  export hne;
  list hnereturn=HNEring;
  setring rettering;
  kill HNEring;
  if (defined(@HNEring))
  {
    def HNEring=@HNEring;  
    export(HNEring);
  }
  dbprint(printlevel-voice+2,"
// 'hnexpansion' created a list containing a ring, which
// contains the Hamburger-Noether expansion as a list hne.
// To see the ring, type (if the name of your list is hn):
     show(hn);
// To access the ring and list, type:
     def hnering = hn[1]; 
     setring hnering;
     hne;
////////////////////////////////////////////////");

  return(hnereturn);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),ls;
  list hn=hnexpansion(x4-y6);
  show(hn);
  
  def hnering=hn[1];
  setring hnering;
  size(hne);           // number of branches
  print(hne[1][1]);    // HN-matrix of 1st branch
  parametrisation(hne);    // parametrization of the two branches
  /////////////////////////////////////////////////////////
  ring s=2,(x,y),ls;
  poly f=(x4+x2y+y2)*(x3+xy2+y3);
  // --------- compute all branches: ---------
  hn=hnexpansion(f);
  hnering=hn[1];
  setring hnering;
  displayHNE(hne[1]);   // HN-matrix of 1st branch
  displayHNE(hne[4]);   // HN-matrix of 4th branch
  setring s;
  // --- compute only one of conjugate branches: ---
  hn=hnexpansion(f,"ess");
  hnering=hn[1];
  setring hnering;
  displayHNE(hne);
  // no. 1 of hnexpansion(f,"ess") represents no. 1 - 3 of hnexpansion(f) and
  // no. 2 of hnexpansion(f,"ess") represents no. 4 + 5 of hnexpansion(f)
}
///////////////////////////////////////////////////////////////////////////////

proc sethnering (list L,list #)
"USAGE:  sethnering(L[,s]); L list, s string (optional)
ASSUME:  L is a list containing a ring (e.g. the output of @code{hnexpansion}).
CREATE:  The procedure creates a ring with name given by the optional parameter
         s resp. with name hnering, if no optional parameter is given, and
         changes your ring to this ring. The new ring will be the ring given 
         as the first entry in the list L.          
RETURN:  nothing. 
SEE ALSO: hnexpansion
EXAMPLE: example sethnering; shows some examples.
"
 
{
  if (typeof(L[1])=="ring")
  {
    if (size(#)>0)
    {
      if (typeof(#[1])=="string")
      {        
        execute("def "+#[1]+"=L[1];");
        execute("export "+#[1]+";");
        execute("setring "+#[1]+";");
        execute("keepring "+#[1]+";");
      }
      else
      {
        ERROR("Optional Input was no string.");
        return();
      }
    }
    else
    {
      def hnering=L[1];
      export hnering;
      setring hnering;
      keepring hnering;
    }
    return();
  }
  else
  {
    ERROR("Input was no hnering.");
    return();
  }
}
example
{
  // -------- prepare for example ---------
  if (defined(hnering))
  {
   def rette@ring=hnering;
   if (nameof(basering)=="hnering")
   {
     int wechsel=1;
   }
   else
   {
     int wechsel;
   }
   kill hnering;
  }
  // ------ the example starts here -------
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),ls;
  nameof(basering);
  sethnering(hnexpansion(x4-y6)); // Creates hnering and changes to it!
  nameof(basering);
  echo = 0;
  // --- restore HNEring if previously defined ---
  kill hnering;
  if (defined(rette@ring)) {
   def hnering=rette@ring;
   export hnering;
   if (wechsel==1)
   {
     setring hnering;
   }
  }
}