source: git/Singular/LIB/hnoether.lib @ c67136

// $Id: hnoether.lib,v 1.9 1998-10-21 10:28:24 krueger Exp $
// author:  Martin Lamm,  email: lamm@mathematik.uni-kl.de
// last change:           26.03.98
///////////////////////////////////////////////////////////////////////////////

version="$Id: hnoether.lib,v 1.9 1998-10-21 10:28:24 krueger Exp $";
info="
LIBRARY:  hnoether.lib   PROCEDURES FOR THE HAMBURGER-NOETHER-DEVELOPMENT

           Important procedures:
 develop(f [,n]);    Hamburger-Noether development from the irred. polynomial f
 reddevelop(f);      Hamburger-Noether development from the (red.) polynomial f
 extdevelop(hne,n);  extension of Hamburger-Noether development hne from f
 param(hne);         returns a parametrization of f (input=output(develop))
 displayHNE(hne);    display Hamburger-Noether development as an ideal
 invariants(hne);    invariants of f, e.g. the characteristic exponents
 displayInvariants(hne);  display invariants of f
 generators(hne);    computes the generators of the semigroup of values
 intersection(hne1,hne2); intersection multiplicity of two curves
 stripHNE(hne);      reduce amount of memory consumed by hne
";
////////////////////////////////////////////////////////////////////////////
//           perhaps useful procedures:
// puiseux2generators(m,n); convert puiseux pairs to generators of semigroup
// multiplicities(hne);     multiplicities of blowed up curves
// newtonpoly(f);      newtonpolygon of polynom f
// newtonhoehne(f);    same as newtonpoly, but uses internal procedure
// getnm(f);           intersection points of the Newton-polygon with the axes
// T_Transform(f,Q,N); returns f(y,xy^Q)/y^NQ (type f(x,y): poly, Q,N: int)
// T1_Transform(f,d,M); returns f(x,y+d*x^M)  (type f(x,y): poly,d:number,M:int)
// T2_Transform(f,d,M,N);   a composition of T1 & T
// koeff(f,I,J);       gets coefficient of indicated monomial of poly f (I,J:int)
// redleit(f,S,E);     restriction of monomials of f to line (S-E)
// squarefree(f);      returns a squarefree divisor of the poly f
// allsquarefree(f,l); returns the maximal squarefree divisor of the poly f
// set_list(L,v,l);    managing interconnected lists
//
//           procedures (more or less) for internal use:
// leit(f,n,m);        special case of redleit (for irred. polynomials)
// testreducible(f,n,m); tests whether f is reducible
// polytest(f);        tests coefficients and exponents of poly f
// 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 in a trivial way before using 'factorize'
// extrafactor(f,M,N); try to factor charPoly(f,M,N) where 'factorize' cannot
// extractHNEs(H,t);   extracts output H of HN to output of reddevelop
// HN(f,grenze);       recursive subroutine for reddevelop
// constructHNEs(...); subroutine for HN
///////////////////////////////////////////////////////////////////////////////
LIB "primitiv.lib";
///////////////////////////////////////////////////////////////////////////////

proc getnm (poly f)
"USAGE:   getnm(f); f polynomial in x,y
ASSUME:  basering = ...,(x,y),ls;  or ds
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
EXAMPLE: example getnm; shows an example
"
{
 // wir gehen davon aus, dass die Prozedur von develop aufgerufen wurde, also
 // die Ringordnung ist ls (ds ginge auch)
 return(ord(subst(f,x,0)),ord(subst(f,y,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)
"// leit(poly f, int n, int m) returns all monomials on the line
// from (0,n) to (m,0) in the Newton diagram
"
{
 return(jet(f,m*n,intvec(n,m))-jet(f,m*n-1,intvec(n,m)))
}
///////////////////////////////////////////////////////////////////////////////
proc testreducible (poly f, int n, int m)
"// testreducible(poly f,int n,int m) returns true if there are points in the
// Newton diagram below the line (0,n)-(m,0)
"
{
 return(size(jet(f,m*n-1,intvec(n,m))) != 0)
}
///////////////////////////////////////////////////////////////////////////////
proc T_Transform (poly f, int Q, int N)
"// returns f(y,xy^Q)/y^NQ
"
{
 map T = basering,y,x*y^Q;
 return(T(f)/y^(N*Q));
}
///////////////////////////////////////////////////////////////////////////////
proc T1_Transform (poly f, number d, int M)
"// returns f(x,y+d*x^M)
"
{
 map T1 = basering,x,y+d*x^M;
 return(T1(f));
}
///////////////////////////////////////////////////////////////////////////////

proc T2_Transform (poly f, number d, int M, int N)
"USAGE: T2_Transform(f,d,M,N); f poly, d number; M,N int
RETURN: list: poly T2(f,d',M,N), number d' in \{ d, 1/d \}
"
{
 //---------------------- 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];
  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,x,y+d*x^Q;
  map Tstar=basering,x^(N-Q*tau)*y^tau,x^(M-sigma*Q)*y^sigma;
  if (defined(Protokoll)) {
   "Trafo. T2: x->x^"+string(N-Q*tau)+"*y^"+string(tau)+", y->x^"
    +string(M-sigma*Q)+"*y^"+string(sigma);
   "delta =",d,"Q =",Q,"tau,sigma =",tau,sigma;
  }
 //------------------- Durchfuehrung der Transformation T2 --------------------
 // man koennte an T2_Transform noch eine Liste mit allen Eckpunkten des
 // Newtonpolys uebergeben und nur die Transformierten der Eckpunkte betrachten
 // um max. Exponenten e,f von x,y zu finden--noch einfacher: ein Referenzpo-
 // lynom, das das gleiche Newtonpolygon hat (z.B. Einschraenkung auf die Mono-
 // me in den Eckpunkten, oder alle Monome haben Koeff. 1, das ist einfacher)
 //----------------------------------------------------------------------------

  f=Tstar(f); // hier dann refpoly=Tstar(refpoly)
  poly hilf;
 // dividiere f so lange durch x, wie die Div. aufgeht:
  for (hilf=f/x; hilf*x==f; hilf=f/x) {f=hilf;}
  for (hilf=f/y; hilf*y==f; hilf=f/y) {f=hilf;} // gleiches fuer y
  return(list(T1(f),d));
}
///////////////////////////////////////////////////////////////////////////////

proc koeff (poly f, int I, int J)
"USAGE:  koeff(f,I,J); f polynomial in x and y, 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)
EXAMPLE: example koeff;  shows an example
"{
  matrix mat = coeffs(coeffs(f,y)[J+1,1],x);
  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 in x and y
RETURN: a squarefree divisor of f
        Normally the return value is a greatest squarefree divisor, but there
        is an exception: factors with a p-th root, p the ring characteristic,
        are lost.
EXAMPLE: example squarefree; shows some examples.
"{
 // es gibt noch ein Problem: syz gibt manchmal in Koerpererweiterung (p,a)
 // nicht f/g zurueck, obwohl die Division aufgeht ... In dem Fall ist die
 // Rueckgabe evtl. nicht quadratfrei! Bsp.: squarefree((1/a*x3+y7)^3);
 // in char 7
 // In Singular 1.1.3 funktioniert syz, dafuer geht das gleiche Bsp. schief,
 // weil gcd auf einmal falsch wird! (ring (7,a),(x,y),dp; minpoly=a2+1; )

 //----------------- 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 (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]; }
     }
 //----------------------------------------------------------------------------
 // Polynomdivision geht, wenn factory eingebunden ist
 // dann also l=f/g; aber: ist ebenfalls fehlerhaft! (s.o. als Bsp.)
 // wenn / zur Polynomdivision einsetzbar ist, kann der Ringwechsel gespart
 // werden, sollte dann also zum Einsatz gebracht werden (Zeitersparnis)
 //----------------------------------------------------------------------------
    }
    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 (deg(l)==deg(f)) { e=1;}
            else {e=0;}
     }
    }
    else { l=syz(ideal(dif,f))[1,1];
           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);
 squarefree(x3);
 squarefree(x2);
 squarefree(xy+x);
 squarefree((x+y)^3*(x-y)^2); // (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,l);  poly f,l
        l is the squarefree divisor of f you get by l=squarefree(f);
