source: git/Singular/LIB/gaussman.lib @ 337919

///////////////////////////////////////////////////////////////////////////////
version="$Id: gaussman.lib,v 1.25 2001-01-11 17:12:58 Singular Exp $";
category="Singularities";

info="
LIBRARY:  gaussman.lib  Gauss-Manin Connection of a Singularity

AUTHOR:   Mathias Schulze, email: mschulze@mathematik.uni-kl.de

OVERVIEW: A library to compute invariants related to the Gauss-Manin connection
          of a an isolated hypersurface singularity.

PROCEDURES:
 monodromy(f[,...]);    monodromy matrix of f, spectrum of monodromy of f
 vfiltration(f[,...]);  V-filtration of f on H''/H', singularity spectrum of f
 vfiltjacalg(...);      V-filtration on Jacobian algebra
 gamma(...);            Hertling's gamma invariant
 gamma4(...);           Hertling's gamma4 invariant

SEE ALSO: mondromy_lib, spectrum_lib, jordan_lib

KEYWORDS: singularities; Gauss-Manin connection; monodromy; spectrum;
          Brieskorn lattice; Hodge filtration; V-filtration;
";

LIB "jordan.lib";

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

static proc maxintdif(ideal e)
{
  dbprint(printlevel-voice+2,"//gaussman::maxintdif");
  int i,j,id;
  int mid=0;
  for(i=ncols(e);i>=1;i--)
  {
    for(j=i-1;j>=1;j--)
    {
      id=int(e[i]-e[j]);
      if(id<0)
      {
        id=-id;
      }
      if(id>mid)
      {
        mid=id;
      }
    }
  }
  return(mid);
}
///////////////////////////////////////////////////////////////////////////////

static proc maxorddif(matrix H)
{
  dbprint(printlevel-voice+2,"//gaussman::maxorddif");
  int i,j,d;
  int d0,d1=-1,-1;
  for(i=nrows(H);i>=1;i--)
  {
    for(j=ncols(H);j>=1;j--)
    {
      d=ord(H[i,j]);
      if(d>=0)
      {
        if(d0<0||d<d0)
        {
          d0=d;
        }
        if(d1<0||d>d1)
        {
          d1=d;
        }
      }
    }
  }
  return(d1-d0);
}
///////////////////////////////////////////////////////////////////////////////

static proc invunit(poly u,int n)
{
  dbprint(printlevel-voice+2,"//gaussman::invunit");
  if(ord(u)>0)
  {
    ERROR("no unit");
  }
  poly u0=jet(u,0);
  u=jet(1-u/u0,n);
  poly w=u;
  poly v=1+u;
  for(int i=n div ord(u);i>1;i--)
  {
    w=jet(w*u,n);
    v=v+w;
  }
  v=jet(v,n)/u0;
  return(v);
}
///////////////////////////////////////////////////////////////////////////////

static proc expand(matrix M,matrix U,int n)
{
  dbprint(printlevel-voice+2,"//gaussman::expand");
  for(int i=ncols(U);i>=1;i--)
  {
    U[i,i]=invunit(U[i,i],n-ord(M[i]));
  }
  return(jet(M*U,n));
}
///////////////////////////////////////////////////////////////////////////////

