///////////////////////////////////////////////////////////////////////////////
version="$Id: gmspoly.lib,v 1.16 2008-03-18 15:50:14 Singular Exp $";
category="Singularities";

info="
LIBRARY:  gmspoly.lib  Gauss-Manin System of Tame Polynomials

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

OVERVIEW: A library to compute invariants related to the Gauss-Manin system
          of a cohomologically tame polynomial

PROCEDURES:
  isTame(f);  test if the polynomial f is tame
  goodBasis(f);  a good basis of the Brieskorn lattice of a cohomologically tame f

SEE ALSO: gmssing_lib

KEYWORDS: tame polynomial; Gauss-Manin system; Brieskorn lattice;
          mixed Hodge structure; V-filtration; weight filtration;
          monodromy; spectrum; spectral pairs; good basis
";

LIB "linalg.lib";
LIB "ring.lib";

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

static proc mindegree(matrix A)
{
  int d=0;

  while(A/var(1)^(d+1)*var(1)^(d+1)==A)
  {
    d++;
  }

  return(d);
}
///////////////////////////////////////////////////////////////////////////////

static proc maxdegree(matrix A)
{
  int d=0;
  matrix N[nrows(A)][ncols(A)];

  while(A/var(1)^(d+1)!=N)
  {
    d++;
  }

  return(d);
}
///////////////////////////////////////////////////////////////////////////////

proc isTame(poly f)
"USAGE:    isTame(f); poly f
ASSUME:   basering has no variables named w(1),w(2),...
RETURN:
@format
int k=
  0;  if f is not tame
  1;  if f is tame
@end format
KEYWORDS: tame polynomial
EXAMPLE:  example isTame; shows examples
"
{
  int d=vdim(std(jacob(f)));
  def @X=basering;
  int n=nvars(@X);
  def @WX=changechar("(0,w(1.."+string(n)+"))");
  setring @WX;
  ideal J=jacob(imap(@X,f));
  int i;
  for(i=1;i<=n;i++)
  {
    J[i]=J[i]+w(i);
  }
  int D=vdim(std(J));

  setring(@X);
  kill @WX;

  return(d>0&&d==D);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),dp;
  isTame(x2y+x);
  isTame(x3+y3+xy);
}
///////////////////////////////////////////////////////////////////////////////

static proc chart(matrix A)
{
  A=ideal(homog(transpose(ideal(A)),var(2)));
  def r=basering;
  map h=r,1,var(1);
  return(h(A));
}
///////////////////////////////////////////////////////////////////////////////

static proc pidbasis(module M0,module M)
{
  int m=nrows(M);
  int n=ncols(M);

  module L,N;
  module T=freemodule(m);
  while(matrix(L)!=matrix(M))
  {
    L=M;

    M=T,M;
    N=transpose(std(transpose(M)));
    T=N[1..m];
    M=N[m+1..ncols(N)];

    M=freemodule(n),transpose(M);
    N=std(transpose(M));
    N=transpose(simplify(N,1));
    M=N[n+1..ncols(N)];
    M=transpose(M);
  }

  if(maxdegree(M)>0)
  {
    print("   ? module not free");
    return(module());
  }

  attrib(M,"isSB",1);
  N=lift(T,simplify(reduce(M0,M),2));

  return(N);
}
///////////////////////////////////////////////////////////////////////////////

static proc vfilt(matrix B,int d)
{
  int mu=ncols(B);

  module V=freemodule(mu);
  module V0=var(1)^(d-1)*freemodule(mu);
  attrib(V0,"isSB",1);
  module V1=B;
  option(redSB);
  while(size(reduce(V1,V0))>0)
  {
    V=std(V0+V1);
    V0=var(1)^(d-1)*V;
    attrib(V0,"isSB",1);
    V1=B*matrix(V1)-var(1)^d*diff(matrix(V1),var(1));
  }
  option("noredSB");

  B=lift(V0,B*matrix(V)-var(1)^d*diff(matrix(V),var(1)));
  list l=eigenvals(B);
  def e0,s0=l[1..2];

  module U;
  int i,j,i0,j0,i1,j1,k;
  for(k=int(e0[ncols(e0)]-e0[1]);k>=1;k--)
  {
    U=0;
    for(i=1;i<=ncols(e0);i++)
    {
      U=U+syz(power(jet(B,0)-e0[i],s0[i]));
    }
    B=lift(U,B*U);
    V=V*U;

    for(i0,i=1,1;i0<=ncols(e0);i0++)
    {
      for(i1=1;i1<=s0[i0];i1,i=i1+1,i+1)
      {
        for(j0,j=1,1;j0<=ncols(e0);j0++)
        {
          for(j1=1;j1<=s0[j0];j1,j=j1+1,j+1)
          {
            if(leadcoef(e0[i0]-e0[1])>=1&&leadcoef(e0[j0]-e0[1])<1)
            {
              B[i,j]=B[i,j]/var(1);
            }
            if(leadcoef(e0[i0]-e0[1])<1&&leadcoef(e0[j0]-e0[1])>=1)
            {
              B[i,j]=B[i,j]*var(1);
            }
          }
        }
      }
    }

    for(i0,i=1,1;i0<=ncols(e0);i0++)
    {
      if(leadcoef(e0[i0]-e0[1])>=1)
      {
        for(i1=1;i1<=s0[i0];i1,i=i1+1,i+1)
        {
          B[i,i]=B[i,i]-1;
          V[i]=V[i]*var(1);
        }
        e0[i0]=e0[i0]-1;
      }
      else
      {
        i=i+s0[i0];
      }
    }

    l=spnf(list(e0,s0));
    e0,s0=l[1..2];
  }

  U=0;
  for(i=1;i<=ncols(e0);i++)
  {
    U=U+syz(power(jet(B,0)-e0[i],s0[i]));
  }
  B=lift(U,B*U);
  V=V*U;

  d=mindegree(V);
  V=V/var(1)^d;
  B=B+d*matrix(freemodule(mu));
  for(i=ncols(e0);i>=1;i--)
  {
    e0[i]=e0[i]+d;
  }

  return(e0,s0,V,B);
}
///////////////////////////////////////////////////////////////////////////////