RETURN: the squarefree divisor of f consisting of all irreducible factors of f
NOTE  : * this proc uses factorize to get the missing factors of f not in l and
          therefore may be slow. Furthermore factorize can still return
          a factor more than once (what should not occur)
        * this proc uses polynomial division of factory (which still has a bug
         at the moment: it can return nonsense, if the basering has parameters)
EXAMPLE: example allsquarefree;  shows an example
"{
 //---------- eliminiere bereits mit squarefree gefundene Faktoren ------------
 poly g=l;                            // g = gcd(f,l);
 while (deg(g)!=0) {
  f=f/g;
  g=gcd(f,l);
 }                                    // jetzt f=h^p
 //--------------- Berechne uebrige Faktoren mit factorize --------------------
 if (deg(f)>0) {
  g=1;
  ideal factf=factorize(f,1);
  for (int i=1; 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 the correct result..."; g=testp;}
  }
  l=l*g;
 }
 return(l);
}
example
{ "EXAMPLE:"; echo = 2;
  ring exring=7,(x,y),dp;
  poly f=(x+y)^7*(x-y)^8;
  poly l=squarefree(f);
  l; // x-y found because 7 does not divide 8
  allsquarefree(f,l); // all factors (x+y)*(x-y) found
}
///////////////////////////////////////////////////////////////////////////////

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 polynomial in two variables, n integer
ASSUME:  f irreducible as a power series (for reducible f use reddevelop)
RETURN:  Hamburger-Noether development of f, in form of a list:
         [1]: Hamburger-Noether matrix
            each row contains the coefficients of the corresponding line of the
            Hamburger-Noether extension (HNE). The end of the line is marked in
            the matrix as the first ring variable (usually x)
         [2]: intvec indicating the length of lines of the HNE
         [3]: int:
           0  if the 1st ring variable was transverse (with respect to f),
           1  if the variables were changed at the beginning of the computation
           -1 if an error has occurred
         [4]: 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)
DISPLAY: The (non zero) elements of the HNE

NOTE:    If the optional parameter n is given, the HN-matrix will have at least
         n columns. Otherwise the number of columns will be chosen minimal s.t.
         the matrix contains all necessary information (i.e. all lines of the
         HNE but the last (which is in general infinite) appear).
         If n is negative, the algorithm is stopped as soon as possible, i.e.
         the information computed is enough for 'invariants', but the HN-matrix
         may contain undetermined elements, which are marked with the
         2nd variable (of the basering).
         In any case, the optional parameter only affects the calculation of
         the LAST line of the HNE; develop(f) gives already all necessary
         information for the procedure 'invariants'. A negative value of n will
         result in a very poor parametrization, but it can make 'develop'
         faster; a positive value will improve the exactness of the
         parametrization.

         For time critical computations it is recommended to use
         \"ring ...,(x,y),ls\" as basering - it increases the algorithm's speed.

