source: git/Singular/LIB/gaussman.lib @ 0fcdf0

///////////////////////////////////////////////////////////////////////////////
version="$Id: gaussman.lib,v 1.43 2001-05-30 18:32:37 mschulze 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:
 gmsring(f,s);              Brieskorn lattice in the Gauss-Manin system
 gmsnf(p,K[,Kmax]);         Gauss-Manin system normal form
 gmscoeffs(p,K[,Kmax]);     Gauss-Manin system basis representation
 monodromy(f[,opt]);        monodromy matrix, spectrum of monodromy of f
 vfiltration(f[,opt]);      V-filtration on H''/H', singularity spectrum of f
 spectrum(f);               singularity spectrum of f
 endfilt(poly f,list V);    endomorphism filtration of filtration V
 spprint(list S);           print spectrum S
 spadd(list S1,list S2);    sum of spectra S1 and S2
 spsub(list S1,list S2);    difference of spectra S1 and S2
 spmul(list S,int k);       product of spectrum S and integer k
 spmul(list S,intvec k);    linear combination of spectra S with coefficients k
 spissemicont(list S[,opt]);        test spectrum S for semicontinuity
 spsemicont(list S0,list S[,opt]);  relative semicontinuity of spectra S0 and S
 spmilnor(list S);          milnor number of spectrum S
 spgeomgenus(list S);       geometrical genus of spectrum S
 spgamma(list S);           gamma invariant of spectrum S

SEE ALSO: mondromy_lib, spectrum_lib

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

LIB "linalg.lib";

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

static proc stdtrans(ideal I)
{
  def R=basering;

  string s=ordstr(R);
  int j=find(s,",C");
  if(j==0)
  {
    j=find(s,"C,");
  }
  if(j==0)
  {
    j=find(s,",c");
  }
  if(j==0)
  {
    j=find(s,"c,");
  }
  if(j>0)
  {
    s[j..j+1]="  ";
  }

  execute("ring @S="+charstr(R)+",(gmspoly,"+varstr(R)+"),(c,dp,"+s+");");

  ideal I=homog(imap(R,I),gmspoly);
  module M=transpose(transpose(I)+freemodule(ncols(I)));
  M=std(M);

  setring(R);
  execute("map h=@S,1,"+varstr(R)+";");

  module M=h(M);

  for(int i=ncols(M);i>=1;i--)
  {
    for(j=ncols(M);j>=1;j--)
    {
      if(M[i][1]==0)
      {
        M[i]=0;
      }
      if(i!=j&&M[j][1]!=0)
      {
        if(lead(M[i][1])/lead(M[j][1])!=0)
        {
          M[i]=0;
        }
      }
    }
  }

  M=transpose(simplify(M,2));
  I=M[1];
  attrib(I,"isSB",1);
  M=M[2..ncols(M)];
  module U=transpose(M);

  return(list(I,U));
}
///////////////////////////////////////////////////////////////////////////////

proc gmsring(poly t,string s)
"USAGE:    gmsring(f,s); poly f, string s;
ASSUME:   basering has characteristic 0 and local degree ordering,
          f has isolated singularity at 0
RETURN:
@format
ring G:
  G=C{{s}}*basering,
  s is the inverse of the Gauss-Manin connection of f,
  G contains:
    poly gmspoly: image of f
    ideal gmsjacob: image of Jacobian ideal
    ideal gmsstd: image of standard basis of Jacobian ideal
    matrix gmsmatrix: matrix(gmsjacob)*gmsmatrix=matrix(gmsstd)
    ideal gmsbasis: image of monomial vector space basis of Jacobian algebra
    int gmsmaxweight: maximal weight of variables of basering
  G projects to H''=C{{s}}*gmsbasis
