Changeset 86c1f0 in git for Singular/LIB/gaussman.lib


Ignore:
Timestamp:
Aug 9, 2001, 10:39:54 AM (22 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
85140660b24ae4a0f2f58450c04b514020b5d8f1
Parents:
1ac173caaa52482213be39480214765a2e70eb5f
Message:
*mschulze: added sppairs


git-svn-id: file:///usr/local/Singular/svn/trunk@5576 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/gaussman.lib

    r1ac173 r86c1f0  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: gaussman.lib,v 1.48 2001-08-01 13:21:36 mschulze Exp $";
     2version="$Id: gaussman.lib,v 1.49 2001-08-09 08:39:54 mschulze Exp $";
    33category="Singularities";
    44
    … …  
    1212
    1313PROCEDURES:
    14  gmsring(f,s);              Brieskorn lattice in the Gauss-Manin system
     14 gmsring(t,s);              Brieskorn lattice in Gauss-Manin system of t
    1515 gmsnf(p,K[,Kmax]);         Gauss-Manin system normal form
    1616 gmscoeffs(p,K[,Kmax]);     Gauss-Manin system basis representation
    17  monodromy(f[,opt]);        monodromy matrix or spectrum of monodromy of f
    18  vfilt(f[,opt]);            V-filtration on H''/H' or spectrum of f
    19  endfilt(poly f,list V);    endomorphism filtration of filtration V
    20  spectrum(f);               spectrum of f
    21  sppairs(f[,opt]);          spectral pairs or spectrum of f
     17 monodromy(t[,opt]);        Jordan data or eigenvalues of monodromy of t
     18 spectrum(t);               spectrum of t
     19 sppairs(t[,opt]);          spectral pairs or spectrum of t
     20 vfilt(t[,opt]);            V-filtration on H''/H' or spectrum of t
     21 endfilt(t,V);              endomorphism filtration of V-filtration V
    2222 spgen(a);                  generate spectrum defined by a
    2323 sppgen(a,w);               generate spectral pairs defined by a and w
    … …  
    106106
    107107proc gmsring(poly t,string s)
    108 "USAGE:    gmsring(f,s); poly f, string s;
     108"USAGE:    gmsring(t,s); poly t, string s;
    109109ASSUME:   basering with characteristic 0 and local degree ordering,
    110           f with isolated singularity at 0
     110          t with isolated singularity at 0
    111111RETURN:
    112112@format
    113 ring G:
    114   G=C{{s}}*basering,
    115   s is the inverse of the Gauss-Manin connection of f,
    116   G contains:
    117     poly gmspoly: image of f
    118     ideal gmsjacob: image of Jacobian ideal
    119     ideal gmsstd: image of standard basis of Jacobian ideal
    120     matrix gmsmatrix: matrix(gmsjacob)*gmsmatrix=matrix(gmsstd)
    121     ideal gmsbasis: image of monomial vector space basis of Jacobian algebra
    122     int gmsmaxweight: maximal weight of variables of basering
    123   G projects to H''=C{{s}}*gmsbasis
     113ring G: C{{s}}*basering,
     114  poly gmspoly: image of t
     115  ideal gmsjacob: image of Jacobian ideal
     116  ideal gmsstd: image of standard basis of Jacobian ideal
     117  matrix gmsmatrix: matrix(gmsjacob)*gmsmatrix=matrix(gmsstd)
     118  ideal gmsbasis: image of monomial vector space basis of Jacobian algebra
     119  int gmsmaxweight: maximal weight of variables of basering
    124120@end format
    125121NOTE:     do not modify gms variables if you want to use gms procedures
    … …  
    203199{ "EXAMPLE:"; echo=2;
    204200  ring R=0,(x,y),ds;
    205   poly f=x5+x2y2+y5;
    206   def G=gmsring(f,"s");
     201  poly t=x5+x2y2+y5;
     202  def G=gmsring(t,"s");
    207203  setring(G);
    208204  gmspoly;
    … …  
    210206  print(gmsstd);
    211207  print(gmsmatrix);
     208  print(gmsbasis);
    212209  gmsmaxweight;
    213210}
    … …  
    216213proc gmsnf(ideal p,int K,list #)
    217214"USAGE:    gmsnf(p,K[,Kmax]); poly p, int K[, int Kmax];
    218 ASSUME:   basering was constructed by gms, K<=Kmax
     215ASSUME:   basering constructed by gmsring, K<=Kmax
    219216RETURN:
    220217@format
    … …  
    305302{ "EXAMPLE:"; echo=2;
    306303  ring R=0,(x,y),ds;
    307   poly f=x5+x2y2+y5;
    308   def G=gmsring(f,"s");
     304  poly t=x5+x2y2+y5;
     305  def G=gmsring(t,"s");
    309306  setring(G);
    310307  list l0=gmsnf(gmspoly,0);
    … …  
    319316proc gmscoeffs(ideal p,int K,list #)
    320317"USAGE:    gmscoeffs(p,K[,Kmax]); poly p, int K[, int Kmax];
    321 ASSUME:   basering was constructed by gms, K<=Kmax
     318ASSUME:   basering constructed by gmsring, K<=Kmax
    322319RETURN:
    323320@format
    … …  
    345342{ "EXAMPLE:"; echo=2;
    346343  ring R=0,(x,y),ds;
    347   poly f=x5+x2y2+y5;
    348   def G=gmsring(f,"s");
     344  poly t=x5+x2y2+y5;
     345  def G=gmsring(t,"s");
    349346  setring(G);
    350347  list l0=gmscoeffs(gmspoly,0);
    … …  
    359356static proc maxintdif(ideal e)
    360357{
    361   int i,j,id;
    362   int mid=0;
     358  int i,j,d;
     359  int d0=0;
    363360  for(i=ncols(e);i>=1;i--)
    364361  {
    365362    for(j=i-1;j>=1;j--)
    366363    {
    367       id=int(e[i]-e[j]);
    368       if(id<0)
    369       {
    370         id=-id;
    371       }
    372       if(id>mid)
    373       {
    374         mid=id;
    375       }
    376     }
    377   }
    378   return(mid);
     364      d=int(e[i]-e[j]);
     365      if(d<0)
     366      {
     367        d=-d;
     368      }
     369      if(d>d0)
     370      {
     371        d0=d;
     372      }
     373    }
     374  }
     375  return(d0);
    379376}
    380377///////////////////////////////////////////////////////////////////////////////
    … …  
    406403///////////////////////////////////////////////////////////////////////////////
    407404
    408 proc monodromy(poly f,list #)
    409 "USAGE:    monodromy(f[,opt]); poly f, int opt
     405proc monodromy(poly t,list #)
     406"USAGE:    monodromy(t[,opt]); poly t[, int opt]
    410407ASSUME:   basering with characteristic 0 and local degree ordering,
    411           f with isolated singularity at 0
     408          t with isolated singularity at 0
    412409RETURN:
    413410@format
    … …  
    435432  def R=basering;
    436433  int n=nvars(R)-1;
    437   def G=gmsring(f,"s");
     434  def G=gmsring(t,"s");
    438435  setring G;
    439436
    440   int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis);
     437  int mu=ncols(gmsbasis);
    441438  ideal r=gmspoly*gmsbasis;
    442439  list l;
    443440  matrix U=freemodule(mu);
    444441  matrix A0[mu][mu],C;
    445   module H,dH=freemodule(mu),freemodule(mu);
     442  module H,H1=freemodule(mu),freemodule(mu);
    446443  module H0;
    447444  int k=-1;
    … …  
    464461    dbprint(printlevel-voice+2,"// compute saturation of H''");
    465462    H0=H;
    466     dH=jet(module(A0*dH+s^2*diff(matrix(dH),s)),k);
    467     H=H*s+dH;
     463    H1=jet(module(A0*H1+s^2*diff(matrix(H1),s)),k);
     464    H=H*s+H1;
    468465  }
    469466  A0=A0-k*s;
    … …  
    471468  dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
    472469  H=std(H0);
    473   int modH=maxorddif(H)/deg(s);
    474   dbprint(printlevel-voice+2,"// k="+string(modH+1));
     470  int d0=maxorddif(H)/deg(s);
     471  dbprint(printlevel-voice+2,"// k="+string(d0+1));
    475472  dbprint(printlevel-voice+2,"// compute matrix A of t");
    476473  if(opt==0)
    477474  {
    478     l=gmscoeffs(r,modH+1,modH+1);
     475    l=gmscoeffs(r,d0+1,d0+1);
    479476  }
    480477  else
    481478  {
    482     l=gmscoeffs(r,modH+1,modH+1+n);
     479    l=gmscoeffs(r,d0+1,d0+1+n);
    483480  }
    484481  C,r=l[1..2];
    … …  
    502499  }
    503500
    504   dbprint(printlevel-voice+2,"// compute maximal integer difference of e");
    505   int mide=maxintdif(e);
    506   dbprint(printlevel-voice+2,"// mide="+string(mide));
    507 
    508   if(mide>0)
    509   {
    510     dbprint(printlevel-voice+2,"// k="+string(modH+1+mide));
     501  dbprint(printlevel-voice+2,"// compute maximal integer difference d1 of e");
     502  int d1=maxintdif(e);
     503  dbprint(printlevel-voice+2,"// d1="+string(d1));
     504
     505  if(d1>0)
     506  {
     507    dbprint(printlevel-voice+2,"// k="+string(d0+d1+1));
    511508    dbprint(printlevel-voice+2,"// compute matrix A of t");
    512     l=gmscoeffs(r,modH+1+mide,modH+1+mide);
     509    l=gmscoeffs(r,d0+d1+1,d0+d1+1);
    513510    C,r=l[1..2];
    514511    A0=A0+C;
    515512    dbprint(printlevel-voice+2,"// transform A to saturation of H''");
    516513    l=division(H*s,A0*H+s^2*diff(matrix(H),s));
    517     A=jet(l[1],l[2],mide);
    518 
    519     intvec ide;
    520     ide[mu]=0;
     514    A=jet(l[1],l[2],d1);
     515
     516    intvec d;
     517    d[mu]=0;
    521518    int i,j;
    522519    for(i=ncols(e);i>=1;i--)
    … …  
    525522      {
    526523        k=int(e[j]-e[i]);
    527         if(k>ide[j])
    528         {
    529           ide[j]=k;
    530         }
    531         if(-k>ide[i])
    532         {
    533           ide[i]=-k;
     524        if(k>d[j])
     525        {
     526          d[j]=k;
     527        }
     528        if(-k>d[i])
     529        {
     530          d[i]=-k;
    534531        }
    535532      }
    … …  
    539536      for(i=m[j];i>=1;i--)
    540537      {
    541         ide[k]=ide[j];
     538        d[k]=d[j];
    542539        k--;
    543540      }
    544541    }
    545542
    546     while(mide>0)
    547     {
    548       dbprint(printlevel-voice+2,"// transform basis to reduce mide by 1");
     543    while(d1>0)
     544    {
     545      dbprint(printlevel-voice+2,"// transform basis to reduce d1 by 1");
    549546
    550547      A0=jet(A,0);
    … …  
    560557        for(j=mu;j>=1;j--)
    561558        {
    562           if(ide[i]==0&&ide[j]>0)
     559          if(d[i]==0&&d[j]>0)
    563560          {
    564561            A[i,j]=A[i,j]/s;
    … …  
    566563          else
    567564          {
    568             if(ide[i]>0&&ide[j]==0)
     565            if(d[i]>0&&d[j]==0)
    569566            {
    570567              A[i,j]=A[i,j]*s;
    … …  
    575572      for(i=mu;i>=1;i--)
    576573      {
    577         if(ide[i]>0)
     574        if(d[i]>0)
    578575        {
    579576          A[i,i]=A[i,i]-1;
    580           ide[i]=ide[i]-1;
    581         }
    582       }
    583 
    584       mide--;
    585       dbprint(printlevel-voice+2,"// mide="+string(mide));
     577          d[i]=d[i]-1;
     578        }
     579      }
     580
     581      d1--;
     582      dbprint(printlevel-voice+2,"// d1="+string(d1));
    586583    }
    587584
    … …  
    601598///////////////////////////////////////////////////////////////////////////////
    602599
    603 proc vfilt(poly f,list #)
    604 "USAGE:    vfilt(f[,opt]); poly f, int opt
     600proc spectrum(poly t)
     601"USAGE:    spectrum(t); poly t
    605602ASSUME:   basering with characteristic 0 and local degree ordering,
    606           f with isolated singularity at 0
     603          t with isolated singularity at 0
    607604RETURN:
    608605@format
    609 list V: V-filtration of f on H''/H'
    610 if opt=0 or opt=1:
    611   intvec V[1]: numerators of spectral numbers
    612   intvec V[2]: denominators of spectral numbers
    613   intvec V[3]:
    614     int V[3][i]: multiplicity of spectral number V[1][i]/V[2][i]
    615 if opt=1:
     606list Sp: spectrum of t
     607  ideal Sp[1]: spectral numbers in increasing order
     608  intvec Sp[2]:
     609    int Sp[2][i]: multiplicity of spectral number Sp[1][i]
     610@end format
     611SEE ALSO: spectrum_lib
     612KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; spectrum
     613EXAMPLE:  example spnumbers; shows examples
     614"
     615{
     616  return(sppairs(t,0));
     617}
     618example
     619{ "EXAMPLE:"; echo=2;
     620  ring R=0,(x,y),ds;
     621  poly t=x5+x2y2+y5;
     622  spprint(spectrum(t));
     623}
     624///////////////////////////////////////////////////////////////////////////////
     625
     626proc sppairs(poly t,list #)
     627"USAGE:    sppairs(t[,opt]); poly t[, int opt]
     628ASSUME:   basering with characteristic 0 and local degree ordering,
     629          t with isolated singularity at 0
     630RETURN: list l:
     631@format
     632  ideal l[1]: spectral numbers in increasing order
     633if opt<=0:
     634  intvec l[2]:
     635    int l[2][i]: multiplicity of spectral number l[1][i]
     636if opt>=1:
     637  intvec l[2]: weight filtration indices in decreasing order
     638  intvec l[3]:
     639    int l[3][i]: multiplicity of spectral pair (l[1][i],l[2][i])
     640default: opt=1
     641@end format
     642SEE ALSO: spectrum_lib
     643KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
     644          spectrum; spectral pairs
     645EXAMPLE:  example sppairs; shows examples
     646"
     647{
     648  int opt=1;
     649  if(size(#)>0)
     650  {
     651    if(typeof(#[1])=="int")
     652    {
     653      opt=#[1];
     654    }
     655  }
     656
     657  def R=basering;
     658  int n=nvars(R)-1;
     659  def G=gmsring(t,"s");
     660  setring(G);
     661
     662  int mu=ncols(gmsbasis);
     663  ideal r=gmspoly*gmsbasis;
     664  list l;
     665  matrix U=freemodule(mu);
     666  matrix A0[mu][mu],C;
     667  module H0;
     668  module H,H1=freemodule(mu),freemodule(mu);
     669  int k=-1;
     670  while(size(reduce(H,std(H0*s)))>0)
     671  {
     672    k++;
     673    dbprint(printlevel-voice+2,"// k="+string(k));
     674    dbprint(printlevel-voice+2,"// compute matrix A of t");
     675    l=gmscoeffs(r,k,mu+n);
     676    C,r=l[1..2];
     677    A0=A0+C;
     678
     679    dbprint(printlevel-voice+2,"// compute saturation of H''");
     680    H0=H;
     681    H1=jet(module(A0*H1+s^2*diff(matrix(H1),s)),k);
     682    H=H*s+H1;
     683  }
     684  A0=A0-k*s;
     685
     686  dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
     687  H=std(H0);
     688  int d0=maxorddif(H)/deg(s);
     689  dbprint(printlevel-voice+2,"// transform H'' to saturation of H''");
     690  l=division(H,freemodule(mu)*s^k);
     691  H0=l[1];
     692
     693  dbprint(printlevel-voice+2,"// k="+string(d0+1));
     694  dbprint(printlevel-voice+2,"// compute matrix A of t");
     695  l=gmscoeffs(r,d0+1,d0+1+n);
     696  C,r=l[1..2];
     697  A0=A0+C;
     698  dbprint(printlevel-voice+2,"// transform A to saturation of H''");
     699  l=division(H*s,A0*H+s^2*diff(matrix(H),s));
     700  matrix A=jet(l[1],l[2],0);
     701
     702  int i,j;
     703  matrix V;
     704  if(opt<=0)
     705  {
     706    dbprint(printlevel-voice+2,"// transform to Jordan basis");
     707    U=jordanbasis(A)[1];
     708    V=lift(U,freemodule(mu));
     709    A=V*A*U;
     710    dbprint(printlevel-voice+2,"// compute normal form of H''");
     711    H0=std(V*H0);
     712
     713    dbprint(printlevel-voice+2,"// compute spectral numbers");
     714    ideal a;
     715    for(i=1;i<=mu;i++)
     716    {
     717      j=leadexp(H0[i])[nvars(basering)+1];
     718      a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)-1;
     719    }
     720
     721    setring(R);
     722    return(spgen(imap(G,a)));
     723  }
     724
     725  dbprint(printlevel-voice+2,
     726    "// compute eigenvalues e with multiplicities m of A");
     727  l=eigenval(A);
     728  def e,m=l[1..2];
     729  dbprint(printlevel-voice+2,"// e="+string(e));
     730  dbprint(printlevel-voice+2,"// m="+string(m));
     731  dbprint(printlevel-voice+2,"// compute maximal integer difference d1 of e");
     732  int d1=maxintdif(e);
     733  dbprint(printlevel-voice+2,"// d1="+string(d1));
     734
     735  if(d1>0)
     736  {
     737    dbprint(printlevel-voice+2,"// k="+string(d0+d1+1));
     738    dbprint(printlevel-voice+2,"// compute matrix A of t");
     739    l=gmscoeffs(r,d0+d1+1,d0+d1+1);
     740    C,r=l[1..2];
     741    A0=A0+C;
     742    dbprint(printlevel-voice+2,"// transform A to saturation of H''");
     743    l=division(H*s,A0*H+s^2*diff(matrix(H),s));
     744    A=jet(l[1],l[2],d1);
     745
     746    intvec d;
     747    d[mu]=0;
     748    for(i=ncols(e);i>=1;i--)
     749    {
     750      for(j=i-1;j>=1;j--)
     751      {
     752        k=int(e[j]-e[i]);
     753        if(k>d[j])
     754        {
     755          d[j]=k;
     756        }
     757        if(-k>d[i])
     758        {
     759          d[i]=-k;
     760        }
     761      }
     762    }
     763    for(j,k=ncols(e),mu;j>=1;j--)
     764    {
     765      for(i=m[j];i>=1;i--)
     766      {
     767        d[k]=d[j];
     768        k--;
     769      }
     770    }
     771
     772    while(d1>0)
     773    {
     774      dbprint(printlevel-voice+2,"// transform to separate eigenvalues");
     775      A0=jet(A,0);
     776      U=0;
     777      for(i=ncols(e);i>=1;i--)
     778      {
     779        U=syz(power(A0-e[i],n+1))+U;
     780      }
     781      V=lift(U,freemodule(mu));
     782      A=V*A*U;
     783      H0=V*H0;
     784
     785      dbprint(printlevel-voice+2,"// transform to reduce d1 by 1");
     786      for(i=mu;i>=1;i--)
     787      {
     788        for(j=mu;j>=1;j--)
     789        {
     790          if(d[i]==0&&d[j]>0)
     791          {
     792            A[i,j]=A[i,j]/s;
     793          }
     794          else
     795          {
     796            if(d[i]>0&&d[j]==0)
     797            {
     798              A[i,j]=A[i,j]*s;
     799            }
     800          }
     801        }
     802      }
     803      H0=transpose(H0);
     804      for(i=mu;i>=1;i--)
     805      {
     806        if(d[i]>0)
     807        {
     808          A[i,i]=A[i,i]-1;
     809          d[i]=d[i]-1;
     810          H0[i]=H0[i]*s;
     811        }
     812      }
     813      H0=transpose(H0);
     814
     815      d1--;
     816      dbprint(printlevel-voice+2,"// d1="+string(d1));
     817    }
     818
     819    A=jet(A,0);
     820  }
     821
     822  dbprint(printlevel-voice+2,"// transform to Jordan basis");
     823  intvec w0;
     824  l=jordanbasis(A);
     825  U,w0=l[1..2];
     826  V=lift(U,freemodule(mu));
     827  A0=V*A*U;
     828
     829  dbprint(printlevel-voice+2,"// compute weight filtration basis");
     830  vector u;
     831  i=1;
     832  while(i<ncols(A))
     833  {
     834    j=i+1;
     835    while(j<ncols(A)&&A[i,i]==A[j,j])
     836    {
     837      if(w0[i]<w0[j])
     838      {
     839        k=w0[i];
     840        w0[i]=w0[j];
     841        w0[i]=k;
     842        u=U[i];
     843        U[i]=U[j];
     844        U[j]=u;
     845      }
     846      j++;
     847    }
     848    if(j==ncols(A)&&A[i,i]==A[j,j]&&w0[i]<w0[j])
     849    {
     850      k=w0[i];
     851      w0[i]=w0[j];
     852      w0[i]=k;
     853      u=U[i];
     854      U[i]=U[j];
     855      U[j]=u;
     856    }
     857    i=j;
     858  }
     859
     860  dbprint(printlevel-voice+2,"// transform to weight filtration basis");
     861  V=lift(U,freemodule(mu));
     862  A=V*A*U;
     863  dbprint(printlevel-voice+2,"// compute normal form of H''");
     864  H0=std(V*H0);
     865
     866  dbprint(printlevel-voice+2,"// compute spectral pairs");
     867  ideal a;
     868  intvec w;
     869  for(i=1;i<=mu;i++)
     870  {
     871    j=leadexp(H0[i])[nvars(basering)+1];
     872    a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)-1;
     873    w[i]=w0[j]+n;
     874  }
     875  setring(R);
     876  return(sppgen(imap(G,a),w));
     877}
     878example
     879{ "EXAMPLE:"; echo=2;
     880  ring R=0,(x,y),ds;
     881  poly t=x5+x2y2+y5;
     882  spprint(sppairs(t));
     883}
     884///////////////////////////////////////////////////////////////////////////////
     885
     886proc vfilt(poly t,list #)
     887"USAGE:    vfilt(t[,opt]); poly t, int opt
     888ASSUME:   basering with characteristic 0 and local degree ordering,
     889          t with isolated singularity at 0
     890RETURN:
     891@format
     892list V: V-filtration of t on H''/H'
     893  intvec V[1]: spectral numbers in increasing order
     894  intvec V[2]:
     895    int V[2][i]: multiplicity of spectral number V[1][i]/V[2][i]
     896if opt>=1:
    616897  list V[4]:
    617     module V[4][i]: vector space basis of V[1][i]/V[2][i]-th graded part
     898    module V[3][i]: vector space basis of V[1][i]/V[2][i]-th graded part
    618899                    in terms of V[5]
    619   ideal V[5]: monomial vector space basis
     900  ideal V[4]: monomial vector space basis of H''/H'
     901  ideal V[5]: standard basis of Jacobian ideal
    620902default: opt=1
    621903@end format
    … …  
    638920  def R=basering;
    639921  int n=nvars(R)-1;
    640   def G=gmsring(f,"s");
     922  def G=gmsring(t,"s");
    641923  setring G;
    642924
    643   int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis);
     925  int mu=ncols(gmsbasis);
    644926  ideal r=gmspoly*gmsbasis;
    645927  list l;
    646928  matrix U=freemodule(mu);
    647929  matrix A[mu][mu],C;
    648   module H,dH=freemodule(mu),freemodule(mu);
     930  module H,H1=freemodule(mu),freemodule(mu);
    649931  module H0;
    650932  int k=-1;
    … …  
    662944    dbprint(printlevel-voice+2,"// compute saturation of H''");
    663945    H0=H;
    664     dH=jet(module(A*dH+s^2*diff(matrix(dH),s)),k);
    665     H=H*s+dH;
     946    H1=jet(module(A*H1+s^2*diff(matrix(H1),s)),k);
     947    H=H*s+H1;
    666948  }
    667949  A=A-k*s;
    … …  
    669951  dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
    670952  H=std(H0);
    671   int modH=maxorddif(H)/deg(s);
    672   dbprint(printlevel-voice+2,"// k="+string(N+modH));
     953  int d0=maxorddif(H)/deg(s);
     954  dbprint(printlevel-voice+2,"// k="+string(d0+N));
    673955  dbprint(printlevel-voice+2,"// compute matrix A of t");
    674   l=gmscoeffs(r,N+modH,N+modH);
     956  l=gmscoeffs(r,d0+N,d0+N);
    675957  C,r=l[1..2];
    676958  A=A+C;
    … …  
    7751057  V[ncols(eM)+1]=interred(V1);
    7761058  intvec d;
    777   if(opt==0)
     1059  if(opt<=0)
    7781060  {
    7791061    for(i=ncols(eM);number(eM[i])-1>number(n-1)/2;i--)
    … …  
    8761158    list v=imap(G,v);
    8771159    ideal m=imap(G,gmsbasis);
    878     return(list(a,d,v,m));
     1160    ideal g=imap(G,gmsstd);
     1161    attrib(g,"isSB",1);
     1162    return(list(a,d,v,m,g));
    8791163  }
    8801164}
    … …  
    8821166{ "EXAMPLE:"; echo=2;
    8831167  ring R=0,(x,y),ds;
    884   poly f=x5+x2y2+y5;
    885   vfilt(f);
    886 }
    887 ///////////////////////////////////////////////////////////////////////////////
    888 
    889 proc endfilt(poly f,list V)
    890 "USAGE:   endfilt(f,V); poly f, list V
    891 ASSUME:  basering with characteristic 0 and local degree ordering,
    892          f with isolated singularity at 0
     1168  poly t=x5+x2y2+y5;
     1169  vfilt(t);
     1170}
     1171///////////////////////////////////////////////////////////////////////////////
     1172
     1173proc endfilt(list V)
     1174"USAGE:   endfilt(V); list V
     1175ASSUME:  V computed by vfilt
    8931176RETURN:
    8941177@format
    895 list V1: endomorphim filtration of V on the Jacobian algebra of f
     1178list V1: endomorphim filtration of V on the Jacobian algebra
    8961179  ideal V1[1]: spectral numbers in increasing order
    8971180  intvec V1[2]:
    … …  
    9081191"
    9091192{
    910   if(charstr(basering)!="0")
    911   {
    912     ERROR("characteristic 0 expected");
    913   }
    914   int n=nvars(basering)-1;
    915   for(int i=n+1;i>=1;i--)
    916   {
    917     if(var(i)>1)
    918     {
    919       ERROR("local ordering expected");
    920     }
    921   }
    922   ideal sJ=std(jacob(f));
    923   if(vdim(sJ)<=0)
    924   {
    925     if(vdim(sJ)==0)
    926     {
    927       ERROR("singularity at 0 expected");
    928     }
    929     else
    930     {
    931       ERROR("isolated singularity at 0 expected");
    932     }
    933   }
    934 
    935   def a,d,v,m=V[1..4];
     1193  def a,d,v,m,g=V[1..5];
    9361194  int mu=ncols(m);
    9371195
    9381196  module V0=v[1];
    939   for(i=2;i<=size(v);i++)
     1197  for(int i=2;i<=size(v);i++)
    9401198  {
    9411199    V0=V0,v[i];
    … …  
    9471205  for(i=ncols(m);i>=1;i--)
    9481206  {
    949     M[i]=lift(V0,coeffs(reduce(m[i]*m,U,sJ),m)*V0);
     1207    M[i]=lift(V0,coeffs(reduce(m[i]*m,U,g),m)*V0);
    9501208  }
    9511209
    … …  
    10621320{ "EXAMPLE:"; echo=2;
    10631321  ring R=0,(x,y),ds;
    1064   poly f=x5+x2y2+y5;
    1065   endfilt(f,vfilt(f));
    1066 }
    1067 ///////////////////////////////////////////////////////////////////////////////
    1068 
    1069 proc spectrum(poly f)
    1070 "USAGE:    spectrum(f); poly f
    1071 ASSUME:   basering with characteristic 0 and local degree ordering,
    1072           f with isolated singularity at 0
    1073 RETURN:
    1074 @format
    1075 list Sp: spectrum of f
    1076   ideal Sp[1]: spectral numbers in increasing order
    1077   intvec Sp[2]:
    1078     int Sp[2][i]: multiplicity of spectral number Sp[1][i]
    1079 @end format
    1080 SEE ALSO: spectrum_lib
    1081 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; spectrum
    1082 EXAMPLE:  example spnumbers; shows examples
    1083 "
    1084 {
    1085   return(sppairs(f,0));
    1086 }
    1087 example
    1088 { "EXAMPLE:"; echo=2;
    1089   ring R=0,(x,y),ds;
    1090   poly f=x5+x2y2+y5;
    1091   spprint(spectrum(f));
    1092 }
    1093 ///////////////////////////////////////////////////////////////////////////////
    1094 
    1095 proc sppairs(poly f,list #)
    1096 "USAGE:    sppairs(f[,opt]); poly f, int opt
    1097 ASSUME:   basering with characteristic 0 and local degree ordering,
    1098           f with isolated singularity at 0
    1099 RETURN:
    1100 @format
    1101 if opt=0
    1102 list Sp: spectrum of f
    1103   ideal Sp[1]: spectral numbers in increasing order
    1104   intvec Sp[2]: corresponding multiplicities of spectral numbers
    1105 if opt=1:
    1106 list Spp: spectral pairs of f
    1107   ideal Spp[1]: spectral numbers in increasing order
    1108   intvec Spp[2]: corresponding weight filtration indices in increasing order
    1109   intvec Spp[3]: corresponding multiplicities of spectral pairs
    1110 default: opt=1
    1111 @end format
    1112 SEE ALSO: spectrum_lib
    1113 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
    1114           spectrum; spectral pairs
    1115 EXAMPLE:  example sppairs; shows examples
    1116 "
    1117 {
    1118   int opt=1;
    1119   if(size(#)>0)
    1120   {
    1121     if(typeof(#[1])=="int")
    1122     {
    1123       opt=#[1];
    1124     }
    1125   }
    1126 
    1127   def R=basering;
    1128   int n=nvars(R)-1;
    1129   def G=gmsring(f,"s");
    1130   setring(G);
    1131 
    1132   int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis);
    1133   ideal r=gmspoly*gmsbasis;
    1134   list l;
    1135   matrix U=freemodule(mu);
    1136   matrix A0[mu][mu],C;
    1137   module H0;
    1138   module H,dH=freemodule(mu),freemodule(mu);
    1139   int k=-1;
    1140   while(size(reduce(H,std(H0*s)))>0)
    1141   {
    1142     k++;
    1143     dbprint(printlevel-voice+2,"// k="+string(k));
    1144     dbprint(printlevel-voice+2,"// compute matrix A of t");
    1145     l=gmscoeffs(r,k,mu+n);
    1146     C,r=l[1..2];
    1147     A0=A0+C;
    1148 
    1149     dbprint(printlevel-voice+2,"// compute saturation of H''");
    1150     H0=H;
    1151     dH=jet(module(A0*dH+s^2*diff(matrix(dH),s)),k);
    1152     H=H*s+dH;
    1153   }
    1154   A0=A0-k*s;
    1155 
    1156   dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
    1157   H=std(H0);
    1158   int modH=maxorddif(H)/deg(s);
    1159   dbprint(printlevel-voice+2,"// transform H'' to saturation of H''");
    1160   l=division(H,freemodule(mu)*s^k);
    1161   H0=l[1];
    1162 
    1163   dbprint(printlevel-voice+2,"// k="+string(modH+1));
    1164   dbprint(printlevel-voice+2,"// compute matrix A of t");
    1165   l=gmscoeffs(r,modH+1,modH+1+n);
    1166   C,r=l[1..2];
    1167   A0=A0+C;
    1168   dbprint(printlevel-voice+2,"// transform A to saturation of H''");
    1169   l=division(H*s,A0*H+s^2*diff(matrix(H),s));
    1170   matrix A=jet(l[1],l[2],0);
    1171 
    1172   int i,j;
    1173   matrix V;
    1174   if(!opt)
    1175   {
    1176     dbprint(printlevel-voice+2,"// transform to Jordan basis");
    1177     U=jordanbasis(A)[1];
    1178     V=lift(U,freemodule(mu));
    1179     A=V*A*U;
    1180     dbprint(printlevel-voice+2,"// compute normal form of H''");
    1181     H0=std(V*H0);
    1182 
    1183     dbprint(printlevel-voice+2,"// compute spectral numbers");
    1184     ideal a;
    1185     for(i=1;i<=mu;i++)
    1186     {
    1187       j=leadexp(H0[i])[nvars(basering)+1];
    1188       a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)-1;
    1189     }
    1190 
    1191     setring(R);
    1192     return(spgen(imap(G,a)));
    1193   }
    1194 
    1195   dbprint(printlevel-voice+2,
    1196     "// compute eigenvalues e with multiplicities m of A");
    1197   l=eigenval(A);
    1198   def e,m=l[1..2];
    1199   dbprint(printlevel-voice+2,"// e="+string(e));
    1200   dbprint(printlevel-voice+2,"// m="+string(m));
    1201   dbprint(printlevel-voice+2,"// compute maximal integer difference of e");
    1202   int mide=maxintdif(e);
    1203   dbprint(printlevel-voice+2,"// mide="+string(mide));
    1204 
    1205   if(mide>0)
    1206   {
    1207     dbprint(printlevel-voice+2,"// k="+string(modH+1+mide));
    1208     dbprint(printlevel-voice+2,"// compute matrix A of t");
    1209     l=gmscoeffs(r,modH+1+mide,modH+1+mide);
    1210     C,r=l[1..2];
    1211     A0=A0+C;
    1212     dbprint(printlevel-voice+2,"// transform A to saturation of H''");
    1213     l=division(H*s,A0*H+s^2*diff(matrix(H),s));
    1214     A=jet(l[1],l[2],mide);
    1215 
    1216     intvec ide;
    1217     ide[mu]=0;
    1218     for(i=ncols(e);i>=1;i--)
    1219     {
    1220       for(j=i-1;j>=1;j--)
    1221       {
    1222         k=int(e[j]-e[i]);
    1223         if(k>ide[j])
    1224         {
    1225           ide[j]=k;
    1226         }
    1227         if(-k>ide[i])
    1228         {
    1229           ide[i]=-k;
    1230         }
    1231       }
    1232     }
    1233     for(j,k=ncols(e),mu;j>=1;j--)
    1234     {
    1235       for(i=m[j];i>=1;i--)
    1236       {
    1237         ide[k]=ide[j];
    1238         k--;
    1239       }
    1240     }
    1241 
    1242     while(mide>0)
    1243     {
    1244       dbprint(printlevel-voice+2,"// transform to separate eigenvalues");
    1245       A0=jet(A,0);
    1246       U=0;
    1247       for(i=ncols(e);i>=1;i--)
    1248       {
    1249         U=syz(power(A0-e[i],n+1))+U;
    1250       }
    1251       V=lift(U,freemodule(mu));
    1252       A=V*A*U;
    1253       H0=V*H0;
    1254 
    1255       dbprint(printlevel-voice+2,"// transform to reduce mide by 1");
    1256       for(i=mu;i>=1;i--)
    1257       {
    1258         for(j=mu;j>=1;j--)
    1259         {
    1260           if(ide[i]==0&&ide[j]>0)
    1261           {
    1262             A[i,j]=A[i,j]/s;
    1263           }
    1264           else
    1265           {
    1266             if(ide[i]>0&&ide[j]==0)
    1267             {
    1268               A[i,j]=A[i,j]*s;
    1269             }
    1270           }
    1271         }
    1272       }
    1273       H0=transpose(H0);
    1274       for(i=mu;i>=1;i--)
    1275       {
    1276         if(ide[i]>0)
    1277         {
    1278           A[i,i]=A[i,i]-1;
    1279           ide[i]=ide[i]-1;
    1280           H0[i]=H0[i]*s;
    1281         }
    1282       }
    1283       H0=transpose(H0);
    1284 
    1285       mide--;
    1286       dbprint(printlevel-voice+2,"// mide="+string(mide));
    1287     }
    1288 
    1289     A=jet(A,0);
    1290   }
    1291 
    1292   dbprint(printlevel-voice+2,"// transform to Jordan basis");
    1293   intvec w0;
    1294   l=jordanbasis(A);
    1295   U,w0=l[1..2];
    1296   V=lift(U,freemodule(mu));
    1297   A0=V*A*U;
    1298 
    1299   dbprint(printlevel-voice+2,"// compute weight filtration basis");
    1300   vector u;
    1301   i=1;
    1302   while(i<ncols(A))
    1303   {
    1304     j=i+1;
    1305     while(j<ncols(A)&&A[i,i]==A[j,j])
    1306     {
    1307       if(w0[i]<w0[j])
    1308       {
    1309         k=w0[i];
    1310         w0[i]=w0[j];
    1311         w0[i]=k;
    1312         u=U[i];
    1313         U[i]=U[j];
    1314         U[j]=u;
    1315       }
    1316       j++;
    1317     }
    1318     if(j==ncols(A)&&A[i,i]==A[j,j]&&w0[i]<w0[j])
    1319     {
    1320       k=w0[i];
    1321       w0[i]=w0[j];
    1322       w0[i]=k;
    1323       u=U[i];
    1324       U[i]=U[j];
    1325       U[j]=u;
    1326     }
    1327     i=j;
    1328   }
    1329 
    1330   dbprint(printlevel-voice+2,"// transform to weight filtration basis");
    1331   V=lift(U,freemodule(mu));
    1332   A=V*A*U;
    1333   dbprint(printlevel-voice+2,"// compute normal form of H''");
    1334   H0=std(V*H0);
    1335 
    1336   dbprint(printlevel-voice+2,"// compute spectral pairs");
    1337   ideal a;
    1338   intvec w;
    1339   for(i=1;i<=mu;i++)
    1340   {
    1341     j=leadexp(H0[i])[nvars(basering)+1];
    1342     a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)-1;
    1343     w[i]=w0[j]+n;
    1344   }
    1345   setring(R);
    1346   return(sppgen(imap(G,a),w));
    1347 }
    1348 example
    1349 { "EXAMPLE:"; echo=2;
    1350   ring R=0,(x,y),ds;
    1351   poly f=x5+x2y2+y5;
    1352   spprint(sppairs(f));
     1322  poly t=x5+x2y2+y5;
     1323  endfilt(vfilt(t));
    13531324}
    13541325///////////////////////////////////////////////////////////////////////////////
    … …  
    13601331list Sp: numbers in a with multiplicities
    13611332  ideal Sp[1]: numbers in a in increasing order
    1362   intvec Sp[2]: corresponding multiplicities
     1333  intvec Sp[2]:
     1334    int Sp[2][i]: multiplicity of number Sp[1][i] in a
    13631335@end format
    13641336EXAMPLE: example spgen; shows examples
    … …  
    14071379
    14081380proc sppgen(ideal a,intvec w)
    1409 "USAGE:   sppgen(a); ideal a
     1381"USAGE:   sppgen(a,w); ideal a, intvec w
    14101382RETURN:
    14111383@format
    14121384list Spp: pairs in a and w with multiplicities
    1413   ideal Spp[1]: numbers in a
    1414   intvec Spp[2]: corresponding integers in w
    1415   intvec Spp[3]: corresponding multiplicities
     1385  ideal Spp[1]: numbers in a in increasing order
     1386  intvec Spp[2]: integers in w in decreasing order
     1387  intvec Spp[3]:
     1388    int Spp[3][i]: multiplicity of pair (Spp[1][i],Spp[2][i]) in a,w
    14161389@end format
    14171390EXAMPLE: example sppgen; shows examples
Note: See TracChangeset for help on using the changeset viewer.