EXAMPLES: example develop; shows an example
          example param;   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)));}
 if (nvars(basering) > 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);

 //------------- 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.";
      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(Protokoll))
 {"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
 intvec hqs;                // erhaelt die Werte von h(zeile)=Q;

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

 number c,delta;
 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)-2]+";";
                            // 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) {"You have given me 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)
   {" The given polynomial is a unit as power series!!"; Abbruch=1; Ausgabe=1;}
   else {
    if (N == -1) {"The HNE is x = 0"; Abbruch=1; Ausgabe=2; getauscht=1;}
    else {
     if (M == -1) {"The HNE is y = 0"; Abbruch=1; Ausgabe=2;}
   }}
   if (N+M==-2) {"Remark: The polynomial is divisible both by x and y,",
                 "so it is not irreducible.";
                 "The result of develop is only one of the existing HNE's.";
   }
 }
 //--------------------- 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) {
    "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];
    setring guenstig;
    map polyhinueber=zweitring,x,y;
    if (str=="r") {
      poly test_sqr=squarefree(f);
      if (test_sqr != f) {
       "The given polynomial is in fact not squarefree.";
       "I'll continue with the radical.";
       pause;
       f=test_sqr;
      }
      else {"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;
          if(system("with","Namespaces")) { kill Top::guenstig; }
          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(Protokoll)) {"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;
    if(system("with","Namespaces")) { kill Top::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.";}
    pause;
    f=test_sqr;
    nm = getnm(f);              // N,M haben sich veraendert
    N = nm[1]; M = nm[2];       // Berechne Schnittpunkte Newtonpoly mit Achsen
    if (defined(Protokoll)) {"I continue with the polynomial",f; }
   }

  }                             // endelse(p==0)

  if (N==0) {
    " Sorry. The remaining polynomial is a unit in the power series ring...";
    if(system("with","Namespaces")) { kill Top::guenstig; }
    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(Protokoll)) {
   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)) {
    " The given polynomial is not irreducible";
    kill guenstig;
    if(system("with","Namespaces")) { kill Top::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) {
      delta = koeff(f,M/ e,N - N/ e) / (-1*e*c);
      if (defined(Protokoll)) {"quasihomogeneous leading form:",
         leit(f,N,M)," = ",c,"* (y -",delta,"* x^"+string(M/ e)+")^",e," ?";}
      if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e)
        { " The given Polynomial is reducible !"; Abbruch=1; Ausgabe=1;}
    }
    else {                     // p!=0
      if (e%p != 0) {
        delta = koeff(f,M/ e,N - N/ e) / (-1*e*c);
        if (defined(Protokoll)) {"quasihomogeneous leading form:",
           leit(f,N,M)," = ",c,"* (y -",delta,"* x^"+string(M/ e)+")^",e," ?";}
        if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e)
         { " 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(Protokoll)) {e," -> ",eps,"*",p,"^",l;}
        delta = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);

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

 //-- 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                                        -
 //----------------------------------------------------------------------------
          }}
          delta = generat^(extgcd(n_elements-1,p^l)[3]*exponent);
        }

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

      nm=getnm(f); N=nm[1]; M=nm[2];        // Neuberechnung des Newtonpolygons
  }
 //--------------------------- Fall R > 0 : -----------------------------------
  else {
    "h("+string(zeile)+") =",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(Protokoll)) {"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) {
     "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);
     f=0;               // transformiertes Polynom wird nicht mehr gebraucht
     Abbruch=1;
  }
  else {if (M<N) {"Something has gone wrong: M<N";}}
  if(defined(Protokoll)) {"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) { f=zurueck(f); }
 }

 if(system("with","Namespaces")) { kill Top::guenstig; }
 kill guenstig;
 if (Ausgabe == 0) { return(list(amat,hqs,getauscht,f));}
 if (Ausgabe == 1) { return(return_error);}             // error has occurred
 if (Ausgabe == 2) { return(list(matrix(ideal(0,x)),intvec(1),getauscht,
                                 poly(0)));}            // 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),   (note that there is  1 zero   in the 2nd row before x)
  // z(1) = z(2)^3*z(3), (note that there are 3 zeroes in the 3rd row before x)
  // 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 matrix indicates that this line is not complete.
  //    It can only occur in the last line of the HNE, and normally does.)
  hne[2]; pause;
  param(hne);
  // returns the parametrisation x(t)= -t14+O(t21), y(t)= -3t98+O(t105)
  // (the term -t109 in y may have a wrong coefficient)
  displayHNE(hne);
}

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

proc param
"USAGE:  param(l) takes the output of develop(f)
        (list l (matrix m, intvec v, int s[,poly g]))
        and gives a parametrisation for f; the first variable of the active
        ring is chosen as indefinite. If the ring contains more than
        two variables, the 3rd variable is chosen (remember that develop takes
        the first two variables and therefore the other vars should be unused)
RETURN: ideal of two polynomials: if the HNE was finite, f(param[1],param[2])
        will be  zero. If not, the real parametrisation will be
        two power series; then param will return a truncation of these series.
EXAMPLE: example param;      shows an example
         example develop;   shows another example

"{
 //-------------------------- Initialisierungen -------------------------------
 matrix m=#[1];
 intvec v=#[2];
 int switch=#[3];
 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)
      { "// ** Warning: result is exact up to order",fehlervor-1,"in",
        maxideal(1)[1],"and",fehler-1,"in",maxideal(1)[2],"!";}
     else { "// ** Warning: result is exact up to order",fehler-1,"in",
        maxideal(1)[2],"!";}}
   return(ideal(z(2),z(1)));}
 else {
   if (fehler>0) {
     if (fehlervor>0)
      { "// ** Warning: result is exact up to order",fehler-1,"in",
        maxideal(1)[1],"and",fehlervor-1,"in",maxideal(1)[2],"!";}
     else { "// ** Warning: result is exact up to order",fehler-1,"in",
        maxideal(1)[1],"!";}}
   return(ideal(z(1),z(2)));}
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y,t),ds;
 poly f=x3+2xy2+y2;
 list hne=develop(f);
 list hne_extended1=develop(f,6);
 list hne_extended2=develop(f,10);
 pause;
                                   //   compare the different matrices ...
 print(hne[1]);
 print(hne_extended1[1]);
 print(hne_extended2[1]);
                                   // ... and the resulting parametrizations:
 param(hne);
 param(hne_extended1);
 param(hne_extended2);
}
///////////////////////////////////////////////////////////////////////////////

proc invariants
"USAGE:  invariants(l) takes the output of develop(f)
        (list l (matrix m, intvec v, int s[,poly g]))
        and computes some invariants:
RETURN:  a list, if l contains a valid HNE:
    - invariants(l)[1]: intvec         (characteristic exponents)
    - invariants(1)[2],[3]: 2 x intvec (puiseux pairs, 1st and 2nd components)
    - invariants(l)[4]: int            (degree of the conductor)
    - invariants(l)[5]: intvec         (sequence of multiplicities)
         an empty list, if an error occurred in the procedure develop
EXAMPLE: example invariants; shows an example
"{
 //-------------------------- 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,multips,multseq,mpuiseux,npuiseux;
 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;
 //------------------- Multiplizitaetensequenz: -------------------------------
 k=1;
 multips=multiplicities(#);
 for (i=1; i<size(v); i++) {
   for (j=1; j<=v[i]; j++) {
     multseq[k]=multips[i];
     k++;
 }}
 multseq[k]=1;
 //------------------------- Rueckgabe ----------------------------------------
 ergebnis=beta,mpuiseux,npuiseux,conductor,multseq;
 return(ergebnis);
}
example
{ "EXAMPLE:"; echo = 2;
 ring exring=0,(x,y),dp;
 list hne=develop(y4+2x3y2+x6+x5y);
 list erg=invariants(hne);
 erg[1];                   // the characteristic exponents
 erg[2],erg[3];            // the puiseux pairs in packed form
 erg[4] / 2;               // the delta-invariant
 erg[5];                   // the sequence of multiplicities
                           // To display the invariants more 'nicely':
 displayInvariants(hne);
}
///////////////////////////////////////////////////////////////////////////////

proc displayInvariants
"USAGE:   displayInvariants(l); takes the output both of
         develop(f) and reddevelop(f)
          ( list l (matrix m, intvec v, int s[,poly g])
            or list of lists in the form l            )
DISPLAY: invariants of the corresponding branch, resp. of all branches,
         in a better readable form
RETURN:  nothing
EXAMPLE: example displayInvariants;  shows an example
"{
 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];
    if (size(ergebnis[1])>1) {
     for (i=1; i<=size(ergebnis[2]); i++) {
       Ausgabe=Ausgabe+"("+string(ergebnis[2][i])+","
       +string(ergebnis[3][i])+")";
    }}
    " puiseux pairs             :",Ausgabe;
    " degree of the conductor   :",ergebnis[4];
    " delta invariant           :",ergebnis[4]/2;
    " generators of semigroup   :",puiseux2generators(ergebnis[2],ergebnis[3]);
    " sequence of multiplicities:",ergebnis[5];
 }}
 //-------------------- Ausgabe aller Zweige ----------------------------------
 else {
  for (j=1; j<=size(#); j++) {
    ergebnis=invariants(#[j]);
    " --- invariants of branch number",j,": ---";
    " characteristic exponents  :",ergebnis[1];
    Ausgabe="";
    if (size(ergebnis[1])>1) {
     for (i=1; i<=size(ergebnis[2]); i++) {
       Ausgabe=Ausgabe+"("+string(ergebnis[2][i])+","
       +string(ergebnis[3][i])+")";
    }}
    " puiseux pairs             :",Ausgabe;
    " degree of the conductor   :",ergebnis[4];
    " delta invariant           :",ergebnis[4]/2;
    " generators of semigroup   :",puiseux2generators(ergebnis[2],ergebnis[3]);
    " sequence of multiplicities:",ergebnis[5];
    "";
  }
  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=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;
  }
 }
 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) takes the output of develop(f)
         (list l (matrix m, intvec v, int s[,poly g]))
RETURN:  intvec of the different multiplicities that occur during the succes-
         sive blowing up of the curve corresponding to f
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;
  ring r=0,(x,y),dp;
  multiplicities(develop(x5+y7));
  // The first value is the multiplicity of the curve itself, here it's 5
}
///////////////////////////////////////////////////////////////////////////////

proc puiseux2generators (intvec m, intvec n)
"USAGE:   puiseux2generators (m,n);
         m,n intvecs with 1st resp. 2nd components of puiseux pairs
RETURN:  intvec of the generators of the semigroup of values
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 generators
"USAGE:   generators(l); takes the output of develop(f)
         (list l (matrix m, intvec v, int s[,poly g]))
RETURN:  intvec of the generators of the semigroup of values
EXAMPLE: example generators;  shows an example
"
{
 list inv=invariants(#);
 return(puiseux2generators(inv[2..3]));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  poly f=x15-21x14+8x13y-6x13-16x12y+20x11y2-1x12+8x11y-36x10y2+24x9y3+4x9y2
         -16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-1y8;
  list hne=develop(f);
  generators(hne);
}
///////////////////////////////////////////////////////////////////////////////

proc intersection (list hn1, list hn2)
"USAGE:   intersection(hne1,hne2);
         hne1,hne2 two lists representing a HNE (normally two entries out of
         the output of reddevelop), i.e. list(matrix,intvec,int[,poly])
RETURN:  intersection multiplicity of the branches corresponding to hne1 & hne2
         (of type int)
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;
 //h1[1],h2[1];size(h1),size(h2);
   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.)
       }
 // print(ma1);print(ma2);
       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
{
  if (nameof(basering)=="HNEring") {
   def rettering=HNEring;
   kill HNEring;
  }
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  list hne=reddevelop((x2-y3)*(x2+y3));
  intersection(hne[1],hne[2]);
  kill HNEring,r;
  echo = 0;
  if (defined(rettering)) { // wenn aktueller Ring vor Aufruf von example
   setring rettering;       // HNEring war, muss dieser erst wieder restauriert
   def HNEring=rettering;   // werden
   export HNEring;
  }
}
///////////////////////////////////////////////////////////////////////////////

