///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Singularities";

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

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

OVERVIEW:
A library for computing the Gauss-Manin system of a cohomologically tame
polynomial f. Schulze's algorithm [Sch05], based on Sabbah's theory [Sab98],
is used to compute a good basis of (the Brieskorn lattice of) the Gauss-Manin system and the differential operation of f in terms of this basis.
In addition, there is a test for tameness in the sense of Broughton.
Tame polynomials can be considered as an affine algebraic analogue of local
analytic isolated hypersurface singularities. They have only finitely many
citical points, and those at infinity do not give rise to atypical values
in a sense depending on the precise notion of tameness considered. Well-known
notions of tameness like tameness, M-tameness, Malgrange-tameness, and
cohomological tameness, and their relations, are reviewed in [Sab98,8].
For ordinary tameness, see Broughton [Bro88,3].
Sabbah [Sab98] showed that the Gauss-Manin system, the D-module direct image
of the structure sheaf, of a cohomologically tame polynomial carries a
similar structure as in the isolated singularity case, coming from a Mixed
Hodge structure on the cohomology of the Milnor (typical) fibre (see
gmssing.lib). The data computed by this library encodes the differential structure of the Gauss-Manin system, and the Mixed Hodge structure of the Milnor fibre over the complex numbers. As a consequence, it yields the Hodge numbers, spectral pairs, and monodromy at infinity.

REFERENCES:
[Bro88] S. Broughton: Milnor numbers and the topology of polynomial
        hypersurfaces. Inv. Math. 92 (1988) 217-241.
[Sab98] C. Sabbah: Hypergeometric periods for a tame polynomial.
        arXiv.org math.AG/9805077.
[Sch05] M. Schulze: Good bases for tame polynomials.
        J. Symb. Comp. 39,1 (2005), 103-126.

PROCEDURES:
  isTame(f);     test whether the polynomial f is tame
  goodBasis(f);  good basis of Brieskorn lattice of cohom. tame polynomial 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=
  1;  if f is tame in the sense of Broughton [Bro88,3]
  0;  if f is not tame
@end format
REMARKS:  procedure implements Proposition 3.1 in [Bro88]
KEYWORDS: tame polynomial
EXAMPLE:  example isTame; shows examples
"
{
  int d=vdim(std(jacob(f)));
  def @X=basering;
  int n=nvars(@X);
  ring ext_ring=0,(w(1..n)),dp;
  setring @X;
  def @WX=changechar(ringlist(ext_ring));
  setring @WX;
  kill ext_ring;
  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 vfiltmat(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 in the sense of Sabbah [Sab98,8]
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
REMARKS:  procedure implements Algorithm 6 in [Sch05]
KEYWORDS: tame polynomial; Gauss-Manin system; Brieskorn lattice;
          mixed Hodge structure; V-filtration; weight filtration;
          monodromy; spectrum; spectral pairs; good basis
SEE ALSO: gmssing_lib
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=vfiltmat(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));
}
///////////////////////////////////////////////////////////////////////////////