source: git/Singular/LIB/gaussman.lib @ 962b1f

///////////////////////////////////////////////////////////////////////////////
version="$Id: gaussman.lib,v 1.70 2002-02-25 18:14:39 mschulze Exp $";
category="Singularities";

info="
LIBRARY:  gaussman.lib  Algorithmic Gauss-Manin Connection

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

OVERVIEW: A library to compute Hodge-theoretic invariants 
          of isolated hypersurface singularities

PROCEDURES:
 gmsring(t,s);              Gauss-Manin connection of t with variable s
 gmsnf(p,K);                Gauss-Manin connection normal form of p
 gmscoeffs(p,K);            Gauss-Manin connection basis representation of p
 monodromy(t);              Jordan data of monodromy of t
 spectrum(t);               singularity spectrum of t
 sppairs(t);                spectral pairs of t
 spnf(a[,m][,V]);           spectrum normal form of (a,m,V)
 sppnf(a,w[,m][,V]);        spectral pairs normal form of (a,w,m,V)
 vfilt(t);                  V-filtration of t on Brieskorn lattice
 vwfilt(t);                 weighted V-filtration of t on Brieskorn lattice
 tmatrix(t);                t-matrix on Brieskorn lattice
 endvfilt(V);               endomorphism V-filtration on Jacobian algebra
 spprint(sp);               print spectrum sp
 sppprint(spp);             print spectral pairs spp
 spadd(sp1,sp2);            sum of spectra sp1 and sp2
 spsub(sp1,sp2);            difference of spectra sp1 and sp2
 spmul(sp0,k);              linear combination of spectra sp
 spissemicont(sp[,opt]);    semicontinuity test of spectrum sp
 spsemicont(sp0,sp[,opt]);  semicontinuous combinations of spectra sp0 in sp
 spmilnor(sp);              milnor number of spectrum sp
 spgeomgenus(sp);           geometrical genus of spectrum sp
 spgamma(sp);               gamma invariant of spectrum sp

SEE ALSO: mondromy_lib, spectrum_lib

KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
          monodromy; spectrum; spectral pairs;
          mixed Hodge structure; V-filtration; weight filtration
";

LIB "linalg2.lib";

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

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

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

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

  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(t,s); poly t, string s
ASSUME:   characteristic 0; local degree ordering;
          isolated critical point 0 of t
RETURN:
@format
ring G;  Gauss-Manin connection of t with variable s
  poly gmspoly=t;
  ideal gmsjacob;  Jacobian ideal of t
  ideal gmsstd;  standard basis of Jacobian ideal
  matrix gmsmatrix;  matrix(gmsjacob)*gmsmatrix==matrix(gmsstd)
  ideal gmsbasis;  monomial vector space basis of Jacobian algebra
  int gmsmaxdeg;  maximal weight of variables
@end format
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice
EXAMPLE:  example gmsring; 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 critical point 0 expected");
    }
  }  
  matrix B=l[2];
  ideal m=kbase(g);

  int gmsmaxdeg;
  for(i=nvars(R);i>=1;i--)
  {
    if(deg(var(i))>gmsmaxdeg)
    {
      gmsmaxdeg=deg(var(i));
    }
  }

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

  execute("ring G="+string(charstr(R))+",("+s+","+varstr(R)+"),(ws("+
    string(deg(highcorner(g))+2*gmsmaxdeg)+"),"+os+",c);");

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

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

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

proc gmsnf(ideal p,int K)
"USAGE:    gmsnf(p,K); poly p, int K
ASSUME:   basering returned by gmsring
RETURN:
@format
list nf; 
  ideal nf[1];  projection of p to gmsbasis mod s^(K+1)
  ideal nf[2];  p=nf[1]+nf[2]