static proc spec(ideal e0,intvec s0,module V,matrix B)
{
  int mu=ncols(B);

  int i,j,k;

  int d=maxdegree(V);
  int d0=d;
  V=chart(V);
  module U=std(V);
  while(size(reduce(var(1)^d*freemodule(mu),U))>0)
  {
    d++;
  }
  if(d>d0)
  {
    k=d-d0;
    B=B-k*freemodule(mu);
    for(i=1;i<=ncols(e0);i++)
    {
      e0[i]=e0[i]-k;
    }
  }
  module G=lift(V,var(1)^d*freemodule(mu));
  G=std(G);
  G=simplify(G,1);

  ideal e;
  intvec s;
  e[mu]=0;
  for(j,k=1,1;j<=ncols(e0);j++)
  {
    for(i=s0[j];i>=1;i,k=i-1,k+1)
    {
      e[k]=e0[j];
      s[k]=j;
    }
  }

  ideal a;
  a[mu]=0;
  for(i=1;i<=mu;i++)
  {
    a[i]=leadcoef(e[leadexp(G[i])[nvars(basering)+1]])+leadexp(G[i])[1];
  }

  return(a,e0,e,s,V,B,G);
}
///////////////////////////////////////////////////////////////////////////////

static proc fsplit(ideal e0,ideal e,intvec s,module V,matrix B,module G)
{
  int mu=ncols(e);

  int i,j,k;

  number n,n0;
  vector v,v0;
  list F;
  for(i=ncols(e0);i>=1;i--)
  {
    F[i]=module(matrix(0,mu,1));
  }
  for(i=mu;i>=1;i--)
  {
    v=G[i];
    v0=lead(v);
    n0=leadcoef(e[leadexp(v0)[nvars(basering)+1]])+leadexp(v0)[1];
    v=v-lead(v);
    while(v!=0)
    {
      n=leadcoef(e[leadexp(v)[nvars(basering)+1]])+leadexp(v)[1];
      if(n==n0)
      {
        v0=v0+lead(v);
        v=v-lead(v);
      }
      else
      {
        v=0;
      }
    }
    j=s[leadexp(v0)[nvars(basering)+1]];
    F[j]=F[j]+v0;
  }

  matrix B0=jet(B,0);
  module U,U0,U1,U2;
  matrix N;
  for(i=size(F);i>=1;i--)
  {
    N=B0-e0[i];
    U0=0;
    while(size(F[i])>0)
    {
      k=0;
      U1=jet(F[i],0);
      while(size(U1)>0)
      {
        for(j=ncols(U1);j>=1;j--)
        {
          if(size(reduce(U1[j],std(U0)))>0)
          {
            U0=U1[j]+U0;
          }
        }
        U1=N*U1;
        k++;
      }
      F[i]=module(F[i]/var(1));
    }
    U=U0+U;
  }

  V=V*U;
  G=lift(U,G);
  B=lift(U,B*U);

  return(e,V,B,G);
}
///////////////////////////////////////////////////////////////////////////////

