version="$Id: alexpoly.lib,v 1.14 2006-08-02 10:14:06 Singular Exp $";
category="Singularities";
info="
LIBRARY:  alexpoly.lib   Resolution Graph and Alexander Polynomial
AUTHOR:   Fernando Hernando Carrillo, hernando@agt.uva.es
          Thomas Keilen,   keilen@mathematik.uni-kl.de

OVERVIEW:
 A library for computing the resolution graph of a plane curve singularity f,
 the total multiplicities of the total transforms of the branches of f along
 the exceptional divisors of a minimal good resolution of f, the Alexander
 polynomial of f, and the zeta function of its monodromy operator.

PROCEDURES:
 resolutiongraph(f);        resolution graph f
 totalmultiplicities(f);    resolution graph, total multiplicities and strict multiplicities of f
 alexanderpolynomial(f);    Alexander polynomial of f
 semigroup(f);              calculates generators for the semigroup of f
 multseq2charexp(v);        converts multiplicity sequence to characteristic exponents
 charexp2multseq(v);        converts characteristic exponents to multiplicity sequence
 charexp2generators(v);     converts characteristic exponents to generators of the semigroup
 charexp2inter(c,e);        converts contact matrix and charact. exp. to intersection matrix
 charexp2conductor(v);      converts characteristic exponents to conductor
 charexp2poly(v,a);         calculates a polynomial f with characteristic exponents v
 tau_es2(f);                equisingular Tjurina number of f

KEYWORDS: Hamburger-Noether expansion; Puiseux expansion; curve singularities;
          topological invariants; Alexander polynomial; resolution graph;
          total multiplicities; equisingular Tjurina number
";

///////////////////////////////////////////////////////////////////////////////////////////
LIB "hnoether.lib";
///////////////////////////////////////////////////////////////////////////////////////////