proc displayHNE(list ldev,list #)
"USAGE:   displayHNE(ldev,[,n]); ldev=list
         (output list of develop(f) or reddevelop(f)), n=int
RETURN:  - if only one argument is given, no return value, but
         display an ideal HNE of the following form:
            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+......
     where x,y are the indeterminates of the basering. The values of [] are
     the coefficients of the Hamburger-Noether matrix, the values of <> are
     represented in the HN-matrix as 'x'
         the 1st line (HNE[1]) means that y==[]*z(0)^1+... ,
         the 2nd line (HNE[2]) means that x==[]*z(1)^1+...
     so you can see which indeterminate corresponds to which line
     (it's also possible that x corresponds to the 1st line and y to the 2nd)
         - if a second argument is given, create and export a new ring with
         name 'displayring' containing an ideal 'HNE' as described above

     If ldev contains the output of reddevelop(f), displayHNE shows the HNE's
     of all branches of f in the form described above.
     The optional parameter is then ignored.
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 in develop, so there is no HNE.";
   return(ideal(0));
 }
 if (parstr(basering)!="") {
   if (charstr(basering)!=string(char(basering))+","+parstr(basering)) {
     "Sorry -- not implemented for rings of type (p^k,a),...";
     return(ideal(0));
   }
 }
 def altring=basering;
 ring dazu=char(altring),z(0..size(v)-1),ls;
 def displayring=dazu+altring;
 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)
"USAGE:   newtonpoly(f);   f poly
RETURN:  list of intvec(x,y) of coordinates of the newtonpolygon of f
NOTE: There are many assumptions: (to improve speed)
      - The basering must be ...,(x,y),ls
      - f(x,0) != 0 != f(0,y), f(0,0) = 0
EXAMPLE: example newtonpoly;  shows an example
"
{
 intvec A=(0,ord(subst(f,x,0)));
 intvec B=(ord(subst(f,y,0)),0);
 intvec C,D,H; list L;
 int abbruch,i;
 poly hilf;

 L[1]=A;
 //-------- wirf alle Monome auf oder oberhalb der Geraden AB raus: -----------
 f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
 map xytausch=basering,y,x;
 for (i=2; f!=0; i++) {
   abbruch=0;
   while (abbruch==0) {
 // finde den Punkt aus {verbliebene Pkte (a,b) mit a minimal} mit b minimal: -

     C=leadexp(f);           // Ordnung ls ist wesentlich!

     if (jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
       { abbruch=1; }        // keine Monome unterhalb der Geraden AC

 // ----- alle Monome auf der Parallelen zur y-Achse durch C wegwerfen: -------
 // ------------------ (links von C gibt es sowieso keine mehr) ---------------
     else { f=jet(f,-C[1]-1,intvec(-1,0)); }
   }
 //- finde alle Monome auf der Geraden durch A und C (unterhalb gibt's keine) -
   hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));

   H=leadexp(xytausch(hilf));
   D=H[2],H[1];

 // die Alternative waere ein Ringwechsel nach ..,(y,x),ds gewesen
 // D ist der naechste Eckpunkt (unterster Punkt auf Geraden durch A,C)

   L[i]=D;
   A=D;
 //----------------- alle Monome auf oder unterhalb DB raus -------------------
   f=jet(f,D[2]*B[1]-1,intvec(D[2],B[1]-D[1]));
 }
 L[i]=B;
 return(L);
}
example
{ "EXAMPLE:";
  ring @exring_Newt=0,(x,y),ls;
  export @exring_Newt;
  "  ring exring=0,(x,y),ls;";
  echo = 2;
  poly f=x5+2x3y-x2y2+3xy5+y6-y7;
  newtonpoly(f);
  echo = 0;
  kill @exring_Newt;
}
///////////////////////////////////////////////////////////////////////////////

proc charPoly(poly f, int M, int N)
"USAGE:  charPoly(f,M,N);  f poly in x,y,  M,N int: length and height
                          of newtonpolygon 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,x,1);
 for(charp=0; f<>0; f=f-lead(f))
  { charp=charp+leadcoef(f)*y^(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;
"
{
 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);
}
///////////////////////////////////////////////////////////////////////////////
proc set_list
"USAGE: set_list(L,v,l); L list, v intvec, l element of L (arbitrary type)
RETURN: L with L[v[1]][v[2]][...][v[n]] set to l
WARNING: Make sure there is no conflict of types!
         size(v) must be bigger than 1! (i.e. v must NOT be simply an integer)
EXAMPLE:                 L=set_list(L,intvec(a,b,c),q);
         does exactly what you would expect when typing
                         L[a][b][c]=q;
         (but the last command will only produce an error message)
"
{
 list L=#[1];
 intvec v=#[2];
 def ersetze=#[3];
 //---------------- Bilde den Verzweigungsbaum nach ---------------------------
 def hilf(1)=L[v[1]];
 for (int tiefe=2; tiefe<size(v); tiefe++) {
   def hilf(tiefe)=hilf(tiefe-1)[v[tiefe]];
 }
 //---------------- Fuege Aenderung in die Liste ein --------------------------
 hilf(size(v)-1)[v[size(v)]]=ersetze;
 for (tiefe=size(v)-1; tiefe>1; tiefe--) {
   hilf(tiefe-1)[v[tiefe]]=hilf(tiefe);
 }
 L[v[1]]=hilf(1);
 return(L);
}
///////////////////////////////////////////////////////////////////////////////
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): two different points on a line in the newtonpolygon
         of f (e.g. the starting and the ending point)
RETURN:  poly g: all monomials of f which lay on that line
NOTE:    It is not checked whether S and E really belong to the newtonpolygon!
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 (matrix m, intvec v, int s, poly g), int n
         takes the output l of develop(f) ( or extdevelop(L,n),
         or one entry of the output of reddevelop(f))  and
