source: git/Singular/LIB/alexpoly.lib @ e009ef

////
version="version alexpoly.lib 4.0.0.0 Jun_2013 "; // $Id$
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
 proximitymatrix(f);        calculates the proximity matrix 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_es2 - 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;
  int inputispoly=0;
  /////////////////////////////////////////////////////////////////////////////////
  // 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);
    inputispoly=1;
  }
  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 polynomial 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_es2.
  int check_tau;
  if (size(#)==1)
  {
    if (#[1]=="tau")
    {
      check_tau=1;
    }
  }
  /////////////////////////////////////////////////////////////////////////////////
  // The procedure tau_es2 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_es2, 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 two 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;
  list newpos;  // when reordering the branches, here we save the new order of the old branches
  for (ii=1;ii<=r;ii++)
  {
    newpos[ii]=ii;
  }
  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.
      newpos=insert(newpos,j,i-1);
      newpos=delete(newpos,j+1);
      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);
    }
  }
  if (inputispoly==1) // if the input was a polynomial the ordering of the branches is fine
  {
    list result=rgraph,rtm,rmt;
    return(result);
  }
  else // reorder the branches according to the ordering of the input
  {
    // reordered total multiplicities, multiplicities and resolution graph
    intmat rtmro[nrows(rtm)][ncols(rtm)];
    intmat rmtro[nrows(rmt)][ncols(rmt)];
    intmat rgraphro=rgraph;
    for (i=1;i<=r;i++)
    {
      rtmro[1..nrows(rtm),newpos[i]]=rtm[1..nrows(rtm),i];
      rmtro[1..nrows(rmt),newpos[i]]=rmt[1..nrows(rmt),i];
      for (j=1;j<=nrows(rgraph);j++)
      {
        if (rgraph[j,j]==-i)
        {
          rgraphro[j,j]=-newpos[i];
        }
      }
    }
    list result=rgraphro,rtmro,rmtro;
    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]) div 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) div 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)) div 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 lengths 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, intvec rtm_var, intvec rmt_fix, intvec rmt_var, def p, def q, def 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,int 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,int 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,int 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_es2.
"
{
  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,int 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_es2.
"
{
  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_es2.
"
{
  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);
}

proc proximitymatrix (def INPUT)
"USAGE:  proximitymatrix(INPUT); INPUT poly or list or intmat
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, or the resolution
         graph of a plane curve singularity (i.e. the output of resolutiongraph or
         the first entry in the output of totalmultiplicities).
RETURN:  list, of three integer matrices. The first one is the proximity matrix of
         the plane curve defined by the INPUT, i.e. the entry i,j is 1 if the
         infinitely near point corresponding to row i is proximate to the infinitely
         near point corresponding to row j. The second integer matrix is the incidence matrix of the
         resolution graph of the plane curve. The entry on the diagonal in row i is -s-1
         if s is the number of points proximate to the infinitely near point corresponding
         to the ith row in the matrix. The third integer matrix is the incidence matrix of
         the Enriques diagram of the plane curve singularity, i.e. each row corresponds to
         an infinitely near point in the minimal standard resolution, including the
         strict transforms of the branches, the diagonal element gives
         the level of the point, and the entry i,j is -1 if row i is proximate to row j.
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.
         @*
         If the input is a smooth curve, then the output will consist of
         three one-by-one zero matrices.
         @*
         For the definitions of the computed objects see e.g. the book
         Eduardo Casas-Alvero, Singularities of Plane Curves.
SEE ALSO: develop, hnexpansion, totalmultiplicities, alexanderpolynomial
EXAMPLE: example proximitymatrix;  shows an example
"
{
  ///////// Decide on the type of input. //////////
  if (typeof(INPUT)=="intmat")
  {
    intmat resgr=INPUT;
  }
  else
  {
    intmat resgr=totalmultiplicities(INPUT)[1];
  }
  ////////  Deal with the case of a smooth curve ////////////////
  if (size(resgr)==1)
  {
    return(list(intmat(intvec(1),1,1),intmat(intvec(-1),1,1),intmat(intvec(0),1,1)));
  }
  ////////  Calculate the proximity resolution graph ////////////
  int i,j;
  int n=nrows(resgr);
  intvec diag; // Diagonal of the Enriques diagram.
  int si; // number of points proximate to the point corresponding to the ith row
  // Adjust the weights of the nodes corresponding to strict transforms so
  // as if there had been one additional blow up.
  for (i=1;i<=n;i++)
  {
    // Check if the row corresponds to an exceptional divisor ...
    if (resgr[i,i]<0)
    {
      j=1;
      while ((resgr[i,j]==0) or (i==j))
      {
        j++;
      }
      resgr[i,i]=resgr[j,j]+1;
    }
  }
  // Set the weights in the resolution graph to the blowing up level in the resolution.
  for (i=1;i<=n;i++)
  {
    resgr[i,i]=resgr[i,i]-1;
    diag[i]=resgr[i,i]; // The level of the corresponding infinitely near point.
  }
  // Replace the weights in the resolution graph by
  // -s-1, where s is the number of points which are proximate to the point.
  for (i=1;i<=n;i++)
  {
    si=-1;
    // Find the points of higher weight which are connected to the ith row.
    for (j=i+1;j<=n;j++)
    {
      // If the point in row j is connected to the ith row, then all the points of
      // weight resgr[i,i]+1 to weight resgr[j,j] in the corresponding subgraph
      // are proximate to the point of the ith row. This number is just resgr[j,j]-resgr[i,i].
      if ((resgr[i,j]!=0) and (resgr[j,j]>0))
      {
        si=si-(resgr[j,j]-resgr[i,i]);
      }
    }
    resgr[i,i]=si;
  }
  ///////////////  Calculate the proximity matrix  ///////////////////
  intmat proximity=proxgauss(resgr);
  intmat enriques=proximity;
  for (i=1;i<=nrows(enriques);i++)
  {
    enriques[i,i]=diag[i];
  }
  return(list(proximity,resgr,enriques));
}
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 proximity=proximitymatrix(f1*f2*f3);
  /// The proximity matrix P ///
  print(proximity[1]);
  /// The proximity resolution graph N ///
  print(proximity[2]);
  /// They satisfy N=-transpose(P)*P ///
  print(-transpose(proximity[1])*proximity[1]);
  /// The incidence matrix of the Enriques diagram ///
  print(proximity[3]);
  /// If M is the matrix of multiplicities and TM the matrix of total
  /// multiplicities of the singularity, then  M=P*TM.
  /// We therefore calculate the (total) multiplicities. Note that
  /// they have to be slightly extended.
  list MULT=extend_multiplicities(totalmultiplicities(f1*f2*f3));
  intmat TM=MULT[1];  // Total multiplicites.
  intmat M=MULT[2];   // Multiplicities.
  /// Check: M-P*TM=0.
  M-proximity[1]*TM;
  /// Check: inverse(P)*M-TM=0.
  intmat_inverse(proximity[1])*M-TM;
}