static proc glift(ideal e,module V,matrix B,module G)
{
  int mu=ncols(e);

  int d=maxdegree(B);
  B=chart(B);
  G=std(G);
  G=simplify(G,1);

  int i,j,k;

  ideal v;
  for(i=mu;i>=1;i--)
  {
    v[i]=e[leadexp(G[i])[nvars(basering)+1]]+leadexp(G[i])[1];
  }

  number c;
  matrix g[mu][1];
  matrix m[mu][1];
  matrix a[mu][1];
  matrix A[mu][mu];
  module M=var(1)^d*G;
  module N=var(1)*B*matrix(G)+var(1)^(d+2)*diff(matrix(G),var(1));
  while(size(N)>0)
  {
    j=mu;
    for(k=mu-1;k>=1;k--)
    {
      if(N[k]>N[j])
      {
        j=k;
      }
    }

    i=mu;
    while(leadexp(M[i])[nvars(basering)+1]!=leadexp(N[j])[nvars(basering)+1])
    {
      i--;
    }

    k=leadexp(N[j])[1]-leadexp(M[i])[1];
    if(k==0||i==j)
    {
      c=leadcoef(N[j])/leadcoef(M[i]);
      A[i,j]=A[i,j]+c*var(1)^k;
      N[j]=N[j]-c*var(1)^k*M[i];
    }
    else
    {
      c=leadcoef(N[j])/leadcoef(M[i])/(1-k-leadcoef(v[i])+leadcoef(v[j]));
      G[j]=G[j]+c*var(1)^(k-1)*G[i];
      M[j]=M[j]+c*var(1)^(k-1)*M[i];
      g=c*var(1)^(k-1)*G[i];
      N[j]=N[j]+(var(1)*B*g+var(1)^(d+2)*diff(g,var(1)))[1];
      m=M[i];
      a=transpose(A)[j];
      N=N-c*var(1)^(k-1)*m*transpose(a);
    }
  }

  G=V*G;
  G=G/var(1)^mindegree(G);

  return(G,A);
}
///////////////////////////////////////////////////////////////////////////////

proc goodBasis(poly f)
"USAGE:    goodBasis(f); poly f
ASSUME:   f is cohomologically tame
RETURN:
@format
ring R;  basering with new variable s
  ideal b;  [matrix(b)] is a good basis of the Brieskorn lattice
  matrix A;  A(s)=A0+s*A1 and t[matrix(b)]=[matrix(b)](A(s)+s^2*(d/ds))
@end format
KEYWORDS: tame polynomial; Gauss-Manin system; Brieskorn lattice;
          mixed Hodge structure; V-filtration; weight filtration;
          monodromy; spectrum; spectral pairs; good basis
EXAMPLE:  example goodBasis; shows examples
"
{
  def @X=basering;
  int n=nvars(@X);
  ideal J=jacob(f);

  if(vdim(std(J))<=0)
  {
    ERROR("input is not cohomologically tame");
  }

  int i,j,k,l;

  ideal X=maxideal(1);
  string c="ring @XS=0,(s,"+varstr(@X)+"),(C,dp(1),dp(n));";
  execute(c);
  poly f=imap(@X,f);
  ideal J=imap(@X,J);
  ideal JS=std(J+var(1));
  ideal b0=kbase(JS);
  int mu=ncols(b0);
  ideal X=imap(@X,X);

  ideal b;
  matrix A;
  module B,B0;
  ideal K,L,M=1,J,1;
  ideal K0,L0,M0=X,X,X;
  module K1,L1,K2,L2;
  module LL1;
  for(i=1;i<=deg(f)-1;i++)
  {
    M=M,M0;
    M0=M0*X;
  }

  ring @S=0,(s,t),(dp,C);
  number a0;
  ideal a;
  int d;
  ideal e,e0;
  intvec s,s0;
  matrix A,B;
  module V,G;

  while(2*a0!=mu*n)
  {
    setring(@XS);

    B=0;
    while(size(B)<mu||size(B0)<mu||maxdegree(b)+deg(f)>k)
    {
      k++;
      K=K,K0;
      K0=K0*X;

      B=0;
      while(size(B)==0||size(B)>mu)
      {
        l++;
        for(i=1;i<=size(L0);i++)
        {
          for(j=1;j<=n;j++)
          {
            L=L,J[j]*L0[i]-var(1)*diff(L0[i],var(j+1));
          }
        }
        L0=L0*X;
        M=M,M0;
        M0=M0*X;

        K1=coeffs(K,K,product(X));
        L1=std(coeffs(L,M,product(X)));
        LL1=jet(lead(L1),0);
        attrib(LL1,"isSB",1);
        K2=simplify(reduce(K1,LL1),2);
        L2=intersect(K2,L1);

        B=pidbasis(K2,L2);
      }
      B0=std(coeffs(reduce(matrix(K,nrows(K),nrows(B))*B,JS),b0));
      b=matrix(K,nrows(K),nrows(B))*B;
    }

    A=lift(B,reduce(coeffs(f*b+var(1)^2*diff(b,var(1)),M,product(X)),L1));
    d=maxdegree(A);
    A=chart(A);

    setring(@S);

    e0,s0,V,B=vfilt(imap(@XS,A),d);
    a,e0,e,s,V,B,G=spec(e0,s0,V,B);

    a0=leadcoef(a[1]);
    for(i=2;i<=mu;i++)
    {
      a0=a0+leadcoef(a[i]);
    }
  }

  G,A=glift(fsplit(e0,e,s,V,B,G));

  setring(@XS);
  b=matrix(b)*imap(@S,G);
  A=imap(@S,A);
  export(b,A);
  kill @S;

  setring(@X);
  return(@XS);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y,z),dp;
  poly f=x+y+z+x2y2z2;
  def Rs=goodBasis(f);
  setring(Rs);
  b;
  print(jet(A,0));
  print(jet(A/var(1),0));
}
///////////////////////////////////////////////////////////////////////////////