RETURN:  an extension of the Hamburger-Noether development of f in the form
         of develop, i.e. a list L (matrix m,intvec v,...)
         The new HN-matrix will have at least n columns
         (if the HNE isn't finite).
         Thus if f is irreducible, extdevelop(develop(f),n);
         (in most cases) will produce the same result as develop(f,n).
         Type  help develop;  for more details
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,lastside;
 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); // Ringcharakteristik
  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;
  poly f=hole(f);
  ideal a=hole(lastrow);
  //----- map hat sich als sehr zeitaufwendig bei f gross herausgestellt ------
  //-------- (fetch,imap nur 5% schneller), daher nach Moeglichkeit -----------
  //-----------------Beibehaltung des alten Ringes: ---------------------------
 }
 else { ideal a=lastrow; }
 list Newton;
 number delta;
 //--------------------- Fortsetzung der HNE ----------------------------------
 while (Q<Exaktheit) {
  Newton=newtonpoly(f);
  lastside=size(Newton)-1;
  if (Newton[lastside][2]!=1) {
    " *** The transformed polynomial was not valid!!";}
  Q=Newton[lastside+1][1]-Newton[lastside][1];      // Laenge der letzten Seite

 //--- quasihomogenes Leitpolynom ist c*x^a*y+d*x^(a+Q) => delta=-d/c: --------
  delta=-koeff(f,Newton[lastside+1][1],0)/koeff(f,Newton[lastside][1],1);
  a[Q]=delta;
  "a("+string(zeile-1)+","+string(Q)+") =",delta;
  if (Q<Exaktheit) {
   f=T1_Transform(f,delta,Q);
   if (defined(Protokoll)) { "transformed polynomial:",f; }
   if (subst(f,y,0)==0) {
     "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);
  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 (nameof(basering)=="HNEring") {
   def rettering=HNEring;
   kill HNEring;
  }
  "EXAMPLE:"; echo = 2;
  ring exring=0,(x,y),dp;
  list hne=reddevelop(x14-3y2x11-y3x10-y2x9+3y4x8+y5x7+3y4x6+x5*(-y6+y5)-3y6x3-y7x2+y8);
  print(hne[1][1]); // finite HNE
  print(extdevelop(hne[1],5)[1]); pause;
  print(hne[2][1]); // HNE that can be extended
  list ehne=extdevelop(hne[2],5);
  print(ehne[1]); // new HN-matrix has 5 columns
  param(hne[2]);
  param(ehne);
  param(extdevelop(ehne,7)); pause;
  //////////////////////////
  print(develop(x-y2-2y3-3y4)[1]);
  print(extdevelop(develop(x-y2-2y3-3y4),4)[1]);
  print(extdevelop(develop(x-y2-2y3-3y4),10)[1]);
  kill HNEring,exring;
  echo = 0;
  if (defined(rettering)) { // wenn aktueller Ring vor Aufruf von example
   setring rettering;       // HNEring war, muss dieser erst wieder restauriert
   def HNEring=rettering;   // werden
   export HNEring;
  }
}
///////////////////////////////////////////////////////////////////////////////

proc stripHNE (list l)
"USAGE:   stripHNE(l);    takes the output both of develop(f) and reddevelop(f)
           ( list l (matrix m, intvec v, int s[,poly g])
             or list of lists in the form l            )
RETURN:  a list in the same format as l, but all polynomials g 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 consumes much memory. It is still possible to compute invariants,
         parametrizations etc. with the stripped HNE(s), but it is not possible
         to use extdevelop with them.
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);
}
///////////////////////////////////////////////////////////////////////////////

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 reddevelop
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(Protokoll)) {"To process:";HNEaktu;}
  if (size(HNEs[hnezaehler])!=size(HNEs[hnezaehler][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 delta;
 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) { delta = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
 else {
   if (e%p != 0) { delta = 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(Protokoll)) {e," -> ",eps,"*",p,"^",l;}
     delta = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);

     if ((charstr(basering) != string(p)) and (delta != 0)) {
 //------ coefficient field is not Z/pZ => (p^l)th root is not identity -------
       delta=0;
       if (defined(Protokoll)) {
         "trivial factorization not implemented for",
         "parameters---I've to use 'factorize'";
       }
     }
   }
 }
 if (defined(Protokoll)) {"quasihomogeneous leading form:",f," = ",c,
        "* (y^"+string(N/ e),"-",delta,"* x^"+string(M/ e)+")^",e," ?";}
 if (f != c*(y^(N/ e) - delta*x^(M/ e))^e) {return(0);}
 else {return(delta);}
}
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);
}
///////////////////////////////////////////////////////////////////////////////


proc reddevelop (poly f)
"USAGE:   reddevelop(f); f poly
RETURN:  Hamburger-Noether development of f :
         A list of lists in the form of develop(f); each entry contains the
         data for one of the branches of f.
         For more details type 'help develop;'
NOTE:    as the Hamburger-Noether development of a reducible curve normally
         exists only in a ring extension, reddevelop automatically changes
         the basering to  ring HNEring=...,(x,y),ls; !