@end format
NOTE:     do not modify gms variables if you want to use gms procedures
KEYWORDS: singularities; Gauss-Manin connection
EXAMPLE:  example gms; shows examples
"
{
  def R=basering;
  if(charstr(R)!="0")
  {
    ERROR("characteristic 0 expected");
  }
  for(int i=nvars(R);i>=1;i--)
  {
    if(var(i)>1)
    {
      ERROR("local ordering expected");
    }
  }

  ideal dt=jacob(t);
  list l=stdtrans(dt);
  ideal g=l[1];
  if(vdim(g)<=0)
  {
    if(vdim(g)==0)
    {
      ERROR("singularity at 0 expected");
    }
    else
    {
      ERROR("isolated singularity at 0 expected");
    }
  }  
  matrix a=l[2];
  ideal m=kbase(g);

  intvec w;
  int gmsmaxweight;
  for(i=nvars(R);i>=1;i--)
  {
    w[i+1]=deg(var(i));
    if(deg(var(i))>gmsmaxweight)
    {
      gmsmaxweight=deg(var(i));
    }
  }
  w[1]=deg(highcorner(g))+2*gmsmaxweight;

  execute("ring G="+string(charstr(R))+",("+s+","+varstr(R)+"),ws("+
    string(w)+");");

  poly gmspoly=imap(R,t);
  ideal gmsjacob=imap(R,dt);
  ideal gmsstd=imap(R,g);
  matrix gmsmatrix=imap(R,a);
  ideal gmsbasis=imap(R,m);

  attrib(gmsstd,"isSB",1);
  export gmspoly,gmsjacob,gmsstd,gmsmatrix,gmsbasis,gmsmaxweight;

  return(G);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  def G=gmsring(f,"s");
  setring(G);
  gmspoly;
  print(gmsjacob);
  print(gmsstd);
  print(gmsmatrix);
  gmsmaxweight;
}
///////////////////////////////////////////////////////////////////////////////

proc gmsnf(ideal p,int K,list #)
"USAGE:    gmsnf(p,K[,Kmax]); poly p, int K[, int Kmax];
ASSUME:   basering was constructed by gms, K<=Kmax
RETURN:
@format
list l:
  ideal l[1]: projection of p to H''=C{{s}}*gmsbasis mod s^{K+1}
  ideal l[2]: p=l[1]+l[2] mod s^(Kmax+1)
@end format
NOTE:     by setting p=l[2] the computation can be continued up to order
          at most Kmax, by default Kmax=infinity
KEYWORDS: singularities; Gauss-Manin connection
EXAMPLE:  example gmsnf; shows examples
"
{
  int Kmax=-1;
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      Kmax=#[1];
      if(K>Kmax)
      {
        Kmax=K;
      }
    }
  }

  intvec v=1;
  v[nvars(basering)]=0;

  int k;
  if(Kmax>=0)
  {
    p=jet(jet(p,K,v),(Kmax+1)*deg(var(1))-2*gmsmaxweight);
  }

  ideal r,q;
  r[ncols(p)]=0;
  q[ncols(p)]=0;

  poly s;
  int i,j;
  for(k=ncols(p);k>=1;k--)
  {
    while(p[k]!=0&&deg(lead(p[k]),v)<=K)
    {
      i=1;
      s=lead(p[k])/lead(gmsstd[i]);
      while(i<ncols(gmsstd)&&s==0)
      {
        i++;
        s=lead(p[k])/lead(gmsstd[i]);
      }
      if(s!=0)
      {
        p[k]=p[k]-s*gmsstd[i];
        for(j=1;j<=nrows(gmsmatrix);j++)
        {
          if(Kmax>=0)
          {
            p[k]=p[k]+
              jet(jet(diff(s*gmsmatrix[j,i],var(j+1))*var(1),Kmax,v),
                (Kmax+1)*deg(var(1))-2*gmsmaxweight);
          }
          else
          {
            p[k]=p[k]+diff(s*gmsmatrix[j,i],var(j+1))*var(1);
          }
        }
      }
      else
      {
        r[k]=r[k]+lead(p[k]);
        p[k]=p[k]-lead(p[k]);
      }
      while(deg(lead(p[k]))>(K+1)*deg(var(1))-2*gmsmaxweight&&
        deg(lead(p[k]),v)<=K)
      {
        q[k]=q[k]+lead(p[k]);
        p[k]=p[k]-lead(p[k]);
      }
    }
    q[k]=q[k]+p[k];
  }

  return(list(r,q));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  def G=gmsring(f,"s");
  setring(G);
  list l0=gmsnf(gmspoly,0);
  print(l0[1]);
  list l1=gmsnf(gmspoly,1);
  print(l1[1]);
  list l=gmsnf(l0[2],1);
  print(l[1]);
}
///////////////////////////////////////////////////////////////////////////////

