///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Singularities";
info="
LIBRARY:  surfacesignature.lib        signature of surface singularity

AUTHORS:  Gerhard Pfister             pfister@mathematik.uni-kl.de
@*        Muhammad Ahsan Banyamin     ahsanbanyamin@gmail.com
@*        Stefan Steidel              steidel@mathematik.uni-kl.de

OVERVIEW:

  A library for computing the signature of irreducible surface singularity.
  The signature of a surface singularity is defined in [3]. The algorithm we
  use has been proposed in [9].
  Let g in C[x,y] define an isolated curve singularity at 0 in C^2 and
  f:=z^N+g(x,y). The zero-set V:=V(f) in C^3 of f has an isolated singularity
  at 0. For a small e>0 let V_e:=V(f-e) in C^3 be the Milnor fibre of (V,0) and
  s: H_2(V_e,R) x H_2(V_e,R) ---> R be the intersection form (cf. [1],[7]).
  H_2(V_e,R) is an m-dimensional R-vector space, m the Milnor number of (V,0)
  (cf. [1],[4],[5],[6]), and s is a symmetric bilinear form.
  Let sigma(f) be the signature of s, called the signature of the surface
  singularity (V,0). Formulaes to compute the signature are given by Nemethi
  (cf. [8],[9]) and van Doorn, Steenbrink (cf. [2]).
  We have implemented three approaches using Puiseux expansions, the resolution
  of singularities resp. the spectral pairs of the singularity.

REFERENCES:

 [1] Arnold, V.I.; Gusein-Zade, S.M.; Varchenko, A.N.: Singularities of
     Differentiable Mappings. Vol. 1,2, Birkh\"auser (1988).
 [2] van Doorn, M.G.M.; Steenbrink, J.H.M.: A supplement to the monodromy
     theorem. Abh. Math. Sem. Univ. Hamburg 59, 225-233 (1989).
 [3] Durfee, A.H.: The Signature of Smoothings of Complex Surface
     Singularities. Mathematische Annalen 232, 85-98 (1978).
 [4] de Jong, T.; Pfister, G.: Local Analytic Geometry. Vieweg (2000).
 [5] Kerner, D.; Nemethi, A.: The Milnor fibre signature is not semi-continous.
     arXiv:0907.5252 (2009).
 [6] Kulikov, V.S.: Mixed Hodge Structures and Singularities. Cambridge Tracts
     in Mathematics 132, Cambridge University Press (1998).
 [7] Nemethi, A.: The real Seifert form and the spectral pairs of isolated
     hypersurface singularities. Compositio Mathematica 98, 23-41 (1995).
 [8] Nemethi, A.: Dedekind sums and the signature of f(x,y)+z^N. Selecta
     Mathematica, New series, Vol. 4, 361-376 (1998).
 [9] Nemethi, A.: The Signature of f(x,y)+z^$. Proceedings of Real and Complex
     Singularities (C.T.C. Wall's 60th birthday meeting, Liverpool (England),
     August 1996), London Math. Soc. Lecture Notes Series 263, 131--149 (1999).

PROCEDURES:
 signatureBrieskorn(a1,a2,a3);  signature of singularity x^a1+y^a2+z^a3
 signaturePuiseux(N,f);         signature of singularity z^N+f(x,y)=0, f irred.
 signatureNemethi(N,f);         signature of singularity z^N+f(x,y)=0
";

LIB "hnoether.lib";
LIB "alexpoly.lib";
LIB "gmssing.lib";

///////////////////////////////////////////////////////////////////////////////
//------- sigma(z^N + f) in terms of Puiseux pairs of f for f irreducible -----

static proc exponentSequence(poly f)
//=== computes the sequence a_1,...,a_s of exponents as described in [Nemethi]
//=== using the Puiseux pairs (m_1, n_1),...,(m_s, n_s) of f:
//===  - a_1 = m_1,
//===  - a_i = m_i - n_i * (m_[i-1] - n_[i-1] * a_[i-1]).
//===
//=== Return: list of two intvecs:
//===         1st entry: A = (a_1,...,a_s)
//===         2nd entry: N = (n_1,...,n_s)
{
   def R = basering;
   ring S = 0,(x,y),dp;
   poly f = fetch(R,f);
   list puiseuxPairs = invariants(f);
   setring R;

   intvec M = puiseuxPairs[1][3];
   intvec N = puiseuxPairs[1][4];

   int i;
   int a = M[1];
   intvec A = a;
   for(i = 2; i <= size(M); i++)
   {
      a = M[i] - N[i] * (M[i-1] - N[i-1] * a);
      A[size(A)+1] = a;
   }

   return(list(A,N));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y),dp;
   exponentSequence(y4+2x3y2+x6+x5y);
}

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