EXAMPLE: example reddevelop;  shows an example
"
{
 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) {
   " Warning: all but the first two variables are 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(Protokoll))
 {"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) {
   "The given polynomial is a unit in the power series ring!";
   keepring HNEring;
   return(list());                   // there are no HNEs
 }
 if (f==0) {
   "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) {
   "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];                     // : liest 1 Zeichen von der Tastatur
   setring HNEring;
   map polyhinueber=zweitring,x,y;
   if (str=="r") {
     poly test_sqr=squarefree(f);
     if (test_sqr != f) {
      "The given polynomial is in fact not squarefree.";
      "I'll continue with the radical.";
      f=test_sqr;
     }
     else {"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") {
         if(system("with","Namespaces")) { kill Top::HNEring; }
         setring altring;kill HNEring;kill zweitring;
         return(list());}
       else { "if the algorithm doesn't terminate, you were wrong...";}
   }}
   kill zweitring;
   if (defined(Protokoll)) {"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 (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)
   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(Protokoll)) {"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);
 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);
 ideal deltais;
 poly delta;                   // nicht number, weil delta von a abhaengen kann
 export f,azeilen,HNEs,aneu,deltais,delta;
 //----- hier steht die Anzahl bisher benoetigter Ringerweiterungen drin: -----
 int EXTHNEnumber=0; export EXTHNEnumber;
 setring HNEring;

 // -------- Berechne HNE von allen Zweigen, fuer die x transvers ist: --------
 if (defined(Protokoll))
   {"1st step: Treat newton polygon until height",grenze1;}
 if (grenze1>0) {
  hilflist=HN(f,grenze1,1);
  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(Protokoll)) {" 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(Protokoll))
    {"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);
  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(Protokoll)) {" 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;
  if(system("with","Namespaces")) { kill Primitiv::EXTHNEring(1..EXTHNEnumber);}
  kill EXTHNEring(1..EXTHNEnumber);
 }
 if(system("with","Namespaces")){kill Top::HNE_noparam;}
 kill HNE_noparam;
 kill EXTHNEnumber;
 keepring basering;
 "  ~~~ result:",size(Ergebnis)," branch(es) successfully computed ~~~";
 return(Ergebnis);
}
example
{
  if (nameof(basering)=="HNEring") {
   def rettering=HNEring;
   kill HNEring;
  }
  "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  list hne=reddevelop(x4-y6);
  size(hne);           // i.e. x4-y6 has two branches
  print(hne[1][1]);    // HN-matrix of 1st branch
  param(hne[1]);
  param(hne[2]);
  displayInvariants(hne);
  kill HNEring,r;
  pause;
  // a more interesting example:
  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);
  print(hne[1][1]);
  print(hne[4][1]);
  // a ring change was necessary, a is a parameter
  HNEring;
  kill HNEring,r;
  echo = 0;
  if (defined(rettering)) { // wenn aktueller Ring vor Aufruf von example
   setring rettering;       // HNEring war, muss dieser erst wieder restauriert
   def HNEring=rettering;   // werden
   export HNEring;
  }
}
///////////////////////////////////////////////////////////////////////////////
static
proc HN (poly f,int grenze, int Aufruf_Ebene)
"NOTE: This procedure is only for internal use, it is called via reddevelop
"
{
 //---------- Variablendefinitionen fuer den unverzweigten Teil: --------------
 if (defined(Protokoll)) {"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 delta;
 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);
  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(Protokoll)) {"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);
    delta=factorfirst(redleit(f,Newton[zeiger],Newton[zeiger+1])
                      /x^Newton[zeiger][1],M,N);
    if (delta==0) {
      if (defined(Protokoll)) {" 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,delta,M/ e); }
        else     { ende=1; }
        if (defined(Protokoll)) {"a("+string(zeile)+","+string(Q)+") =",delta;}
        azeilen=set_list(azeilen,intvec(zeile+1,Q),delta);
      }
      else {
 //------------- R > 0 : transformiere f mit T2 -------------------------------
        erg=T2_Transform(f,delta,M,N);
        f=erg[1];delta=erg[2];
 //------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ---------
        while (R!=0) {
         if (defined(Protokoll)) { "h("+string(zeile)+") =",Q; }
         hqs[zeile+1]=Q;  //denn zeile beginnt mit dem Wert 0
 //------------------ markiere das Zeilenende der HNE: ------------------------
         azeilen=set_list(azeilen,intvec(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(Protokoll)) {"a("+string(zeile)+","+string(Q)+") =",delta;}
        azeilen=set_list(azeilen,intvec(zeile+1,Q),delta);
      }
      if (defined(Protokoll)) {"transformed polynomial: ",f;}
      grenze=e;
 //----------------------- teste Abbruchbedingungen: --------------------------
      if (subst(f,y,0)==0) {              // <==> y|f
        "finite HNE of one branch found";
        azeilen=set_list(azeilen,intvec(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) {
          "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 faktoren=ideal(0); list aneu=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(Protokoll)) {
    "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(Protokoll)) { "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;

 //-- 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]);
    delta=factorfirst(leitf,M,N);
    needext=0;
    if (delta==0) {

 //---------- Sonderbehandlung: faktorisere einige Polynome ueber Q(a): -------
     if (charstr(basering)=="0,a") {
       faktoren,flag=extrafactor(leitf,M,N);  // damit funktion. Bsp. Baladi 5
       if (flag==0)
         { " ** I cannot proceed, will produce some error messages...";}
     }
     else
      {
 //------------------ faktorisiere das charakt. Polynom: ----------------------
      faktoren=factorize(charPoly(leitf,M,N));
     }

     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(Protokoll)) {"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 {
          " Change of basering necessary!!";
          if (defined(Protokoll)) { teiler,"is not properly factored!"; }
          if (needext==0) { poly zerlege=teiler; }
          needext=1;
        }
      }
     }                             // end(for j)
    }
    else { deltais=ideal(delta); eis=e;}
    if (defined(Protokoll)) {"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));
         poly transf=translist[1]; poly transfproc=translist[2];
        }
        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 faktoren; list aneu=ideal(0);
      ideal deltais;
      number delta;
      setring HNE_noparam;
      if (defined(letztring)) { kill letztring; }
      if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
      else                  { def letztring=HNEring; }
      f=imap(letztring,f);
      setring EXTHNEring(EXTHNEnumber);
      map hole=HNE_noparam,transf,x,y;
      poly f=hole(f);
    }
   }    // end (while needext==1) bzw. for (numberofRingchanges)

   if (eis==0) { i++; continue; }
   if (ringischanged==1) {
    list erg,hilflist;              // dienen nur zum Sp. von Zwi.erg.
    ideal hilfid;
    erg=ideal(0); hilflist=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(Protokoll)) {
        "a("+string(zeile)+","+string(Q)+") =",deltais[j];
        "transformed polynomial: ",transformiert;
      }
      if (subst(transformiert,y,0)==0) {
       "finite HNE found";
       hnezaehler++;
 //------------ trage deltais[j],x ein in letzte Zeile, fertig: ---------------
       HNEs[hnezaehler]=azeilen+list(poly(0));
       hilfid=HNEs[hnezaehler][zeile+2]; hilfid[Q]=deltais[j]; hilfid[Q+1]=x;
       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+2),hilfid);
       HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+1),Q);
                                  // aktualisiere Vektor mit den hqs
       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);
      transformiert=erg[1];delta=erg[2];
      if (defined(Protokoll)) {"transformed polynomial: ",transformiert;}
      if (subst(transformiert,y,0)==0) {
       "finite HNE found";
       hnezaehler++;
 //---------------- trage endliche HNE in HNEs ein: ---------------------------
       HNEs[hnezaehler]=azeilen;  // dupliziere den gemeins. Anfang der HNE's
       zl=2;
 //----------------------------------------------------------------------------
 // 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(Protokoll)) { "h("+string(zeile+zl-2)+") =",Q1; }
        HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+zl-1),Q1);
        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl,Q1+1),x);
                                  // markiere das Zeilenende der HNE
        zl=zl+1;
 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),ideal(0));

        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
       }
       if (defined(Protokoll)) {
         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delta;
       }
       hilfid=ideal(0);           // = HNEs[hnezaehler][zeile+zl];
       hilfid[Q1]=delta;hilfid[Q1+1]=x;
       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),hilfid);
       HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+zl-1),Q1);
                                  // aktualisiere Vektor mit hqs
       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl+1),poly(0));
 //-------------------- 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);
     if (HNerg[2]==-1) {          // kein Ringwechsel in HN aufgetreten
       aneu=HNerg[1];  }
     else {
       if (defined(Protokoll))
          {" 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;
       ideal hilfid;
       erg=ideal(0); hilflist=erg; faktoren=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);
       delta=imap(letztring,delta);
       f=imap(letztring,f);

       setring EXTHNEring(EXTHNEnumber);
       list HNEs=hole(HNEs);
       list azeilen=hole(azeilen);
       ideal deltais=hole(deltais);
       number delta=number(hole(delta));
       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++;
       HNEs[hnezaehler]=azeilen; // dupliziere den gemeinsamen Anfang der HNE's
       zl=2;
 //---------------- 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(Protokoll)) { "h("+string(zeile+zl-2)+") =",Q1; }
        HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+zl-1),Q1);
        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl,Q1+1),x);
                                 // Markierung des Zeilenendes der HNE
        zl=zl+1;
 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),ideal(0));
        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
       }
       if (defined(Protokoll)) {
         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delta;
       }
       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl,Q1),delta);

 //--- Daten aus T2_Transform sind eingetragen; haenge Daten von HN an: -------
       hilfid=HNEs[hnezaehler][zeile+zl];
       for (zl1=Q1+1; zl1<=ncols(aneu[zaehler][2]); zl1++) {
        hilfid[zl1]=aneu[zaehler][2][zl1];
       }
       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),hilfid);
 //--- vorher HNEs[.][zeile+zl]<-aneu[.][2], jetzt [zeile+zl+1] <- [3] usw.: --
       for (zl1=3; zl1<=size(aneu[zaehler]); zl1++) {
        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl+zl1-2),
             aneu[zaehler][zl1]);
       }
 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
       hilfvec=HNEs[hnezaehler][1],aneu[zaehler][1];
       HNEs=set_list(HNEs,intvec(hnezaehler,1),hilfvec);
 //----------- weil nicht geht: liste[1]=liste[1],aneu[zaehler][1] ------------
      }                     // 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=HNEs[hnezaehler][zeile+2];
     hilfid[Q]=deltais[j];
     for (zl=Q+1; zl<=ncols(aneu[zaehler][2]); zl++) {
      hilfid[zl]=aneu[zaehler][2][zl];
     }
     hilflist=HNEs[hnezaehler];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 extrafactor (poly leitf, int M, int N)