static proc addmultiplrows (intmat M, int i, int j, int ki, int kj)
"USAGE:   addmultiplrows(M,i,j,ki,kj);  intmat M, int i,j,ki,kj
RETURN:   intmat, replaces the j-th row in M by ki-times the i-th row plus
                  kj times the j-th
NOTE:     This procedure is only for internal use; it is called in intmat_inverse
          and proxgauss.
"
{
  int k=ncols(M);
  M[j,1..k]=kj*M[j,1..k]+ki*M[i,1..k];
  return(M);
}


static proc proxgauss (intmat M)
"USAGE:   proxgauss(M);  intmat M
ASSUME:   M is the output of proximity_resgr
RETURN:   intmat, replaces the j-th row in M by ki-times the i-th row plus
                  kj times the j-th
NOTE:     This procedure is only for internal use; it is called in intmat_inverse.
"
{
  int i;
  int n=nrows(M);
  if (n==1)
  {
    M[1,1]=1;
    return(M);
  }
  else
  {
    if (M[n,n]<0)
    {
      M=addmultiplrows(M,n,n,-1,0);
    }
    for (i=n-1;i>=1;i--)
    {
      if (M[i,n]!=0)
      {
        M=addmultiplrows(M,n,i,-M[i,n],M[n,n]);
      }
    }
    intmat N[n-1][n-1]=M[1..n-1,1..n-1];
    N=proxgauss(N);
    M[1..n-1,1..n-1]=N[1..n-1,1..n-1];
    return(M);
  }
}

proc extend_multiplicities (list rg)
"USAGE:      extend_multiplicities(rg); list rg
ASSUME:      rg is the output of the procedure totalmultiplicities
RETURN:      list, the first entry is an integer matrix containing the total
             multiplicities and the second entry is an integer matrix containing
             the multiplicies of the resolution given by rg, where the zeros
             corresponding to the strict transforms of the branches of the curve
             have been replaced by the (total) multiplicities of the infinitely near
             point corresponding to one further blow up for each branch.
KEYWORDS:    total multiplicities; multiplicity sequence; resolution graph
EXAMPLE:     example extend_multiplicities;   shows an example
"
{
  intmat resgr,tm,mt=rg[1],rg[2],rg[3];
  int i,j;
  for (i=1;i<=nrows(resgr);i++)
  {
    if (resgr[i,i]<0)
    {
      j=1;
      while ((resgr[i,j]==0) or (i==j))
      {
        j++;
      }
      tm[i,1..ncols(tm)]=tm[j,1..ncols(tm)];
      tm[i,-resgr[i,i]]=tm[i,-resgr[i,i]]+1;
      mt[i,-resgr[i,i]]=1;
    }
  }
  return(list(tm,mt));
}
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;
  // Calculate the resolution graph and the (total) multiplicities of f1*f2*f3.
  list RESGR=totalmultiplicities(f1*f2*f3);
  // Extend the (total) multiplicities.
  list MULT=extend_multiplicities(RESGR);
  // Compare the total multiplicities.
  RESGR[2];
  MULT[1];
  // Compare the multiplicities.
  RESGR[3];
  MULT[2];
}

proc intmat_inverse (intmat M)
"USAGE:      intmat_inverse(M); intmat M
ASSUME:      M is a lower triangular integer matrix with diagonal entries 1 or -1
RETURN:      intmat, the inverse of M
KEYWORDS:    integer matrix, inverse
EXAMPLE:     example intmat_inverse;   shows an example
"
{
  int i,j,k;
  int n=nrows(M);
  intmat U[n][n];
  U=U+1;
  for (i=1;i<=n;i++)
  {
    if (M[i,i]==-1)
    {
      M=addmultiplrows(M,i,i,-1,0);
      U=addmultiplrows(U,i,i,-1,0);
    }
    for (j=i+1;j<=n;j++)
    {
      U=addmultiplrows(U,i,j,-M[j,i],M[i,i]);
      M=addmultiplrows(M,i,j,-M[j,i],M[i,i]);
    }
  }
  return(U);
}
example
{
  "EXAMPLE:";echo=2;
  intmat M[5][5]=1,0,0,0,0,1,1,0,0,0,2,1,1,0,0,3,1,1,1,0,4,1,1,1,1 ;
  intmat U=intmat_inverse(M);
  print(U);
  print(U*M);
}