proc resolutiongraph (def INPUT,list #)
"USAGE:  resolutiongraph(INPUT); INPUT poly or list
ASSUME:  INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity,
         or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in
         the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of
         @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing
         the contact matrix and a list of integer vectors with the characteristic exponents
         of the branches of a plane curve singularity, or an integer vector containing the
         characteristic exponents of an irreducible plane curve singularity.
RETURN:  intmat, the incidence matrix of the resolution graph of the plane curve
         defined by INPUT, where the entries on the diagonal are the weights of the
         vertices of the graph and a negative entry corresponds to the strict transform
         of a branch of the curve.
NOTE:    In case the Hamburger-Noether expansion of the curve f is needed
         for other purposes as well it is better to calculate this first
         with the aid of @code{hnexpansion} and use it as input instead of
         the polynomial itself.
         @*
         If you are not sure whether the INPUT polynomial is reduced or not, use
         @code{squarefree(INPUT)} as input instead.
SEE ALSO: develop, hnexpansion, totalmultiplicities, alexanderpolynomial
EXAMPLE: example resolutiongraph;  shows an example
"
{
  return(totalmultiplicities(INPUT,#)[1]);
}
example
{
  "EXAMPLE:";
  echo=2;
  ring r=0,(x,y),ls;
  poly f1=(y2-x3)^2-4x5y-x7;
  poly f2=y2-x3;
  poly f3=y3-x2;
  resolutiongraph(f1*f2*f3);
}

proc totalmultiplicities (def INPUT, list #)
"USAGE:  totalmultiplicities(INPUT); INPUT poly or list
ASSUME:  INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity,
         or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in
         the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of
         @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing
         the contact matrix and a list of integer vectors with the characteristic exponents
         of the branches of a plane curve singularity, or an integer vector containing the
         characteristic exponents of an irreducible plane curve singularity.
RETURN:  list @code{L} of three integer matrices. @code{L[1]} is the incidence matrix of
         the resolution graph of the plane curve defined by INPUT, where the entries on the
         diagonal are the weights of the vertices of the graph and a negative entry corresponds
         to the strict transform of a branch of the curve. @code{L[2]} is an integer matrix,
         which has for each vertex in the graph a row and for each branch of the curve a column.
         The entry in position [i,j] contains the total multiplicity of the j-th branch (i.e. the
         branch with weight -j in @code{L[1]}) along the exceptional divisor corresponding
         to the i-th row in @code{L[1]}. In particular, the i-th row contains
         the total multiplicities of the branches of the plane curve (defined by INPUT) along
         the exceptional divisor which corresponds to the i-th row in the incidence matrix
         @code{L[1]}. @code{L[3]} is an integer matrix which contains the (strict) multiplicities
         of the branches of the curve along the exceptional divisors in the same way as @code{L[2]}
         contains the total multiplicities.
NOTE:    The total multiplicty of a branch along an exceptional divisor is the multiplicity
         with which this exceptional divisor occurs in the total transform of this branch
         under the resolution corresponding to the resolution graph.
         @*
         In case the Hamburger-Noether expansion of the curve f is needed
         for other purposes as well it is better to calculate this first
         with the aid of @code{hnexpansion} and use it as input instead of
         the polynomial itself.
         @*
         If you are not sure whether the INPUT polynomial is reduced or not, use
         @code{squarefree(INPUT)} as input instead.
SEE ALSO: develop, hnexpansion, alexanderpolynomial, resolutiongraph
EXAMPLE: example totalmultiplicities;  shows an example
"
{
  ///////////////////////////////////////////////////////////////////////////////////////////////
  /// The algorithm described in [JP] DeJong, Pfister, Local Analytic Geometry, Vieweg 2000,
  /// is used for the implementation -- the missing case is included by appropriate sorting.
  ///////////////////////////////////////////////////////////////////////////////////////////////
  /// If # is empty, then an additional sorting of the branches will be made, so that
  /// the branches can be split with split_graph in order to get the graphs of the curves
  /// separated on the first exceptional divisor. Otherwise split_graph should not be applied!
  ///////////////////////////////////////////////////////////////////////////////////////////////
  /// If #[1]="tau", then totalmultiplicities is called for the use in tau_es - thus it returns
  /// the prolonged resolutiongraphs of the branches, their multiplicity sequences and their
  /// contact matrix - sorted appropriately.
  ///////////////////////////////////////////////////////////////////////////////////////////////
  int i,j,s,e,ii;
  intmat helpmati,helpmatii;
  intvec helpvec;
  /////////////////////////////////////////////////////////////////////////////////
  // Decide, which kind of input we have, and define the contact matrix
  /////////////////////////////////////////////////////////////////////////////////
  // Input: a polynomial defining a plane curve singularity.
  //////////////////////////////////////////////////////////////////////////////
  if (typeof(INPUT)=="poly")
  {
    /*
    poly f=squarefree(INPUT);
    if ( deg(f)!=deg(INPUT) )
    {
      dbprint(printlevel-voice+4,"// input polynomial was not reduced");
      dbprint(printlevel-voice+4,"// we continue with its reduction");
    }
    list I@N@V=invariants(f);
    */
    list I@N@V=invariants(INPUT);
  }
  else
  {
    ///////////////////////////////////////////////////////////////////////////////////
    // Input: the vector of characteristic exponents of an irreducible plane curve
    ///////////////////////////////////////////////////////////////////////////////////
    if (typeof(INPUT)=="intvec")
    {
      list charexp;
      charexp[1]=INPUT;
      intmat contact[1][1]=0;
    }
    else
    {
      if (typeof(INPUT)=="list")
      {
        /////////////////////////////////////////////////////////////////////////////////
        // Input: intersection-matrix and characteristic exponents.
        //////////////////////////////////////////////////////////////////////////////
        if (typeof(INPUT[1])=="intmat")
        {
          intmat contact=INPUT[1];
          list charexp=INPUT[2];
        }
        else
        {
          /////////////////////////////////////////////////////////////////////////////////
          // Input: output of hnexpansion or hne in the ring created by hnexpansion
          //////////////////////////////////////////////////////////////////////////////
          if ((typeof(INPUT[1])=="ring") or (typeof(INPUT[1])=="list"))
          {
            list I@N@V=invariants(INPUT);
          }
          else
          {
            ////////////////////////////////////////////////////////////////////////////////////
            // Input: output of reddevelop or extdevelop -- irreducible plane curve singularity
            ////////////////////////////////////////////////////////////////////////////////////
            if (typeof(INPUT[1])=="matrix")
            {
              list charexp=invariants(INPUT)[1];
              intmat contact[1][1]=0;
            }
            else
            {
              ERROR("The input is invalid!");
            }
          }
        }
      }
      else
      {
        ERROR("The input is invalid!");
      }
    }
  }
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  // If the input was a poly or a HN-Expansion, then calculate the contact matrix and char.exponents.
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  if (defined(I@N@V))
  {
    intmat contact=I@N@V[size(I@N@V)][1];   // contact numbers
    list charexp;                       // characteristic exponents
    for (i=1;i<=size(I@N@V)-1;i++)
    {
      charexp[i]=I@N@V[i][1];
    }
  }
  //////////////////////////////////////////////////////////////////////////////
  // Find the maximal contact numbers of the branches
  //////////////////////////////////////////////////////////////////////////////
  int r=ncols(contact);   // number of branches
  intvec maxcontact;      // maximal contactnumber of branch i with other branches
  for (i=1;i<=r;i++)
  {
    maxcontact[i]=max_in_intvec(intvec(contact[i,1..r]));
  }
  ///////////////////////////////////////////////////////////////////////////////
  // Define the graphs of the branches and prolong them if necessary.
  ///////////////////////////////////////////////////////////////////////////////
  intmat gr,gr_red,grp;      // a subgraph, a reduced subgraph, a prolonged subgraph
  int omega;              // point of highest weight in subgraph
  list graphs;            // contains the subgraphs of the C_i
  list totmult;           // contains the total multiplicities of the subgraphs
  list multipl;           // contains the multiplicities of the subgraphs
  list gr_tm;             // takes a single subgraph and tot.mult.
  intvec tm,mt;           // total multiplicities and multiplicities
  for (i=1;i<=r;i++)
  {
    gr_tm=irred_resgraph_totmult(charexp[i]); // graph, total multipl. and multipl. of the ith branch
    gr=gr_tm[1];  // the graph of the ith branch
    tm=gr_tm[2];  // the total multiplicities of the ith branch
    mt=gr_tm[3];  // the multiplicities of the ith branch
    omega=nrows(gr)-1;  // maximal weight in the graph of the ith branch
    // If the maximal contact of the ith branch with some other branch is larger
    // than the maximal weight omega, then the graph has to be polonged.
    if (omega<maxcontact[i])
    {
      grp=intmat(intvec(0),maxcontact[i]+1,maxcontact[i]+1);
      // save the graph without the point of the strict transform
      if (omega>=1) // otherwise the graph consists only of the strict transform
      {
        grp[1..omega,1..omega]=gr[1..omega,1..omega];
      }
      // add the points of multiplicity 1 up to weight maxcontact[i] and the strict transform
      // and connect them
      for (j=omega+1;j<=maxcontact[i]+1;j++)
      {
        // adding the vertex to the graph and adding the total multiplicities
        if (j<=maxcontact[i])
        {
          grp[j,j]=j;
          if (j>1)
          {
            tm[j]=tm[j-1]+1;
            mt[j]=1;
          }
          else  // then there is no previous point in the graph
          {
            tm[j]=1;
            mt[j]=1;
          }
        }
        // connecting the vertex with the previous part of the graph
        if (j>1)  // otherwise there is nothing to connect to
        {
          grp[j-1,j]=1;
          grp[j,j-1]=1;
        }
      }
      gr=grp;  // replace the subgraph by the prolonged subgraph
    }
    gr[nrows(gr),ncols(gr)]=-i;  // give the strict transform in the ith branch weight -i
    graphs[i]=gr;
    totmult[i]=tm;
    multipl[i]=mt;
  }
  ///////////////////////////////////////////////////////////////////////////////
  // Check, if the procedure is called from inside tau_es.
  int check_tau;
  if (size(#)==1)
  {
    if (#[1]=="tau")
    {
      check_tau=1;
    }
  }
  /////////////////////////////////////////////////////////////////////////////////
  // The procedure tau_es expects the branches to be ordered according to their
  // contact during the resolution process.
  /////////////////////////////////////////////////////////////////////////////////
  if ((size(#)==0) or (check_tau==1))
  {
    list SORT;
    for (i=1;i<=ncols(contact)-2;i++)
    {
      SORT=sort_branches(contact,graphs,totmult,multipl,i,ncols(contact));
      contact=SORT[1];
      graphs=SORT[2];
      totmult=SORT[3];
      multipl=SORT[4];
    }
  }
  ///////////////////////////////////////////////////////////////////////////////////
  /// If the procedure is called from inside tau_es, the output will be the prolonged
  /// resolutiongraphs of the branches, their multiplicity sequences and their contact matrix.
  if (check_tau==1)
  {
    return(list(graphs,multipl,contact));
  }
  /////////////////////////////////////////////////////////////////////////////////////
  // Sort the branches in such a way, that in the blowing up procedure whenever to branches
  // C_i and C_j have contact k (i.e. after k steps they are separated) and one of the following
  // situations occurs:
  // A) some exceptional divisor E_l (l<k) still intersects C_i after k steps of blowing up
  //    and no such E_l intersects C_j
  // OR
  // B) some exceptional divisor E_l (l<k-1) still intersects C_i after k steps of blowing up
  //    and another one (necessarily E_k-1) intersects C_j,
  // THEN: i<j!
  /////////////////////////////////////////////////////////////////////////////////////////
  // This means in the algorithm for glueing graphs together described in [JP] p.217-19
  // the Case 3 never occurs -- and Case 1 only occurs, if it is unavoidable, that is when
  // C_1 there meets an E_l with l<k.
  // (Otherwise we say that (C_i,C_j) has "bad contact".)
  ////////////////////////////////////////////////////////////////////////////////////////
  // We can detect this by looking at the graphs of C_i and C_j. If E_l stays together with
  // C_i and not with C_j in the k-th step of blowing up, then in the graph of C_i
  // k and l are NOT connected, while in the graph of C_j they are!
  //////////////////////////////////////////////////////////////////////////////////////////
  // Note also, that for a fixed pair (i,j) this situation can only occur for ONE l (in Case A: l<k
  // and possibly l=k-1;  in Case B: l<k-1).
  //////////////////////////////////////////////////////////////////////////////////////////
  // Our algorithm now works as follows. Start with i=1 and j=2.
  // While (i<r=number of branches) do as follows: Check C_i with C_j.
  // If a bad contact (C_i,C_j) occurs, then MOVE C_j to the position of C_i (and move all C_t
  // with t>=i, t!=j  one position further). Else do nothing particular.
  // Then, if j<r set j=j+1. Else set j=i+2 and i=i+1.  End of While.
  // To see that this works we note that if (C_i,C_j) has bad contact for some j>i, then
  // (C_j,C_t) does NOT have bad contact for t=i+1,...,j-1. For this we consider the 3 cases:
  // 1)  contact(C_t,C_j)=contact(C_i,C_j)=k:  E_l stays with C_j in k-th step and C_t and C_j
  //                                           separate there, so (C_t,C_j) has bad contact
  //                                           and hence (C_j,C_t) does NOT have bad contact
  // 2)  contact(C_t,C_j)>contact(C_i,C_j)=k:  CANNOT HAPPEN - C_i had no bad contact with C_t,
  //                                           hence E_l did not stay with C_t in the k-th step;
  //                                           however, E_l stays with C_j there, so C_t and C_j
  //                                           cannot have a contact higher than k
  // 3)  contact(C_t,C_j)<contact(C_i,C_j)=k:  when C_t and C_j separate, C_j and C_i are still
  //                                           together; and since (C_i,C_t) did not have bad
  //                                           contact also (C_j,C_t) cannot have bad contact.
  // This ensures that when we do the moving of the graphs, we do NOT create any new bad contacts.
  //////////////////////////////////////////////////////////////////////////////
  i=1;
  j=2;
  intvec connections_i,connections_j;
  int n_i,n_j;
  int bad_contact;
  while (i<r)
  {
    // Test graphs[i] with graphs[j] for bad contact.
    connections_i=find_lower_connection_points(graphs[i],contact[i,j]);
    if (connections_i[1]==0) // no connection point of lower weight there
    {
      n_i=0;
    }
    else
    {
      n_i=size(connections_i);
    }
    connections_j=find_lower_connection_points(graphs[j],contact[i,j]);
    if (connections_j[1]==0) // no connection point of lower weight there
    {
      n_j=0;
    }
    else
    {
      n_j=size(connections_j);
    }
    bad_contact=0;
    ii=1;
    // The points of weight k=contact(C_i,C_j) in C_i and C_j are connectd to n_i respectively
    // n_j points of weight less than k, where n_i and n_j might take values among 0, 1 and 2.
    // n_j=2:        Then E_k is intersected by E_k-1 and E_l (l<k-1) in the k-th step, however
    //               C_j does not stay with anyone of those, so (C_i,C_j) does not have bad contact.
    // n_i=0:        Then C_i stays with E_k-1 and there is no E_l with l<k-1 in the k-th step,
    //               hence (C_i,C_j) does not have bad contact.
    // n_i=1, n_j=0: Analogously, then C_j stays with E_k-1 and there is no E_l (l<k-1). This means
    //               (C_i,C_j) has bad contact.
    // n_i=1, n_j=1: EITHER there is no E_l (l<k-1) and C_i and C_j both separate from E_k-1
    //               in the k-th step of blowing up (hence (C_i,C_j) does not have bad contact).
    //               OR there is an E_l (l<k-1) with which one stays and the other one stays
    //               with E_k-1 (bad contact, if E_l stays with C_j; good contact otherwise).
    // n_i=2,        E_k is intersected by E_k-1 and E_l in the k-th step and C_i does not
    //               stay together with any of those.
    //        n_j=0: This CANNOT occur, since then C_j would have to stay together with E_l
    //               and E_k-1, which is impossible.
    //        n_j=1: Then C_j stays together with either E_l or with E_k-1, however both cases
    //               imply that (C_i,C_j) has bad contact
    // ALTERNATIVE CONSIDERATION FOR DISTINGUISHING THE CASES:
    // n_i<n_j   : C_i stays with "more" E's than C_j does (only one is possible of course!),
    //             hence (C_i,C_j) does not have bad contact.
    // n_i>n_j   : C_i stays with "less" E's than C_j does (i.e. with none) hence (C_i,C_j) has bad contact.
    // n_i=n_j=0 : Not possible, since then both would stay with E_k-1.
    // n_i=n_j=2 : Both separate from E_k-1 as well as from some E_l (l<k) - hence NO bad contact.
    // n_i=n_j=1 : One stays with E_k-1 and one with E_l (l<k). Bad contact, if C_j stays with E_l,
    //             i.e. if in the graph of C_i the point k is connected to l, and no bad contact otherwise.
    if ((n_i>n_j) or ((n_i==1) and (n_j==1) and (connections_i[1]<contact[i,j]-1)))
    {
      bad_contact=1;
      // Move the graph (etc.) of C_j into the position i.
      graphs=insert(graphs,graphs[j],i-1);
      graphs=delete(graphs,j+1);
      totmult=insert(totmult,totmult[j],i-1);
      totmult=delete(totmult,j+1);
      multipl=insert(multipl,multipl[j],i-1);
      multipl=delete(multipl,j+1);
      contact=move_row_col(contact,i,j);
    }
    if (j<r)  // There are still some C_j against which C_i has to be tested.
    {
      j++;
    }
    else      // Now C_i has been tested against all C_j, and we may continue with C_i+1.
    {
      j=i+2;
      i=i+1;
    }
  }
  /////////////////////////////////////////////////////////////////////////////////////
  // Glue the graphs together.
  //////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////////
  intmat rgraph=graphs[1];                 // keeps the resolution graph
  intmat rtm=intmat(totmult[1],nrows(rgraph),1); // keeps the tot.mult. at the vertices of the graph as rows
  intmat rmt=intmat(multipl[1],nrows(rgraph),1); // keeps the mult. at the vertices of the graph as rows
  intvec stricttransforms=0,ncols(rgraph);   // keeps the position of the ith strict
                                             // transform in rgraph at position i+1 !!!
  intvec k,kp,p,q,o;   // highest contact numbers, num. of branch with hgt contact, separation points
  int maxc,maxcp;
  for (i=2;i<=r;i++)
  {
    //////////////////////////////////////////////////////////////////////////////////
    // Find j<i minimal s.t. contact[i,j] is maximal.
    maxcp=i;
    for (j=1;j<i;j++)
    {
      if (contact[i,j]>contact[i,maxcp]){maxcp=j;}
    }
    kp[i]=maxcp;            // the j such that C_i has its maximal contact to C_j
    k[i]=contact[i,maxcp];  // the maximal contact of C_i with C_1,...,C_i-1
    ///////////////////////////////////////////////////////////////////////////////////
    // Find in the graph of C_1,...,C_i-1 the points p of wgt k and q of wgt k-1
    // connected to C_maxcp.
    // Since non of C_j for j<maxcp has contact k with C_i, the point p lies in
    // the remaining part of the graph of C_maxcp.
    s=rgraph[stricttransforms[maxcp]+1,stricttransforms[maxcp]+1];
    p[i]=stricttransforms[maxcp]+1+k[i]-s; // pt. to which reduced subgraph of C_i is glued
    // If s<k[i], then q also lies in this part, otherwise it lies in the remaining part
    // of the graph of the C_j to which C_maxcp is connected, i.e. j=kp[maxcp], since
    // the contact of C_i and of C_maxcp to this C_j is strictly less than k.
    // If s=k[i]=1, then p=1 and there is no q! We may thus set q=0.
    if ((s<k[i]) or ((s==k[i]) and (s==1)))  // i.e. q is on the same subgraph as p, or q does not exist
    {
      q[i]=p[i]-1;
    }
    else               // i.e. q is on the subgraph to which the subgraph of p has been glued
    {
      s=rgraph[stricttransforms[kp[maxcp]]+1,stricttransforms[kp[maxcp]]+1];
      q[i]=stricttransforms[kp[maxcp]]+k[i]-s;
    }
    //////////////////////////////////////////////////////////////////////////////////////
    // Reduce the graph of C_i and add it to the graph of C_1,...,C_i-1.
    gr=graphs[i];
    s=nrows(rgraph);
    // Delete in the graph of C_i everything of weight <=k.
    gr_red=intmat(intvec(gr[k[i]+1..nrows(gr),k[i]+1..ncols(gr)]),nrows(gr)-k[i],ncols(gr)-k[i]);
    // Add this part to the graph of C_1,...,C_i-1.
    rgraph=addmat(rgraph,gr_red);
    /////////////////////////////////////////////////////////////////////////////////////
    // Connect the two parts of the graph.
    /////////////////////////////////////////////////////////////////////////////////////
    // Connect the points connected to the point of wgt k in the graph of C_i to p[i].
    for (j=k[i]+1;j<=ncols(gr);j++)
    {
      if(gr[k[i],j]==1)
      {
        rgraph[s+j-k[i],p[i]]=1;
        rgraph[p[i],s+j-k[i]]=1;
      }
    }
    // If pt. of wgt k is not connected to pt of wgt k-1 in graph of C_i, then points
    // connected to the one of wgt k-1 have to be connected to q[i].
    if (k[i]>1)
    {
      if (gr[k[i],k[i]-1]!=1)
      {
        for (j=k[i]+1;j<=ncols(gr);j++)
        {
          if(gr[k[i]-1,j]==1)
          {
            rgraph[s+j-k[i],q[i]]=1;
            rgraph[q[i],s+j-k[i]]=1;
          }
        }
        // Cut the connection from p[i] to q[i].
        rgraph[p[i],q[i]]=0;
        rgraph[q[i],p[i]]=0;
      }
    }
    stricttransforms[i+1]=ncols(rgraph);
    ////////////////////////////////////////////////////////////////////////////////
    // Adjust the total multiplicities
    ////////////////////////////////////////////////////////////////////////////////
    // Add the total multiplicities for the added subgraph to rtm
    tm=totmult[i];
    mt=multipl[i];
    if (k[i]<size(tm)) // if the reduced subgraph of C_i has more than one point
    {
      rtm=addmat(rtm,intmat(intvec(tm[k[i]+1..size(tm)]),nrows(gr_red),1));
      rmt=addmat(rmt,intmat(intvec(mt[k[i]+1..size(mt)]),nrows(gr_red),1));
    }
    else // if the reduced subgraph of C_i consists only of the strict transform
    {
      rtm=addmat(rtm,0);
      rmt=addmat(rmt,0);
    }
    // Adjust the total multiplicities at the places where the subgraph has been glued.
    e=k[i];    // the highest weight of a point that has not yet been assigned its tot. mult.
    while(e>=1)
    {
      s=stricttransforms[maxcp]+1;  // Line in the graph of the start. pt. of the subgraph of C_maxcp.
      for (j=rgraph[s,s];j<=e;j++)  // Adjust the multiplicities.
      {
        rtm[s+j-rgraph[s,s],i]=tm[j];
        rmt[s+j-rgraph[s,s],i]=mt[j];
      }
      maxcp=kp[maxcp];  // To which subgraph has the subgraph of C_maxcp been glued?
      e=rgraph[s,s]-1;  // What is the highest weight of a pt. that has not yet been assigned its tot.mult.?
    }
    e=nrows(rtm);  // Number of rows in the matrix of totalmultiplicities.
    // The total multiplicities of the C_s for s=1,...,i-1 along the exceptional divisors
    // which are introduced after the strict transform of C_s has separated (i.e. the entries
    // in rows stricttransform[i]+1,...,stricttransform[i+1]-1 in the s-th column of the matrix
    // of total multiplicities still have to be calculated.
    for (s=1;s<i;s++)
    {
      rtm[1..e,s]=adjust_tot_mult(intvec(rtm[1..e,i]),intvec(rtm[1..e,s]),intvec(rmt[1..e,i]),intvec(rmt[1..e,s]),p,q,stricttransforms,i);
    }
    // The total multiplicities of the C_i along the exceptional divisors
    // which are introduced for the sake of C_s, s=1,...,i-1, after the strict transform
    // of C_i has separated (i.e. the entries in rows stricttransform[s]+1,...,stricttransform[s+1]-1
    // in the i-th column of the matrix of total multiplicities still have to be calculated.
    for (s=1;s<i;s++)
    {
      rtm[1..e,i]=adjust_tot_mult(intvec(rtm[1..e,s]),intvec(rtm[1..e,i]),intvec(rmt[1..e,s]),intvec(rmt[1..e,i]),p,q,stricttransforms,s);
    }
  }
  list result=rgraph,rtm,rmt;
  return(result);
}
example
{
  "EXAMPLE:";
  echo=2;
  ring r=0,(x,y),ls;
  poly f1=(y2-x3)^2-4x5y-x7;
  poly f2=y2-x3;
  poly f3=y3-x2;
  totalmultiplicities(f1*f2*f3);
}



proc alexanderpolynomial (def INPUT)
"USAGE:  alexanderpolynomial(INPUT); INPUT poly or list
ASSUME:  INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity,
         or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in
         the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of
         @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing
         the contact matrix and a list of integer vectors with the characteristic exponents
         of the branches of a plane curve singularity, or an integer vector containing the
         characteristic exponents of an irreducible plane curve singularity.
CREATE:  a ring with variables t, t(1), ..., t(r) (where r is the number of branches of
         the plane curve singularity f defined by INPUT) and ordering ls over the
         ground field of the basering. @*
         Moreover, the ring contains the Alexander polynomial of f in variables t(1), ..., t(r)
         (@code{alexpoly}), the zeta function of the monodromy operator of f in the variable t
         (@code{zeta_monodromy}), and a list containing the factors of the Alexander
         polynomial with multiplicities (@code{alexfactors}).
RETURN:  a list, say @code{ALEX}, where @code{ALEX[1]} is the created ring
NOTE:    to use the ring type: @code{def ALEXring=ALEX[i]; setring ALEXring;}.
         @*
         Alternatively you may use the procedure sethnering and type: sethnering(ALEX,\"ALEXring\");
         @*
         To access the Alexander polynomial resp. the zeta function resp. the
         factors of the Alexander polynomial type: @code{alexpoly} resp. @code{zeta_monodromy}
         resp. @code{alexfactors}.@*
         In case the Hamburger-Noether expansion of the curve f is needed
         for other purposes as well it is better to calculate this first
         with the aid of @code{hnexpansion} and use it as input instead of
         the polynomial itself.
         @*
         If you are not sure whether the INPUT polynomial is reduced or not, use
         @code{squarefree(INPUT)} as input instead.
SEE ALSO: resolutiongraph, totalmultiplicities
EXAMPLE: example alexanderpolynomial;  shows an example
"
{
  // Get the resolution graph and the total multiplicities.
  list gr_tm=totalmultiplicities(INPUT);
  intmat gr=gr_tm[1];
  intmat tm=gr_tm[2];
  int r=ncols(tm);
  int e=ncols(gr);
  // Define the Ring for the Alexander Polynomial and the Zeta Function of the Monodromy.
  execute("ring ALEXring=("+charstr(basering)+"),(t,t(1..r)),dp;");
  poly hilfspoly=1;
  poly alexnumerator=1;    // numerator of the Alexander polynomial
  poly alexdenominator=1;  // denominator of the Alexander polynomial
  list alexfactors;        // the factors of the Alexanderpolynomial with multiplicities
  list alexfactor;
  int chi=2;
  int i,j,k;
  int s=1;
  // Consider every vertex of the resolution graph.
  for (i=1;i<=e;i++)
  {
    if (gr[i,i]>0)  // If it belongs to an exceptional curve.
    {
      for (j=1;j<=e;j++)  // Calculate the Euler charateristik of the smooth locus of the exc. curve.
      {
        if ((gr[i,j]==1) and (i!=j))
        {
          chi=chi-1;
        }
      }
      if (chi!=0)         // If the Euler characteristik is not zero, then it gives a factor in the AP.
      {
        for (k=1;k<=r;k++)
        {
          hilfspoly=hilfspoly*t(k)^tm[i,k];
        }
        hilfspoly=1-hilfspoly;
        if (chi<0)       // ... either in the numerator ...
        {
          alexnumerator=alexnumerator * hilfspoly^-chi;
        }
        else             // ... or in the denominator.
        {
          alexdenominator=alexdenominator * hilfspoly^chi;
        }
        alexfactor=hilfspoly,-chi;
        alexfactors[s]=alexfactor;
        s++;
      }
      chi=2;
      hilfspoly=1;
    }
  }
  // Calculate the Alexander Polynomial.
  if (ncols(tm)==1)  // If we have only one branch, then the numerator has to be multiplied by 1-t.
  {
    alexnumerator=alexnumerator*(1-t(1));
    alexfactor=1-t(1),1;
    alexfactors[size(alexfactors)+1]=alexfactor;
  }
  poly alexpoly=alexnumerator / alexdenominator;
  // Calculate the Zeta Function of the Monodromy Operator.
  poly zeta_monodromy=alexpoly;
  for (i=1;i<=r;i++)
  {
    zeta_monodromy=subst(zeta_monodromy,t(i),t);
  }
  export zeta_monodromy;
  export alexnumerator;
  export alexdenominator;
  export alexfactors;
  export alexpoly;
  list ALEX=ALEXring;
  return(ALEX);
}
example
{
  "EXAMPLE:";
  echo=2;
  ring r=0,(x,y),ls;
  poly f1=(y2-x3)^2-4x5y-x7;
  poly f2=y2-x3;
  poly f3=y3-x2;
  list ALEX=alexanderpolynomial(f1*f2*f3);
  def ALEXring=ALEX[1];
  setring ALEXring;
  alexfactors;
  alexpoly;
  zeta_monodromy;
}

proc semigroup (def INPUT)
"USAGE:    semigroup(INPUT); INPUT poly or list
ASSUME:   INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity,
          or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in
          the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of
          @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing
          the contact matrix and a list of integer vectors with the characteristic exponents
          of the branches of a plane curve singularity, or an integer vector containing
          the characteristic exponents of an irreducible plane curve singularity.
RETURN:   a list with three entries. The first and the second are lists @code{v_1,...,v_s}
          and @code{w_1,...,w_r} respectively of integer vectors such that the semigroup
          of the plane curve defined by the INPUT is generated by the vectors
          @code{v_1,...,v_s,w_1+k*e_1,...,w_r+k*e_r}, where e_i denotes the i-th standard
          basis vector and k runs through all non-negative integers. The thrid entry is the conductor
          of the plane curve singularity. Note that r is the number of branches of the plane curve
          singularity and integer vectors thus have size r.
NOTE:     If the output is zero this means that the curve has one branch and is regular.
          In the reducible case the set of generators may not be minimal.
          @*
          If you are not sure whether the INPUT polynomial is reduced or not, use
          @code{squarefree(INPUT)} as input instead.
SEE ALSO: resolutiongraph, totalmultiplicities
EXAMPLE:  example semigroup;  shows an example
"
{
  ////////////////////////////////////////////////////////////////////////////////////////
  // Uses the algorithm in [CDG99] A. Campillo, F. Delgado, S.M. Gusein-Zade, On the
  // generators of the semigroup of a plane curve singularity,
  // J. London Math. Soc. (2) 60 (1999), 420-430.
  ////////////////////////////////////////////////////////////////////////////////////////
  intvec conductor; // conductor of the singularity
  list charexp;     // characteristic exponents of the singularity
  int i,j;
  /////////////////////////////////////////////////////////////////////////////////
  // Decide, which kind of input we have, and define the contact matrix
  /////////////////////////////////////////////////////////////////////////////////
  // Input: a polynomial defining a plane curve singularity.
  //////////////////////////////////////////////////////////////////////////////
  if (typeof(INPUT)=="poly")
  {
    /*
    poly FFF=squarefree(INPUT);
    if ( deg(FFF)!=deg(INPUT) )
    {
      dbprint(printlevel-voice+3,"// input polynomial was not reduced");
      dbprint(printlevel-voice+3,"// we continue with its reduction");
    }
    list I@N@V=invariants(FFF);
    */
    list I@N@V=invariants(INPUT);
  }
  else
  {
    ///////////////////////////////////////////////////////////////////////////////////
    // Input: the vector of characteristic exponents of an irreducible plane curve
    ///////////////////////////////////////////////////////////////////////////////////
    if (typeof(INPUT)=="intvec")
    {
      // Calculate the semigroup and the conductor directly.
      conductor[1]=charexp2conductor(INPUT);
      intvec gener=charexp2generators(INPUT);
      list genera;
      for (i=1;i<=size(gener);i++)
      {
        genera[i]=gener[i];
      }
      return(list(genera,list(),conductor));
    }
    else
    {
      if (typeof(INPUT)=="list")
      {
        /////////////////////////////////////////////////////////////////////////////////
        // Input: intersection-matrix and characteristic exponents.
        //////////////////////////////////////////////////////////////////////////////
        if (typeof(INPUT[1])=="intmat")
        {
          intmat contact=INPUT[1];
          charexp=INPUT[2];
          int n=ncols(contact); // to know how many branches we have
          if (n==1) // Only one branch!
          {
            return(semigroup(charexp[1]));
          }
          intmat intersecmult=charexp2inter(contact,charexp);
          for(i=1;i<=ncols(contact);i++)
          {
            conductor[i]=charexp2conductor(charexp[i]);//list with the characteristic exponents
          }
        }
        else
        {
          /////////////////////////////////////////////////////////////////////////////////
          // Input: output of hnexpansion or hne in the ring created by hnexpansion
          //////////////////////////////////////////////////////////////////////////////
          if ((typeof(INPUT[1])=="ring") or (typeof(INPUT[1])=="list"))
          {
            list I@N@V=invariants(INPUT);
          }
          else
          {
            ////////////////////////////////////////////////////////////////////////////////////
            // Input: output of develop or extdevelop -- irreducible plane curve singularity
            ////////////////////////////////////////////////////////////////////////////////////
            if (typeof(INPUT[1])=="matrix")
            {
              return(semigroup(invariants(INPUT)[1]));
            }
            else
            {
              ERROR("The input is invalid!");
            }
          }
        }
      }
      else
      {
        ERROR("The input is invalid!");
      }
    }
  }
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  // If the input was a poly or an HN-Expansion, then calculate the contact matrix and char.exponents.
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  if (defined(I@N@V))
  {
    int n =size(I@N@V)-1;// number of branches
    // If the singularity is irreducible, then we calculate the semigroup directly.
    if (n==1)
    {
      return(semigroup(I@N@V[1][1]));
    }
    // If the singularity is not irreducible, then we go on.
    intmat contact=I@N@V[size(I@N@V)][1];        // contact matrix
    intmat intersecmult=I@N@V[size(I@N@V)][2];   // intersection multiplicities
    for(i=1;i<=n;i++)
        {
          conductor[i]=I@N@V[i][5];
      charexp[i]=I@N@V[i][1];
        }
  }
  /////////////////////////////////////////////////////////////////////////////////////
  // If we have come so far, the curve is reducible!
  /////////////////////////////////////////////////////////////////////////////////////
  // Compute the conductor of the curve.
  /////////////////////////////////////////////////////////////////////////////////////
  for(i=1;i<=size(conductor);i++)
  {
    for(j=1;j<=size(conductor);j++)
        {
          conductor[i]=conductor[i]+intersecmult[i,j];
        }
  }
  /////////////////////////////////////////////////////////////////////////////////////
  /// The total multiplicity vectors of the exceptional divisors generate the semigroup,
  /// however, this system of generators is not minimal. Theorem 2 in [CDG99] leads
  /// to the following algorithm for minimizing the set of generators.
  /////////////////////////////////////////////////////////////////////////////////////
  list resgr_totmult=totalmultiplicities(list(contact,charexp)); // resolution graph and tot. mult.
  intmat resgr=resgr_totmult[1];    // resolution graph
  intmat totmult=resgr_totmult[2];  // total multiplicities
  list conpts;                      // ith entry = points to which i is connected and their weights
  intvec series;                    // ith entry = row in totmult which corresponds to w_i
  intvec deadarc;                   // Star point and the point connected to the star point of a dead arc.
  list deadarcs;                    // Saves deadarc for all dead arcs.
  int stop,ctp,ctpstop;
  list ordinary_generators, series_generators; // the v_i and the w_i from the description
  // Find for each branch of the curve all the points in the graph to which it is connected
  // and save the numbers of the corresponding rows in the graph as well as the weight of the point.
  for (i=1;i<=ncols(resgr);i++)
  {
    conpts[i]=find_connection_points(resgr,i);
  }
  //////////////////////////////////////////////////////////////////////////////////
  // Check for each point in the graph, whether it contributes to the semigroup.
  // If it does not contribute to the semigroup, then set the weight of the point
  // to zero, i.e. resgr[i,i]=0. Note that the vertex 1 always contributes!
  //////////////////////////////////////////////////////////////////////////////////
  // Find first all the points corresponding to the branches and the end points of dead arcs.
  // The points to which the branch k is connected contributes to w_k. The end points of
  // dead arcs contribute to the v_i, while the remaining points of the dead arcs are to be removed.
  j=1;      // Counter for dead arcs - the star points of dead arcs have to be saved for future use.
  for (i=2;i<=ncols(resgr);i++)
  {
    // The end points of dead arcs and the end points corresponding to the branches are
    // recognized by the fact that they are only connected to one other point.
    if (size(conpts[i][1])==1)  // row i in the graph corresponds to  an end point
    {
      // For an end point corresponding to a branch k the last point E_alpha(k), to which it is
      // connected, contributes to the generator w_k.
      if (resgr[i,i]<0)
      {
        series[-resgr[i,i]]=conpts[i][1][1];  // Find E_alpha(k), where k=-resgr[i,i]
        ctp=conpts[i][1][1];
        resgr[ctp,ctp]=0;                     // So E_alph(k) does not contribute to the v_i.
      }
      // For an end point of a dead arc the end point contributes, but all the other points
      // of the dead arc - including the star point = separation pt of the dead arc - do not contribute.
      if (resgr[i,i]>0)
      {
        stop=0;  // Stop checks whether the star point of the dead arc has been reached.
        ctp=conpts[i][1][1]; // The point to which the end point is connected.
        // Set the weights of all the other points in the dead arc to zero.
        while ((stop==0) and (ctp!=1)) // If the star point or the vertex 1 has been reached, stop.
        {
          deadarc[2]=i;          // i might be the point connected to the star point of the dead arc.
          resgr[ctp,ctp]=0;
          if (size(conpts[ctp][1])==2)  // If ctp was an ordinary vertex, we go on.
          {
            deadarc[2]=ctp; // ctp might be the point connectd to the star point of the dead arc.
            ctp=conpts[ctp][1][2];  // This is the second point (of higher weight) to which ctp is connected.
          }
          else  // If ctp is the star point, we stop.
          {
            deadarc[1]=ctp;       // Star point of this dead arc.
            deadarcs[j]=deadarc;  // Save the j-th dead arc.
            j++;
            stop=1;
          }
        }
      }
    }
  }
  //////////////////////////////////////////////////////////////////////////////////
  // Points (!=1) which are on the geodesics of every branch don't contribute.
  // (The geodesic of a branch is the shortest connection of the strict transform to 1 in the graph.)
  stop=0;              // Checks whether the points on all geodesics have been found.
  ctp=1;               // Point (=row in resgr) to be considered.
  int prev_ctp;        // Previous point in the graph, to which ctp was connected.
  int dead_arc_ctp;    // If ctp is the star pt of a dead arc, this is the connection pt of ctp in the d.a.
  int dastarptcheck;   // Checks whether a point is a star point of a dead arc.
  if (size(conpts[1][1])>=2) // The graphs separate already at point 1.
  {
    stop=1;
  }
  else  // The graphs do not separate at point 1.
  {
    prev_ctp=1;
    ctp=conpts[1][1][1];  // Next point to be considered.
  }
  // Pass on along the graph until we reach the first point where some branches separate.
  while (stop==0)
  {
    if (size(conpts[ctp][1])==2)  // Point ctp is connected to 2 other points, hence is a normal vertex.
    {
      resgr[ctp,ctp]=0;       // Point ctp is a normal vertex.
      prev_ctp=ctp;           // Save the position of ctp for future use.
      ctp=conpts[ctp][1][2];  // Next point to which ctp is connected.
    }
    if (size(conpts[ctp][1])>3)  // More than three points are connected to ctp.
    {
      resgr[ctp,ctp]=0;
      stop=1;                    // The graphs separate at point ctp.
    }
    if (size(conpts[ctp][1])==3)  // At ctp a dead arc might depart or some branch(es)!
    {                             // If a dead arc departs, then the branches stay together.
      resgr[ctp,ctp]=0;
      // Check if a dead arc departs at point ctp (i.e. if ctp is the star point of a dead arc),
      // then the branches do not separate at ctp.
      dastarptcheck=0;
      i=1;
      while ((i<=size(deadarcs)) and (dastarptcheck==0))
      {
        if (ctp==deadarcs[i][1])  // ctp is the star point of a dead arc.
        {
          dastarptcheck=1;
          dead_arc_ctp=deadarcs[i][2];  // The point in the dead arc to which ctp is connected.
        }
        i++;
      }
      if (dastarptcheck==0)  // ctp is not the star point of a dead arc, hence the graphs separate at ctp.
      {
        stop=1;
      }
      else
      {
        // Set ctp to the point connected to ctp which is not in the dead arc and is not prev_ctp.
        i=1;
        ctpstop=0;
        while ((i<=3) and (ctpstop==0))
        {
          if ((conpts[ctp][1][i]!=prev_ctp) and (conpts[ctp][1][i]!=dead_arc_ctp))
          {
            prev_ctp=ctp;
            ctp=conpts[ctp][1][i];
            ctpstop=1;
          }
          i++;
        }
      }
    }
  }
  /////////////////////////////////////////////////////////////////////////////////////////////
  // Collect the generators v_j by checking which points in the graph still have
  // a positive weight! These points contribute their total multiplicity vector as
  // generator v_j.
  j=1;
  for (i=1;i<=ncols(resgr);i++)
  {
    if (resgr[i,i]>0)
    {
      ordinary_generators[j]=intvec(totmult[i,1..ncols(totmult)]);
      j++;
    }
  }
  // The "exceptional" generators w_i, for which we have to include w_i+ke_i (for all k)
  // are the total multiplicity vectors of the points saved in series.
  for (i=1;i<=ncols(totmult);i++)
  {
    series_generators[i]=intvec(totmult[series[i],1..ncols(totmult)]);
  }
  return(list(ordinary_generators,series_generators,conductor));
}
example
{
  "EXAMPLE:";
  echo=2;
  ring r=0,(x,y),ls;
  // Irreducible Case
  semigroup((x2-y3)^2-4x5y-x7);
  // In the irreducible case, invariants() also calculates a minimal set of
  // generators of the semigroup.
  invariants((x2-y3)^2-4x5y-x7)[1][2];
  // Reducible Case
  poly f=(y2-x3)*(y2+x3)*(y4-2x3y2-4x5y+x6-x7);
  semigroup(f);
}



proc charexp2generators (intvec charexp)
"USAGE:      charexp2generators(v),  v intvec
ASSUME:      v contains the characteristic exponents of an irreducible plane
             curve singularity
RETURN:      intvec, the minimal set of generators of the semigroup of the plane curve singularity
SEE ALSO:    invariants, resolutiongraph, totalmultiplicities, alexanderpolynomial
KEYWORDS:    generators; semigroup; characteristic exponents; topological invariants
EXAMPLE:     example charexp2generators;   shows an example"
{
  int end=size(charexp);
  // If the singularity is smooth!
  if (end==1)
  {
    return(1);
  }
  int i,j;
  intvec gener;
  intvec GGT;
  for (i=1;i<=end;i++)
  {
    // Calculate the sequence of gcd's of the characteristic exponents.
    if (i==1)
    {
      GGT[1]=charexp[1];
    }
    else
    {
      GGT[i]=gcd(GGT[i-1],charexp[i]);
    }
    // Calculate the generators.
    gener[i]=charexp[i];
    for (j=2;j<=i-1;j++)
    {
      gener[i]=gener[i]+((GGT[j-1]-GGT[j])/GGT[i-1])*charexp[j];
    }
  }
  return(gener);
}
example
{
  "EXAMPLE:";
  echo=2;
  charexp2generators(intvec(28,64,66,77));
}


proc charexp2multseq (intvec charexp)
"USAGE:      charexp2multseq(v),  v intvec
ASSUME:      v contains the characteristic exponents of an irreducible plane
             curve singularity
RETURN:      intvec, the multiplicity sequence of the plane curve singularity
NOTE:        If the curve singularity is smooth, then the multiplicity sequence is empty.
             This is expressed by returning zero.
SEE ALSO:    invariants, resolutiongraph, totalmultiplicities, alexanderpolynomial
KEYWORDS:    characteristic exponents; multiplicity sequence; topological invariants
EXAMPLE:     example charexp2multseq;   shows an example"
{
  int end=size(charexp);
  // If the singularity is smooth!
  if (end==1)
  {
    return(1); // ERROR: Should be 0, but for the time being, Singular returns 1.
  }
  intvec multseq=euclidseq(charexp[2],charexp[1]);
  for (int i=3;i<=end;i++)
  {
    multseq=multseq,euclidseq(charexp[i]-charexp[i-1],multseq[size(multseq)]);
  }
  return(multseq);
}
example
{
  "EXAMPLE:";
  echo=2;
  charexp2multseq(intvec(28,64,66,77));
}

proc  multseq2charexp(def v)   // Procedure written by Fernando.
"USAGE:  multseq2charexp(v), v intvec
ASSUME:  The input is an intvec, which contains the mutiplicity sequence
         of an irreducible plane curve singularity .
RETURN:  An intvec, which contains the sequence of characteristic
         exponents of the irreducible plane curve singularity defined by v.
EXAMPLE: example multseq2charexp; shows an example.
"
{
  ///////  Preamble which reduces the input of an intvec to  /////////////
  ///////  the originally assumed input of a list of intvecs /////////////
  ///////  and tests the input.                              /////////////
  if (typeof(v)=="intvec")
  {
    list RESULT=multseq2charexp(list(v));
    return(RESULT[1]);
  }
  if (typeof(v)!="list")
  {
    ERROR("Invalid Input!");
  }
  if (typeof(v)=="list")
  {
    int TESTV;
    for (int III=1;III<=size(v);III++)
    {
      if (typeof(v[III])!="intvec")
      {
        TESTV=1;
      }
    }
    if (TESTV==1)
    {
      ERROR("Invalid Input!");
    }
  }
  ///////////////////////////////////////////////////////////
  list L=v;
  int n =size(L);
  // ///////////////////////////////////////////////////////
  // we look the size of each vector
  intvec mm;
  for(int j=1;j<=n;j++)
  {
    mm[j]=size(L[j]);
  }
  // ///////////////////////////////////////////////////////
  // we define some variables
  list LL;
  int temp, temp1,temp2;
  int ind,r,l,boolean;
  int old,new;
  int contador;
  list EUCLI,EUCLI1;
  intvec exponent,exponentes1;
  int new1,old1;
  int contador1;
  int a,b,f;
  //with the for we round each branch.
  for(int k=1;k<=n;k++)
  {
    l=1;
    old=L[k][l];
    //if the vertor has more than one element
    if(mm[k]<>1)
    {
      // ///////////////////////////////////////////////////////////////////////////////
      // the first step is special because is easy to know the two first characteristic exponents
      new=L[k][l+1];
      contador=1;
      while(old==new)//we check how many consecutives elements are equal, starting in the first.
      {
        contador=contador+1;
        old=new;
        new=L[k][contador+1];
      }
      exponent=L[k][l],contador*L[k][l]+L[k][l+contador];// those are the first two characteristic exponents.
      LL[k]=exponent;// we insert them in the final matrix
      EUCLI=euclides(LL[k][2],LL[k][1]);// compute the euclides algorithm for the two first characteristic exponents.
      temp=size(EUCLI[1]);
      // measure how many elements of the multiplicity sequence belong to the first euclidean algorithm.
      for(ind=1;ind<=temp;ind=ind+1)
      {
        l=l+EUCLI[1][ind];
      }
      l=l-1;//this is the number we are looking for.
      ///////////////////////////////////////////////////////////////
      r=1;
      //repeat the same process until the end of the multiplicity sequence.
      if(mm[k]-1>l)
      {
        while( l<mm[k]-1)
        {
          r=r+1;
          old1=L[k][l];
          new1=L[k][l+1];
          contador1=0;
          boolean=1;
          if(old1==new1)
          {
            while(old1==new1 and boolean==1)
            {
              contador1=contador1+1;
              old1=new1;
              new1=L[k][l+contador1+1];
              if(size(L[k])<=l+contador1+1)
              {
                boolean =0;
              }
            }
          }
          temp1=size(LL[k]);
          exponentes1=LL[k],LL[k][temp1]+(contador1*L[k][l])+L[k][contador1+l+1];
          LL[k]=exponentes1;
          EUCLI1=euclides(LL[k][temp1+1]-LL[k][temp1],L[k][l]);
          temp2=size(EUCLI1[1]);
          for(ind=1;ind<=temp2;ind=ind+1)
          {
            l=l+EUCLI1[1][ind];
          }
        }
      }
    }
    // if the vector has only one element then the charexp is only 1.
    else
    {
      LL[k]=1;
    }
  }
  return(LL);
}
example
{
  "EXAMPLE:";echo=2;
  intvec v=2,1,1;
  multseq2charexp(v);
  intvec v1=4,2,2,1,1;
  multseq2charexp(v1);
}

proc charexp2inter (intmat contact, list charexp)
"USAGE:      charexp2inter(contact,charexp),  contact matrix, charexp list
ASSUME:      charexp contains the integer vectors of characteristic exponents
             of the branches of a plane curve singularity, and contact is their
             contact matrix
RETURN:      intmat, the matrix intersection multiplicities of the branches
SEE ALSO:    invariants, resolutiongraph, totalmultiplicities, semigroup
KEYWORDS:    contact matrix; characteristic exponents; intersection multiplicity; topological invariants
EXAMPLE:     example charexp2inter;   shows an example"
{
  int n=ncols(contact);
  int i,j,k;
  list multpl;
  int max=0;
  intvec helpvect;
  intmat inters[n][n];
  // Calculate the multiplicity sequences of the branches.
  for (i=1;i<=n;i++)
  {
    multpl[i]=charexp2multseq(charexp[i]);
    // Find the maximal length of a multiplicity sequence.
    if (max<size(multpl[i]))
    {
      max=size(multpl[i]);
    }
  }
  // If the contact of certain branches is higher than the maximal length of the
  // multiplicity sequence, then max should be the maximal contact!
  helpvect=max,intvec(contact);
  max=max_in_intvec(helpvect)[1];
  // Prolong them by 1s, in order to take care of higher contact.
  for (i=1;i<=n;i++)
  {
    helpvect=multpl[i];
    for (j=size(multpl[i]+1);j<=max;j++)
    {
      helpvect[j]=1;
    }
    multpl[i]=helpvect;
  }
  // Calculate the intersection numbers of the branches: for two branches f_i and f_j
  // this is the sum over mult_k(f_i)*mult_k(f_j), where k runs over all infinitely
  // near points which f_i and f_j share, i.e. from 1 to contact[i,j].
  for (i=1;i<=n;i++)
  {
    for (j=i+1;j<=n;j++)
    {
      for (k=1;k<=contact[i,j];k++)
      {
        inters[i,j]=inters[i,j]+multpl[i][k]*multpl[j][k];
      }
      inters[j,i]=inters[i,j];
    }
  }
  return(inters);
}
example
{
  "EXAMPLE:";echo=2;
  ring r=0,(x,y),ds;
  list INV=invariants((x2-y3)*(x3-y2)*((x2-y3)^2-4x5y-x7));
  intmat contact=INV[4][1];
  list charexp=INV[1][1],INV[2][1],INV[3][1];
  // The intersection matrix is INV[4][2].
  print(INV[4][2]);
  // And it is calulated as ...
  print(charexp2inter(contact,charexp));
}


proc charexp2conductor(intvec B)  // Procedure written by Fernando
"USAGE:      charexp2conductor(v),  v intvec
ASSUME:      v contains the characteristic exponents of an irreducible plane
             curve singularity
RETURN:      int, the conductor of the plane curve singularity
NOTE:        If the curve singularity is smooth, the conductor is zero.
SEE ALSO:    invariants, resolutiongraph, totalmultiplicities, semigroup
KEYWORDS:    conductor; characteristic exponents; multiplicity sequence; topological invariants
EXAMPLE:     example charexp2conductor;   shows an example"
{
  intvec E;
  int i,conductor;
  E[1]=B[1];
  for(i=2;i<=size(B);i++)
    {
      E[i]=gcd(E[i-1],B[i]);
    }
  conductor=1-E[1];
  for(i=2;i<=size(B);i++)
    {
      conductor=conductor+(B[i]*(E[i-1]-E[i]));
    }
  return(conductor);
}
example
{
  "EXAMPLE:";
  echo=2;
  charexp2conductor(intvec(2,3));  // A1-Singularity
  charexp2conductor(intvec(28,64,66,77));
}


proc charexp2poly(intvec v, vector a)  // Procedure written by Fernando.
"USAGE:  charexp2poly(v,a); intvec v, vector a.
ASSUME:  v an intvec containing the characterictic exponents of an irreducible plane curve singularity.
         a a vector containing the coefficients of a parametrization given by x(t)=x^v[1],
         y(t)=a(1)t^v[2]+...+a[n-1]t^v[n], i.e. the entries of a are of type number.
RETURN:  A polynomial f in the first two variables of the basering, such that f defines an
         irreducible plane curve singularity with characteristic exponents v.
NOTE:    The entries in a should be of type number and the vector v should be the sequence of
         characteristic exponents of an irreducible plane curve singularity in order to
         get a sensible result,
SEE ALSO: charexp2multseq, multseq2charexp.
EXAMPLE: example charexp2poly;  shows an example
"
{
  int n=size(v);
  vector expo;
  int i,j,s;
  for(i=1;i<=v[1];i++)
  {
    expo[i]=0;//initialize to 0.
  }
  for(i=2;i<=n;i++)
  {
    s=v[i] mod v[1];//calculate the position.
    expo=expo-a[i-1]*var(1)^((v[i]-s)/v[1])*gen(s+1);//save in expo -var(1) to the power the corresponding
                                                     //but only in the right positions
  }
  matrix M[v[1]][v[1]];
  //construct the matrix that generates the polynomial
  for(i=1;i<=v[1];i++)
  {
    M[i,i]=var(2)+expo[1];//The  diagonal
    for(j=1;j<=v[1]-i;j++)
         {
           M[j,j+i]=expo[i+1];//over diagonal
         }
    for(j=1;j<=v[1]-i;j++)
         {
           M[j+i,j]=var(1)*expo[1+v[1]-i];//under diagonal
         }
  }
  //the poynomial is the determinant of the matrix
  poly irredpoly=det(M);
  return(irredpoly)
}
example
{
  "EXAMPLE:";echo=2;
  ring r=0,(x,y),dp;
  intvec v=8,12,14,17;
  vector a=[1,1,1];
  poly f=charexp2poly(v,a);
  f;
  invariants(f)[1][1];  // The characteristic exponents of f.
}

proc tau_es2 (def INPUT, list #)
"USAGE:  tau_es2(INPUT); INPUT poly or list
ASSUME:  INPUT is either a REDUCED bivariate polynomial defining a plane curve singularity,
         or the output of @code{hnexpansion(f[,\"ess\"])}, or the list @code{hne} in
         the ring created by @code{hnexpansion(f[,\"ess\"])}, or the output of
         @code{develop(f)} resp. of @code{extdevelop(f,n)}, or a list containing
         the contact matrix and a list of integer vectors with the characteristic exponents
         of the branches of a plane curve singularity, or an integer vector containing the
         characteristic exponents of an irreducible plane curve singularity.
RETURN:  int, the equisingular Tjurina number of f, i. e. the codimension of the mu-constant
         stratum in the semiuniversal deformation of f, where mu is the Milnor number of f.
NOTE:    The equisingular Tjurina number is calculated with the aid of a Hamburger-Noether
         expansion, which is the hard part of the calculation.
         In case the Hamburger-Noether expansion of the curve f is needed
         for other purposes as well it is better to calculate this first
         with the aid of @code{hnexpansion} and use it as input instead of
         the polynomial itself.
         @*
         If you are not sure whether the INPUT polynomial is reduced or not, use
         @code{squarefree(INPUT)} as input instead.
SEE ALSO: tau_es, develop, hnexpansion, totalmultiplicities, equising_lib
EXAMPLE: example tau_es2;  shows an example
"
{
  // If the input is a weighted homogeneous polynomial, then use a direct algorithm to
  // calculate the equisingular Tjurina number, by caluclating a K-basis of the
  // Tjurina algebra and omitting those elements with weighted degree at least the
  // weighted degree of the polynomial. -- If an additional input # is given (which is
  // not just the number 1 !!!), the procedure always uses the recursive algorithm.
  if ((typeof(INPUT)=="poly") and (size(#)==0))
  {
    if (qhweight(INPUT)[1]!=0)
    {
      /*
      poly f=squarefree(INPUT);
      if ( deg(f)!=deg(INPUT) )
      {
        dbprint(printlevel-voice+3,"// input polynomial was not reduced");
        dbprint(printlevel-voice+3,"// we continue with its reduction");
      }
      return(tau_es_qh(f));
      */
      return(tau_es_qh(INPUT));
    }
  }
  // Else apply the recursive algorithm from Eugenii Shustin, On Manifolds of Singular
  // Algebraic Curves, Selecta Math. Sov. Vol 10, No. 1 (1991), p. 31.
  int i,j,k;                  // Laufvariablen
  int tau,tau_i;              // Variable for es-Tjurina no., multiplicitiy, and es-Tjurina at blow ups
  intmat contact;             // contact matrix
  list graphs, multipl;       // resolution graphs and multiplicity sequences
  list graphs_ed, multipl_ed; // res.graphs and mult.seq. of one curve on 1st ex. div.
  intmat contact_ed;          // contact matrix of one curve on 1st exceptional divisor
  int d_i;                    // correction term
  intvec connections;
  intmat graph;
  intvec multseq;
  int s,e;
  int multi;                  // multiplicity of the curve defined by the input.
  /////////////////////////////////////////////////////////
  // Check what type the input is, and act accordingly.
  if (typeof(INPUT)=="list")
  {
    if (size(INPUT)==3)
    {
      // If the INPUT is the output of totalmultiplicities(INPUT,"tau")
      if ((typeof(INPUT[1])=="list") and (typeof(INPUT[2])=="list") and (typeof(INPUT[3])=="intmat"))
      {
        graphs=INPUT[1];
        multipl=INPUT[2];
        contact=INPUT[3];
      }
      // Otherwise call the procedure with the output of totalmultiplicities(INPUT,"tau").
      else
      {
        return(tau_es2(totalmultiplicities(INPUT,"tau")));
      }
    }
    else
    {
      return(tau_es2(totalmultiplicities(INPUT,"tau")));
    }
  }
  // Otherwise call the procedure with the output of totalmultiplicities(INPUT,"tau").
  else
  {
    return(tau_es2(totalmultiplicities(INPUT,"tau")));
  }
  /// End of checking the input
  ///////////////////////////////////////////////////////
  // If the singularity is smooth, the equisingular Tjurina number is zero.
  if ((ncols(contact)==1) and (multipl[1][1]<=1))
  {
    return(0);
  }
  // Otherwise calculate the multiplicity of the singularity.
  for (i=1;i<=size(multipl);i++)
  {
    multi=multi+multipl[i][1];
  }
  // Find the branches which stay together after blowing up once, and group them together.
  k=0;
  intvec curves_on_first_ex_div=k;  // If this is 0=i_0,i_1,...i_s=ncols(contact), then the branches
  while (k<ncols(contact))          // (i_j)+1,...,i_(j+1) stay together after the first blowing up,
  {                                 // and s is the number of infinitely near points of first order.
    k=find_last_non_one(intvec(contact[k+1,1..ncols(contact)]),k+2);
    curves_on_first_ex_div=curves_on_first_ex_div,k;
  }
  // And then calculate the equisingular Tjurina numbers in the infinitely near points of
  // first order plus some correction term (= the number of blow ups needed to separate the
  // strict transform from the exceptional divisor in the corresponding infinitely near point)
  // and add them.
  for (i=1;i<=size(curves_on_first_ex_div)-1;i++)
  {
    s=curves_on_first_ex_div[i]+1;
    e=curves_on_first_ex_div[i+1];
    graphs_ed=list(graphs[s..e]);
    multipl_ed=list(multipl[s..e]);
    contact_ed=intmat_minus_one(intmat(intvec(contact[s..e,s..e]),e-s+1,e-s+1));
    connections=0;
    // Adjust the graphs and multiplicity sequences by cutting the first level.
    // And find in each graph the point to which the point 1 is connected = 1 + number of
    // blow ups needed to separate the strict transform from the first exceptional divisor.
    for (j=1;j<=size(graphs_ed);j++)
    {
      // Adjust the graphs and find the connection points.
      graph=graphs_ed[j];
      connections[j]=find_connection_point(intvec(graph[1,1..ncols(graph)]),1);
      graph=delete_row(delete_col(graph,1),1);
      graphs_ed[j]=graph;
      // Adjust the multiplicity sequences.
      multseq=multipl_ed[j];
      if (size(multseq)==1)  // If multseq has only one entry, then their is only
      {                      // one branch left and it is smooth! So set multipl_ed[j]=1.
        multipl_ed[j]=intvec(1);
      }
      else  // Otherwise just cut the first entry.
      {
        multipl_ed[j]=intvec(multseq[2..size(multseq)]);
      }
    }
    // Calculate the equisingular Tjurina number of the strict transform at the i-th
    // intersection point with the first exceptional divisor.
    tau_i=tau_es2(list(graphs_ed,multipl_ed,contact_ed));
    // Calculate the number no. of blow ups needed to separate the strict transform of
    // the curve defined by the input from the first exceptional divisor at their i-th
    // intersection point.
    d_i=max_in_intvec(connections)-1;
    // If d_i is negative, then there was only one branch left and it was smooth, so d_i should be 1.
    if (d_i<0)
    {
      d_i=1;
    }
    // If the curve defined by graphs_ed was smooth, then the summand has to be reduced by 1.
    if (tau_i==0)
    {
      tau_i=-1;
    }
    tau=tau+tau_i+d_i;
  }
  // The equisingular Tjurina number is then calculated by adding the following term.
  tau=tau+((multi*(multi+1))/2)-2;
  return(tau);
}
example
{
  "EXAMPLE:";
  echo=2;
  ring r=0,(x,y),ls;
  poly f1=y2-x3;
  poly f2=(y2-x3)^2-4x5y-x7;
  poly f3=y3-x2;
  tau_es2(f1);
  tau_es2(f2);
  tau_es2(f1*f2*f3);
}


///////////////////////////////////////////////////////////////////////////////////////////
// Static procedures.
///////////////////////////////////////////////////////////////////////////////////////////

static proc euclidseq(int a,int b)
"USAGE:      euclidseq(a,b),  a,b integers
RETURN:      intvec,  the divisors in the euclidean alogrithm with multiplicities
KEYWORDS:    Euclidean algorithm; multiplicity sequence
NOTE :       This procedure is for internal use only; it is called by charexp2multseq.
"
{
  int i;
  // multseq saves in each step of the Euclidean algorithm q-times the divisor b
  intvec multseq;
  int q=a div b;
  int r=a mod b;
  for (i=1;i<=q;i++)
  {
    multseq[i]=b;
  }
  int s=q;        // size of multseq
  a=b;
  b=r;
  while(r<>0)
  {
        q=a div b;
        r=a mod b;
    for (i=1;i<=q;i++)
    {
      multseq[s+i]=b;
    }
    s=s+q;        // size of multseq
        a=b;
        b=r;
  }
  return(multseq);
}

static proc puiseuxchainpart (int piij, int muij, intvec multpl, int initial_tm, int end_tm, int j)
"USAGE:   puiseuxchainpart(piij,muiij,multpl,initial_tm,end_tm,j);
RETURN:   list L, L[1] incidence matrix of a part of the Puiseux chain, L[2] the total
          multiplicities of this part of the Puiseux chain.
NOTE:     This procedure is only for internal use; it is called by puiseuxchain.
"
{
  int delta=1;
  if (j==1){delta=0;}                 // Delta measures whether j is 1 or not.
  intvec totalmultiplicity;           // Keeps the total multiplicities.
  intmat pcp[muij][muij];             // Keeps the incidence matrix of the Puiseuxchainpart S_i,j.
  // Calculate the total multiplicities and the Puiseuxchainpart S_i,j.
  totalmultiplicity[1]=initial_tm+end_tm+multpl[1];
  pcp[1,1]=piij+1;
  for (int k=2;k<=muij;k++)
  {
    pcp[k,k]=piij+k;
    pcp[k-1,k]=1;
    pcp[k,k-1]=1;
    totalmultiplicity[k]=totalmultiplicity[k-1]+delta*initial_tm+multpl[k];
  }
  list result=pcp,totalmultiplicity;
  return(result);
}

static proc puiseuxchain (int initial, intvec divseq, intvec multpl, int initial_tm)
"USAGE:   puiseuxchain(initial,divseq,multpl,initial_tm); int initial, initial_tm, intvec divseq, multpl
RETURN:   list L, L[1] incidence matrix of a Puiseux chain, L[2] the weight of the point to which the
          previous Puiseux chain has to be connected, L[3] the sequence of total multiplicities of
          the points in this Puiseux chain.
NOTE:     This procedure is only for internal use; it is called by irred_resgraph_totmult.
"
{
  int j,k,l,connectpoint;
  intvec multpli;
  int pc_tm=initial_tm;       // Keeps the total multipl. of the endpoint of P_i-1.
  int end_tm=0;
  int start=1;
  int omega=size(divseq);
  // Keep the endpoints of the puiseuxchainparts (minus initial)  s_i,j in divseq_sum.
  intvec divseq_sum=divseq[1];
  for (j=2;j<=omega;j++)
    {
      divseq_sum[j]=divseq_sum[j-1]+divseq[j];
    }
  // Define the connecting point of the Puiseuxchain P_i-1 with P_i.
  if (divseq[1]==0)
  {
    // If divseq[1]=mu_i,1=0, then the Puiseuxchainpart S_i,1 is empty.
    // We may start building the Puiseuxchain with part 2, i.e. set start=2.
    start=2;
    if (omega>=3)
    {
      connectpoint=initial+divseq_sum[2]+1;  // startpoint of s_i+1,3
    }
    else
    {
      connectpoint=initial+divseq_sum[2];    // endpoint of s_i+1,2
    }
  }
  else
  {
    connectpoint=initial+1;
  }
  // Build the Puiseuxchainparts s_i,j and put them as blocks into pc=P_i.
  multpli=multpl[initial+1..initial+divseq_sum[start]];
  list pcp=puiseuxchainpart(initial,divseq[start],multpli,initial_tm,end_tm,start);
  intmat pc=pcp[1];
  intvec tm=pcp[2];
  for (j=start+1;j<=omega;j++)
    {
      // Keep the final total multipl. of the puiseuxchainpart S_i,j-2 resp. P_i-1, if S_i,1 empty.
      if (j>2){end_tm=initial_tm;}
      // Calculate the endpoint of S_i,j-1.
      initial=initial+divseq[j-1];
      // Keep the total multiplicity of the endpoint of S_i,j-1
      initial_tm=pcp[2][size(pcp[2])];
      // Build the new puiseuxchainpart S_i,j
      multpli=multpl[initial+1..initial+divseq[j]];
      pcp=puiseuxchainpart(initial,divseq[j],multpli,initial_tm,end_tm,j);
      pc=addmat(pc,pcp[1]);
      tm=tm,pcp[2];
    }
  // Connect the Puiseuxchainparts s_i,j.
  for (j=start;j<=omega-2;j++)
    {
      k=divseq_sum[j];          // endpoint of s_i,j
      l=divseq_sum[j+1]+1;      // startpoint of s_i,j+2
      pc[k,l]=1;            // connecting these
      pc[l,k]=1;
    }
  if (omega>=start+1)   // If there are at least two non-empty s_i,j.
    {
      k=divseq_sum[omega-1];    // endpoint of s_i,omega-1
      l=divseq_sum[omega];      // endpoint of s_i,omega
      pc[k,l]=1;            // connecting these
      pc[l,k]=1;
    }
  list ergebnis=pc,connectpoint,tm;
  return(ergebnis);
}

static proc irred_resgraph_totmult (intvec charexp)
"USAGE:   irred_resgraph_totmult(charexp); charexp intvec
ASSUME:   charexp is an integer vector containg the characteristic exponents of
          an irreducible plane curve singularity
RETURN:   list L, L[1] is the incidence matrix of the resolution graph of the plane curve
          singularity defined by INPUT, and L[2] is its sequence of total multiplicities
NOTE:     This procedure is only for internal use; it is called by resgraph.
"
{
  int k,l;
  intvec multpl=charexp2multseq(charexp);  // multiplicity sequence of the singularity
  // Do first the case where the singularity is actually smooth.
  if (size(charexp)==1)
  {
    intmat resgraph[1][1]=0;
    intvec tm=1;   // there is no exceptional divisor in the resolution - ERROR: should be 0,
                   // but for the time being, Singular returns 1 as multiplicity of smooth curves
    list result=resgraph,tm,multpl;
    return(result);
  }
  // End of the smooth case
  int initial_tm=0;                  // total multipl. of the endpoint of Puiseux chain P_i-1
  int g=size(charexp);
  list es=divsequence(charexp[2],charexp[1]);   // keeps the lenghts of the Puiseuxchainparts s_i,j
  intvec divseq=es[1];
  int r=es[2];
  int initial=0;
  // Build the Puiseuxchains P_i and put them as blocks into a matrix.
  list pc=puiseuxchain(initial,divseq,multpl,initial_tm);
  intmat resgraph=pc[1];
  intvec endpoints=resgraph[nrows(resgraph),ncols(resgraph)];
  intvec connectpoints=pc[2];
  intvec tm=pc[3];
  for (int i=3;i<=g;i++)
    {
      initial_tm=tm[size(tm)];
      es=divsequence(charexp[i]-charexp[i-1],r);
      divseq=es[1];
      r=es[2];
      initial=endpoints[size(endpoints)];
      pc=puiseuxchain(initial,divseq,multpl,initial_tm);
      resgraph=addmat(resgraph,pc[1]);
      endpoints=endpoints,resgraph[nrows(resgraph),ncols(resgraph)];
      connectpoints=connectpoints,pc[2];
      tm=tm,pc[3];
    }
  // Adding the * for the strict transform to the Graph.
  resgraph=addmat(resgraph,0);
  // The connecting point of the * with the graph.
  connectpoints=connectpoints,nrows(resgraph);
  // Connect the P_i with each other and with *.
  for (i=2;i<=g;i++)
  {
    k=endpoints[i-1];          // endpoint of P_i-1
    l=connectpoints[i];        // conncting point of P_i resp. of *
    resgraph[k,l]=1;          // connecting these
    resgraph[l,k]=1;
  }
  list result=resgraph,tm,multpl; //HIER GEAENDERT!!!!
  return(result);
}


static proc max_in_intvec (intvec v, list #)
"USAGE:   max_in_intvec(v); v intvec
RETURN:   int m, maximum of the integers in v
USAGE:    max_in_intvec(v,1); v intvec
RETURN:   intvec m, m[1] maximum of the integers in v, m[2] position of the
          last occurence of the maximum in v
NOTE:     This procedure is only for internal use; this procedure is called by
          totalmultiplicities and semigroup.
"
{
  int max=v[1];
  int maxpos=1;
  for (int i=2;i<=size(v);i++)
  {
    if (v[i]>max)
    {
      max=v[i];
      maxpos=i;
    }
  }
  if (size(#)==0)
  {
    return(max);
  }
  else
  {
    return(intvec(max,maxpos));
  }
}

static proc addmat (intmat A,intmat B)
"USAGE:   max_in_intvec(A,B); A, B integer matrices
RETURN:   intmat C, block matrix with left-upper block A, right-lower block B
NOTE:     This procedure is only for internal use; this procedure is called several times.
"
{
  int nc=ncols(A);
  int nr=nrows(A);
  int mc=ncols(B);
  int mr=nrows(B);
  int i,j;
  intmat AB[nr+mr][nc+mc];
  for (i=1;i<=nr;i++)
    {
      for (j=1;j<=nc;j++)
        {
          AB[i,j]=A[i,j];
        }
    }
  for (i=1;i<=mr;i++)
    {
      for (j=1;j<=mc;j++)
        {
          AB[i+nr,j+nc]=B[i,j];
        }
    }
  return(AB);
}

static proc divsequence(int a,int b)
"USAGE:   divsequence(a,b); a,b integers
RETURN:   list l, l[1] the multiplicities of the divisors in the Euclidean algorithm
          and l[2] the last non-zero remainder in the Euclidean algorithm
NOTE:     This procedure is only for internal use; it is called in irred_res_graph.
"
{
  int q=a div b;
  int r=a mod b;
  intvec divseq=q;
  while(r<>0)
  {
    a=b;
    b=r;
    q=a div b;
    r=a mod b;
    divseq = divseq,q;
  }
  list result=divseq,b;
  return(result);
}



static proc adjust_tot_mult (intvec rtm_fix, rtm_var, rmt_fix, rmt_var, p, q, stricttransforms, int k)
"USAGE:   adjust_tot_mult(v1,v2,v3,v4,p,q,st,k); v1,...,st intvecs and k an integer
RETURN:   intvec rtm_var, adjusted total multiplicities
NOTE:     This procedure is only for internal use; it is called in totalmultiplicities.
"
{
  int j,l,store;  // Help variables.
  // Recalculate the entries in rtm_var from stricttransforms[k]+1,...,stricttransforms[k+1]-1.
  for (j=stricttransforms[k]+1;j<stricttransforms[k+1];j++)
  {
    if (rtm_var[j]==0) // If the entry is non-zero, we know that it is already correct.
    {
      // Check if the vertex in the fixed part is connected to one or to two vertices of lower weight.
      if (j==stricttransforms[k]+1)  // The vertex of weight 1 less is p[k], to which the subgraph is glued.
      {
        store=rtm_fix[j]-rmt_fix[j]-rtm_fix[p[k]];
      }
      else                           // The vertex of weight 1 less belongs to the subgraph.
      {
        store=rtm_fix[j]-rmt_fix[j]-rtm_fix[j-1];
      }
      // It is connected to two vertices V (which has weight one less) and W.
      if (store>0)
      {
        if (j==stricttransforms[k]+1)  // V is p[k] (to which the subgraph was glued) and W is q[k], the
        {                              // vertex of weight one less, to which p[k] is connected.
          rtm_var[j]=rtm_var[p[k]]+rtm_var[q[k]];  // In this case the subgraph separates p[k] and q[k]!
        }
        if (j==stricttransforms[k]+2)  // V belongs to the subgraph (it is the vertex considerd in the
        {  // previous case!), and W is p[k] or q[k]. In this case the subgraph separates p[k] and q[k].
          if (store==rtm_fix[p[k]])  // If W=p[k], ...
          {
            rtm_var[j]=rtm_var[j-1]+rtm_var[p[k]];
          }
          else                       // else W=q[k] ... .
          {
            rtm_var[j]=rtm_var[j-1]+rtm_var[q[k]];  // separates p[k] and q[k].
          }
        }
        if (j>stricttransforms[k]+2)  // V belongs to the subgraph and W either does as well or is p[k].
        {
          l=j-2;
          while (store<rtm_fix[l]) // Find the second vertex W which is connected to which it is
          {                        // connected. It has total multipl. = store!
            if (l>stricttransforms[k]+1) // If l-1 belongs to the graph, then reduce l.
            {
              l=l-1;
            }
            else // If l-1 is stricttransform[k], hence isn't in the graph, then reducing l gives p[k].
            {    // If l=p[k] and still store<rtm_fix[l], then j must be connected to q[k]!
              if (l==stricttransforms[k]+1)
              {
                l=p[k];
              }
              else
              {
                l=q[k];
              }
            }
          }
          rtm_var[j]=rtm_var[j-1]+rtm_var[l];
        }
      }
      // It is only connected to one vertex V, which then must be the one of weight one less.
      else
      {
        if (j==stricttransforms[k]+1) // V is p[k], the vertex, to which the subgraph was glued.
        {
          rtm_var[j]=rtm_var[p[k]];
        }
        else
        {
          rtm_var[j]=rtm_var[j-1];  // V is belongs already to the subgraph.
        }
      }
    }
  }
  return(rtm_var);
}


static proc find_connection_point (intvec v, int k)
"USAGE:   find_connection_point(v,c); where v is an intvec, and k is an integer
RETURN:   The largest integer i>k, such that v[i]=1, or 0 if no such i exists.
NOTE:     This procedure is only for internal use; it is called in totalmultiplicities.
"
{
  for (int i=size(v)-1;i>=k+1;i--)
  {
    if (v[i]==1)
    {
      return(i);
    }
  }
  return(0);
}

static proc find_connection_points (intmat resgr, int k)
"USAGE:   find_connection_points(resgr,k); where resgr is an intmat, and k is an integer
RETURN:   list of two intvecs ctps and ctpswts, where ctps contains all integers i!=k, such
          that resgr[k,i]=1, and ctpswts contains for each such i the weight resgr[i,i].
NOTE:     This procedure is only for internal use; it is called in totalmultiplicities.
"
{
  intvec ctps;
  intvec ctpswts;
  int j=1;
  for (int i=1;i<=ncols(resgr);i++)
  {
    if ((resgr[k,i]==1) and (i!=k))
    {
      ctps[j]=i;
      ctpswts[j]=resgr[i,i];
      j++;
    }
  }
  return(list(ctps,ctpswts));
}

static proc find_lower_connection_points (intmat resgr, int k)
"USAGE:   find_lower_connection_points(resgr,k); where resgr is an intmat, and k is an integer
ASSUME:   resgr is the resolutiongraph of an IRREDUCIBLE curve singularity and k<ncols(resgr).
RETURN:   intvec, which contains the weights of points of lower weight than k, to which
          the point of weight k in resgr is connected, and 0 if no such point exists.
NOTE:     This procedure is only for internal use; it is called in totalmultiplicities.
"
{
  intvec ctps=find_connection_points(resgr,k)[2];
  intvec lower_ctps;
  int i=1;
  while ((ctps[i]<k) and (ctps[i]>0))
  {
    lower_ctps[i]=ctps[i];
    i++;
  }
  return(lower_ctps);
}


static proc euclides(int a,int b)  // Procedure of Fernando.
"USAGE:   euclides(m,n);where m,n are integers.
RETURN:   The divisor, the quotients and the remainders of the euclidean algorithm performing for m and n.
NOTE:     This procedure is only for internal use; it is called in multseq2charexp.
"
{
 int c=a div b;//we compute in c the integer division between a and b.
 int r=a mod b;//in r the remainer of the division between a and b
 intvec dividendo=c;// we define the intvec of the dividens and we initialize it to c
 intvec remain=r;// we define the intvec of the remainders and we initialize it to r
 a=b;//change a to b
 b=r;// and b to r

 while(r<>0)// while the current remainder is diferent to 0 we make the same af before
  {
    c=a div b;
    r=a mod b;
    dividendo =dividendo,c;
    if(r<>0)
      {
        remain=remain,r;
      }
       a=b;
       b=r;
     }
   list L=dividendo,remain;//we put in a list all the dividens and all the remainders
   return(L);// and return the list
}



static proc tau_es_qh (poly f)
"USAGE:    tau_es_qh(f), poly f
RETURN:   int, equisingular Tjurina number
NOTE:     This procedure is only for internal use; it is called in Tau_es.
"
{
  intvec qh=qhweight(f);
  if (qh[1]==0)
  {
    ERROR("Input is not quasi-homogenous.");
  }
  else
  {
    int d_f = deg(f,qh);
    list Tj=Tjurina(f,1);
    ideal tj=Tj[2];
    int Taues=size(tj);
    for (int i=1;i<=size(tj);i++)
    {
      if (deg(tj[i],qh)>=d_f)
      {
        Taues--;
      }
    }
  }
  return(Taues);
}


static proc move_row (intmat M, int i,j)
"USAGE:    move_row(M,i,j), intmat M, int i,j
RETURN:   intmat, the matrix M with j-th row now i-th row and the remaining rows moved accordingly.
NOTE:     This procedure is only for internal use; it is called in move_row_col.
"
{
  if ((i>nrows(M)) or (j>nrows(M)))
  {
    ERROR("The matrix has not enough rows.");
  }
  if (i==j)
  {
    return(M);
  }
  if (i>1)
  {
    intmat N[nrows(M)+1][ncols(M)]=intvec(M[1..i-1,1..ncols(M)]),intvec(M[j,1..ncols(M)]),intvec(M[i..nrows(M),1..ncols(M)]);
  }
  if (i==1)
  {
    intmat N[nrows(M)+1][ncols(M)]=intvec(M[j,1..ncols(M)]),intvec(M[i..nrows(M),1..ncols(M)]);
  }
  if (i<j)
  {
    N=delete_row(N,j+1);
  }
  if (i>j)
  {
    N=delete_row(N,j);
  }
  return(N);
}

static proc move_col (intmat M, int i,j)
"USAGE:    move_col(M,i,j), intmat M, int i,j
RETURN:   intmat, the matrix M with j-th column now i-th column and the remaining columns moved accordingly.
NOTE:     This procedure is only for internal use; it is called in move_row_col.
"
{
  return(transpose(move_row(transpose(M),i,j)));
}

static proc move_row_col (intmat M,int i,j)
"USAGE:    move_row(M,i,j), intmat M, int i,j
RETURN:   intmat, the matrix M with j-th row/column now i-th row/column and the remaining
          rows and columns moved accordingly.
NOTE:     This procedure is only for internal use; it is called in totalmultiplicities.
"
{
  return(move_col(move_row(M,i,j),i,j));
}


static proc delete_row (intmat M,int i)
"USAGE:    delete_row(M,i); M intmat, i int
RETURN:   intmat, which is derived from M by removing the ith row
NOTE:     This procedure is only for internal use; it is called in move_row and tau_es.
"
{
  if ((i!=1) and (i!=nrows(M)))
  {
    return(intmat(intvec(M[1..i-1,1..ncols(M)],M[i+1..nrows(M),1..ncols(M)]),nrows(M)-1,ncols(M)));
  }
  if (i==1)
  {
    return(intmat(intvec(M[2..nrows(M),1..ncols(M)]),nrows(M)-1,ncols(M)));
  }
  else
  {
    return(intmat(intvec(M[1..nrows(M)-1,1..ncols(M)]),nrows(M)-1,ncols(M)));
  }
}

static proc delete_col (intmat M, int i)
"USAGE:    delete_col(M,i); M intmat, i int
RETURN:   intmat, which is derived from M by removing the ith column
NOTE:     This procedure is only for internal use; it is called in tau_es.
"
{
  return(transpose(delete_row(transpose(M),i)));
}



static proc sort_branches (intmat contact, list graphs, list totmult, list multpl, int k,l)
"USAGE:    sort_branches(K,L,M,N,k,l); K intmat, L,M,N lists, k,l integers
ASSUME:   K = contact matrix of the branches of a curve, L = their resolutiongraphs,
          M = their totalmultiplicities, N = their multiplicities
RETURN:   list LL, LL[1] transformed K, LL[2..4] transformed L,M,N.
          The procedure sorts the branches from k+1 to l according to descending contact with
          with the branch k.
NOTE:     This procedure is only for internal use; it is called in totalmultiplicities.
"
{
  intvec max;
  for (int i=k+1;i<=l;i++)
  {
    // Find the last graph max between i and l which has maximal contact with k.
    max=max_in_intvec(intvec(contact[k,i..l]),1);
    max[2]=max[2]+i-1;
    if (i!=max[2])  // If max is not i, then move max to position i!
    {
      graphs=insert(graphs,graphs[max[2]],i-1);
      graphs=delete(graphs,max[2]+1);
      totmult=insert(totmult,totmult[max[2]],i-1);
      totmult=delete(totmult,max[2]+1);
      multpl=insert(multpl,multpl[max[2]],i-1);
      multpl=delete(multpl,max[2]+1);
      contact=move_row_col(contact,i,max[2]);
    }
  }
  return(list(contact,graphs,totmult,multpl));
}


static proc find_last_non_one (intvec v,int k)
"USAGE:    find_last_non_one (v,k); intvec v, int k
RETURN:   int i, the last index i>=k such that v[i]!=1, or k-1 if no such i exists.
NOTE:     This procedure is only for internal use; it is called in tau_es.
"
{
  int i=size(v);
  while (i>=k)
  {
    if (v[i]!=1)
    {
      return(i);
    }
    else
    {
      i--;
    }
  }
  return(i);
}

static proc intmat_minus_one (intmat M)
"USAGE:    intmat_minus_one(M);  intmat M
RETURN:   intmat, all non-zero entries of M deminuished by 1.
NOTE:     This procedure is only for internal use; it is called in tau_es.
"
{
  int i,j;
  for (i=1;i<=nrows(M);i++)
  {
    for (j=1;j<=ncols(M);j++)
    {
      if (M[i,j]!=0)
      {
        M[i,j]=M[i,j]-1;
      }
    }
  }
  return(M);
}