static proc redNF(module M,module N,list #)
{
  matrix U=freemodule(ncols(M));
  if(size(#)>0)
  {
    if(typeof(#[1])=="matrix")
    {
      U=#[1];
    }
  }

  for(int i=ncols(M);i>=1;i--)
  {
    M[i]=M[i]/lead(U[i,i]);
    U[i,i]=U[i,i]/lead(U[i,i]);
  }
  N=std(N);
  module redNFM=matrix(0,nrows(M),ncols(M));
  dbprint(printlevel-voice+2,"//gaussman::redNF: reduce");
  module weakNFM=reduce(M,N);
  while(size(weakNFM)>0)
  {
    redNFM=matrix(redNFM)+matrix(lead(weakNFM));
    M=matrix(M)-matrix(lead(weakNFM))*U;
    dbprint(printlevel-voice+2,"//gaussman::redNF: reduce");
    weakNFM=reduce(M,N);
  }
  return(redNFM);
}
///////////////////////////////////////////////////////////////////////////////

proc monodromy(poly f,list #)
"USAGE:    monodromy(f[,mode]); poly f, int mode[=1]
ASSUME:   basering has local ordering, f has isolated singularity at 0
RETURN:
@format
          if mode=0:
            matrix M: exp(-2*pi*i*M) is a monodromy matrix of f
          if mode=1:
            ideal e: exp(-2*pi*i*e) is the spectrum of the monodromy of f
@end format
SEE ALSO: monodromy.lib, jordan.lib
KEYWORDS: singularities; Gauss-Manin connection; monodromy;
          Brieskorn lattice
EXAMPLE:  example monodromy; shows an example
"
{
  int mode=1;
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      mode=#[1];
    }
  }

  int i,j;
  int n=nvars(basering)-1;
  for(i=n+1;i>=1;i--)
  {
    if(var(i)>1)
    {
      ERROR("no local ordering");
    }
  }
  ideal J=jacob(f);
  ideal sJ=std(J);
  if(vdim(sJ)<=0)
  {
    if(vdim(sJ)==0)
    {
      ERROR("no singularity at 0");
    }
    else
    {
      ERROR("no isolated singularity at 0");
    }
  }
  ideal m=kbase(sJ);
  int mu,modm=ncols(m),maxorddif(m);

  ideal w=f*m;
  matrix U=freemodule(mu);
  matrix A0[mu][mu],A,C,D;
  list l;
  module H,dH=freemodule(mu),freemodule(mu);
  module H0;
  int sdH=1;
  int k=-1;
  int K,N,mide;

  while(k<K||sdH>0)
  {
    k++;
    dbprint(printlevel-voice+2,"//gaussman::monodromy: k="+string(k));

    dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
    C=coeffs(redNF(w,sJ,U),m);
    A0=A0+C*var(1)^k;

    if(sdH>0)
    {
      H0=H;
      dbprint(printlevel-voice+2,"//gaussman::monodromy: compute dH");
      dH=jet(module(A0*dH+var(1)^2*diff(matrix(dH),var(1))),k);
      H=H*var(1)+dH;

      dbprint(printlevel-voice+2,"//gaussman::monodromy: test dH==0");
      sdH=size(reduce(H,std(H0*var(1))));
      if(sdH>0)
      {
        A0=A0-var(1);
      }
      else
      {
        dbprint(printlevel-voice+2,
          "//gaussman::monodromy: compute basis of saturation");
        H=minbase(H0);
        int modH=maxorddif(H);
        K=modH+1;
      }
    }

    if(k==K&&sdH==0)
    {
      N=k-modH;
      dbprint(printlevel-voice+2,
        "//gaussman::monodromy: compute A on saturation");
      l=division(H*var(1),A0*H+var(1)^2*diff(matrix(H),var(1)));
      A=expand(l[1],l[2],N-1);
      if(mide==0)
      {
        dbprint(printlevel-voice+2,
          "//gaussman::monodromy: compute eigenvalues e and"+
          "multiplicities b of A");
        l=jordan(A);
        ideal e=l[1];
        intvec b;
        for(i=ncols(e);i>=1;i--)
        {
          b[i]=sum(l[2][i]);
        }
        dbprint(printlevel-voice+2,"//gaussman::monodromy: e="+string(e));
        dbprint(printlevel-voice+2,"//gaussman::monodromy: b="+string(b));
        if(mode==1)
        {
          mide=maxintdif(e);
          K=K+mide;
        }
      }
    }

    if(k<K||sdH>0)
    {
      dbprint(printlevel-voice+2,"//gaussman::monodromy: divide by J");
      l=division(J,ideal(matrix(w)-matrix(m)*C*U));
      D=l[1];

      dbprint(printlevel-voice+2,"//gaussman::monodromy: compute w/U");
      for(j=mu;j>=1;j--)
      {
        if(l[2][j,j]!=0)
        {
          dbprint(printlevel-voice+2,
            "//gaussman::monodromy: compute U["+string(j)+"]");
          U[j,j]=U[j,j]*l[2][j,j];
        }
        dbprint(printlevel-voice+2,
          "//gaussman::monodromy: compute w["+string(j)+"]");
        w[j]=0;
        for(i=n+1;i>=1;i--)
        {
          w[j]=w[j]+U[j,j]*diff(D[i,j],var(i))-diff(U[j,j],var(i))*D[i,j];
        }
      }
      U=U*U;
    }
  }

  if(mide>0)
  {
    intvec ide;
    ide[mu]=0;
    module dU;
    matrix A0e;
    for(i=ncols(e);i>=1;i--)
    {
      for(j=i-1;j>=1;j--)
      {
        k=int(e[j]-e[i]);
        if(k>ide[i])
        {
          ide[i]=k;
        }
        if(-k>ide[j])
        {
          ide[j]=-k;
        }
      }
    }
    for(j,k=ncols(e),mu;j>=1;j--)
    {
      for(i=b[j];i>=1;i--)
      {
        ide[k]=ide[j];
        k--;
      }
    }
  }
  while(mide>0)
  {
    dbprint(printlevel-voice+2,"//gaussman::monodromy: mide="+string(mide));

    U=0;
    A0=jet(A,0);
    for(i=ncols(e);i>=1;i--)
    {
      A0e=freemodule(mu);
      for(j=n;j>=0;j--) // Potenzen von Matrizen?
      {
        A0e=A0e*(A0-e[i]);
      }
      dU=syz(A0e);
      U=dU+U;
    }
    A=division(U,A*U)[1];

    for(i=mu;i>=1;i--)
    {
      for(j=mu;j>=1;j--)
      {
        if(ide[i]==0&&ide[j]>0)
        {
          A[i,j]=A[i,j]*var(1);
        }
        else
        {
          if(ide[i]>0&&ide[j]==0)
          {
            A[i,j]=A[i,j]/var(1);
          }
        }
      }
    }
    for(i=mu;i>=1;i--)
    {
      if(ide[i]>0)
      {
        A[i,i]=A[i,i]+1;
        e[i]=e[i]+1;
        ide[i]=ide[i]-1;
      }
    }
    mide--;
  }

  if(mode==1)
  {
    return(jet(A,0));
  }
  else
  {
    return(e);
  }
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  matrix M=monodromy(f);
  print(M);
}
///////////////////////////////////////////////////////////////////////////////

proc vfiltration(poly f,list #)
"USAGE:    vfiltration(f[,mode]); poly f, int mode[=1]
ASSUME:   basering has local ordering, f has isolated singularity at 0
RETURN:
@format
          list l:
          if mode=0 or mode=1:
            ideal l[1]: spectral numbers in increasing order
            intvec l[2]:
              int l[2][i]: multiplicity of spectral number l[1][i]
          if mode=1:
            list l[3]:
              module l[3][i]: vector space basis of l[1][i]-th graded part
                              of the V-filtration on H''/H' in terms of l[4]
            ideal l[4]: monomial vector space basis of H''/H'
            ideal l[5]: standard basis of the Jacobian ideal
@end format
NOTE:     H' and H'' denote the Brieskorn lattices
SEE ALSO: spectrum.lib
KEYWORDS: singularities; Gauss-Manin connection; spectrum;
          Brieskorn lattice; Hodge filtration; V-filtration;
EXAMPLE:  example vfiltration; shows an example
"
{
  int mode=1;
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      mode=#[1];
    }
  }

  int i,j;
  int n=nvars(basering)-1;
  for(i=n+1;i>=1;i--)
  {
    if(var(i)>1)
    {
      ERROR("no local ordering");
    }
  }
  ideal J=jacob(f);
  ideal sJ=std(J);
  if(vdim(sJ)<=0)
  {
    if(vdim(sJ)==0)
    {
      ERROR("no singularity at 0");
    }
    else
    {
      ERROR("no isolated singularity at 0");
    }
  }
  ideal m=kbase(sJ);
  int mu,modm=ncols(m),maxorddif(m);

  ideal w=f*m;
  matrix U=freemodule(mu);
  matrix A[mu][mu],C,D;
  list l;
  module H,dH=freemodule(mu),freemodule(mu);
  module H0;
  int sdH=1;
  int k=-1;
  int K;

  while(k<K||sdH>0)
  {
    k++;
    dbprint(printlevel-voice+2,"//gaussman::vfiltration: k="+string(k));

    dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute C");
    C=coeffs(redNF(w,sJ,U),m);
    A=A+C*var(1)^k;

    if(sdH>0)
    {
      H0=H;
      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute dH");
      dH=jet(module(A*dH+var(1)^2*diff(matrix(dH),var(1))),k);
      H=H*var(1)+dH;

      dbprint(printlevel-voice+2,"//gaussman::vfiltration: test dH==0");
      sdH=size(reduce(H,std(H0*var(1))));
      if(sdH>0)
      {
        A=A-var(1);
      }
      else
      {
        dbprint(printlevel-voice+2,
          "//gaussman::vfiltration: compute basis of saturation");
        H=minbase(H0);
        int modH=maxorddif(H);
        if(k<n)
        {
          K=modH+n+1;
        }
        else
        {
          K=modH+k+1;
        }
        H0=freemodule(mu)*var(1)^k;
      }
    }

    if(k<K||sdH>0)
    {
      dbprint(printlevel-voice+2,"//gaussman::vfiltration: divide by J");
      l=division(J,ideal(matrix(w)-matrix(m)*C*U));
      D=l[1];

      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute w/U");
      for(j=mu;j>=1;j--)
      {
        if(l[2][j,j]!=0)
        {
          dbprint(printlevel-voice+2,
            "//gaussman::vfiltration: compute U["+string(j)+"]");
          U[j,j]=U[j,j]*l[2][j,j];
        }
        dbprint(printlevel-voice+2,
          "//gaussman::vfiltration: compute w["+string(j)+"]");
        w[j]=0;
        for(i=n+1;i>=1;i--)
        {
          w[j]=w[j]+U[j,j]*diff(D[i,j],var(i))-diff(U[j,j],var(i))*D[i,j];
        }
      }
      U=U*U;
    }
  }
  int N=k-modH;

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: transform H0 to saturation");
  l=division(H,H0);
  H0=expand(l[1],l[2],N-1);

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: compute H0 as vector space V0");
  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: compute H1 as vector space V1");
  poly p;
  int i0,j0,i1,j1;
  matrix V0[mu*N][mu*N];
  matrix V1[mu*N][mu*(N-1)];
  for(i0=mu;i0>=1;i0--)
  {
    for(i1=mu;i1>=1;i1--)
    {
      p=H0[i1,i0];
      while(p!=0)
      {
        j1=leadexp(p)[1];
        for(j0=N-j1-1;j0>=0;j0--)
        {
          V0[i1+(j1+j0)*mu,i0+j0*mu]=V0[i1+(j1+j0)*mu,i0+j0*mu]+leadcoef(p);
          if(j1+j0+1<N)
          {
            V1[i1+(j1+j0+1)*mu,i0+j0*mu]=
            V1[i1+(j1+j0+1)*mu,i0+j0*mu]+leadcoef(p);
          }
        }
        p=p-lead(p);
      }
    }
  }

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: compute A on saturation");
  l=division(H*var(1),A*H+var(1)^2*diff(matrix(H),var(1)));
  A=expand(l[1],l[2],N-1);

  dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute matrix M of A");
  matrix M[mu*N][mu*N];
  for(i0=mu;i0>=1;i0--)
  {
    for(i1=mu;i1>=1;i1--)
    {
      p=A[i1,i0];
      while(p!=0)
      {
        j1=leadexp(p)[1];
        for(j0=N-j1-1;j0>=0;j0--)
        {
          M[i1+(j0+j1)*mu,i0+j0*mu]=leadcoef(p);
        }
        p=p-lead(p);
      }
    }
  }
  for(i0=mu;i0>=1;i0--)
  {
    for(j0=N-1;j0>=0;j0--)
    {
      M[i0+j0*mu,i0+j0*mu]=M[i0+j0*mu,i0+j0*mu]+j0;
    }
  }

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: compute eigenvalues eA of A");
  ideal eA=jordan(A,-1)[1];
  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eA="+string(eA));

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: compute eigenvalues eM of M");
  ideal eM;
  k=0;
  intvec u;
  for(i=N;i>=1;i--)
  {
    u[i]=1;
  }
  i0=1;
  while(u[N]<=ncols(eA))
  {
    for(i,i1=i0+1,i0;i<=N;i++)
    {
      if(eA[u[i]]+i<eA[u[i1]]+i1)
      {
        i1=i;
      }
    }
    k++;
    eM[k]=eA[u[i1]]+i1-1;
    u[i1]=u[i1]+1;
    if(u[i1]>ncols(eA))
    {
      i0=i1+1;
    }
  }
  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eM="+string(eM));

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: compute V-filtration on H0/H1");
  ideal s;
  k=0;
  list V;
  matrix Me;
  V[ncols(eM)+1]=std(V1);
  intvec d;
  if(mode==0)
  {
    for(i=ncols(eM);number(eM[i])-1>number(n-1)/2;i--)
    {
      Me=freemodule(mu*N);
      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
      {
        Me=Me*(M-eM[i]);
      }

      dbprint(printlevel-voice+2,
        "//gaussman::vfiltration: compute V["+string(i)+"]");
      V1=module(V1)+syz(Me);
      V[i]=std(intersect(V1,V0));

      if(size(V[i])>size(V[i+1]))
      {
        k++;
        s[k]=eM[i]-1;
        d[k]=size(V[i])-size(V[i+1]);
      }
    }

    dbprint(printlevel-voice+2,
      "//gaussman::vfiltration: symmetry index found");
    j=k;
    if(number(eM[i])-1==number(n-1)/2)
    {
      Me=freemodule(mu*N);
      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
      {
        Me=Me*(M-eM[i]);
      }

      dbprint(printlevel-voice+2,
        "//gaussman::vfiltration: compute V["+string(i)+"]");
      V1=module(V1)+syz(Me);
      V[i]=std(intersect(V1,V0));

      if(size(V[i])>size(V[i+1]))
      {
        k++;
        s[k]=eM[i]-1;
        d[k]=size(V[i])-size(V[i+1]);
      }
    }

    dbprint(printlevel-voice+2,"//gaussman::vfiltration: apply symmetry");
    while(j>=1)
    {
      k++;
      s[k]=s[j];
      s[j]=n-1-s[k];
      d[k]=d[j];
      j--;
    }

    return(list(s,d));
  }
  else
  {
    list v;
    j=-1;
    for(i=ncols(eM);i>=1;i--)
    {
      Me=freemodule(mu*N);
      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
      {
        Me=Me*(M-eM[i]);
      }

      dbprint(printlevel-voice+2,
        "//gaussman::vfiltration: compute V["+string(i)+"]");
      V1=module(V1)+syz(Me);
      V[i]=std(intersect(V1,V0));

      if(size(V[i])>size(V[i+1]))
      {
        if(number(eM[i]-1)>=number(n-1)/2)
        {
          k++;
          s[k]=eM[i]-1;
          dbprint(printlevel-voice+2,
            "//gaussman::vfiltration: transform to V0");
          v[k]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
        }
        else
        {
          if(j<0)
          {
            if(s[k]==number(n-1)/2)
            {
              j=k-1;
            }
            else
            {
              j=k;
            }
          }
          k++;
          s[k]=s[j];
          s[j]=eM[i]-1;
          v[k]=v[j];
          dbprint(printlevel-voice+2,
            "//gaussman::vfiltration: transform to V0");
          v[j]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
          j--;
        }
      }
    }

    dbprint(printlevel-voice+2,
      "//gaussman::vfiltration: compute graded parts");
    option(redSB);
    for(k=1;k<size(v);k++)
    {
      v[k]=std(reduce(v[k],std(v[k+1])));
      d[k]=size(v[k]);
    }
    v[k]=std(v[k]);
    d[k]=size(v[k]);

    return(list(s,d,v,m,sJ));
  }
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  list l=vfiltration(f);
  print(l);
}
///////////////////////////////////////////////////////////////////////////////

proc vfiltjacalg(list l)
"USAGE:   vfiltjacalg(vfiltration(f)); poly f
ASSUME:  basering has local ordering, f has isolated singularity at 0
RETURN:
@format
         list l:
           ideal l[1]: spectral numbers of the V-filtration
                       on the Jacobian algebra in increasing order
           intvec l[2]:
             int l[2][i]: multiplicity of spectral number l[1][i]
           list l[3]:
             module l[3][i]: vector space basis of the l[1][i]-th graded part
                             of the V-filtration on the Jacobian algebra
                             in terms of l[4]
           ideal l[4]: monomial vector space basis of the Jacobian algebra
           ideal l[5]: standard basis of the Jacobian ideal
@end format
EXAMPLE: example vfiltjacalg; shows an example
"
{
  def s,d,v,m,sJ=l[1..5];
  int mu=ncols(m);

  int i,j,k;
  module V=v[1];
  for(i=2;i<=size(v);i++)
  {
    V=V,v[i];
  }

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltjacalg: compute multiplication in Jacobian algebra");
  list M;
  for(i=ncols(m);i>=1;i--)
  {
    M[i]=lift(V,coeffs(redNF(m[i]*m,sJ),m)*V);
  }

  int i0,j0,i1,j1;
  number r0=number(s[1]-s[ncols(s)]);
  number r1,r2;
  matrix M0;
  module L;
  list v0=freemodule(ncols(m));
  ideal s0=r0;

  while(r0<number(s[ncols(s)]-s[1]))
  {
    dbprint(printlevel-voice+2,
      "//gaussman::vfiltjacalg: find next possible index");
    r1=number(s[ncols(s)]-s[1]);
    for(j=ncols(s);j>=1;j--)
    {
      for(i=ncols(s);i>=1;i--)
      {
        r2=number(s[i]-s[j]);
        if(r2>r0&&r2<r1)
        {
          r1=r2;
        }
        else
        {
          if(r2<=r0)
          {
            i=0;
          }
        }
      }
    }
    r0=r1;

    l=ideal();
    for(k=ncols(m);k>=2;k--)
    {
      l=l+list(ideal());
    }

    dbprint(printlevel-voice+2,
      "//gaussman::vfiltjacalg: collect conditions for V["+string(r0)+"]");
    j=ncols(s);
    j0=mu;
    while(j>=1)
    {
      i0=1;
      i=1;
      while(i<ncols(s)&&s[i]<s[j]+r0)
      {
        i0=i0+d[i];
        i++;
      }
      if(s[i]<s[j]+r0)
      {
        i0=i0+d[i];
        i++;
      }
      for(k=1;k<=ncols(m);k++)
      {
        M0=M[k];
        if(i0>1)
        {
          l[k]=l[k],M0[1..i0-1,j0-d[j]+1..j0];
        }
      }
      j0=j0-d[j];
      j--;
    }

    dbprint(printlevel-voice+2,
      "//gaussman::vfiltjacalg: compose condition matrix");
    L=transpose(module(l[1]));
    for(k=2;k<=ncols(m);k++)
    {
      L=L,transpose(module(l[k]));
    }

    dbprint(printlevel-voice+2,
      "//gaussman::vfiltjacalg: compute kernel of condition matrix");
    v0=v0+list(syz(L));
    s0=s0,r0;
  }

  dbprint(printlevel-voice+2,"//gaussman::vfiltjacalg: compute graded parts");
  option(redSB);
  for(i=1;i<size(v0);i++)
  {
    v0[i+1]=std(v0[i+1]);
    v0[i]=std(reduce(v0[i],v0[i+1]));
  }

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltjacalg: remove trivial graded parts");
  i=1;
  while(size(v0[i])==0)
  {
    i++;
  }
  list v1=v0[i];
  intvec d1=size(v0[i]);
  ideal s1=s0[i];
  i++;
  while(i<=size(v0))
  {
    if(size(v0[i])>0)
    {
      v1=v1+list(v0[i]);
      d1=d1,size(v0[i]);
      s1=s1,s0[i];
    }
    i++;
  }
  return(list(s1,d1,v1,m));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  vfiltjacalg(vfiltration(f));
}
///////////////////////////////////////////////////////////////////////////////

proc gamma(list l)
"USAGE:   gamma(vfiltration(f,0)); poly f
ASSUME:  basering has local ordering, f has isolated singularity at 0
RETURN:  number g: Hertling's gamma invariant
EXAMPLE: example gamma; shows an example
"
{
  ideal s=l[1];
  intvec d=l[2];
  int n=nvars(basering)-1;
  number g=0;
  int i,j;
  for(i=1;i<=ncols(s);i++)
  {
    for(j=1;j<=d[i];j++)
    {
      g=g+(number(s[i])-number(n-1)/2)^2;
    }
  }
  g=-g/4+sum(d)*number(s[ncols(s)]-s[1])/48;
  return(g);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  gamma(vfiltration(f,0));
}
///////////////////////////////////////////////////////////////////////////////

proc gamma4(list l)
"USAGE:   gamma4(vfiltration(f,0)); poly f
ASSUME:  basering has local ordering, f has isolated singularity at 0
RETURN:  number g4: Hertling's gamma4 invariant
EXAMPLE: example gamma4; shows an example
"
{
  ideal s=l[1];
  intvec d=l[2];
  int n=nvars(basering)-1;
  number g4=0;
  int i,j;
  for(i=1;i<=ncols(s);i++)
  {
    for(j=1;j<=d[i];j++)
    {
      g4=g4+(number(s[i])-number(n-1)/2)^4;
    }
  }
  g4=g4-(number(s[ncols(s)]-s[1])/12-1/30)*
    (sum(d)*number(s[ncols(s)]-s[1])/4-24*gamma(l));
 return(g4);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  gamma4(vfiltration(f,0));
}
///////////////////////////////////////////////////////////////////////////////

proc tst_gaussm(poly f)
{
  echo=2;
  basering;
  f;
  print(monodromy(f));
  list l=vfiltration(f);
  l;
  vfiltjacalg(l);
  gamma(l);
  gamma4(l);
}
///////////////////////////////////////////////////////////////////////////////