"USAGE:   factors,flag=extrafactor(f,M,N);
         list factors, int flag,M,N, poly f in x,y
ASSUME:  basering is (0,a),(x,y),ds or ls
         Newtonpolygon of f is one side with height N, length M
RETURN:  for some special f, factors is a list of the factors of
         charPoly(f,M,N), same format as factorize
         (see help charPoly; help factorize)
         In this case, flag!=0 (TRUE). If extrafactor cannot factorize f,
         flag will be 0 (FALSE), factors will be the empty list.
NOTE:    This procedure is designed to make reddevelop able to compute some
         special cases that need a ring extension in char 0.
         It becomes obsolete as soon as factorize works also in such rings.
EXAMPLE: example extrafactor;  shows an example
"
{
 list faktoren;

      if (a2==-1) {
       poly testpol=charPoly(leitf,M,N);
       if (testpol==1+2y2+y4) {
         faktoren=list(ideal(y+a,y-a),intvec(2,2));
         testpol=0;
       }
 //------------ ist Poly von der Form q*i+r*y^n, n in N, q,r in Q ?: ----------
       if ((size(testpol)==2) && (find(string(lead(testpol)),"a")>0)
           && (find(string(testpol-lead(testpol)),"a")==0)) {
        faktoren=list(ideal(testpol),intvec(1));
        testpol=0;
       }
      }
      if (a4==-625) {
        poly testpol=charPoly(leitf,M,N);
        if (testpol==4a2-4y2)
         { faktoren=list(ideal(-4,y+a,y-a),intvec(1,1,1)); testpol=0;}
        if (testpol==-4a2-4y2)
         { faktoren=list(ideal(-4,y+1/25*a3,y-1/25*a3),intvec(1,1,1));
           testpol=0;}
      }
      if (defined(testpol)==0) { poly testpol=1; }
      if (testpol!=0) {
        "factorize not implemented in char (0,a)!";
        "could not factorize:",charPoly(leitf,M,N); pause;
      }
 return(faktoren,(testpol==0)); // Test: faktoren==list() geht leider nicht
}
example
{ "EXAMPLE:"; echo=2;
  ring r=(0,a),(x,y),ds;
  minpoly=a2+1;
  poly f=x4+2x2y2+y4;
  charPoly(f,4,4);
  list factors;
  int flag;
  factors,flag=extrafactor(f,4,4);
  if (flag!=0) {factors;}
}