proc signatureBrieskorn(int a1, int a2, int a3)
"USAGE:  signatureBrieskorn(a1,a2,a3); a1,a2,a3 = integers
RETURN:  signature of Brieskorn singularity x^a1+y^a2+z^a3
EXAMPLE: example signatureBrieskorn; shows an example
"
{
   int a_temp, t, k1, k2, k3, s_t, sigma;
   number s;

   if(a1 > a2) { a_temp = a1; a1 = a2; a2 = a_temp; }
   if(a2 > a3) { a_temp = a2; a2 = a3; a3 = a_temp; }
   if(a1 > a2) { a_temp = a1; a1 = a2; a2 = a_temp; }

   for(t = 0; t <= 2; t++)
   {
      s_t = 0;
      for(k1 = 1; k1 <= a1-1; k1++)
      {
         for(k2 = 1; k2 <= a2-1; k2++)
         {
            for(k3 = 1; k3 <= a3-1; k3++)
            {
               s = number(k1)/a1 + number(k2)/a2 + number(k3)/a3;
               if(t < s)
               {
                  if(s < t+1)
                  {
                     s_t = s_t + 1;
                  }
                  else
                  {
                     break;
                  }
               }
            }
            if(k3 == 1) { break; }
         }
         if(k2 == 1) { break; }
      }
      sigma = sigma + (-1)^t * s_t;
   }
   return(sigma);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,x,dp;
   signatureBrieskorn(11,3,5);
}

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

proc signaturePuiseux(int N, poly f)
"USAGE:  signaturePuiseux(N,f); N = int, f = irreducible poly in 2 variables
RETURN:  signature of surface singularity defined by z^N + f(x,y) = 0
EXAMPLE: example signaturePuiseux; shows an example
"
{
   int i, d, prod, sigma;
   list L = exponentSequence(f);
   int s = size(L[2]);

   if(s == 1)
   {
      return(signatureBrieskorn(L[1][1], L[2][1], N));
   }

   prod = 1;
   sigma = signatureBrieskorn(L[1][s], L[2][s], N);
   for(i = s - 1; i >= 1; i--)
   {
      prod = prod * L[2][i+1];
      d = gcd(N, prod);
      sigma = sigma + d * signatureBrieskorn(L[1][i], L[2][i], N div d);
   }

   return(sigma);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y),dp;
   int N  = 3;
   poly f = x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2
            +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;
   signaturePuiseux(N,f);
}

///////////////////////////////////////////////////////////////////////////////
//------- sigma(z^N + f) in terms of the imbedded resolution graph of f -------

static proc dedekindSum(number b, number c, int a)
{
   number s,d,e;
   int k;
   for(k=1;k<=a-1;k++)
   {
      d=k*b mod a;
      e=k*c mod a;
      if(d*e!=0)
      {
         s=s+(d/a-1/2)*(e/a-1/2);
      }
   }
   return(s);
}

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

static proc isRupture(intvec v)
//=== decides whether the exceptional divisor given by the row v in the
//=== incidence matrix of the resolution graph intersects at least 3 other
//=== divisors
{
   int i,j;
   for(i=1;i<=size(v);i++)
   {
       if(v[i]<0){return(0);}
       if(v[i]!=0){j++;}
   }
   return(j>=4);
}

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

static proc sumExcepDiv(intmat N, list M, int K, int n)
//=== computes part of the formulae for eta(g,K), g defining an
//=== isolated curve singularity
//=== N the incidence matrix of the resolution graph of g
//=== M list of total multiplicities
//=== n = nrows(N)
{
   int i,j,m,d;
   for(i=1;i<=n;i++)
   {
      if(N[i,i]>0)
      {
         m=gcd(K,M[i]);
         for(j=1;j<=n;j++)
         {
            if((i!=j)&&(N[i,j]!=0))
            {
               if(m==1){break;}
               m=gcd(m,M[j]);
            }
         }
         d=d+m-1;
      }
   }
   return(d);
}

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