proc gmscoeffs(ideal p,int K,list #)
"USAGE:    gmscoeffs(p,K[,Kmax]); poly p, int K[, int Kmax];
ASSUME:   basering was constructed by gms, K<=Kmax
RETURN:
@format
list l:
  matrix l[1]: projection of p to H''=C{{s}}*gmsbasis=C{{s}}^mu mod s^(K+1)
  ideal l[2]: p=matrix(gmsbasis)*l[1]+l[2] mod s^(Kmax+1)
@end format
NOTE:     by setting p=l[2] the computation can be continued up to order 
          at most Kmax, by default Kmax=infinity
KEYWORDS: singularities; Gauss-Manin connection
EXAMPLE:  example gmscoeffs; shows examples
"
{
  list l=gmsnf(p,K,#);
  ideal r,q=l[1..2];
  poly v=1;
  for(int i=2;i<=nvars(basering);i++)
  {
    v=v*var(i);
  }
  matrix C=coeffs(r,gmsbasis,v);
  return(C,q);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  def G=gmsring(f,"s");
  setring(G);
  list l0=gmscoeffs(gmspoly,0);
  print(l0[1]);
  list l1=gmscoeffs(gmspoly,1);
  print(l1[1]);
  list l=gmscoeffs(l0[2],1);
  print(l[1]);
}
///////////////////////////////////////////////////////////////////////////////

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);
}
///////////////////////////////////////////////////////////////////////////////