@end format
NOTE:     by setting p=nf[2] the computation can be continued
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice
EXAMPLE:  example gmsnf; shows examples
"
{
  return(system("gmsnf",p,gmsstd,gmsmatrix,(K+1)*deg(var(1))-2*gmsmaxdeg,K));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  def G=gmsring(t,"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)
"USAGE:    gmscoeffs(p,K); poly p, int K
ASSUME:   basering constructed by gmsring
RETURN:
@format
list l; 
  matrix l[1];  gmsbasis representation of p mod s^(K+1)
  ideal l[2];  p=matrix(gmsbasis)*l[1]+l[2]
@end format
NOTE:     by setting p=l[2] the computation can be continued
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice
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 t=x5+x2y2+y5;
  def G=gmsring(t,"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 min(ideal e)
{
  int i;
  number m=number(e[1]);
  for(i=2;i<=ncols(e);i++)
  {
    if(number(e[i])<m)
    {
      m=number(e[i]);
    }
  }
  return(m);
}
///////////////////////////////////////////////////////////////////////////////

static proc max(ideal e)
{
  int i;
  number m=number(e[1]);
  for(i=2;i<=ncols(e);i++)
  {
    if(number(e[i])>m)
    {
      m=number(e[i]);
    }
  }
  return(m);
}
///////////////////////////////////////////////////////////////////////////////

static proc saturate()
{
  int mu=ncols(gmsbasis);
  ideal r=gmspoly*gmsbasis;
  matrix A0[mu][mu],C;
  module H0;
  module H,H1=freemodule(mu),freemodule(mu);
  int k=-1;
  list l;

  while(size(reduce(H,std(H0*s)))>0)
  {
    dbprint(printlevel-voice+2,"// compute matrix A of t");
    k++;
    dbprint(printlevel-voice+2,"// k="+string(k));
    l=gmscoeffs(r,k);
    C,r=l[1..2];
    A0=A0+C;

    dbprint(printlevel-voice+2,"// compute saturation of H''");
    H0=H;
    H1=jet(module(A0*H1+s^2*diff(matrix(H1),s)),k);
    H=H*s+H1;
  }

  A0=A0-k*s;
  dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
  H=std(H0);

  dbprint(printlevel-voice+2,"// transform H'' to saturation of H''");
  H0=division(freemodule(mu)*s^k,H,k,intvec(1))[1];

  return(A0,r,H,H0,k);
}
///////////////////////////////////////////////////////////////////////////////

static proc basisrep(matrix A0,ideal r,module H,int k0,int K)
{
  dbprint(printlevel-voice+2,"// compute matrix A of t");
  dbprint(printlevel-voice+2,"// k="+string(K+k0+1));
  list l=gmscoeffs(r,K+k0+1);
  matrix C;
  C,r=l[1..2];
  A0=A0+C;

  dbprint(printlevel-voice+2,"// transform A to saturation of H''");
  matrix A=division(module(A0*H+s^2*diff(matrix(H),s)),H,K+1,intvec(1))[1]/s;

  return(A,A0,r);
}
///////////////////////////////////////////////////////////////////////////////

static proc eigvals(matrix A0,ideal r,module H,int k0)
{
  dbprint(printlevel-voice+2,
    "// compute eigenvalues e with multiplicities m of A0");
  matrix A;
  A,A0,r=basisrep(A0,r,H,k0,0);
  list l=eigenvals(A);
  def e,m=l[1..2];
  dbprint(printlevel-voice+2,"// e="+string(e));
  dbprint(printlevel-voice+2,"// m="+string(m));

  return(e,m,A0,r);
}
///////////////////////////////////////////////////////////////////////////////

static proc transf(matrix A,matrix A0,ideal r,module H,module H0,ideal e,intvec m,int k0,int K,int opt)
{
  int mu=ncols(gmsbasis);

  number e0,e1=min(e),max(e);

  int i,j,k;
  int k1;
  intvec d;
  d[ncols(e)]=0;
  if(opt)
  {
    dbprint(printlevel-voice+2,
      "// compute maximal differences d of e");
    for(i=1;i<=ncols(e);i++)
    {
      d[i]=int(e[i]-e0);
    }
  }
  else
  {
    dbprint(printlevel-voice+2,
      "// compute maximal integer differences d of e");
    for(i=1;i<ncols(e);i++)
    {
      for(j=i+1;i<=ncols(e);i++)
      {
        k=int(e[i]-e[j]);
        if(number(e[i]-e[j])==k)
        {
          if(k>d[i])
          {
            d[i]=k;
          }
          if(-k>d[j])
          {
            d[j]=-k;
          }
        }
      }
    }
  }
  dbprint(printlevel-voice+2,"// d="+string(d));

  for(i,k=1,0;i<=size(d);i++)
  {
    if(k<d[i])
    {
      k=d[i];
    }
  }

  A,A0,r=basisrep(A0,r,H,k0,K+k1);
  module U0=s^k0*freemodule(mu);

  if(k>0)
  {
    int i0,j0,i1,j1;
    module U,V;
    list l;

    while(k>0)
    {
      dbprint(printlevel-voice+2,"// transform to separate eigenvalues");
      U=0;
      for(i=1;i<=ncols(e);i++)
      {
        U=U+syz(power(jet(A,0)-e[i],m[i]));
      }
      V=inverse(U);
      A=V*A*U;
      H0=V*H0;
      U0=U0*U;

      dbprint(printlevel-voice+2,
        "// transform to reduce d by 1");
      for(i0,i=1,1;i0<=ncols(e);i0++)
      {
        for(i1=1;i1<=m[i0];i1,i=i1+1,i+1)
        {
          for(j0,j=1,1;j0<=ncols(e);j0++)
          {
            for(j1=1;j1<=m[j0];j1,j=j1+1,j+1)
            {
              if(d[i0]==0&&d[j0]>0)
              {
                A[i,j]=A[i,j]/s;
              }
              if(d[i0]>0&&d[j0]==0)
              {
                A[i,j]=A[i,j]*s;
              }
            }
          }
        }
      }

      H0=transpose(H0);
      for(i0,i=1,1;i0<=ncols(e);i0++)
      {
        if(d[i0]>0)
        {
          for(i1=1;i1<=m[i0];i1,i=i1+1,i+1)
          {
            A[i,i]=A[i,i]-1;
            H0[i]=H0[i]*s;
            U0[i]=U0[i]/s;
          }
          e[i0]=e[i0]-1;
          d[i0]=d[i0]-1;
        }
        else
        {
          i=i+m[i0];
        }
      }
      H0=transpose(H0);

      l=sppnf(e,d,m);
      e,d,m=l[1..3];

      k--;
    }

    A=jet(A,K);
  }

  return(A,A0,r,H0,U0,e,m);
}
///////////////////////////////////////////////////////////////////////////////

proc monodromy(poly t)
"USAGE:    monodromy(t); poly t
ASSUME:   characteristic 0; local degree ordering;
          isolated critical point 0 of t
RETURN:   list l;  Jordan data jordan(M) of monodromy matrix exp(-2*pi*i*M)
SEE ALSO: mondromy_lib, linalg.lib
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; monodromy
EXAMPLE:  example monodromy; shows examples
"
{
  def R=basering;
  int n=nvars(R)-1;
  def G=gmsring(t,"s");
  setring(G);

  matrix A;
  module U0;
  ideal e;
  intvec m;

  def A0,r,H,H0,k0=saturate();
  e,m,A0,r=eigvals(A0,r,H,k0);
  A,A0,r,H0,U0,e,m=transf(A,A0,r,H,H0,e,m,k0,0,0);

  list l=jordan(A,e,m);
  setring(R);
  list l=imap(G,l);
  kill G,gmsmaxdeg;

  return(l);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  monodromy(t);
}
///////////////////////////////////////////////////////////////////////////////

proc spectrum(poly t)
"USAGE:    spectrum(t); poly t
ASSUME:   characteristic 0; local degree ordering;
          isolated critical point 0 of t
RETURN:
@format
list sp;  singularity spectrum of t
  ideal sp[1];
    number sp[1][i];  i-th spectral number
  intvec sp[2];
    int sp[2][i];  multiplicity of i-th spectral number
@end format
SEE ALSO: spectrum_lib
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; 
          mixed Hodge structure; V-filtration; spectrum
EXAMPLE:  example spectrum; shows examples
"
{
  list l=vwfilt(t);
  return(spnf(l[1],l[3]));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  spprint(spectrum(t));
}
///////////////////////////////////////////////////////////////////////////////

proc sppairs(poly t)
"USAGE:    sppairs(t); poly t
ASSUME:   characteristic 0; local degree ordering;
          isolated critical point 0 of t
RETURN:
@format
list spp;  spectral pairs of t
  ideal spp[1];
    number spp[1][i];  V-filtration index of i-th spectral pair
  intvec spp[2];
    int spp[2][i];  weight filtration index of i-th spectral pair
  intvec spp[3]; 
    int spp[3][i];  multiplicity of i-th spectral pair
@end format
SEE ALSO: spectrum_lib
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
          mixed Hodge structure; V-filtration; weight filtration;
          spectrum; spectral pairs
EXAMPLE:  example sppairs; shows examples
"
{
  list l=vwfilt(t);
  return(list(l[1],l[2],l[3]));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  sppprint(sppairs(t));
}
///////////////////////////////////////////////////////////////////////////////

proc spnf(ideal e,list #)
"USAGE:   spnf(e[,m][,V]); ideal e, intvec m, list V
ASSUME:  ncols(e)==size(m)==size(V); typeof(V[i])=="int"
RETURN:
@format
list sp;  spectrum normal form of (e,m,V)
  ideal sp[1];  numbers in e in increasing order
  intvec sp[2]; 
    int sp[2][i];  multiplicity of number sp[1][i] in e
  list sp[3]; 
    module sp[3][i];  module associated to number sp[1][i]
@end format
EXAMPLE: example spnf; shows examples
"
{
  int n=ncols(e);
  intvec m;
  module v;
  list V;
  int i,j;
  while(i<size(#))
  {
    i++;
    if(typeof(#[i])=="intvec")
    {
      m=#[i];
    }
    if(typeof(#[i])=="module")
    {
      v=#[i];
      for(j=n;j>=1;j--)
      {
        V[j]=module(v[j]);
      }
    }
    if(typeof(#[i])=="list")
    {
      V=#[i];
    }
  }
  if(m==0)
  {
    for(i=n;i>=1;i--)
    {
      m[i]=1;
    }
  }

  int k;
  ideal e0;
  intvec m0;
  list V0;
  number e1;
  int m1;
  for(i=n;i>=1;i--)
  {
    if(m[i]!=0)
    {
      for(j=i-1;j>=1;j--)
      {
        if(m[j]!=0)
        {
          if(number(e[i])>number(e[j]))
          {
            e1=number(e[i]);
            e[i]=e[j];
            e[j]=e1;
            m1=m[i];
            m[i]=m[j];
            m[j]=m1;
            if(size(V)>0)
            {
              v=V[i];
              V[i]=V[j];
              V[j]=v;
            }
          }
          if(number(e[i])==number(e[j]))
          {
            m[i]=m[i]+m[j];
            m[j]=0;
            if(size(V)>0)
            {
              V[i]=V[i]+V[j];
            }
          }
        }
      }
      k++;
      e0[k]=e[i];
      m0[k]=m[i];
      if(size(V)>0)
      {
        V0[k]=V[i];
      }
    }
  }

  if(size(V0)>0)
  {
    n=size(V0);
    module U=std(V0[n]);
    for(i=n-1;i>=1;i--)
    {
      V0[i]=simplify(reduce(V0[i],U),1);
      if(i>=2)
      {
        U=std(U+V0[i]);
      }
    }
  }

  list l;
  if(k>0)
  {
    l=e0,m0;
    if(size(V0)>0)
    {
      l[3]=V0;
    }
  }
  return(l);
}
example
{ "EXAMPLE:"; echo=2;
}
///////////////////////////////////////////////////////////////////////////////

proc sppnf(ideal a,intvec w,list #)
"USAGE:    sppnf(a,w[,m][,V]); ideal a, intvec w, intvec m, list V
ASSUME:   ncols(e)=size(w)=size(m)=size(V); typeof(V[i])=="module"
RETURN:
@format
list spp;  spectral pairs normal form of (a,w,m,V)
  ideal spp[1];
    number spp[1][i];  V-filtration index of i-th spectral pair
  intvec spp[2];
    int spp[2][i];  weight filtration index of i-th spectral pair
  intvec spp[3]; 
    int spp[3][i];  multiplicity of i-th spectral pair
  list spp[4]; 
    module spp[4][i];  vector space of i-th spectral pair
@end format
SEE ALSO: spectrum_lib
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
          mixed Hodge structure; V-filtration; weight filtration;
          spectrum; spectral pairs
EXAMPLE:  example sppnorm; shows examples
"
{
  int n=ncols(a);
  intvec m;
  module v;
  list V;
  int i,j;
  while(i<size(#))
  {
    i++;
    if(typeof(#[i])=="intvec")
    {
      m=#[i];
    }
    if(typeof(#[i])=="module")
    {
      v=#[i];
      for(j=n;j>=1;j--)
      {
        V[j]=module(v[j]);
      }
    }
    if(typeof(#[i])=="list")
    {
      V=#[i];
    }
  }
  if(m==0)
  {
    for(i=n;i>=1;i--)
    {
      m[i]=1;
    }
  }

  int k;
  ideal a0;
  intvec w0,m0;
  list V0;
  number a1;
  int w1,m1;
  for(i=n;i>=1;i--)
  {
    if(m[i]!=0)
    {
      for(j=i-1;j>=1;j--)
      {
        if(m[j]!=0)
        {
          if(number(a[i])>number(a[j])||
            (number(a[i])==number(a[j])&&w[i]<w[j]))
          {
            a1=number(a[i]);
            a[i]=a[j];
            a[j]=a1;
            w1=w[i];
            w[i]=w[j];
            w[j]=w1;
            m1=m[i];
            m[i]=m[j];
            m[j]=m1;
            if(size(V)>0)
            {
              v=V[i];
              V[i]=V[j];
              V[j]=v;
            }
          }
          if(number(a[i])==number(a[j])&&w[i]==w[j])
          {
            m[i]=m[i]+m[j];
            m[j]=0;
            if(size(V)>0)
            {
              V[i]=V[i]+V[j];
            }
          }
        }
      }
      k++;
      a0[k]=a[i];
      w0[k]=w[i];
      m0[k]=m[i];
      if(size(V)>0)
      {
        V0[k]=V[i];
      }
    }
  }

  if(size(V0)>0)
  {
    n=size(V0);
    module U=std(V0[n]);
    for(i=n-1;i>=1;i--)
    {
      V0[i]=simplify(reduce(V0[i],U),1);
      if(i>=2)
      {
        U=std(U+V0[i]);
      }
    }
  }

  list l;
  if(k>0)
  {
    l=a0,w0,m0;
    if(size(V0)>0)
    {
      l[4]=V0;
    }
  }
  return(l);
}
example
{ "EXAMPLE:"; echo=2;
}
///////////////////////////////////////////////////////////////////////////////

proc vfilt(poly t)
"USAGE:    vfilt(t); poly t
ASSUME:   characteristic 0; local degree ordering;
          isolated critical point 0 of t
RETURN:
@format
list v;  V-filtration on H''/s*H''
  ideal v[1];
    number v[1][i];  V-filtration index of i-th spectral pair
  intvec v[2]; 
    int v[2][i];  multiplicity of i-th spectral pair
  list v[3]; 
    module v[3][i];  vector space of i-th graded part in terms of v[4]
  ideal v[4];  monomial vector space basis of H''/s*H''
  ideal v[5];  standard basis of Jacobian ideal
@end format
SEE ALSO: spectrum_lib
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
          mixed Hodge structure; V-filtration; spectrum
EXAMPLE:  example vfilt; shows examples
"
{
  list l=vwfilt(t);
  return(spnf(l[1],l[3],l[4])+list(l[5],l[6]));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  vfilt(t);
}
///////////////////////////////////////////////////////////////////////////////

proc vwfilt(poly t)
"USAGE:    vwfilt(t); poly t
ASSUME:   characteristic 0; local degree ordering;
          isolated critical point 0 of t
RETURN:
@format
list vw;  weighted V-filtration on H''/s*H''
  ideal vw[1];
    number vw[1][i];  V-filtration index of i-th spectral pair
  intvec vw[2];
    int vw[2][i];  weight filtration index of i-th spectral pair
  intvec vw[3]; 
    int vw[3][i];  multiplicity of i-th spectral pair
  list vw[4]; 
    module vw[4][i];  vector space of i-th graded part in terms of vw[5]
  ideal vw[5];  monomial vector space basis of H''/s*H''
  ideal vw[6];  standard basis of Jacobian ideal
@end format
SEE ALSO: spectrum_lib
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
          mixed Hodge structure; V-filtration; weight filtration;
          spectrum; spectral pairs
EXAMPLE:  example vwfilt; shows examples
"
{
  def R=basering;
  int n=nvars(R)-1;
  def G=gmsring(t,"s");
  setring(G);

  int mu=ncols(gmsbasis);
  matrix A;
  module U0;
  ideal e;
  intvec m;

  def A0,r,H,H0,k0=saturate();
  e,m,A0,r=eigvals(A0,r,H,k0);
  A,A0,r,H0,U0,e,m=transf(A,A0,r,H,H0,e,m,k0,0,1);

  dbprint(printlevel-voice+2,"// compute weight filtration basis");
  list l=jordanbasis(A,e,m);
  def U,v=l[1..2];
  kill l;
  vector u0;
  int v0;
  int i,j,k,l;
  for(k,l=1,1;l<=ncols(e);k,l=k+m[l],l+1)
  {
    for(i=k+m[l]-1;i>=k+1;i--)
    {
      for(j=i-1;j>=k;j--)
      {
        if(v[i]>v[j])
        {
          v0=v[i];v[i]=v[j];v[j]=v0;
          u0=U[i];U[i]=U[j];U[j]=u0;
        }
      }
    }
  }

  dbprint(printlevel-voice+2,"// transform to weight filtration basis");
  matrix V=inverse(U);
  A=V*A*U;
  dbprint(printlevel-voice+2,"// compute normal form of H''");
  H0=std(V*H0);
  U0=U0*U;

  dbprint(printlevel-voice+2,"// compute spectral pairs");
  ideal a;
  intvec w;
  for(i=1;i<=mu;i++)
  {
    j=leadexp(H0[i])[nvars(basering)+1];
    a[i]=A[j,j]+ord(H0[i])/deg(s)-1;
    w[i]=v[j]+n;
  }
  kill v;
  module v=simplify(jet(H*U0*H0,2*k0)/s^(2*k0),1);

  kill l;
  list l=sppnf(a,w,v)+list(gmsbasis,gmsstd);
  setring(R);
  list l=imap(G,l);
  kill G,gmsmaxdeg;
  attrib(l[5],"isSB",1);

  return(l);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  vwfilt(t);
}
///////////////////////////////////////////////////////////////////////////////

static proc commutator(matrix A)
{
  int n=ncols(A);
  int i,j,k;
  matrix C[n^2][n^2];
  for(i=0;i<n;i++)
  {
    for(j=0;j<n;j++)
    {
      for(k=0;k<n;k++)
      {
        C[i*n+j+1,k*n+j+1]=C[i*n+j+1,k*n+j+1]+A[i+1,k+1];
        C[i*n+j+1,i*n+k+1]=C[i*n+j+1,i*n+k+1]-A[k+1,j+1];
      }
    }
  }
  return(C);
}

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

proc tmatrix(poly t,list #)
"USAGE:    tmatrix(t); poly t
ASSUME:   characteristic 0; local degree ordering;
          isolated critical point 0 of t
RETURN:   list A;  t-matrix A[1]+s*A[2] on H''
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
          mixed Hodge structure; opposite Hodge filtration; V-filtration;
EXAMPLE:  example tmatrix; shows examples
"
{
  def R=basering;
  int n=nvars(R)-1;
  def G=gmsring(t,"s");
  setring(G);

  int mu=ncols(gmsbasis);
  matrix A;
  module U0;
  ideal e;
  intvec m;

  def A0,r,H,H0,k0=saturate();
  e,m,A0,r=eigvals(A0,r,H,k0);
  int k1=int(max(e)-min(e));
  A,A0,r,H0,U0,e,m=transf(A,A0,r,H,H0,e,m,k0,k0+k1,1);

  ring S=0,s,(ds,c);
  matrix A=imap(G,A);
  module H0=imap(G,H0);
  ideal e=imap(G,e);
  kill G,gmsmaxdeg;

  dbprint(printlevel-voice+2,"// transform to Jordan basis");
  module U=jordanbasis(A,e,m)[1];
  matrix V=inverse(U);
  A=V*A*U;
  module H=V*H0;

  dbprint(printlevel-voice+2,"// compute splitting of V-filtration");
  int i,j,k;
  U=freemodule(mu);
  V=matrix(0,mu,mu);
  matrix v[mu^2][1];
  matrix A0=commutator(jet(A,0));
  for(k=1;k<=k0+k1;k++)
  {
    for(j=0;j<k;j++)
    {
      V=matrix(V)-(jet(A,k-j)/s^(k-j))*(jet(U,j)/s^j);
    }
    v=V[1..mu,1..mu];
    v=inverse(A0+k)*v;
    V=v[1..mu^2,1];
    U=matrix(U)+s^k*V;
  }
  attrib(U,"isSB",1);

  dbprint(printlevel-voice+2,"// transform to V-splitting basis");
  A=jet(A,0);
  H=std(division(H,U,k0+k1,intvec(1))[1]);

  dbprint(printlevel-voice+2,"// compute V-leading terms of H''");
  int i0,j0;
  module H1=H;
  for(k=ncols(H1);k>=1;k--)
  {
    i0=leadexp(H1[k])[nvars(basering)+1];
    j0=ord(H1[k]);//deg(s);
    H0[k]=lead(H1[k]);
    H1[k]=H1[k]-lead(H1[k]);
    if(H1[k]!=0)
    {
      i=leadexp(H1[k])[nvars(basering)+1];
      j=ord(H1[k]);//deg(s);
      while(A[i,i]+j==A[i0,i0]+j0)
      {
        H0[k]=H0[k]+lead(H1[k]);
        H1[k]=H1[k]-lead(H1[k]);
        i=leadexp(H1[k])[nvars(basering)+1];
        j=ord(H1[k]);//deg(s);
      }
    }
  }
  H0=simplify(H0,1);

  dbprint(printlevel-voice+2,"// compute N");
  matrix N=A;
  for(i=1;i<=ncols(N);i++)
  {
    N[i,i]=0;
  }

  dbprint(printlevel-voice+2,"// compute splitting of Hodge filtration");
  U=0;
  module U1;
  module C;
  list F,I;
  module F0,I0,U0;
  for(i0,j0=1,1;i0<=ncols(e);i0++)
  {
    C=matrix(0,mu,1);
    for(j=m[i0];j>=1;j,j0=j-1,j0+1)
    {
      C=C+gen(j0);
    }
    F0=intersect(C,H0);

    F=list();
    j=0;
    while(size(F0)>0)
    {
      j++;
      F[j]=matrix(0,mu,1);
      if(size(jet(F0,0))>0)
      {
        for(i=ncols(F0);i>=1;i--)
        {
          if(ord(F0[i])==0)
          {
            F[j]=F[j]+F0[i];
          }
        }
      }
      for(i=ncols(F0);i>=1;i--)
      {
        F0[i]=F0[i]/s;
      }
    }

    I=list();
    I0=module();
    U0=std(module());
    for(i=size(F);i>=1;i--)
    {
      I[i]=module();
    }
    for(i=1;i<=size(F);i++)
    {
      I0=reduce(F[i],U0);
      j=i;
      while(size(I0)>0)
      {
        U0=std(U0+I0);
        I[j]=I[j]+I0;
        I0=reduce(N*I0,U0);
        j++;
      }
    }

    for(i=1;i<=size(I);i++)
    {
      U=U+I[i];
    }
  }

  dbprint(printlevel-voice+2,"// transform to Hodge splitting basis");
  V=inverse(U);
  A=V*A*U;
  H=V*H;

  dbprint(printlevel-voice+2,"// compute reduced standard basis of H''");
  degBound=k0+k1+2;
  option(redSB);
  H=std(H);
  option(noredSB);
  degBound=0;
  H=simplify(jet(H,k0+k1),1);
  attrib(H,"isSB",1);

  dbprint(printlevel-voice+2,"// compute matrix A0+sA1 of t");
  A=division(module(s*A*H+s^2*diff(matrix(H),s)),H,1,intvec(1))[1];
  A0=jet(A,0);
  A=jet(A,1)/s;

  setring(R);
  matrix A0=imap(S,A0);
  matrix A1=imap(S,A);
  kill S;
  return(list(A0,A1));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  list A=tmatrix(t);
  print(A[1]);
  print(A[2]);
}
///////////////////////////////////////////////////////////////////////////////

proc endvfilt(list v)
"USAGE:   endvfilt(v); list v
ASSUME:  v returned by vfilt
RETURN:
@format
list ev;  V-filtration on Jacobian algebra
  ideal ev[1];
    number ev[1][i];  V-filtration index of i-th spectral pair
  intvec ev[2]; 
    int ev[2][i];  multiplicity of i-th spectral pair
  list ev[3]; 
    module ev[3][i];  vector space of i-th graded part in terms of ev[4]
  ideal ev[4];  monomial vector space basis of Jacobian algebra
  ideal ev[5];  standard basis of Jacobian ideal
@end format
KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
          mixed Hodge structure; V-filtration; endomorphism filtration
EXAMPLE:  example endvfilt; shows examples
"
{
  def a,d,V,m,g=v[1..5];
  attrib(g,"isSB",1);
  int mu=ncols(m);

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

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

  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;
  list l;

  while(b0<number(a[ncols(a)]-a[1]))
  {
    dbprint(printlevel-voice+2,"// 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;

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

    dbprint(printlevel-voice+2,"// collect conditions for EV["+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,"// 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,"// compute kernel of condition matrix");
    v0=v0+list(syz(L));
    a0=a0,b0;
  }

  dbprint(printlevel-voice+2,"// 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]));
  }
  option(noredSB);

  dbprint(printlevel-voice+2,"// 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,g));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  poly t=x5+x2y2+y5;
  endvfilt(vfilt(t));
}
///////////////////////////////////////////////////////////////////////////////

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

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

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

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

proc spmul(list #)
"USAGE:   spmul(sp0,k); list sp0, int[vec] k
RETURN:  list sp;  linear combination of spectra sp0 with coefficients k
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 sp0=list(ideal(),intvec(0));
        for(int i=size(#[2]);i>=1;i--)
        {
          sp0=spadd(sp0,spmul(#[1][i],#[2][i]));
        }
        return(sp0);
      }
    }
  }
  return(list(ideal(),intvec(0)));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list sp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(sp);
  spprint(spmul(sp,2));
  list sp1=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(sp1);
  list sp2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
  spprint(sp2);
  spprint(spmul(list(sp1,sp2),intvec(1,2)));
}
///////////////////////////////////////////////////////////////////////////////

proc spissemicont(list sp,list #)
"USAGE:   spissemicont(sp[,1]); list sp, int opt
RETURN:
@format
int k=
  1;  if sum of sp is positive on all intervals [a,a+1) [and (a,a+1)]
  0;  if sum of sp is negative on some interval [a,a+1) [or (a,a+1)]
@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(sp[2]))
  {
    while(j+1<=size(sp[2])&&sp[1][j]<sp[1][i]+1)
    {
      k=k+sp[2][j];
      j++;
    }
    if(j==size(sp[2])&&sp[1][j]<sp[1][i]+1)
    {
      k=k+sp[2][j];
      j++;
    }
    if(k<0)
    {
      return(0);
    }
    k=k-sp[2][i];
    if(k<0&&opt==1)
    {
      return(0);
    }
    i++;
  }
  return(1);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list sp1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(sp1);
  list sp2=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(sp2);
  spissemicont(spsub(sp1,spmul(sp2,5)));
  spissemicont(spsub(sp1,spmul(sp2,5)),1);
  spissemicont(spsub(sp1,spmul(sp2,6)));
}
///////////////////////////////////////////////////////////////////////////////

proc spsemicont(list sp0,list sp,list #)
"USAGE:   spsemicont(sp0,sp,k[,1]); list sp0, list sp
RETURN:
@format
list l;
  intvec l[i];  if the spectra sp0 occur with multiplicities k 
                in a deformation of a [quasihomogeneous] singularity 
                with spectrum sp then k<=l[i]
@end format
EXAMPLE: example spsemicont; shows examples
"
{
  list l,l0;
  int i,j,k;
  while(spissemicont(sp0,#))
  {
    if(size(sp)>1)
    {
      l0=spsemicont(sp0,list(sp[1..size(sp)-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(sp)]=k;
          }
          else
          {
            l0[i][size(sp)]=k;
            l=l+list(l0[i]);
          }
        }
        else
        {
          l=l0;
        }
      }
    }
    sp0=spsub(sp0,sp[size(sp)]);
    k++;
  }
  if(size(sp)>1)
  {
    return(l);
  }
  else
  {
    return(list(intvec(k-1)));
  }
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list sp0=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(sp0);
  list sp1=list(ideal(-1/6,1/6),intvec(1,1));
  spprint(sp1);
  list sp2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
  spprint(sp2);
  list sp=sp1,sp2;
  list l=spsemicont(sp0,sp);
  l;
  spissemicont(spsub(sp0,spmul(sp,l[1])));
  spissemicont(spsub(sp0,spmul(sp,l[1]-1)));
  spissemicont(spsub(sp0,spmul(sp,l[1]+1)));
}
///////////////////////////////////////////////////////////////////////////////

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

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

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