Changeset 61549b in git


Ignore:
Timestamp:
Nov 5, 2001, 10:17:05 AM (22 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
1f8111cd2acad6f5a631540b24272cf5aca44000
Parents:
5f403d5fffa49541688757a3ff48dc8d052361f2
Message:
*** empty log message ***


git-svn-id: file:///usr/local/Singular/svn/trunk@5663 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular/LIB
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/gaussman.lib

    r5f403d5 r61549b  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: gaussman.lib,v 1.54 2001-08-31 18:18:08 mschulze Exp $";
     2version="$Id: gaussman.lib,v 1.55 2001-11-05 09:16:48 mschulze Exp $";
    33category="Singularities";
    44
    … …  
    1717 monodromy(t);              Jordan data of monodromy of t
    1818 spectrum(t);               spectrum of t
    19  sppnormalize(a,w[,m]);     normalize spectral pairs
    2019 sppairs(t);                spectral pairs of t
     20 sppnf(a,w[,m][,v|V]);      normalize spectral pairs
     21 vfilt(t);                  V-filtration of t on H''/H'
     22 vwfilt(t);                 weight refined V-filtration of t on H''/H'
    2123 saito(t);                  matrix A0+A1*s of t on H''
    22  vfilt(t[,opt]);            V-filtration on H''/H' or spectrum of t
    23  endfilt(t,V);              endomorphism filtration of V-filtration V
     24 endvfilt(V);               endomorphism V-filtration
    2425 spprint(Sp);               print spectrum Sp
    2526 sppprint(Spp);             print spectral pairs Spp
    … …  
    413414  dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
    414415  H=std(H0);
    415   int d0=maxdeg1(H);
    416416
    417417  dbprint(printlevel-voice+2,"// transform H'' to saturation of H''");
    … …  
    419419  H0=l[1];
    420420
    421   return(A0,r,H,H0);
    422 }
    423 ///////////////////////////////////////////////////////////////////////////////
    424 
    425 static proc tmatrix(matrix A0,ideal r,module H,int K,int K0)
     421  return(A0,r,H,H0,k);
     422}
     423///////////////////////////////////////////////////////////////////////////////
     424
     425static proc tmatrix(matrix A0,ideal r,module H,int k0,int K,int K0)
    426426{
    427427  dbprint(printlevel-voice+2,"// compute matrix A of t");
    428   int d0=maxdeg1(H);
    429   dbprint(printlevel-voice+2,"// k="+string(K+d0+1));
    430   list l=gmscoeffs(r,K+d0+1,K0+d0+1);
     428  dbprint(printlevel-voice+2,"// k="+string(K+k0+1));
     429  list l=gmscoeffs(r,K+k0+1,K0+k0+1);
    431430  matrix C;
    432431  C,r=l[1..2];
    … …  
    441440///////////////////////////////////////////////////////////////////////////////
    442441
    443 static proc eigenvals(matrix A0,ideal r,module H,int K0)
     442static proc eigenvals(matrix A0,ideal r,module H,int k0,int K0)
    444443{
    445444  dbprint(printlevel-voice+2,
    446445    "// compute eigenvalues e with multiplicities m of A");
    447446  matrix A;
    448   A,A0,r=tmatrix(A0,r,H,0,K0);
     447  A,A0,r=tmatrix(A0,r,H,k0,0,K0);
    449448  list l=eigenvalues(A);
    450449  def e,m=l[1..2];
    … …  
    452451  dbprint(printlevel-voice+2,"// m="+string(m));
    453452
    454   return(e,m,A0,r);
    455 }
    456 ///////////////////////////////////////////////////////////////////////////////
    457 
    458 static proc transform(ideal e,intvec m,matrix A,matrix A0,ideal r,module H,module H0,int K,int K0)
    459 {
     453  return(e,m,A0,r,int(max(e)-min(e)));
     454}
     455///////////////////////////////////////////////////////////////////////////////
     456
     457static proc transform(matrix A,matrix A0,ideal r,module H,module H0,ideal e,intvec m,int k0,int k1,int K,int K0)
     458{
     459  int mu=ncols(gmsbasis);
     460
    460461  dbprint(printlevel-voice+2,"// compute minimum e0 and maximum e1 of e");
    461462  number e0,e1=min(e),max(e);
    462463  dbprint(printlevel-voice+2,"// e0="+string(e0));
    463464  dbprint(printlevel-voice+2,"// e1="+string(e1));
    464   int d1=int(e1-e0);
    465   A,A0,r=tmatrix(A0,r,H,K+d1,K0+d1);
     465  A,A0,r=tmatrix(A0,r,H,k0,K+k1,K0+k1);
     466  module U0=s^k0*freemodule(mu);
    466467
    467468  if(e1>=e0+1)
    … …  
    482483      A=V*A*U;
    483484      H0=V*H0;
     485      U0=U0*U;
    484486
    485487      dbprint(printlevel-voice+2,"// transform to reduce e1 by 1");
    … …  
    514516            A[i,i]=A[i,i]-1;
    515517            H0[i]=H0[i]*s;
     518            U0[i]=U0[i]/s;
    516519          }
    517520          e[i0]=e[i0]-1;
    … …  
    524527      H0=transpose(H0);
    525528
    526       l=spnormalize(e,m);
     529      l=spnf(e,m);
    527530      e,m=l[1..2];
    528531
    … …  
    534537  }
    535538
    536   return(e,m,A,A0,r,H0);
     539  return(A,A0,r,H0,U0,e,m);
    537540}
    538541///////////////////////////////////////////////////////////////////////////////
    … …  
    553556  setring(G);
    554557
    555   int mu=ncols(gmsbasis);
    556558  matrix A;
     559  module U0;
    557560  ideal e;
    558561  intvec m;
    559 
    560   def A0,r,H,H0=saturate(n);
    561   e,m,A0,r=eigenvals(A0,r,H,n);
    562   e,m,A,A0,r,H0=transform(e,m,A,A0,r,H,H0,0,0);
     562  int k1;
     563
     564  def A0,r,H,H0,k0=saturate(n);
     565  e,m,A0,r,k1=eigenvals(A0,r,H,k0,n);
     566  A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,k1,0,0);
    563567
    564568  setring(R);
    565   matrix A=imap(G,A);
    566   ideal e=imap(G,e);
    567 
    568   return(jordan(A));
     569  return(jordan(imap(G,A),imap(G,e),m));
    569570}
    570571example
    … …  
    583584@format
    584585list Sp: spectrum of t
    585   ideal Sp[1]: spectral numbers in increasing order
    586   intvec Sp[2]:
    587     int Sp[2][i]: multiplicity of spectral number Sp[1][i]
     586  ideal Sp[1]: V-filtration indices in increasing order
     587  intvec Sp[2]: weight filtration indices in decreasing order
     588  intvec Sp[3]:
     589    int Sp[3][i]: multiplicity of spectral number Sp[1][i]
    588590@end format
    589591SEE ALSO: spectrum_lib
    590 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; spectrum
    591 EXAMPLE:  example spnumbers; shows examples
    592 "
    593 {
    594   list l=sppairs(t);
    595   return(spnormalize(l[1],l[3]));
     592KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
     593          mixed Hodge structure; spectrum
     594EXAMPLE:  example spectrum; shows examples
     595"
     596{
     597  list l=vwfilt(t);
     598  return(spnf(l[1],l[3]));
    596599}
    597600example
    … …  
    600603  poly t=x5+x2y2+y5;
    601604  spprint(spectrum(t));
    602 }
    603 ///////////////////////////////////////////////////////////////////////////////
    604 
    605 static proc sppappend(list l,number a,int w,int m)
    606 {
    607   if(size(l)==0)
    608   {
    609     l=list(ideal(a),intvec(w),intvec(m));
    610   }
    611   else
    612   {
    613     int n=ncols(l[1]);
    614     l[1][n+1]=a;
    615     l[2][n+1]=w;
    616     l[3][n+1]=m;
    617   }
    618   return(l);
    619 }
    620 ///////////////////////////////////////////////////////////////////////////////
    621 
    622 proc sppnormalize(ideal a,intvec w,list #)
    623 "USAGE:    sppnormalize(a,w[,m]);
    624 RETURN:
    625 @format
    626 list Spp: normalized spectral pairs (a,w,m)
    627   ideal Spp[1]: numbers in a in increasing order
    628   intvec Spp[2]: integers in w in decreasing order
    629   intvec Spp[3]:
    630     int Spp[3][i]: multiplicity of pair (Spp[1][i],Spp[2][i]) in a,w
    631 @end format
    632 EXAMPLE:  example sppnormalize; shows examples
    633 "
    634 {
    635   intvec m;
    636   int i,j;
    637   if(size(#)==0)
    638   {
    639     for(i=ncols(a);i>=1;i--)
    640     {
    641       m[i]=1;
    642     }
    643   }
    644   else
    645   {
    646     m=#[1];
    647   }
    648 
    649   list l;
    650   number a0;
    651   int w0,m0;
    652   for(i=ncols(a);i>=1;i--)
    653   {
    654     if(m[i]!=0)
    655     {
    656       for(j=i-1;j>=1;j--)
    657       {
    658         if(m[j]!=0)
    659         {
    660           if(number(a[i])>number(a[j])||
    661             (number(a[i])==number(a[j])&&w[i]<w[j]))
    662           {
    663             a0=number(a[i]);
    664             a[i]=a[j];
    665             a[j]=a0;
    666             w0=w[i];
    667             w[i]=w[j];
    668             w[j]=w0;
    669             m0=m[i];
    670             m[i]=m[j];
    671             m[j]=m0;
    672           }
    673           if(number(a[i])==number(a[j])&&w[i]==w[j])
    674           {
    675             m[i]=m[i]+m[j];
    676             m[j]=0;
    677           }
    678         }
    679       }
    680       l=sppappend(l,number(a[i]),w[i],m[i]);
    681     }
    682   }
    683 
    684   return(l);
    685 }
    686 example
    687 { "EXAMPLE:"; echo=2;
    688605}
    689606///////////////////////////////////////////////////////////////////////////////
    … …  
    695612RETURN:
    696613@format
    697 list Spp:
    698   ideal Spp[1]: spectral numbers in increasing order
    699   intvec Spp[2]: weight filtration indices in decreasing order
     614list Spp: spectrum of t
     615  ideal Spp[1],intvec Spp[2]: spectral pairs in in-/decreasing lex. order
    700616  intvec Spp[3]:
    701617    int Spp[3][i]: multiplicity of spectral pair (Spp[1][i],Spp[2][i])
    … …  
    703619SEE ALSO: spectrum_lib
    704620KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
     621          mixed Hodge structure; spectrum; spectral pairs
     622EXAMPLE:  example sppairs; shows examples
     623"
     624{
     625  list l=vwfilt(t);
     626  return(list(l[1],l[2],l[3]));
     627}
     628example
     629{ "EXAMPLE:"; echo=2;
     630  ring R=0,(x,y),ds;
     631  poly t=x5+x2y2+y5;
     632  sppprint(sppairs(t));
     633}
     634///////////////////////////////////////////////////////////////////////////////
     635
     636proc sppnf(ideal a,intvec w,list #)
     637"USAGE:    sppnf(a,w[,m][,v|V]); ideal a, intvec w, intvec m, module v, list V
     638RETURN:
     639@format
     640list Spp: normalized spectral pairs (a,w,m[,V])
     641  ideal Spp[1]: numbers in a in increasing order
     642  intvec Spp[2]: integers in w in decreasing order
     643  intvec Spp[3]:
     644    int Spp[3][i]: multiplicity of pair (Spp[1][i],Spp[2][i]) in (a,w)
     645  list Spp[4]:
     646    module Spp[4][i]: module associated to pair (Spp[1][i],Spp[2][i])
     647@end format
     648EXAMPLE:  example sppnorm; shows examples
     649"
     650{
     651  int n=ncols(a);
     652  intvec m;
     653  module v;
     654  list V;
     655  int i,j;
     656  while(i<size(#))
     657  {
     658    i++;
     659    if(typeof(#[i])=="intvec")
     660    {
     661      m=#[i];
     662    }
     663    if(typeof(#[i])=="module")
     664    {
     665      v=#[i];
     666      for(j=n;j>=1;j--)
     667      {
     668        V[j]=module(v[j]);
     669      }
     670    }
     671    if(typeof(#[i])=="list")
     672    {
     673      V=#[i];
     674    }
     675  }
     676  if(m==0)
     677  {
     678    for(i=n;i>=1;i--)
     679    {
     680      m[i]=1;
     681    }
     682  }
     683
     684  int k;
     685  ideal a0;
     686  intvec w0,m0;
     687  list V0;
     688  number a1;
     689  int w1,m1;
     690  for(i=n;i>=1;i--)
     691  {
     692    if(m[i]!=0)
     693    {
     694      for(j=i-1;j>=1;j--)
     695      {
     696        if(m[j]!=0)
     697        {
     698          if(number(a[i])>number(a[j])||
     699            (number(a[i])==number(a[j])&&w[i]<w[j]))
     700          {
     701            a1=number(a[i]);
     702            a[i]=a[j];
     703            a[j]=a1;
     704            w1=w[i];
     705            w[i]=w[j];
     706            w[j]=w1;
     707            m1=m[i];
     708            m[i]=m[j];
     709            m[j]=m1;
     710            if(size(V)>0)
     711            {
     712              v=V[i];
     713              V[i]=V[j];
     714              V[j]=v;
     715            }
     716          }
     717          if(number(a[i])==number(a[j])&&w[i]==w[j])
     718          {
     719            m[i]=m[i]+m[j];
     720            m[j]=0;
     721            if(size(V)>0)
     722            {
     723              V[i]=V[i]+V[j];
     724            }
     725          }
     726        }
     727      }
     728      k++;
     729      a0[k]=a[i];
     730      w0[k]=w[i];
     731      m0[k]=m[i];
     732      if(size(V)>0)
     733      {
     734        V0[k]=V[i];
     735      }
     736    }
     737  }
     738
     739  if(size(V0)>0)
     740  {
     741    n=size(V0);
     742    module U=std(V0[n]);
     743    for(i=n-1;i>=1;i--)
     744    {
     745      V0[i]=simplify(reduce(V0[i],U),1);
     746      if(i>=2)
     747      {
     748        U=std(U+V0[i]);
     749      }
     750    }
     751  }
     752
     753  list l;
     754  if(k>0)
     755  {
     756    l=a0,w0,m0;
     757    if(size(V0)>0)
     758    {
     759      l[4]=V0;
     760    }
     761  }
     762  return(l);
     763}
     764example
     765{ "EXAMPLE:"; echo=2;
     766}
     767///////////////////////////////////////////////////////////////////////////////
     768
     769proc vfilt(poly t)
     770"USAGE:    vfilt(t); poly t
     771ASSUME:   basering with characteristic 0 and local degree ordering,
     772          t with isolated citical point 0
     773RETURN:
     774@format
     775list V: V-filtration on H''/H'
     776  ideal V[1]: spectral numbers in increasing order
     777  intvec V[2]:
     778    int V[2][i]: multiplicity of spectral number V[1][i]
     779  list V[4]:
     780    module V[4][i]: vector space basis of V[1][i]-th graded part
     781                    in terms of V[4]
     782  ideal V[4]: monomial vector space basis of H''/H'
     783  ideal V[5]: standard basis of Jacobian ideal
     784@end format
     785SEE ALSO: spectrum_lib
     786KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
     787          mixed Hodge structure; V-filtration; spectrum
     788EXAMPLE:  example vfilt; shows examples
     789"
     790{
     791  list l=vwfilt(t);
     792  return(spnf(l[1],l[3],l[4])+list(l[5],l[6]));
     793}
     794example
     795{ "EXAMPLE:"; echo=2;
     796  ring R=0,(x,y),ds;
     797  poly t=x5+x2y2+y5;
     798  vfilt(t);
     799}
     800///////////////////////////////////////////////////////////////////////////////
     801
     802proc vwfilt(poly t)
     803"USAGE:    vwfilt(t); poly t
     804ASSUME:   basering with characteristic 0 and local degree ordering,
     805          t with isolated citical point 0
     806RETURN:
     807@format
     808list VW: weight refined V-filtration on H''/H'
     809  ideal VW[1]: spectral numbers in increasing order
     810  intvec VW[2]: weights in decreasing order
     811  intvec VW[3]:
     812    int VW[3][i]: multiplicity of spectral pair (VW[1][i],VW[2][i])
     813  list VW[4]:
     814    module VW[4][i]: vector space basis of (VW[1][i],VW[2][i])-th graded part
     815                     in terms of VW[5]
     816  ideal VW[5]: monomial vector space basis of H''/H'
     817  ideal VW[6]: standard basis of Jacobian ideal
     818@end format
     819SEE ALSO: spectrum_lib
     820KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
     821          mixed Hodge structure; V-filtration; weight filtration;
    705822          spectrum; spectral pairs
    706 EXAMPLE:  example sppairs; shows examples
     823EXAMPLE:  example vwfilt; shows examples
    707824"
    708825{
    … …  
    714831  int mu=ncols(gmsbasis);
    715832  matrix A;
     833  module U0;
    716834  ideal e;
    717835  intvec m;
    718 
    719   def A0,r,H,H0=saturate(n);
    720   e,m,A0,r=eigenvals(A0,r,H,n);
    721   e,m,A,A0,r,H0=transform(e,m,A,A0,r,H,H0,0,0);
     836  int k1;
     837
     838  def A0,r,H,H0,k0=saturate(n);
     839  e,m,A0,r,k1=eigenvals(A0,r,H,k0,n);
     840  A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,k1,0,0);
    722841
    723842  dbprint(printlevel-voice+2,"// compute weight filtration basis");
    724843  list l=jordanbasis(A,e,m);
    725844  def U,v=l[1..2];
     845  kill l;
    726846  vector u0;
    727847  int v0;
    728   int i,j,k,k0;
    729   for(k,k0=1,1;k0<=ncols(e);k,k0=k+m[k0],k0+1)
    730   {
    731     for(i=k+m[k0]-1;i>=k+1;i--)
     848  int i,j,k,l;
     849  for(k,l=1,1;l<=ncols(e);k,l=k+m[l],l+1)
     850  {
     851    for(i=k+m[l]-1;i>=k+1;i--)
    732852    {
    733853      for(j=i-1;j>=k;j--)
    … …  
    747867  dbprint(printlevel-voice+2,"// compute normal form of H''");
    748868  H0=std(V*H0);
     869  U0=U0*U;
    749870
    750871  dbprint(printlevel-voice+2,"// compute spectral pairs");
    … …  
    754875  {
    755876    j=leadexp(H0[i])[nvars(basering)+1];
    756     a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)-1;
     877    a[i]=A[j,j]+ord(H0[i])/deg(s)-1;
    757878    w[i]=v[j]+n;
    758879  }
     880  kill v;
     881  module v=simplify(grmat(H*U0*H0,2*k0),1);
    759882
    760883  setring(R);
    761 
    762   return(sppnormalize(imap(G,a),w));
     884  ideal g=imap(G,gmsstd);
     885  attrib(g,"isSB",1);
     886  return(sppnf(imap(G,a),w,imap(G,v))+list(imap(G,gmsbasis),g));
    763887}
    764888example
    … …  
    766890  ring R=0,(x,y),ds;
    767891  poly t=x5+x2y2+y5;
    768   sppprint(sppairs(t));
    769 }
    770 ///////////////////////////////////////////////////////////////////////////////
    771 
    772 static proc grmat(matrix A,int k)
    773 {
    774   int i,j;
    775   for(i=1;i<=ncols(A);i++)
    776   {
    777     for(j=1;j<=nrows(A);j++)
    778     {
    779       A[i,j]=jet(A[i,j]/var(1)^k,0);
    780     }
    781   }
    782   return(A);
    783 }
    784 ///////////////////////////////////////////////////////////////////////////////
    785 
    786 proc saito(poly t,list #)
    787 "USAGE:    saito(t); poly t
    788 ASSUME:   basering with characteristic 0 and local degree ordering,
    789           t with isolated citical point 0
    790 RETURN:   list A: matrix A[1]+A[2]*s of t on H''
    791 SEE ALSO: spectrum_lib
    792 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
    793           monodromy; spectrum; spectral pairs
    794 EXAMPLE:  example saito; shows examples
    795 "
    796 {
    797   def R=basering;
    798   int n=nvars(R)-1;
    799   def G=gmsring(t,"s");
    800   setring(G);
    801 
    802   int mu=ncols(gmsbasis);
    803   matrix A;
    804   ideal e;
    805   intvec m;
    806 
    807   def A0,r,H,H0=saturate(2*n+mu+1);
    808   e,m,A0,r=eigenvals(A0,r,H,n+1);
    809   int d0,d1=maxdeg1(H),int(max(e)-min(e));
    810   e,m,A,A0,r,H0=transform(e,m,A,A0,r,H,H0,d0+d1+1,d0+d1+1);
    811   d0=maxdeg1(H0);
    812 
    813   dbprint(printlevel-voice+2,"// transform to Jordan basis");
    814   module U=jordanbasis(A,e,m)[1];
    815   module V=inverse(U);
    816   A=V*A*U;
    817   H=V*H0;
    818 
    819   dbprint(printlevel-voice+2,"// compute splitting of V filtration");
    820   int i,j,k;
    821   U=freemodule(mu);
    822   matrix U0[mu][mu];
    823   matrix u0[mu^2][1];
    824   A0=commutator(jet(A,0));
    825   for(k=1;k<=d0+1;k++)
    826   {
    827     V=0;
    828     for(j=0;j<k;j++)
    829     {
    830       U0=matrix(U0)-grmat(A,k-j)*grmat(U,j);
    831     }
    832     u0=U0[1..mu,1..mu];
    833     u0=inverse(A0+k)*u0;
    834     U0=u0[1..mu^2,1];
    835     U=matrix(U)+s^k*U0;
    836   }
    837 
    838   dbprint(printlevel-voice+2,"// transform to splitting basis");
    839   A=jet(A,0);
    840   H=std(jet(division(U,H)[1],d0));
    841 
    842   dbprint(printlevel-voice+2,"// compute N");
    843   matrix N=A;
    844   for(i=1;i<=ncols(N);i++)
    845   {
    846     N[i,i]=0;
    847   }
    848 
    849   dbprint(printlevel-voice+2,"// compute N splitting of Hodge filtration");
    850   module F;
    851   matrix K[mu][1];
    852   U=0;
    853   i=1;
    854   for(k=1;k<=ncols(e);k++)
    855   {
    856     j=0;
    857     K=matrix(0,mu,0);
    858     while(j<m[k])
    859     {
    860       K=K+gen(i);
    861       i++;
    862       j++;
    863     }
    864     j=0;
    865     F=grmat(H,0);
    866     V=intersect(F,K);
    867     while(size(V)==0)
    868     {
    869       j++;
    870       V=V+grmat(H,j);
    871       V=intersect(F,K)
    872     }
    873     while(size(V)>0)
    874     {
    875       U=U+V;
    876       V=N*V;
    877     }
    878   }
    879 
    880   dbprint(printlevel-voice+2,"// transform to N splitting basis");
    881   V=inverse(U);
    882   A=V*A*U;
    883   H=V*H*U;
    884 
    885   dbprint(printlevel-voice+2,"// compute matrix A0+sA1 of t");
    886 //  degBound=d0;
    887 //  option("redSB");
    888 //  H=std(H);
    889 H=simplify(lead(std(H)),1);
    890 //  degBound=0;
    891   list l=division(H,s*A*H+s^2*diff(matrix(H),s));
    892   A=jet(l[1],l[2],d0+1);
    893   A0=grmat(A,0);
    894   A=grmat(A,1);
    895 
    896   setring(R);
    897 
    898   return(list(imap(G,A0),imap(G,A)));
    899 }
    900 example
    901 { "EXAMPLE:"; echo=2;
    902   ring R=0,(x,y),ds;
    903   poly t=x5+x2y2+y5;
    904   list A=saito(t);
    905   print(A[1]);
    906   print(A[2]);
    907 }
    908 ///////////////////////////////////////////////////////////////////////////////
    909 
    910 proc vfilt(poly t,list #)
    911 "USAGE:    vfilt(t[,opt]); poly t, int opt
     892  vwfilt(t);
     893}
     894///////////////////////////////////////////////////////////////////////////////
     895
     896proc vfilt1(poly t,list #)
     897"USAGE:    vfilt1(t[,opt]); poly t, int opt
    912898ASSUME:   basering with characteristic 0 and local degree ordering,
    913899          t with isolated citical point 0
    … …  
    930916KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
    931917          Hodge filtration; V-filtration; spectrum
    932 EXAMPLE:  example vfilt; shows examples
     918EXAMPLE:  example vfilt1; shows examples
    933919"
    934920{
    … …  
    11901176  ring R=0,(x,y),ds;
    11911177  poly t=x5+x2y2+y5;
    1192   vfilt(t);
    1193 }
    1194 ///////////////////////////////////////////////////////////////////////////////
    1195 
    1196 proc endfilt(list V)
    1197 "USAGE:   endfilt(V); list V
     1178  vfilt1(t);
     1179}
     1180///////////////////////////////////////////////////////////////////////////////
     1181
     1182static proc grmat(matrix A,int k)
     1183{
     1184  int i,j;
     1185  for(i=1;i<=ncols(A);i++)
     1186  {
     1187    for(j=1;j<=nrows(A);j++)
     1188    {
     1189      A[i,j]=jet(A[i,j]/var(1)^k,0);
     1190    }
     1191  }
     1192  return(A);
     1193}
     1194///////////////////////////////////////////////////////////////////////////////
     1195
     1196proc saito(poly t,list #)
     1197"USAGE:    saito(t); poly t
     1198ASSUME:   basering with characteristic 0 and local degree ordering,
     1199          t with isolated citical point 0
     1200RETURN:   list A: matrix A[1]+A[2]*s of t on H''
     1201SEE ALSO: spectrum_lib
     1202KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
     1203          mixed Hodge structure; V-filtration; weight filtration;
     1204          spectrum; spectral pairs
     1205EXAMPLE:  example saito; shows examples
     1206"
     1207{
     1208  def R=basering;
     1209  int n=nvars(R)-1;
     1210  def G=gmsring(t,"s");
     1211  setring(G);
     1212
     1213  int mu=ncols(gmsbasis);
     1214  matrix A;
     1215  module U0;
     1216  ideal e;
     1217  intvec m;
     1218  int k1;
     1219
     1220  def A0,r,H,H0,k0=saturate(2*n+mu-1);
     1221  e,m,A0,r,k1=eigenvals(A0,r,H,k0,n);
     1222  A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,k1,k0+k1,k0+k1);
     1223
     1224  dbprint(printlevel-voice+2,"// transform to Jordan basis");
     1225  module U=jordanbasis(A,e,m)[1];
     1226  matrix V=inverse(U);
     1227  A=V*A*U;
     1228  H=V*H0;
     1229
     1230  dbprint(printlevel-voice+2,"// compute splitting of V-filtration");
     1231  int i,j,k;
     1232  U=freemodule(mu);
     1233  V=matrix(0,mu,mu);
     1234  matrix v[mu^2][1];
     1235  A0=commutator(jet(A,0));
     1236  for(k=1;k<=k0+k1;k++)
     1237  {
     1238    for(j=0;j<k;j++)
     1239    {
     1240      V=matrix(V)-grmat(A,k-j)*grmat(U,j);
     1241    }
     1242    v=V[1..mu,1..mu];
     1243    v=inverse(A0+k)*v;
     1244    V=v[1..mu^2,1];
     1245    U=matrix(U)+s^k*V;
     1246  }
     1247
     1248  dbprint(printlevel-voice+2,"// transform to V-splitting basis");
     1249  A=jet(A,0);
     1250  H=jet(division(U,H)[1],k0+k1);
     1251  H=std(H);
     1252
     1253  dbprint(printlevel-voice+2,"// compute V-leading terms of H''");
     1254  int i0,j0;
     1255  module H1=H;
     1256  for(k=ncols(H1);k>=1;k--)
     1257  {
     1258    i0=leadexp(H1[k])[nvars(basering)+1];
     1259    j0=ord(H1[k])/deg(s);
     1260    H0[k]=lead(H1[k]);
     1261    H1[k]=H1[k]-lead(H1[k]);
     1262    if(H1[k]!=0)
     1263    {
     1264      i=leadexp(H1[k])[nvars(basering)+1];
     1265      j=ord(H1[k])/deg(s);
     1266      while(A[i,i]+j==A[i0,i0]+j0)
     1267      {
     1268        H0[k]=H0[k]+lead(H1[k]);
     1269        H1[k]=H1[k]-lead(H1[k]);
     1270        i=leadexp(H1[k])[nvars(basering)+1];
     1271        j=ord(H1[k])/deg(s);
     1272      }
     1273    }
     1274  }
     1275  H0=simplify(H0,1);
     1276
     1277  dbprint(printlevel-voice+2,"// compute N");
     1278  matrix N=A;
     1279  for(i=1;i<=ncols(N);i++)
     1280  {
     1281    N[i,i]=0;
     1282  }
     1283
     1284  dbprint(printlevel-voice+2,"// compute splitting of Hodge filtration");
     1285  U=0;
     1286  module U1;
     1287  module C;
     1288  list F,I;
     1289  module F0,I0;
     1290  for(i0,j0=1,1;i0<=ncols(e);i0++)
     1291  {
     1292    C=matrix(0,mu,1);
     1293    for(j=m[i0];j>=1;j,j0=j-1,j0+1)
     1294    {
     1295      C=C+gen(j0);
     1296    }
     1297    F0=intersect(C,H0);
     1298    F=list();
     1299    j=0;
     1300    while(size(F0)>0)
     1301    {
     1302      j++;
     1303      F[j]=matrix(0,mu,1);
     1304      if(size(jet(F0,0))>0)
     1305      {
     1306        for(i=ncols(F0);i>=1;i--)
     1307        {
     1308          if(ord(F0[i])==0)
     1309          {
     1310            F[j]=F[j]+F0[i];
     1311          }
     1312        }
     1313      }
     1314      for(i=ncols(F0);i>=1;i--)
     1315      {
     1316        F0[i]=F0[i]/s;
     1317      }
     1318    }
     1319
     1320    I=list();
     1321    I0=module();
     1322    U0=std(0);
     1323    for(i=size(F);i>=1;i--)
     1324    {
     1325      I[i]=module();
     1326    }
     1327    for(i=1;i<=size(F);i++)
     1328    {
     1329      I0=reduce(F[i],U0);
     1330      j=i;
     1331      while(size(I0)>0)
     1332      {
     1333        U0=std(U0+I0);
     1334        I[j]=I[j]+I0;
     1335        I0=reduce(N*I0,U0);
     1336        j++;
     1337      }
     1338    }
     1339
     1340    for(i=1;i<=size(I);i++)
     1341    {
     1342      U=U+I[i];
     1343    }
     1344  }
     1345
     1346  dbprint(printlevel-voice+2,"// transform to Hodge splitting basis");
     1347  V=inverse(U);
     1348  A=V*A*U;
     1349  H=V*H;
     1350
     1351  dbprint(printlevel-voice+2,"// compute reduced standard basis of H''");
     1352  ring S=0,s,ds;
     1353  module H=imap(G,H);
     1354  degBound=k0+k1+1;
     1355  option("redSB");
     1356  H=std(H);
     1357  degBound=0;
     1358  H=simplify(jet(H,k0+k1),1);
     1359  setring(G);
     1360  H=imap(S,H);
     1361
     1362  dbprint(printlevel-voice+2,"// compute matrix A0+sA1 of t");
     1363  list l=division(H,s*A*H+s^2*diff(matrix(H),s));
     1364  A=jet(l[1],l[2],k0+k1+1);
     1365  A0=grmat(A,0);
     1366  A=grmat(A,1);
     1367
     1368  setring(R);
     1369  return(list(imap(G,A0),imap(G,A)));
     1370}
     1371example
     1372{ "EXAMPLE:"; echo=2;
     1373  ring R=0,(x,y),ds;
     1374  poly t=x5+x2y2+y5;
     1375  list A=saito(t);
     1376  print(A[1]);
     1377  print(A[2]);
     1378}
     1379///////////////////////////////////////////////////////////////////////////////
     1380
     1381proc endvfilt(list V)
     1382"USAGE:   endvwfilt(V); list V
    11981383ASSUME:  V computed by vfilt
    11991384RETURN:
    12001385@format
    1201 list V1: endomorphim filtration of V on the Jacobian algebra
    1202   ideal V1[1]: spectral numbers in increasing order
    1203   intvec V1[2]:
    1204     int V1[2][i]: multiplicity of spectral number V1[1][i]
    1205   list V1[3]:
    1206     module V1[3][i]: vector space basis of the V1[1][i]-th graded part
    1207                      in terms of V1[4]
    1208   ideal V1[4]: monomial vector space basis
     1386list EV: endomorphism V-filtration on the Jacobian algebra
     1387  ideal EV[1]: spectral numbers in increasing order
     1388  intvec EV[2]:
     1389    int EV[2][i]: multiplicity of spectral pair (EV[1][i],EV[2][i])
     1390  list EV[3]:
     1391    module EV[3][i]: vector space basis of the (EV[1][i],EV[2][i])-th
     1392                      graded part in terms of EV[4]
     1393  ideal EV[4]: monomial vector space basis
     1394  ideal EV[5]: standard basis of Jacobian ideal
    12091395@end format
    12101396SEE ALSO: spectrum_lib
    1211 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; spectrum;
    1212           Hodge filtration; V-filtration
    1213 EXAMPLE: example endfilt; shows examples
     1397KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
     1398          mixed Hodge structure; V-filtration; weight filtration;
     1399          endomorphism filtration
     1400EXAMPLE: example endvwfilt; shows examples
    12141401"
    12151402{
    … …  
    12691456    }
    12701457
    1271     dbprint(printlevel-voice+2,"// collect conditions for V1["+string(b0)+"]");
     1458    dbprint(printlevel-voice+2,"// collect conditions for EV["+string(b0)+"]");
    12721459    j=ncols(a);
    12731460    j0=mu;
    … …  
    13381525    i++;
    13391526  }
    1340   return(list(a1,d1,v1,m));
     1527  return(list(a1,d1,v1,m,g));
    13411528}
    13421529example
    … …  
    13441531  ring R=0,(x,y),ds;
    13451532  poly t=x5+x2y2+y5;
    1346   endfilt(vfilt(t));
     1533  endvfilt(vfilt(t));
    13471534}
    13481535///////////////////////////////////////////////////////////////////////////////

  • Singular/LIB/linalg.lib

    r5f403d5 r61549b  
    11//GMG last modified: 04/25/2000
    22//////////////////////////////////////////////////////////////////////////////
    3 version="$Id: linalg.lib,v 1.19 2001-08-25 14:56:41 mschulze Exp $";
     3version="$Id: linalg.lib,v 1.20 2001-11-05 09:17:05 mschulze Exp $";
    44category="Linear Algebra";
    55info="
    … …  
    2626 pos_def(A,i);        test symmetric matrix for positive definiteness
    2727 hessenberg(M);       transforms M to Hessenberg form
    28  spnormalize(a[,m]);  normalize spectrum
    29  eigenval(M);         eigenvalues of M with multiplicities
     28 spnf(a[,m][,v|V]);   normalize spectrum
     29 eigenvalues(M);      eigenvalues of M with multiplicities
    3030 jordan(M);           Jordan data of constant square matrix M
    3131 jordanbasis(M);      Jordan basis of constant square matrix M
    … …  
    13351335///////////////////////////////////////////////////////////////////////////////
    13361336
    1337 static proc spappend(list l,number e,int m)
    1338 {
    1339   if(size(l)==0)
    1340   {
    1341     l=list(ideal(e),intvec(m));
    1342   }
    1343   else
    1344   {
    1345     int n=ncols(l[1]);
    1346     l[1][n+1]=e;
    1347     l[2][n+1]=m;
    1348   }
    1349   return(l);
    1350 }
    1351 ///////////////////////////////////////////////////////////////////////////////
    1352 
    1353 proc spnormalize(ideal e,list #)
    1354 "USAGE:    spnormalize(a,w[,m]);
     1337proc spnf(ideal e,list #)
     1338"USAGE:   spnf(e[,m][,v|V]); ideal e, intvec m, module v, list V
    13551339RETURN:
    13561340@format
    1357 list Sp: normalized spectrum (a,m)
     1341list Sp: normalized spectrum (e,m[,V])
    13581342  ideal Sp[1]: numbers in a in increasing order
    13591343  intvec Sp[2]:
    1360     int Sp[2][i]: multiplicity of number Sp[1][i] in a
     1344    int Sp[2][i]: multiplicity of number Sp[1][i] in e
     1345  list Spp[3]:
     1346    module Spp[3][i]: module associated to number Spp[1][i]
    13611347@end format
    1362 EXAMPLE:  example spnormalize; shows examples
     1348EXAMPLE: example spnf; shows examples
    13631349"
    13641350{
     1351  int n=ncols(e);
    13651352  intvec m;
     1353  module v;
     1354  list V;
    13661355  int i,j;
    1367   if(size(#)==0)
    1368   {
    1369     for(i=ncols(e);i>=1;i--)
     1356  while(i<size(#))
     1357  {
     1358    i++;
     1359    if(typeof(#[i])=="intvec")
     1360    {
     1361      m=#[i];
     1362    }
     1363    if(typeof(#[i])=="module")
     1364    {
     1365      v=#[i];
     1366      for(j=n;j>=1;j--)
     1367      {
     1368        V[j]=module(v[j]);
     1369      }
     1370    }
     1371    if(typeof(#[i])=="list")
     1372    {
     1373      V=#[i];
     1374    }
     1375  }
     1376  if(m==0)
     1377  {
     1378    for(i=n;i>=1;i--)
    13701379    {
    13711380      m[i]=1;
    13721381    }
    13731382  }
    1374   else
    1375   {
    1376     m=#[1];
    1377   }
    1378 
    1379   list l;
    1380   number e0;
    1381   int m0;
    1382   for(i=ncols(e);i>=1;i--)
     1383
     1384  int k;
     1385  ideal e0;
     1386  intvec m0;
     1387  list V0;
     1388  number e1;
     1389  int m1;
     1390  for(i=n;i>=1;i--)
    13831391  {
    13841392    if(m[i]!=0)
    … …  
    13901398          if(number(e[i])>number(e[j]))
    13911399          {
    1392             e0=number(e[i]);
     1400            e1=number(e[i]);
    13931401            e[i]=e[j];
    1394             e[j]=e0;
    1395             m0=m[i];
     1402            e[j]=e1;
     1403            m1=m[i];
    13961404            m[i]=m[j];
    1397             m[j]=m0;
     1405            m[j]=m1;
     1406            if(size(V)>0)
     1407            {
     1408              v=V[i];
     1409              V[i]=V[j];
     1410              V[j]=v;
     1411            }
    13981412          }
    13991413          if(number(e[i])==number(e[j]))
    … …  
    14011415            m[i]=m[i]+m[j];
    14021416            m[j]=0;
     1417            if(size(V)>0)
     1418            {
     1419              V[i]=V[i]+V[j];
     1420            }
    14031421          }
    14041422        }
    14051423      }
    1406       l=spappend(l,number(e[i]),m[i]);
    1407     }
    1408   }
    1409 
     1424      k++;
     1425      e0[k]=e[i];
     1426      m0[k]=m[i];
     1427      if(size(V)>0)
     1428      {
     1429        V0[k]=V[i];
     1430      }
     1431    }
     1432  }
     1433
     1434  if(size(V0)>0)
     1435  {
     1436    n=size(V0);
     1437    module U=std(V0[n]);
     1438    for(i=n-1;i>=1;i--)
     1439    {
     1440      V0[i]=simplify(reduce(V0[i],U),1);
     1441      if(i>=2)
     1442      {
     1443        U=std(U+V0[i]);
     1444      }
     1445    }
     1446  }
     1447
     1448  list l;
     1449  if(k>0)
     1450  {
     1451    l=e0,m0;
     1452    if(size(V0)>0)
     1453    {
     1454      l[3]=V0;
     1455    }
     1456  }
    14101457  return(l);
    14111458}
    … …  
    14251472    int l[2][i]: multiplicity of eigenvalue l[1][i]
    14261473@end format
    1427 EXAMPLE: example eigenval; shows examples
     1474EXAMPLE: example eigenvalues; shows examples
    14281475"
    14291476{
    14301477  M=jet(hessenberg(M),0);
    14311478  int n=ncols(M);
     1479  int k;
     1480  ideal e;
     1481  intvec m;
    14321482  number e0;
    14331483  intvec v;
    1434   list l,f;
    1435   int i,j,k;
     1484  list l;
     1485  int i,j;
    14361486  j=1;
    14371487  while(j<=n)
    … …  
    14541504    if(size(v)==1)
    14551505    {
    1456       l=spappend(l,number(M[v,v]),1);
     1506      k++;
     1507      e[k]=M[v,v];
     1508      m[k]=1;
    14571509    }
    14581510    else
    14591511    {
    1460       f=factorize(det(submat(M,v,v)-var(1)));
    1461       for(i=size(f[1]);i>=1;i--)
     1512      l=factorize(det(submat(M,v,v)-var(1)));
     1513      for(i=size(l[1]);i>=1;i--)
    14621514      {
    1463         e0=number(jet(f[1][i]/var(1),0));
     1515        e0=number(jet(l[1][i]/var(1),0));
    14641516        if(e0!=0)
    14651517        {
    1466           l=spappend(l,number((e0*var(1)-f[1][i])/e0),f[2][i]);
     1518          k++;
     1519          e[k]=(e0*var(1)-l[1][i])/e0;
     1520          m[k]=l[2][i];
    14671521        }
    14681522      }
    14691523    }
    14701524  }
    1471   return(spnormalize(l[1],l[2]));
     1525  return(spnf(e,m));
    14721526}
    14731527example
    … …  
    16001654  def e,m=#[1..2];
    16011655
    1602   int i;
    1603   for(i=1;i<=ncols(e);i++)
     1656  for(int i=1;i<=ncols(e);i++)
    16041657  {
    16051658    if(deg(e[i]>0))
    … …  
    16101663  }
    16111664
    1612   int j,k,l;
     1665  int j,k,l,n;
    16131666  matrix N0,N1;
    16141667  module K0,K1;
    … …  
    16181671  intvec w;
    16191672
    1620   for(i=ncols(e);i>=1;i--)
     1673  for(i=1;i<=ncols(e);i++)
    16211674  {
    16221675    N0=M-e[i]*freemodule(ncols(M));
    … …  
    16481701        for(j=l;j>=1;j--)
    16491702        {
    1650           U=module(u)+U;
    1651           w=2*j-l-1,w;
     1703          U=U+module(u);
     1704          n++;
     1705          w[n]=2*j-l-1;
    16521706          u=N0*u;
    16531707        }
    … …  
    16551709    }
    16561710  }
    1657   w=w[1..size(w)-1];
     1711
    16581712  return(list(U,w));
    16591713}
    … …  
    16961750      for(i=s[k];i>=2;i--)
    16971751      {
    1698         J[j,j+1]=1;
     1752        J[j+1,j]=1;
    16991753        j++;
    17001754        J[j,j]=e[k];
Note: See TracChangeset for help on using the changeset viewer.