proc monodromy(poly f,list #)
"USAGE:    monodromy(f[,opt]); poly f, int opt
ASSUME:   basering has characteristic 0 and local degree ordering,
          f has isolated singularity at 0
RETURN:
@format
if opt==0:
  matrix M: exp(-2*pi*i*M) is a monodromy matrix of f
if opt==1:
  ideal e: exp(-2*pi*i*e) are the eigenvalues of the monodromy of f
default: opt=1
@end format
SEE ALSO: mondromy_lib
KEYWORDS: singularities; Gauss-Manin connection; monodromy
EXAMPLE:  example monodromy; shows examples
"
{
  int opt=1;
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      opt=#[1];
    }
  }

  def R=basering;
  def G=gmsring(f,"s");
  setring G;

  int n=nvars(R)-1;
  int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis);
  ideal w=gmspoly*gmsbasis;
  list l;
  matrix U=freemodule(mu);
  matrix A0[mu][mu],A,C;
  module H,dH=freemodule(mu),freemodule(mu);
  module H0;
  int sdH=1;
  int k=-1;

  while(sdH>0)
  {
    k++;
    dbprint(printlevel-voice+2,"//gaussman::monodromy: k="+string(k));
    dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
    if(opt==0)
    {
      l=gmscoeffs(w,k,mu);
    }
    else
    {
      l=gmscoeffs(w,k,mu+n);
    }
    C,w=l[1..2];
    A0=A0+C;

    H0=H;
    dbprint(printlevel-voice+2,"//gaussman::monodromy: compute dH");
    dH=jet(module(A0*dH+s^2*diff(matrix(dH),s)),k);
    H=H*s+dH;

    dbprint(printlevel-voice+2,"//gaussman::monodromy: test dH==0");
    sdH=size(reduce(H,std(H0*s)));
  }

  A0=A0-s^k;

  dbprint(printlevel-voice+2,
    "//gaussman::monodromy: compute basis of saturation");
  H=minbase(H0);
  int modH=maxorddif(H)/deg(s);
  dbprint(printlevel-voice+2,"//gaussman::monodromy: k="+string(modH+1));
  dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
  if(opt==0)
  {
    l=gmscoeffs(w,modH+1,modH+1);
  }
  else
  {
    l=gmscoeffs(w,modH+1,modH+1+n);
  }
  C,w=l[1..2];
  A0=A0+C;

  dbprint(printlevel-voice+2,
    "//gaussman::monodromy: compute A on saturation");
  l=division(H*s,A0*H+s^2*diff(matrix(H),s));
  A=jet(l[1],l[2],0);

  dbprint(printlevel-voice+2,
    "//gaussman::monodromy: compute eigenvalues e and "+
    "multiplicities b of A");
  l=eigenval(jet(A,0));
  ideal e=l[1];
  intvec b=l[2];
  dbprint(printlevel-voice+2,"//gaussman::monodromy: e="+string(e));
  dbprint(printlevel-voice+2,"//gaussman::monodromy: b="+string(b));

  if(opt==0)
  {
    setring(R);
    ideal e=imap(G,e);
    return(e);
  }

  int mide=maxintdif(e);

  if(mide>0)
  {
    dbprint(printlevel-voice+2,
      "//gaussman::monodromy: k="+string(modH+1+mide));
    dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
    l=gmscoeffs(w,modH+1+mide,modH+1+mide);
    C,w=l[1..2];
    A0=A0+C;

    intvec ide;
    ide[mu]=0;
    int i,j;
    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));

    matrix U=0;
    A0=jet(A,0);
    for(i=ncols(e);i>=1;i--)
    {
      U=syz(power(A0-e[i],n+1))+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]*d;
        }
        else
        {
          if(ide[i]>0&&ide[j]==0)
          {
            A[i,j]=A[i,j]/d;
          }
        }
      }
    }
    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--;
  }
  A=jet(A,0);

  setring(R);
  matrix A=imap(G,A);
  return(A);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  print(monodromy(f));
}
///////////////////////////////////////////////////////////////////////////////

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

  def R=basering;
  def G=gmsring(f,"s");
  setring G;

  int n=nvars(R)-1;
  int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis);
  ideal w=gmspoly*gmsbasis;
  list l;
  matrix U=freemodule(mu);
  matrix A[mu][mu],C;
  module H,dH=freemodule(mu),freemodule(mu);
  module H0;
  int sdH=1;
  int k=-1;
  int N=n+1;

  while(sdH>0)
  {
    k++;
    dbprint(printlevel-voice+2,"//gaussman::vfiltration: k="+string(k));
    dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute C");
    l=gmscoeffs(w,k);
    C,w=l[1..2];
    A=A+C;

    H0=H;
    dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute dH");
    dH=jet(module(A*dH+s^2*diff(matrix(dH),s)),k);
    H=H*s+dH;

    dbprint(printlevel-voice+2,"//gaussman::vfiltration: test dH==0");
    sdH=size(reduce(H,std(H0*s)));
  }

  A=A-s^k;

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: compute basis of saturation");
  H=minbase(H0);
  int modH=maxorddif(H)/deg(s);
  dbprint(printlevel-voice+2,"//gaussman::monodromy: k="+string(N+modH));
  dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
  l=gmscoeffs(w,N+modH,N+modH);
  C,w=l[1..2];
  A=A+C;

  dbprint(printlevel-voice+2,
    "//gaussman::vfiltration: transform H0 to saturation");
  l=division(H,freemodule(mu)*s^k);
  H0=jet(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*s,A*H+s^2*diff(matrix(H),s));
  A=jet(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=eigenval(jet(A,0))[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(int 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 a;
  k=0;
  list V;
  V[ncols(eM)+1]=interred(V1);
  intvec d;
  if(opt==0)
  {
    for(i=ncols(eM);number(eM[i])-1>number(n-1)/2;i--)
    {
      dbprint(printlevel-voice+2,
        "//gaussman::vfiltration: compute V["+string(i)+"]");
      V1=module(V1)+syz(power(M-eM[i],n+1));
      V[i]=interred(intersect(V1,V0));

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

    dbprint(printlevel-voice+2,
      "//gaussman::vfiltration: symmetry index found");
    int j=k;

    if(number(eM[i])-1==number(n-1)/2)
    {
      dbprint(printlevel-voice+2,
        "//gaussman::vfiltration: compute V["+string(i)+"]");
      V1=module(V1)+syz(power(M-eM[i],n+1));
      V[i]=interred(intersect(V1,V0));

      if(size(V[i])>size(V[i+1]))
      {
        k++;
        a[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++;
      a[k]=a[j];
      a[j]=n-1-a[k];
      d[k]=d[j];
      j--;
    }

    setring(R);
    ideal a=imap(G,a);
    return(list(a,d));
  }
  else
  {
    list v;
    int j=-1;
    for(i=ncols(eM);i>=1;i--)
    {
      dbprint(printlevel-voice+2,
        "//gaussman::vfiltration: compute V["+string(i)+"]");
      V1=module(V1)+syz(power(M-eM[i],n+1));
      V[i]=interred(intersect(V1,V0));

      if(size(V[i])>size(V[i+1]))
      {
        if(number(eM[i]-1)>=number(n-1)/2)
        {
          k++;
          a[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(a[k]==number(n-1)/2)
            {
              j=k-1;
            }
            else
            {
              j=k;
            }
          }
          k++;
          a[k]=a[j];
          a[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");
    for(k=1;k<size(v);k++)
    {
      v[k]=interred(reduce(v[k],std(v[k+1])));
      d[k]=size(v[k]);
    }
    v[k]=interred(v[k]);
    d[k]=size(v[k]);

    setring(R);
    ideal a=imap(G,a);
    list v=imap(G,v);
    ideal m=imap(G,gmsbasis);
    return(list(a,d,v,m));
  }
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  vfiltration(f);
}
///////////////////////////////////////////////////////////////////////////////

proc spectrum(poly f)
"USAGE:    spectrum(f); poly f
ASSUME:   basering has characteristic 0 and local degree ordering,
          f has isolated singularity at 0
RETURN:
@format
list S: singularity spectrum of f
  ideal S[1]: spectral numbers in increasing order
  intvec S[2]:
    int S[2][i]: multiplicity of spectral number S[1][i]
@end format
SEE ALSO: spectrum_lib
KEYWORDS: singularities; Gauss-Manin connection; spectrum
EXAMPLE:  example spectrum; shows examples
"
{
  return(vfiltration(f,0));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly f=x5+x2y2+y5;
  spprint(spectrum(f));
}
///////////////////////////////////////////////////////////////////////////////

proc endfilt(poly f,list V)
"USAGE:   endfilt(f,V); poly f, list V
ASSUME:  basering has characteristic 0 and local degree ordering,
         f has isolated singularity at 0
RETURN:
@format
list V1: endomorphim filtration of V on the Jacobian algebra of f
  ideal V1[1]: spectral numbers in increasing order
  intvec V1[2]:
    int V1[2][i]: multiplicity of spectral number V1[1][i]
  list V1[3]:
    module V1[3][i]: vector space basis of the V1[1][i]-th graded part
                     in terms of V1[4]
  ideal V1[4]: monomial vector space basis
@end format
SEE ALSO: spectrum_lib
KEYWORDS: singularities; Gauss-Manin connection; spectrum;
          Brieskorn lattice; Hodge filtration; V-filtration
EXAMPLE: example endfilt; shows examples
"
{
  if(charstr(basering)!="0")
  {
    ERROR("characteristic 0 expected");
  }
  int n=nvars(basering)-1;
  for(int i=n+1;i>=1;i--)
  {
    if(var(i)>1)
    {
      ERROR("local ordering expected");
    }
  }
  ideal sJ=std(jacob(f));
  if(vdim(sJ)<=0)
  {
    if(vdim(sJ)==0)
    {
      ERROR("singularity at 0 expected");
    }
    else
    {
      ERROR("isolated singularity at 0 expected");
    }
  }

  def a,d,v,m=V[1..4];
  int mu=ncols(m);

  module V0=v[1];
  for(i=2;i<=size(v);i++)
  {
    V0=V0,v[i];
  }

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

  int j,k,i0,j0,i1,j1;
  number b0=number(a[1]-a[ncols(a)]);
  number b1,b2;
  matrix M0;
  module L;
  list v0=freemodule(ncols(m));
  ideal a0=b0;

  while(b0<number(a[ncols(a)]-a[1]))
  {
    dbprint(printlevel-voice+2,
      "//gaussman::endfilt: find next possible index");
    b1=number(a[ncols(a)]-a[1]);
    for(j=ncols(a);j>=1;j--)
    {
      for(i=ncols(a);i>=1;i--)
      {
        b2=number(a[i]-a[j]);
        if(b2>b0&&b2<b1)
        {
          b1=b2;
        }
        else
        {
          if(b2<=b0)
          {
            i=0;
          }
        }
      }
    }
    b0=b1;

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

    dbprint(printlevel-voice+2,
      "//gaussman::endfilt: collect conditions for V1["+string(b0)+"]");
    j=ncols(a);
    j0=mu;
    while(j>=1)
    {
      i0=1;
      i=1;
      while(i<ncols(a)&&a[i]<a[j]+b0)
      {
        i0=i0+d[i];
        i++;
      }
      if(a[i]<a[j]+b0)
      {
        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::endfilt: 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::endfilt: compute kernel of condition matrix");
    v0=v0+list(syz(L));
    a0=a0,b0;
  }

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

proc spprint(list S)
"USAGE:   spprint(S); list S
RETURN:  string: spectrum S
EXAMPLE: example spprint; shows examples
"
{
  string s;
  for(int i=1;i<size(S[2]);i++)
  {
    s=s+"("+string(S[1][i])+","+string(S[2][i])+"),";
  }
  s=s+"("+string(S[1][i])+","+string(S[2][i])+")";
  return(s);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S);
}
///////////////////////////////////////////////////////////////////////////////

proc spadd(list S1,list S2)
"USAGE:   spadd(S1,S2); list S1,S2
RETURN:  list: sum of spectra S1 and S2
EXAMPLE: example spadd; shows examples
"
{
  ideal s;
  intvec m;
  int i,i1,i2=1,1,1;
  while(i1<=size(S1[2])||i2<=size(S2[2]))
  {
    if(i1<=size(S1[2]))
    {
      if(i2<=size(S2[2]))
      {
        if(number(S1[1][i1])<number(S2[1][i2]))
        {
          s[i]=S1[1][i1];
          m[i]=S1[2][i1];
          i++;
          i1++;
        }
        else
        {
          if(number(S1[1][i1])>number(S2[1][i2]))
          {
            s[i]=S2[1][i2];
            m[i]=S2[2][i2];
            i++;
            i2++;
          }
          else
          {
            if(S1[2][i1]+S2[2][i2]!=0)
            {
              s[i]=S1[1][i1];
              m[i]=S1[2][i1]+S2[2][i2];
              i++;
            }
            i1++;
            i2++;
          }
        }
      }
      else
      {
        s[i]=S1[1][i1];
        m[i]=S1[2][i1];
        i++;
        i1++;
      }
    }
    else
    {
      s[i]=S2[1][i2];
      m[i]=S2[2][i2];
      i++;
      i2++;
    }
  }
  return(list(s,m));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S1);
  list S2=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(S2);
  spprint(spadd(S1,S2));
}
///////////////////////////////////////////////////////////////////////////////

proc spsub(list S1,list S2)
"USAGE:   spsub(S1,S2); list S1,S2
RETURN:  list: difference of spectra S1 and S2
EXAMPLE: example spsub; shows examples
"
{
  return(spadd(S1,spmul(S2,-1)));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S1);
  list S2=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(S2);
  spprint(spsub(S1,S2));
}
///////////////////////////////////////////////////////////////////////////////

proc spmul(list #)
"USAGE:
@format
1) spmul(S,k); list S, int k
2) spmul(S,k); list S, intvec k
@end format
RETURN:
@format
1) list: product of spectrum S and integer k
2) list: linear combination of spectra S with coefficients k
@end format
EXAMPLE: example spmul; shows examples
"
{
  if(size(#)==2)
  {
    if(typeof(#[1])=="list")
    {
      if(typeof(#[2])=="int")
      {
        return(list(#[1][1],#[1][2]*#[2]));
      }
      if(typeof(#[2])=="intvec")
      {
        list S0=list(ideal(),intvec(0));
        for(int i=size(#[2]);i>=1;i--)
        {
          S0=spadd(S0,spmul(#[1][i],#[2][i]));
        }
        return(S0);
      }
    }
  }
  return(list(ideal(),intvec(0)));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S);
  spprint(spmul(S,2));
  list S1=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(S1);
  list S2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
  spprint(S2);
  spprint(spmul(list(S1,S2),intvec(1,2)));
}
///////////////////////////////////////////////////////////////////////////////

proc spissemicont(list S,list #)
"USAGE:   spissemicont(S[,opt]); list S, int opt
RETURN:
@format
int k=
if opt==0:
  1, if sum of spectrum S over all intervals [a,a+1) is positive
  0, if sum of spectrum S over some interval [a,a+1) is negative
if opt==1:
  1, if sum of spectrum S over all intervals [a,a+1) and (a,a+1) is positive
  0, if sum of spectrum S over some interval [a,a+1) or (a,a+1) is negative
default: opt=0
@end format
EXAMPLE: example spissemicont; shows examples
"
{
  int opt=0;
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      opt=1;
    }
  }
  int i,j,k=1,1,0;
  while(j<=size(S[2]))
  {
    while(j+1<=size(S[2])&&S[1][j]<S[1][i]+1)
    {
      k=k+S[2][j];
      j++;
    }
    if(j==size(S[2])&&S[1][j]<S[1][i]+1)
    {
      k=k+S[2][j];
      j++;
    }
    if(k<0)
    {
      return(0);
    }
    k=k-S[2][i];
    if(k<0&&opt==1)
    {
      return(0);
    }
    i++;
  }
  return(1);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S1);
  list S2=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(S2);
  spissemicont(spsub(S1,spmul(S2,5)));
  spissemicont(spsub(S1,spmul(S2,5)),1);
  spissemicont(spsub(S1,spmul(S2,6)));
}
///////////////////////////////////////////////////////////////////////////////

proc spsemicont(list S0,list S,list #)
"USAGE:   spsemicont(S,k[,opt]); list S0, list S, int opt
RETURN:  list of intvecs l:
         spissemicont(sub(S0,spmul(S,k)),opt)==1 iff k<=l[i] for some i
NOTE:    if the spectra S occur with multiplicities k in a deformation
         of the [quasihomogeneous] spectrum S0 then 
         spissemicont(sub(S0,spmul(S,k))[,1])==1
EXAMPLE: example spsemicont; shows examples
"
{
  list l,l0;
  int i,j,k;
  while(spissemicont(S0,#))
  {
    if(size(S)>1)
    {
      l0=spsemicont(S0,list(S[1..size(S)-1]));
      for(i=1;i<=size(l0);i++)
      {
        if(size(l)>0)
        {
          j=1;
          while(j<size(l)&&l[j]!=l0[i])
          {
            j++;
          }
          if(l[j]==l0[i])
          {
            l[j][size(S)]=k;
          }
          else
          {
            l0[i][size(S)]=k;
            l=l+list(l0[i]);
          }
        }
        else
        {
          l=l0;
        }
      }
    }
    S0=spsub(S0,S[size(S)]);
    k++;
  }
  if(size(S)>1)
  {
    return(l);
  }
  else
  {
    return(list(intvec(k-1)));
  }
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S0=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S0);
  list S1=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(S1);
  list S2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
  spprint(S2);
  list S=S1,S2;
  list l=spsemicont(S0,S);
  l;
  spissemicont(spsub(S0,spmul(S,l[1])));
  spissemicont(spsub(S0,spmul(S,l[1]-1)));
  spissemicont(spsub(S0,spmul(S,l[1]+1)));
}
///////////////////////////////////////////////////////////////////////////////

proc spmilnor(list S)
"USAGE:   spmilnor(S); list S
RETURN:  int: Milnor number of spectrum S
EXAMPLE: example spmilnor; shows examples
"
{
  return(sum(S[2]));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S);
  spmilnor(S);
}
///////////////////////////////////////////////////////////////////////////////

proc spgeomgenus(list S)
"USAGE:   spgeomgenus(S); list S
RETURN:  int: geometrical genus of spectrum S
EXAMPLE: example spgeomgenus; shows examples
"
{
  int g=0;
  int i=1;
  while(i+1<=size(S[2])&&number(S[1][i])<=number(0))
  {
    g=g+S[2][i];
    i++;
  }
  if(i==size(S[2])&&number(S[1][i])<=number(0))
  {
    g=g+S[2][i];
  }
  return(g);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S);
  spgeomgenus(S);
}
///////////////////////////////////////////////////////////////////////////////

proc spgamma(list S)
"USAGE:   spgamma(S); list S
RETURN:  number: gamma invariant of spectrum S
EXAMPLE: example spgamma; shows examples
"
{
  int i,j;
  number g=0;
  for(i=1;i<=ncols(S[1]);i++)
  {
    for(j=1;j<=S[2][i];j++)
    {
      g=g+(number(S[1][i])-number(nvars(basering)-2)/2)^2;
    }
  }
  g=-g/4+sum(S[2])*number(S[1][ncols(S[1])]-S[1][1])/48;
  return(g);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(S);
  spgamma(S);
}
///////////////////////////////////////////////////////////////////////////////