static proc sumEdges(intmat N, list M, int K, int n)
//=== computes part of the formulae for eta(g,K), g defining an
//=== isolated curve singularity
//=== N the incidence matrix of the resolution graph of g
//=== M list of total multiplicities
//=== n = nrows(N)
{
   int i,j,d;
   for(i=1;i<=n-1;i++)
   {
      for(j=i+1;j<=n;j++)
      {
         if(N[i,j]==1)
         {
            d=d+gcd(K,gcd(M[i],M[j]))-1;
         }
      }
   }
   return(d);
}

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

static proc etaRes(list L, int K)
//=== L total multiplicities
//=== eta-invariant in terms of the imbedded resolution graph of f
{
   int i,j,d;
   intvec v;
   number e;
   intmat N = L[1];         // incidence matrix of the resolution graph
   int n = ncols(L[1]);     // number of vertices in the resolution graph
   int a = ncols(L[2]);     // number of branches
   list M;                  // total multiplicities
   for(i=1;i<=n;i++)
   {
      d=L[2][i,1];
      for(j=2;j<=a;j++)
      {
         d=d+L[2][i,j];
      }
      if(d==0){d=1;}
      M[i]=d;
   }
   for(i=1;i<=n;i++)
   {
      v=N[i,1..n];
      if(isRupture(v))      // the divisor intersects more then two others
      {
         for(j=1;j<=n;j++)
         {
            if((i!=j)&&(v[j]!=0))
            {
               e=e+dedekindSum(M[j],K,M[i]);
            }
         }
      }
   }
   if(a==1)
   {
      //the irreducible case
      return(4*e);
   }
   return(a-1+4*e+sumEdges(N,M,K,n)-sumExcepDiv(N,M,K,n));
}

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

static proc signatureRes(int N, poly f)
//=== computes signature of surface singularity defined by z^N + f(x,y) = 0
//=== in terms of the imbedded resolution graph of f
{
   list L = totalmultiplicities(f);
   return(etaRes(L,N) - N*etaRes(L,1));
}

///////////////////////////////////////////////////////////////////////////////
//------------ sigma(z^N + f) in terms of the spectral pairs of f -------------

static proc fracPart(number n)
//=== computes the fractional part n2 of n
//=== i.e. n2 is not in Z but n-n2 is in Z
{
   number a,b;
   int r;
   a = numerator(n);
   b = denominator(n);
   int z = int(number(a));
   int y = int(number(b));
   r = z mod y;
   int q = (z-r) div y;
   number n1 = q;
   number n2 = n-n1;
   return(n2);
}

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

static proc etaSpec(list L, int N)
//=== L spectral numbers
//=== eta-invariant in terms of the spectral pairs of f
{
   int i;
   number e, h;

   int n = ncols(L[1]);

   if((n mod 2) == 0)
   // 0 is not a spectral number, thus f is irreducible
   {
      for(i = n div 2+1; i <= n; i++)
      {
         e = e + (1 - 2 * fracPart(N * number(L[1][i]))) * L[3][i];
      }
      return(2*e);
   }
   else
   // 0 is a spectral number, thus f is reducible
   {
      // sum of Hodge numbers in eta function
      for(i = 1; i <= n; i++)
      {
         if((L[2][i] == 2) && ((denominator(leadcoef(N*L[1][i]))==1)
                              ||(denominator(leadcoef(N*L[1][i]))==-1)))
         {
            h = h + L[3][i];
         }
      }

      // summand coming from spectral number 0 in eta function
      h = h + L[3][(n+1) div 2];

      // sum coming from non-zero spectral numbers in eta function
      for(i = (n+3) div 2; i <= n; i++)
      {
         if(!((denominator(leadcoef(N*L[1][i]))==1)
             ||(denominator(leadcoef(N*L[1][i]))==-1)))
         {
            e = e + (1 - 2 * fracPart(N * number(L[1][i]))) * L[3][i];
         }
      }
      return(h + 2*e);
   }
}

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

