Changeset 0ff6b5 in git


Ignore:
Timestamp:
Aug 25, 2001, 12:52:01 PM (22 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
Children:
a0c62d0d8c47cf16eb73078941175d95d2ab38e7
Parents:
7a7499d2766dbca4b77d47c2cc7f4f8b23e4b382
Message:
*mschulze: new procedures saturate, tmatrix, eigenval, transform


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/gaussman.lib

    r7a7499 r0ff6b5  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: gaussman.lib,v 1.50 2001-08-13 11:40:58 mschulze Exp $";
     2version="$Id: gaussman.lib,v 1.51 2001-08-25 10:52:01 mschulze Exp $";
    33category="Singularities";
    44
    … …  
    108108"USAGE:    gmsring(t,s); poly t, string s;
    109109ASSUME:   basering with characteristic 0 and local degree ordering,
    110           t with isolated singularity at 0
     110          t with isolated citical point 0
    111111RETURN:
    112112@format
    … …  
    148148    else
    149149    {
    150       ERROR("isolated singularity at 0 expected");
     150      ERROR("isolated citical point 0 expected");
    151151    }
    152152  } 
    … …  
    354354///////////////////////////////////////////////////////////////////////////////
    355355
    356 static proc maxintdif(ideal e)
    357 {
    358   int i,j,d;
    359   int d0=0;
    360   for(i=ncols(e);i>=1;i--)
    361   {
    362     for(j=i-1;j>=1;j--)
    363     {
    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);
    376 }
    377 ///////////////////////////////////////////////////////////////////////////////
    378 
    379 proc monodromy(poly t,list #)
    380 "USAGE:    monodromy(t[,opt]); poly t, int opt
    381 ASSUME:   basering with characteristic 0 and local degree ordering,
    382           t with isolated singularity at 0
    383 RETURN:
    384 @format
    385 if opt<=0:
    386 list l:
    387   ideal l[1]: exp(-2*pi*i*l[1]) are the eigenvalues of the monodromy
    388 if opt>=1:
    389 list l: Jordan data jordan(M) of a monodromy matrix exp(-2*pi*i*M)
    390 default: opt=1
    391 @end format
    392 SEE ALSO: mondromy_lib
    393 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; monodromy
    394 EXAMPLE:  example monodromy; shows examples
    395 "
    396 {
    397   int opt=1;
    398   if(size(#)>0)
    399   {
    400     if(typeof(#[1])=="int")
    401     {
    402       opt=#[1];
    403     }
    404   }
    405 
    406   def R=basering;
    407   int n=nvars(R)-1;
    408   def G=gmsring(t,"s");
    409   setring G;
    410 
     356static proc min(ideal e)
     357{
     358  int i;
     359  number m=number(e[1]);
     360  for(i=2;i<=ncols(e);i++)
     361  {
     362    if(number(e[i])<m)
     363    {
     364      m=number(e[i]);
     365    }
     366  }
     367  return(m);
     368}
     369///////////////////////////////////////////////////////////////////////////////
     370
     371static proc max(ideal e)
     372{
     373  int i;
     374  number m=number(e[1]);
     375  for(i=2;i<=ncols(e);i++)
     376  {
     377    if(number(e[i])>m)
     378    {
     379      m=number(e[i]);
     380    }
     381  }
     382  return(m);
     383}
     384///////////////////////////////////////////////////////////////////////////////
     385
     386static proc saturate(int K0)
     387{
    411388  int mu=ncols(gmsbasis);
    412389  ideal r=gmspoly*gmsbasis;
     390  matrix A0[mu][mu],C;
     391  module H0;
     392  module H,H1=freemodule(mu),freemodule(mu);
     393  int k=-1;
    413394  list l;
    414   matrix A0[mu][mu],C;
    415   module H,H1=freemodule(mu),freemodule(mu);
    416   module H0;
    417   int k=-1;
     395
    418396  while(size(reduce(H,std(H0*s)))>0)
    419397  {
     398    dbprint(printlevel-voice+2,"// compute matrix A of t");
    420399    k++;
    421400    dbprint(printlevel-voice+2,"// k="+string(k));
    422     dbprint(printlevel-voice+2,"// compute matrix A of t");
    423     if(opt<=0)
    424     {
    425       l=gmscoeffs(r,k,mu);
    426     }
    427     else
    428     {
    429       l=gmscoeffs(r,k,mu+n);
    430     }
     401    l=gmscoeffs(r,k,mu+K0);
    431402    C,r=l[1..2];
    432403    A0=A0+C;
    … …  
    437408    H=H*s+H1;
    438409  }
     410
    439411  A0=A0-k*s;
    440 
    441412  dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
    442413  H=std(H0);
    443414  int d0=maxdeg1(H);
    444   dbprint(printlevel-voice+2,"// k="+string(d0+1));
     415
     416  dbprint(printlevel-voice+2,"// transform H'' to saturation of H''");
     417  l=division(H,freemodule(mu)*s^k);
     418  H0=l[1];
     419
     420  return(A0,r,H,H0);
     421}
     422///////////////////////////////////////////////////////////////////////////////
     423
     424static proc tmatrix(matrix A0,ideal r,module H,int K,int K0)
     425{
    445426  dbprint(printlevel-voice+2,"// compute matrix A of t");
    446   if(opt<=0)
    447   {
    448     l=gmscoeffs(r,d0+1,d0+1);
    449   }
    450   else
    451   {
    452     l=gmscoeffs(r,d0+1,d0+n+1);
    453   }
     427  int d0=maxdeg1(H);
     428  dbprint(printlevel-voice+2,"// k="+string(K+d0+1));
     429  list l=gmscoeffs(r,K+d0+1,K0+d0+1);
     430  matrix C;
    454431  C,r=l[1..2];
    455432  A0=A0+C;
     433
    456434  dbprint(printlevel-voice+2,"// transform A to saturation of H''");
    457435  l=division(H*s,A0*H+s^2*diff(matrix(H),s));
    458   matrix A=jet(l[1],l[2],0);
    459 
     436  matrix A=jet(l[1],l[2],K);
     437
     438  return(A,A0,r);
     439}
     440///////////////////////////////////////////////////////////////////////////////
     441
     442static proc eigenvals(matrix A0,ideal r,module H,int K0)
     443{
    460444  dbprint(printlevel-voice+2,
    461445    "// compute eigenvalues e with multiplicities m of A");
    462   l=eigenval(jet(A,0));
     446  matrix A;
     447  A,A0,r=tmatrix(A0,r,H,0,K0);
     448  list l=eigenval(A);
    463449  def e,m=l[1..2];
    464450  dbprint(printlevel-voice+2,"// e="+string(e));
    465451  dbprint(printlevel-voice+2,"// m="+string(m));
    466452
    467   if(opt<=0)
    468   {
    469     setring(R);
    470     ideal e=imap(G,e);
    471     return(list(e));
    472   }
    473 
    474   dbprint(printlevel-voice+2,"// compute maximal integer difference d1 of e");
    475   int d1=maxintdif(e);
    476   dbprint(printlevel-voice+2,"// d1="+string(d1));
    477 
    478   module U;
    479   if(d1>0)
    480   {
    481     dbprint(printlevel-voice+2,"// k="+string(d0+d1+1));
    482     dbprint(printlevel-voice+2,"// compute matrix A of t");
    483     l=gmscoeffs(r,d0+d1+1,d0+d1+1);
    484     C,r=l[1..2];
    485     A0=A0+C;
    486     dbprint(printlevel-voice+2,"// transform A to saturation of H''");
    487     l=division(H*s,A0*H+s^2*diff(matrix(H),s));
    488     A=jet(l[1],l[2],d1);
    489 
    490     intvec d;
    491     d[mu]=0;
    492     int i,j;
    493     for(i=ncols(e);i>=1;i--)
    494     {
    495       for(j=i-1;j>=1;j--)
    496       {
    497         k=int(e[j]-e[i]);
    498         if(k>d[j])
    499         {
    500           d[j]=k;
    501         }
    502         if(-k>d[i])
    503         {
    504           d[i]=-k;
    505         }
    506       }
    507     }
    508     for(j,k=ncols(e),mu;j>=1;j--)
    509     {
    510       for(i=m[j];i>=1;i--)
    511       {
    512         d[k]=d[j];
    513         k--;
    514       }
    515     }
    516 
    517     while(d1>0)
    518     {
    519       dbprint(printlevel-voice+2,"// transform basis to reduce d1 by 1");
    520 
     453  return(e,m,A0,r);
     454}
     455///////////////////////////////////////////////////////////////////////////////
     456
     457static proc transform(ideal e,intvec m,matrix A,matrix A0,ideal r,module H,module H0,int K,int K0)
     458{
     459  dbprint(printlevel-voice+2,"// compute minimum e0 and maximum e1 of e");
     460  number e0,e1=min(e),max(e);
     461  dbprint(printlevel-voice+2,"// e0="+string(e0));
     462  dbprint(printlevel-voice+2,"// e1="+string(e1));
     463  int d1=int(e1-e0);
     464  A,A0,r=tmatrix(A0,r,H,K+d1,K0+d1);
     465
     466  if(e1>=e0+1)
     467  {
     468    int i,j,i0,j0,i1,j1;
     469    module U,V;
     470
     471    while(e1>=e0+1)
     472    {
     473      dbprint(printlevel-voice+2,"// transform to separate eigenvalues");
    521474      A0=jet(A,0);
    522475      U=0;
    523476      for(i=ncols(e);i>=1;i--)
    524477      {
    525         U=syz(power(A0-e[i],n+1))+U;
    526       }
    527       A=lift(U,A*U);
    528 
    529       for(i=mu;i>=1;i--)
    530       {
    531         for(j=mu;j>=1;j--)
     478        U=U+syz(power(A0-e[i],m[i]));
     479      }
     480      V=inverse(U);
     481      A=V*A*U;
     482      H0=V*H0;
     483
     484      dbprint(printlevel-voice+2,"// transform to reduce e1 by 1");
     485      for(i0,i=1,1;i0<=ncols(e);i0++)
     486      {
     487        for(i1=1;i1<=m[i0];i1,i=i1+1,i+1)
    532488        {
    533           if(d[i]==0&&d[j]>0)
     489          for(j0,j=1,1;j0<=ncols(e);j0++)
    534490          {
    535             A[i,j]=A[i,j]/s;
     491            for(j1=1;j1<=m[j0];j1,j=j1+1,j+1)
     492            {
     493              if(e[i0]<e0+1&&e[j0]>=e0+1)
     494              {
     495                A[i,j]=A[i,j]/s;
     496              }
     497              if(e[i0]>=e0+1&&e[j0]<e0+1)
     498              {
     499                A[i,j]=A[i,j]*s;
     500              }
     501            }
     502          }
     503        }
     504      }
     505
     506      H0=transpose(H0);
     507      for(i0,i=1,1;i0<=ncols(e);i0++)
     508      {
     509        if(e[i0]>=e0+1)
     510        {
     511          for(i1=1;i1<=m[i0];i1,i=i1+1,i+1)
     512          {
     513            A[i,i]=A[i,i]-1;
     514            H0[i]=H0[i]*s;
    536515          }
    537           else
    538           {
    539             if(d[i]>0&&d[j]==0)
    540             {
    541               A[i,j]=A[i,j]*s;
    542             }
    543           }
     516          e[i0]=e[i0]-1;
    544517        }
    545518      }
    546       for(i=mu;i>=1;i--)
    547       {
    548         if(d[i]>0)
    549         {
    550           A[i,i]=A[i,i]-1;
    551           d[i]=d[i]-1;
    552         }
    553       }
    554 
    555       d1--;
    556       dbprint(printlevel-voice+2,"// d1="+string(d1));
    557     }
    558 
    559     A=jet(A,0);
    560   }
     519      H0=transpose(H0);
     520
     521      e1=e1-1;
     522      dbprint(printlevel-voice+2,"// e1="+string(e1));
     523    }
     524
     525    A=jet(A,K);
     526  }
     527
     528  return(A,A0,r,H0);
     529}
     530///////////////////////////////////////////////////////////////////////////////
     531
     532proc monodromy(poly t,list #)
     533"USAGE:    monodromy(t); poly t
     534ASSUME:   basering with characteristic 0 and local degree ordering,
     535          t with isolated citical point 0
     536RETURN:   list l: Jordan data jordan(M) of a monodromy matrix exp(-2*pi*i*M)
     537SEE ALSO: mondromy_lib
     538KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; monodromy
     539EXAMPLE:  example monodromy; shows examples
     540"
     541{
     542  def R=basering;
     543  int n=nvars(R)-1;
     544  def G=gmsring(t,"s");
     545  setring(G);
     546
     547  int mu=ncols(gmsbasis);
     548  matrix A;
     549  ideal e;
     550  intvec m;
     551
     552  def A0,r,H,H0=saturate(n);
     553  e,m,A0,r=eigenvals(A0,r,H,n);
     554  A,A0,r,H0=transform(e,m,A,A0,r,H,H0,0,0);
    561555
    562556  setring(R);
    563557  matrix A=imap(G,A);
     558
    564559  return(jordan(A));
    565560}
    … …  
    575570"USAGE:    spectrum(t); poly t
    576571ASSUME:   basering with characteristic 0 and local degree ordering,
    577           t with isolated singularity at 0
     572          t with isolated citical point 0
    578573RETURN:
    579574@format
    … …  
    588583"
    589584{
    590   return(sppairs(t,0));
     585  return(spgen(sppairs(t)));
    591586}
    592587example
    … …  
    598593///////////////////////////////////////////////////////////////////////////////
    599594
    600 proc sppairs(poly t,list #)
    601 "USAGE:    sppairs(t[,opt]); poly t, int opt
     595proc sppairs(poly t)
     596"USAGE:    sppairs(t); poly t
    602597ASSUME:   basering with characteristic 0 and local degree ordering,
    603           t with isolated singularity at 0
     598          t with isolated citical point 0
    604599RETURN: list l:
    605600@format
    606601  ideal l[1]: spectral numbers in increasing order
    607 if opt<=0:
    608   intvec l[2]:
    609     int l[2][i]: multiplicity of spectral number l[1][i]
    610 if opt>=1:
    611602  intvec l[2]: weight filtration indices in decreasing order
    612603  intvec l[3]:
    613604    int l[3][i]: multiplicity of spectral pair (l[1][i],l[2][i])
    614 default: opt=1
    615605@end format
    616606SEE ALSO: spectrum_lib
    … …  
    620610"
    621611{
    622   int opt=1;
    623   if(size(#)>0)
    624   {
    625     if(typeof(#[1])=="int")
    626     {
    627       opt=#[1];
    628     }
    629   }
    630 
    631612  def R=basering;
    632613  int n=nvars(R)-1;
    … …  
    635616
    636617  int mu=ncols(gmsbasis);
    637   ideal r=gmspoly*gmsbasis;
    638   list l;
    639   matrix A0[mu][mu],C;
    640   module H0;
    641   module H,H1=freemodule(mu),freemodule(mu);
    642   int k=-1;
    643   while(size(reduce(H,std(H0*s)))>0)
    644   {
    645     k++;
    646     dbprint(printlevel-voice+2,"// k="+string(k));
    647     dbprint(printlevel-voice+2,"// compute matrix A of t");
    648     if(opt<=0)
    649     {
    650       l=gmscoeffs(r,k,mu);
    651     }
    652     else
    653     {
    654       l=gmscoeffs(r,k,mu+n);
    655     }
    656     C,r=l[1..2];
    657     A0=A0+C;
    658 
    659     dbprint(printlevel-voice+2,"// compute saturation of H''");
    660     H0=H;
    661     H1=jet(module(A0*H1+s^2*diff(matrix(H1),s)),k);
    662     H=H*s+H1;
    663   }
    664   A0=A0-k*s;
    665 
    666   dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
    667   H=std(H0);
    668   int d0=maxdeg1(H);
    669   dbprint(printlevel-voice+2,"// transform H'' to saturation of H''");
    670   l=division(H,freemodule(mu)*s^k);
    671   H0=l[1];
    672 
    673   dbprint(printlevel-voice+2,"// k="+string(d0+1));
    674   dbprint(printlevel-voice+2,"// compute matrix A of t");
    675   if(opt<=0)
    676   {
    677     l=gmscoeffs(r,d0+1,d0+1);
    678   }
    679   else
    680   {
    681     l=gmscoeffs(r,d0+1,d0+n+1);
    682   }
    683   C,r=l[1..2];
    684   A0=A0+C;
    685   dbprint(printlevel-voice+2,"// transform A to saturation of H''");
    686   l=division(H*s,A0*H+s^2*diff(matrix(H),s));
    687   matrix A=jet(l[1],l[2],0);
    688 
    689   int i,j;
    690   module U,V;
    691   if(opt<=0)
    692   {
    693     dbprint(printlevel-voice+2,"// transform to Jordan basis");
    694     U=jordanbasis(A)[1];
    695     V=lift(U,freemodule(mu));
    696     A=V*A*U;
    697     dbprint(printlevel-voice+2,"// compute normal form of H''");
    698     H0=std(V*H0);
    699 
    700     dbprint(printlevel-voice+2,"// compute spectral numbers");
    701     ideal a;
    702     for(i=1;i<=mu;i++)
    703     {
    704       j=leadexp(H0[i])[nvars(basering)+1];
    705       a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)-1;
    706     }
    707 
    708     setring(R);
    709     return(spgen(imap(G,a)));
    710   }
    711 
    712   dbprint(printlevel-voice+2,
    713     "// compute eigenvalues e with multiplicities m of A");
    714   l=eigenval(A);
    715   def e,m=l[1..2];
    716   dbprint(printlevel-voice+2,"// e="+string(e));
    717   dbprint(printlevel-voice+2,"// m="+string(m));
    718   dbprint(printlevel-voice+2,"// compute maximal integer difference d1 of e");
    719   int d1=maxintdif(e);
    720   dbprint(printlevel-voice+2,"// d1="+string(d1));
    721 
    722   if(d1>0)
    723   {
    724     dbprint(printlevel-voice+2,"// k="+string(d0+d1+1));
    725     dbprint(printlevel-voice+2,"// compute matrix A of t");
    726     l=gmscoeffs(r,d0+d1+1,d0+d1+1);
    727     C,r=l[1..2];
    728     A0=A0+C;
    729     dbprint(printlevel-voice+2,"// transform A to saturation of H''");
    730     l=division(H*s,A0*H+s^2*diff(matrix(H),s));
    731     A=jet(l[1],l[2],d1);
    732 
    733     intvec d;
    734     d[mu]=0;
    735     for(i=ncols(e);i>=1;i--)
    736     {
    737       for(j=i-1;j>=1;j--)
    738       {
    739         k=int(e[j]-e[i]);
    740         if(k>d[j])
    741         {
    742           d[j]=k;
    743         }
    744         if(-k>d[i])
    745         {
    746           d[i]=-k;
    747         }
    748       }
    749     }
    750     for(j,k=ncols(e),mu;j>=1;j--)
    751     {
    752       for(i=m[j];i>=1;i--)
    753       {
    754         d[k]=d[j];
    755         k--;
    756       }
    757     }
    758 
    759     while(d1>0)
    760     {
    761       dbprint(printlevel-voice+2,"// transform to separate eigenvalues");
    762       A0=jet(A,0);
    763       U=0;
    764       for(i=ncols(e);i>=1;i--)
    765       {
    766         U=syz(power(A0-e[i],n+1))+U;
    767       }
    768       V=lift(U,freemodule(mu));
    769       A=V*A*U;
    770       H0=V*H0;
    771 
    772       dbprint(printlevel-voice+2,"// transform to reduce d1 by 1");
    773       for(i=mu;i>=1;i--)
    774       {
    775         for(j=mu;j>=1;j--)
    776         {
    777           if(d[i]==0&&d[j]>0)
    778           {
    779             A[i,j]=A[i,j]/s;
    780           }
    781           else
    782           {
    783             if(d[i]>0&&d[j]==0)
    784             {
    785               A[i,j]=A[i,j]*s;
    786             }
    787           }
    788         }
    789       }
    790       H0=transpose(H0);
    791       for(i=mu;i>=1;i--)
    792       {
    793         if(d[i]>0)
    794         {
    795           A[i,i]=A[i,i]-1;
    796           d[i]=d[i]-1;
    797           H0[i]=H0[i]*s;
    798         }
    799       }
    800       H0=transpose(H0);
    801 
    802       d1--;
    803       dbprint(printlevel-voice+2,"// d1="+string(d1));
    804     }
    805 
    806     A=jet(A,0);
    807   }
     618  matrix A;
     619  ideal e;
     620  intvec m;
     621
     622  def A0,r,H,H0=saturate(n);
     623  e,m,A0,r=eigenvals(A0,r,H,n);
     624  A,A0,r,H0=transform(e,m,A,A0,r,H,H0,0,0);
    808625
    809626  dbprint(printlevel-voice+2,"// compute weight filtration basis");
    810   intvec w0;
    811   l=jordanbasis(A);
    812   U,w0=l[1..2];
    813   V=lift(U,freemodule(mu));
     627  list l=jordanbasis(A);
     628  def U,v=l[1..2];
     629  module V=inverse(U);
    814630  A0=jet(V*A*U,0);
    815631  vector u;
     632  int i,j,k;
    816633  i=1;
    817634  while(i<ncols(A0))
    … …  
    820637    while(j<ncols(A0)&&A0[i,i]==A0[j,j])
    821638    {
    822       if(w0[i]<w0[j])
    823       {
    824         k=w0[i];
    825         w0[i]=w0[j];
    826         w0[i]=k;
     639      if(v[i]<v[j])
     640      {
     641        k=v[i];
     642        v[i]=v[j];
     643        v[i]=k;
    827644        u=U[i];
    828645        U[i]=U[j];
    … …  
    831648      j++;
    832649    }
    833     if(j==ncols(A0)&&A0[i,i]==A0[j,j]&&w0[i]<w0[j])
    834     {
    835       k=w0[i];
    836       w0[i]=w0[j];
    837       w0[i]=k;
     650    if(j==ncols(A0)&&A0[i,i]==A0[j,j]&&v[i]<v[j])
     651    {
     652      k=v[i];
     653      v[i]=v[j];
     654      v[i]=k;
    838655      u=U[i];
    839656      U[i]=U[j];
    … …  
    844661
    845662  dbprint(printlevel-voice+2,"// transform to weight filtration basis");
    846   V=lift(U,freemodule(mu));
     663  V=inverse(U);
    847664  A=V*A*U;
    848665  dbprint(printlevel-voice+2,"// compute normal form of H''");
    … …  
    856673    j=leadexp(H0[i])[nvars(basering)+1];
    857674    a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)-1;
    858     w[i]=w0[j]+n;
    859   }
     675    w[i]=v[j]+n;
     676  }
     677
    860678  setring(R);
     679
    861680  return(sppgen(imap(G,a),w));
    862681}
    … …  
    872691"USAGE:    vfilt(t[,opt]); poly t, int opt
    873692ASSUME:   basering with characteristic 0 and local degree ordering,
    874           t with isolated singularity at 0
     693          t with isolated citical point 0
    875694RETURN:
    876695@format
    … …  
    11891008  for(i=ncols(m);i>=1;i--)
    11901009  {
    1191     M[i]=lift(V0,coeffs(reduce(m[i]*m,U,g),m)*V0);
     1010    M[i]=division(V0,coeffs(reduce(m[i]*m,U,g),m)*V0)[1];
    11921011  }
    11931012
Note: See TracChangeset for help on using the changeset viewer.