static proc signatureSpec(int N, poly f)
//=== computes signature of surface singularity defined by z^N + f(x,y) = 0
//=== in terms of the spectral pairs of f
{
   def R = basering;
   def Rds = changeord("ds");
   setring Rds;
   poly f = imap(R,f);
   list L = sppairs(f);
   setring R;
   list L = imap(Rds,L);
   return(etaSpec(L,N) - N*etaSpec(L,1));
}

///////////////////////////////////////////////////////////////////////////////
//----------------- Consolidation of the two recent variants ------------------

proc signatureNemethi(int N, poly f, list #)
"USAGE:  signatureNemethi(N,f); N = integer, f = reduced poly in 2 variables,
                                # empty or 1,2,3
@*       - if #[1] = 1 then resolution of singularity is used
@*       - if #[1] = 2 then spectral pairs are used
@*       - if # is empty then both upper variants are used in parallel and the
@*                 fastest returns the result
RETURN:  signature of surface singularity defined by z^N + f(x,y) = 0
REMARK:  computes the signature of some special surface singularities
EXAMPLE: example signatureNemethi; shows an example
"
{
   if(size(#) == 0)
   {
      link l(1) = "ssi:fork"; open(l(1));
      link l(2) = "ssi:fork"; open(l(2));
      list l = list(l(1),l(2));
      write(l(1), quote(signatureRes(N,f)));
      write(l(2), quote(signatureSpec(N,f)));
      int winner = waitfirst(l);
      number sigma = read(l(winner));
      close(l(1));
      close(l(2));
      if(printlevel >= 1)
      {
         if(winner == 1) { "Resolution of singularity has been used."; }
         if(winner == 2) { "Spectral pairs have been used."; }
      }
      return(sigma);
   }

   if(#[1] == 1)
   {
      return(signatureRes(N,f));
   }

   if(#[1] == 2)
   {
      return(signatureSpec(N,f));
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y),dp;
   int N  = 3;
   poly f = x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2
            +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;
   signatureNemethi(N,f,1);
   printlevel = 1;
   signatureNemethi(N,f);
}

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

/*
Further examples

ring r = 0,(x,y),dp;
int N;
poly f,g,g1,g2,g3;


// irreducible polynomials

N = 5;
f = x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2
    +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;
g = f^3 + x17y17;

N = 6;
f = y4+2x3y2+x6+x5y;
g1 = f^2 + x5y5;
g2 = f^3 + x11y11;
g3 = f^3 + x17y17;

N = 7;
f = x5+y11;
g1 = f^3 + x11y11;
g2 = f^3 + x17y17;

N = 6;
// k0 = 30, k1 = 35, k2 = 71
f = x71+6x65+15x59-630x52y6+20x53+6230x46y6+910x39y12+15x47
    -7530x40y6+14955x33y12-285x26y18+6x41+1230x34y6+4680x27y12
    +1830x20y18+30x13y24+x35-5x28y6+10x21y12-10x14y18+5x7y24-y30;

// k0 = 16, k1 = 24, k2 = 28, k3 = 30, k4 = 31
f = x31-781x30+16x29y-3010x29-2464x28y+104x27y2-2805x28-7024x27y
    -5352x26y2+368x25y3+366x27-7136x26y-984x25y2-8000x24y3
    +836x23y4+34x26-320x25y-6464x24y2+6560x23y3-8812x22y4+1392x21y5
    -12x25+256x24y-1296x23y2-1536x22y3+4416x21y4-8864x20y5+1752x19y6
    -x24+16x23y-88x22y2-16x21y3-404x20y4+3056x19y5-6872x18y6+1648x17y7
    +8x21y2-96x20y3+524x19y4-1472x18y5+3464x17y6-3808x16y7+1290x15y8
    -28x18y4+240x17y5-976x16y6+2208x15y7-2494x14y8+816x13y9+56x15y6
    -320x14y7+844x13y8-1216x12y9+440x11y10-70x12y8+240x11y9-344x10y10
    +240x9y11+56x9y10-96x8y11+52x7y12-28x6y12+16x5y13+8x3y14-y16;


// reducible polynomials

N = 12;
f = ((y2-x3)^2 - 4x5y - x7)*(x2-y3);

f = 2x3y3-2y5+x4-xy2;

f = -x3y3+x6y+xy6-x4y4;
*/