source: git/Singular/LIB/gmssing.lib @ 590cbe

spielwiese
Last change on this file since 590cbe was 590cbe, checked in by Hans Schönemann <hannes@…>, 14 years ago
*hannes: fix export git-svn-id: file:///usr/local/Singular/svn/trunk@12016 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.6 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: gmssing.lib,v 1.14 2009-07-28 08:13:38 Singular Exp $";
3category="Singularities";
4
5info="
6LIBRARY:  gaussman.lib  Gauss-Manin System of Isolated Singularities
7
8AUTHOR:   Mathias Schulze, email: mschulze@mathematik.uni-kl.de
9
10OVERVIEW: A library to compute invariants related to the the Gauss-Manin system
11          of an isolated hypersurface singularity
12
13PROCEDURES:
14 gmsring(t,s);              Gauss-Manin system of t with variable s
15 gmsnf(p,K);                Gauss-Manin normal form of p
16 gmscoeffs(p,K);            Gauss-Manin basis representation of p
17 bernstein(t);              Bernstein-Sato polynomial of t
18 monodromy(t);              Jordan data of complex monodromy of t
19 spectrum(t);               singularity spectrum of t
20 sppairs(t);                spectral pairs of t
21 vfilt(t);                  V-filtration of t on Brieskorn lattice
22 vwfilt(t);                 weighted V-filtration of t on Brieskorn lattice
23 tmatrix(t);                matrix of t w.r.t. good basis of Brieskorn lattice
24 endvfilt(V);               endomorphism V-filtration on Jacobian algebra
25 sppnf(a,w[,m]);            spectral pairs normal form of (a,w[,m])
26 sppprint(spp);             print spectral pairs spp
27 spadd(sp1,sp2);            sum of spectra sp1 and sp2
28 spsub(sp1,sp2);            difference of spectra sp1 and sp2
29 spmul(sp0,k);              linear combination of spectra sp
30 spissemicont(sp[,opt]);    semicontinuity test of spectrum sp
31 spsemicont(sp0,sp[,opt]);  semicontinuous combinations of spectra sp0 in sp
32 spmilnor(sp);              Milnor number of spectrum sp
33 spgeomgenus(sp);           geometrical genus of spectrum sp
34 spgamma(sp);               gamma invariant of spectrum sp
35
36SEE ALSO: mondromy_lib, spectrum_lib, gmspoly_lib
37
38KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
39          mixed Hodge structure; V-filtration; weight filtration
40          Bernstein-Sato polynomial; monodromy; spectrum; spectral pairs;
41          good basis
42";
43
44LIB "linalg.lib";
45
46///////////////////////////////////////////////////////////////////////////////
47
48static proc stdtrans(ideal I)
49{
50  def @R=basering;
51
52  string os=ordstr(@R);
53  int j=find(os,",C");
54  if(j==0)
55  {
56    j=find(os,"C,");
57  }
58  if(j==0)
59  {
60    j=find(os,",c");
61  }
62  if(j==0)
63  {
64    j=find(os,"c,");
65  }
66  if(j>0)
67  {
68    os[j..j+1]="  ";
69  }
70
71  execute("ring @S="+charstr(@R)+",(gmspoly,"+varstr(@R)+"),(c,dp,"+os+");");
72
73  ideal I=homog(imap(@R,I),gmspoly);
74  module M=transpose(transpose(I)+freemodule(ncols(I)));
75  M=std(M);
76
77  setring(@R);
78  execute("map h=@S,1,"+varstr(@R)+";");
79  module M=h(M);
80
81  for(int i=ncols(M);i>=1;i--)
82  {
83    for(j=ncols(M);j>=1;j--)
84    {
85      if(M[i][1]==0)
86      {
87        M[i]=0;
88      }
89      if(i!=j&&M[j][1]!=0)
90      {
91        if(lead(M[i][1])/lead(M[j][1])!=0)
92        {
93          M[i]=0;
94        }
95      }
96    }
97  }
98
99  M=transpose(simplify(M,2));
100  I=M[1];
101  attrib(I,"isSB",1);
102  M=M[2..ncols(M)];
103  module U=transpose(M);
104
105  return(list(I,U));
106}
107///////////////////////////////////////////////////////////////////////////////
108
109proc gmsring(poly t,string s)
110"USAGE:    gmsring(t,s); poly t, string s
111ASSUME:   characteristic 0; local degree ordering;
112          isolated critical point 0 of t
113RETURN:
114@format
115ring G;  Gauss-Manin system of t with variable s
116  poly gmspoly=t;
117  ideal gmsjacob;  Jacobian ideal of t
118  ideal gmsstd;  standard basis of Jacobian ideal
119  matrix gmsmatrix;  matrix(gmsjacob)*gmsmatrix==matrix(gmsstd)
120  ideal gmsbasis;  monomial vector space basis of Jacobian algebra
121  int gmsmaxdeg;  maximal weight of variables
122@end format
123NOTE:     gmsbasis is a C[[s]]-basis of H'' and [t,s]=s^2
124KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice
125EXAMPLE:  example gmsring; shows examples
126"
127{
128  def @R=basering;
129  if(charstr(@R)!="0")
130  {
131    ERROR("characteristic 0 expected");
132  }
133  for(int i=nvars(@R);i>=1;i--)
134  {
135    if(var(i)>1)
136    {
137      ERROR("local ordering expected");
138    }
139  }
140
141  ideal dt=jacob(t);
142  list l=stdtrans(dt);
143  ideal g=l[1];
144  if(vdim(g)<=0)
145  {
146    if(vdim(g)==0)
147    {
148      ERROR("singularity at 0 expected");
149    }
150    else
151    {
152      ERROR("isolated critical point 0 expected");
153    }
154  }
155  matrix B=l[2];
156  ideal m=kbase(g);
157
158  int gmsmaxdeg;
159  for(i=nvars(@R);i>=1;i--)
160  {
161    if(deg(var(i))>gmsmaxdeg)
162    {
163      gmsmaxdeg=deg(var(i));
164    }
165  }
166
167  string os=ordstr(@R);
168  int j=find(os,",C");
169  if(j==0)
170  {
171    j=find(os,"C,");
172  }
173  if(j==0)
174  {
175    j=find(os,",c");
176  }
177  if(j==0)
178  {
179    j=find(os,"c,");
180  }
181  if(j>0)
182  {
183    os[j..j+1]="  ";
184  }
185
186  execute("ring G="+string(charstr(@R))+",("+s+","+varstr(@R)+"),(ws("+
187    string(deg(highcorner(g))+2*gmsmaxdeg)+"),"+os+",c);");
188
189  poly gmspoly=imap(@R,t);
190  ideal gmsjacob=imap(@R,dt);
191  ideal gmsstd=imap(@R,g);
192  matrix gmsmatrix=imap(@R,B);
193  ideal gmsbasis=imap(@R,m);
194
195  attrib(gmsstd,"isSB",1);
196  export gmspoly,gmsjacob,gmsstd,gmsmatrix,gmsbasis,gmsmaxdeg;
197  if (system("with","Namespaces")) { exportto(Top, gmsmaxdeg);}
198  else { export(gmsmaxdeg);}
199
200  return(G);
201}
202example
203{ "EXAMPLE:"; echo=2;
204  ring @R=0,(x,y),ds;
205  poly t=x5+x2y2+y5;
206  def G=gmsring(t,"s");
207  setring(G);
208  gmspoly;
209  print(gmsjacob);
210  print(gmsstd);
211  print(gmsmatrix);
212  print(gmsbasis);
213  gmsmaxdeg;
214}
215///////////////////////////////////////////////////////////////////////////////
216
217proc gmsnf(ideal p,int K)
218"USAGE:    gmsnf(p,K); poly p, int K
219ASSUME:   basering returned by gmsring
220RETURN:
221  list nf;
222  ideal nf[1];  projection of p to <gmsbasis>C{{s}} mod s^(K+1)
223  ideal nf[2];  p==nf[1]+nf[2]
224NOTE:     computation can be continued by setting p=nf[2]
225KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice
226EXAMPLE:  example gmsnf; shows examples
227"
228{
229  if(system("with","gms"))
230  {
231    return(system("gmsnf",p,gmsstd,gmsmatrix,(K+1)*deg(var(1))-2*gmsmaxdeg,K));
232  }
233
234  intvec v=1;
235  v[nvars(basering)]=0;
236
237  int k;
238  ideal r,q;
239  r[ncols(p)]=0;
240  q[ncols(p)]=0;
241
242  poly s;
243  int i,j;
244  for(k=ncols(p);k>=1;k--)
245  {
246    while(p[k]!=0&&deg(lead(p[k]),v)<=K)
247    {
248      i=1;
249      s=lead(p[k])/lead(gmsstd[i]);
250      while(i<ncols(gmsstd)&&s==0)
251      {
252        i++;
253        s=lead(p[k])/lead(gmsstd[i]);
254      }
255      if(s!=0)
256      {
257        p[k]=p[k]-s*gmsstd[i];
258        for(j=1;j<=nrows(gmsmatrix);j++)
259        {
260          p[k]=p[k]+diff(s*gmsmatrix[j,i],var(j+1))*var(1);
261        }
262      }
263      else
264      {
265        r[k]=r[k]+lead(p[k]);
266        p[k]=p[k]-lead(p[k]);
267      }
268      while(deg(lead(p[k]))>(K+1)*deg(var(1))-2*gmsmaxdeg&&
269        deg(lead(p[k]),v)<=K)
270      {
271        q[k]=q[k]+lead(p[k]);
272        p[k]=p[k]-lead(p[k]);
273      }
274    }
275    q[k]=q[k]+p[k];
276  }
277
278  return(list(r,q));
279}
280example
281{ "EXAMPLE:"; echo=2;
282  ring R=0,(x,y),ds;
283  poly t=x5+x2y2+y5;
284  def G=gmsring(t,"s");
285  setring(G);
286  list l0=gmsnf(gmspoly,0);
287  print(l0[1]);
288  list l1=gmsnf(gmspoly,1);
289  print(l1[1]);
290  list l=gmsnf(l0[2],1);
291  print(l[1]);
292}
293///////////////////////////////////////////////////////////////////////////////
294
295proc gmscoeffs(ideal p,int K)
296"USAGE:    gmscoeffs(p,K); poly p, int K
297ASSUME:   basering constructed by gmsring
298RETURN:
299@format
300list l;
301  matrix l[1];  C{{s}}-basis representation of p mod s^(K+1)
302  ideal l[2];  p==matrix(gmsbasis)*l[1]+l[2]
303@end format
304NOTE:     computation can be continued by setting p=l[2]
305KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice
306EXAMPLE:  example gmscoeffs; shows examples
307"
308{
309  list l=gmsnf(p,K);
310  ideal r,q=l[1..2];
311  poly v=1;
312  for(int i=2;i<=nvars(basering);i++)
313  {
314    v=v*var(i);
315  }
316  matrix C=coeffs(r,gmsbasis,v);
317  return(list(C,q));
318}
319example
320{ "EXAMPLE:"; echo=2;
321  ring R=0,(x,y),ds;
322  poly t=x5+x2y2+y5;
323  def G=gmsring(t,"s");
324  setring(G);
325  list l0=gmscoeffs(gmspoly,0);
326  print(l0[1]);
327  list l1=gmscoeffs(gmspoly,1);
328  print(l1[1]);
329  list l=gmscoeffs(l0[2],1);
330  print(l[1]);
331}
332///////////////////////////////////////////////////////////////////////////////
333
334static proc mindegree(matrix A)
335{
336  int d=-1;
337
338  int i,j;
339  for(i=nrows(A);i>=1;i--)
340  {
341    for(j=ncols(A);j>=1;j--)
342    {
343      if(d==-1||(ord(A[i,j])<d&&ord(A[i,j])>-1))
344      {
345        d=ord(A[i,j]);
346      }
347    }
348  }
349
350  return(d);
351}
352///////////////////////////////////////////////////////////////////////////////
353
354static proc maxdegree(matrix A)
355{
356  int d=-1;
357
358  int i,j;
359  for(i=nrows(A);i>=1;i--)
360  {
361    for(j=ncols(A);j>=1;j--)
362    {
363      if(deg(A[i,j])>d)
364      {
365        d=deg(A[i,j]);
366      }
367    }
368  }
369
370  return(d);
371}
372///////////////////////////////////////////////////////////////////////////////
373
374static proc saturate()
375{
376  int mu=ncols(gmsbasis);
377  ideal r=gmspoly*gmsbasis;
378  matrix A0[mu][mu],C;
379  module H0;
380  module H,H1=freemodule(mu),freemodule(mu);
381  int k=-1;
382  list l;
383
384  dbprint(printlevel-voice+2,"// compute saturation of H''");
385  while(size(reduce(H,std(H0*var(1))))>0)
386  {
387    dbprint(printlevel-voice+2,"// compute matrix A of t");
388    k++;
389    dbprint(printlevel-voice+2,"// k="+string(k));
390    l=gmscoeffs(r,k);
391    C,r=l[1..2];
392    A0=A0+C;
393
394    dbprint(printlevel-voice+2,"// compute saturation step");
395    H0=H;
396    H1=jet(module(A0*H1+var(1)^2*diff(matrix(H1),var(1))),k);
397    H=H*var(1)+H1;
398  }
399
400  A0=A0-k*var(1);
401  dbprint(printlevel-voice+2,"// compute basis of saturation of H''");
402  H=std(H0);
403
404  dbprint(printlevel-voice+2,"// transform H'' to saturation of H''");
405  H0=division(freemodule(mu)*var(1)^k,H,k*deg(var(1)))[1];
406
407  return(A0,r,H,H0,k);
408}
409///////////////////////////////////////////////////////////////////////////////
410
411static proc tjet(matrix A0,ideal r,module H,int K0,int K)
412{
413  dbprint(printlevel-voice+2,"// compute matrix A of t");
414  dbprint(printlevel-voice+2,"// k="+string(K0+K+1));
415  list l=gmscoeffs(r,K0+K+1);
416  matrix C;
417  C,r=l[1..2];
418  A0=A0+C;
419  dbprint(printlevel-voice+2,"// transform A to saturation of H''");
420  matrix A=division(A0*H+var(1)^2*diff(matrix(H),var(1)),H,
421    (K+1)*deg(var(1)))[1]/var(1);
422  return(A,A0,r);
423}
424///////////////////////////////////////////////////////////////////////////////
425
426static proc eigenval(matrix A0,ideal r,module H,int K0)
427{
428  dbprint(printlevel-voice+2,
429    "// compute eigenvalues e with multiplicities m of A1");
430  matrix A;
431  A,A0,r=tjet(A0,r,H,K0,0);
432  list l=eigenvals(A);
433  def e,m=l[1..2];
434  dbprint(printlevel-voice+2,"// e="+string(e));
435  dbprint(printlevel-voice+2,"// m="+string(m));
436  return(e,m,A0,r);
437}
438///////////////////////////////////////////////////////////////////////////////
439
440static proc transform(matrix A,matrix A0,ideal r,module H,module H0,ideal e,
441  intvec m,int K0,int K,int opt)
442{
443  int mu=ncols(gmsbasis);
444
445  int i,j,k;
446  intvec d;
447  d[ncols(e)]=0;
448  if(opt)
449  {
450    dbprint(printlevel-voice+2,
451      "// compute rounded maximal differences d of e");
452    for(i=1;i<=ncols(e);i++)
453    {
454      d[i]=int(e[ncols(e)]-e[i]);
455    }
456  }
457  else
458  {
459    dbprint(printlevel-voice+2,
460      "// compute maximal integer differences d of e");
461    for(i=1;i<ncols(e);i++)
462    {
463      for(j=i+1;j<=ncols(e);j++)
464      {
465        k=int(e[j]-e[i]);
466        if(number(e[j]-e[i])==k)
467        {
468          if(k>d[i])
469          {
470            d[i]=k;
471          }
472          if(-k>d[j])
473          {
474            d[j]=-k;
475          }
476        }
477      }
478    }
479  }
480  dbprint(printlevel-voice+2,"// d="+string(d));
481
482  for(i,k=1,0;i<=size(d);i++)
483  {
484    if(k<d[i])
485    {
486      k=d[i];
487    }
488  }
489  A,A0,r=tjet(A0,r,H,K0,K+k);
490
491  module U,V;
492  if(k>0)
493  {
494    int i0,j0,i1,j1;
495    list l;
496
497    while(k>0)
498    {
499      dbprint(printlevel-voice+2,"// transform to separate eigenvalues");
500      U=0;
501      for(i=1;i<=ncols(e);i++)
502      {
503        U=U+syz(power(jet(A,0)-e[i],m[i]));
504      }
505      V=inverse(U);
506      A=V*A*U;
507      H=H*U;
508      H0=V*H0;
509
510      dbprint(printlevel-voice+2,"// transform to reduce maximum of d by 1");
511      for(i0,i=1,1;i0<=ncols(e);i0++)
512      {
513        for(i1=1;i1<=m[i0];i1,i=i1+1,i+1)
514        {
515          for(j0,j=1,1;j0<=ncols(e);j0++)
516          {
517            for(j1=1;j1<=m[j0];j1,j=j1+1,j+1)
518            {
519              if(d[i0]==0&&d[j0]>=1)
520              {
521                A[i,j]=A[i,j]*var(1);
522              }
523              if(d[i0]>=1&&d[j0]==0)
524              {
525                A[i,j]=A[i,j]/var(1);
526              }
527            }
528          }
529        }
530      }
531
532      H0=transpose(H0);
533      for(i0,i=1,1;i0<=ncols(e);i0++)
534      {
535        if(d[i0]>=1)
536        {
537          for(i1=1;i1<=m[i0];i1,i=i1+1,i+1)
538          {
539            H[i]=H[i]*var(1);
540          }
541          d[i0]=d[i0]-1;
542        }
543        else
544        {
545          for(i1=1;i1<=m[i0];i1,i=i1+1,i+1)
546          {
547            A[i,i]=A[i,i]-1;
548            H0[i]=H0[i]*var(1);
549          }
550          e[i0]=e[i0]-1;
551        }
552      }
553      H0=transpose(H0);
554
555      l=sppnf(list(e,d,m));
556      e,d,m=l[1..3];
557
558      k--;
559      K0++;
560    }
561
562    A=jet(A,K);
563  }
564
565  dbprint(printlevel-voice+2,"// transform to separate eigenvalues");
566  U=0;
567  for(i=1;i<=ncols(e);i++)
568  {
569    U=U+syz(power(jet(A,0)-e[i],m[i]));
570  }
571  V=inverse(U);
572  A=V*A*U;
573  H=H*U;
574  H0=V*H0;
575
576  return(A,A0,r,H,H0,e,m,K0);
577}
578///////////////////////////////////////////////////////////////////////////////
579
580proc bernstein(poly t)
581"USAGE:    bernstein(t); poly t
582ASSUME:   characteristic 0; local degree ordering;
583          isolated critical point 0 of t
584RETURN:
585@format
586list bs;  Bernstein-Sato polynomial b(s) of t
587  ideal bs[1];
588    number bs[1][i];  i-th root of b(s)
589  intvec bs[2];
590    int bs[2][i];  multiplicity of i-th root of b(s)
591@end format
592KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
593          Bernstein-Sato polynomial
594EXAMPLE:  example bernstein; shows examples
595"
596{
597  def @R=basering;
598  int n=nvars(@R)-1;
599  def @G=gmsring(t,"s");
600  setring(@G);
601
602  matrix A;
603  module U0;
604  ideal e;
605  intvec m;
606
607  def A0,r,H,H0,K0=saturate();
608  A,A0,r=tjet(A0,r,H,K0,0);
609  list l=minipoly(A);
610  e,m=l[1..2];
611  e=-e;
612  l=spnf(spadd(list(e,m),list(ideal(-1),intvec(1))));
613
614  setring(@R);
615  list l=imap(@G,l);
616  kill @G,gmsmaxdeg;
617
618  return(l);
619}
620example
621{ "EXAMPLE:"; echo=2;
622  ring R=0,(x,y),ds;
623  poly t=x5+x2y2+y5;
624  bernstein(t);
625}
626///////////////////////////////////////////////////////////////////////////////
627
628proc monodromy(poly t)
629"USAGE:    monodromy(t); poly t
630ASSUME:   characteristic 0; local degree ordering;
631          isolated critical point 0 of t
632RETURN:
633@format
634list l;  Jordan data jordan(M) of monodromy matrix exp(-2*pi*i*M)
635  ideal l[1];
636    number l[1][i];  eigenvalue of i-th Jordan block of M
637  intvec l[2];
638    int l[2][i];  size of i-th Jordan block of M
639  intvec l[3];
640    int l[3][i];  multiplicity of i-th Jordan block of M
641@end format
642SEE ALSO: mondromy_lib, linalg_lib
643KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice; monodromy
644EXAMPLE:  example monodromy; shows examples
645"
646{
647  def @R=basering;
648  int n=nvars(@R)-1;
649  def @G=gmsring(t,"s");
650  setring(@G);
651
652  matrix A;
653  module U0;
654  ideal e;
655  intvec m;
656
657  def A0,r,H,H0,K0=saturate();
658  e,m,A0,r=eigenval(A0,r,H,K0);
659  A,A0,r,H,H0,e,m,K0=transform(A,A0,r,H,H0,e,m,K0,0,0);
660
661  list l=jordan(A,e,m);
662  setring(@R);
663  list l=imap(@G,l);
664  kill @G,gmsmaxdeg;
665
666  return(l);
667}
668example
669{ "EXAMPLE:"; echo=2;
670  ring R=0,(x,y),ds;
671  poly t=x5+x2y2+y5;
672  monodromy(t);
673}
674///////////////////////////////////////////////////////////////////////////////
675
676proc spectrum(poly t)
677"USAGE:    spectrum(t); poly t
678ASSUME:   characteristic 0; local degree ordering;
679          isolated critical point 0 of t
680RETURN:
681@format
682list sp;  singularity spectrum of t
683  ideal sp[1];
684    number sp[1][i];  i-th spectral number
685  intvec sp[2];
686    int sp[2][i];  multiplicity of i-th spectral number
687@end format
688SEE ALSO: spectrum_lib
689KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
690          mixed Hodge structure; V-filtration; spectrum
691EXAMPLE:  example spectrum; shows examples
692"
693{
694  list l=vwfilt(t);
695  return(spnf(list(l[1],l[3])));
696}
697example
698{ "EXAMPLE:"; echo=2;
699  ring R=0,(x,y),ds;
700  poly t=x5+x2y2+y5;
701  spprint(spectrum(t));
702}
703///////////////////////////////////////////////////////////////////////////////
704
705proc sppairs(poly t)
706"USAGE:    sppairs(t); poly t
707ASSUME:   characteristic 0; local degree ordering;
708          isolated critical point 0 of t
709RETURN:
710@format
711list spp;  spectral pairs of t
712  ideal spp[1];
713    number spp[1][i];  V-filtration index of i-th spectral pair
714  intvec spp[2];
715    int spp[2][i];  weight filtration index of i-th spectral pair
716  intvec spp[3];
717    int spp[3][i];  multiplicity of i-th spectral pair
718@end format
719SEE ALSO: spectrum_lib
720KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
721          mixed Hodge structure; V-filtration; weight filtration;
722          spectrum; spectral pairs
723EXAMPLE:  example sppairs; shows examples
724"
725{
726  list l=vwfilt(t);
727  return(list(l[1],l[2],l[3]));
728}
729example
730{ "EXAMPLE:"; echo=2;
731  ring R=0,(x,y),ds;
732  poly t=x5+x2y2+y5;
733  sppprint(sppairs(t));
734}
735///////////////////////////////////////////////////////////////////////////////
736
737proc vfilt(poly t)
738"USAGE:    vfilt(t); poly t
739ASSUME:   characteristic 0; local degree ordering;
740          isolated critical point 0 of t
741RETURN:
742@format
743list v;  V-filtration on H''/s*H''
744  ideal v[1];
745    number v[1][i];  V-filtration index of i-th spectral number
746  intvec v[2];
747    int v[2][i];  multiplicity of i-th spectral number
748  list v[3];
749    module v[3][i];  vector space of i-th graded part in terms of v[4]
750  ideal v[4];  monomial vector space basis of H''/s*H''
751  ideal v[5];  standard basis of Jacobian ideal
752@end format
753SEE ALSO: spectrum_lib
754KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
755          mixed Hodge structure; V-filtration; spectrum
756EXAMPLE:  example vfilt; shows examples
757"
758{
759  list l=vwfilt(t);
760  return(spnf(list(l[1],l[3],l[4]))+list(l[5],l[6]));
761}
762example
763{ "EXAMPLE:"; echo=2;
764  ring R=0,(x,y),ds;
765  poly t=x5+x2y2+y5;
766  vfilt(t);
767}
768///////////////////////////////////////////////////////////////////////////////
769
770proc vwfilt(poly t)
771"USAGE:    vwfilt(t); poly t
772ASSUME:   characteristic 0; local degree ordering;
773          isolated critical point 0 of t
774RETURN:
775@format
776list vw;  weighted V-filtration on H''/s*H''
777  ideal vw[1];
778    number vw[1][i];  V-filtration index of i-th spectral pair
779  intvec vw[2];
780    int vw[2][i];  weight filtration index of i-th spectral pair
781  intvec vw[3];
782    int vw[3][i];  multiplicity of i-th spectral pair
783  list vw[4];
784    module vw[4][i];  vector space of i-th graded part in terms of vw[5]
785  ideal vw[5];  monomial vector space basis of H''/s*H''
786  ideal vw[6];  standard basis of Jacobian ideal
787@end format
788SEE ALSO: spectrum_lib
789KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
790          mixed Hodge structure; V-filtration; weight filtration;
791          spectrum; spectral pairs
792EXAMPLE:  example vwfilt; shows examples
793"
794{
795  def @R=basering;
796  int n=nvars(@R)-1;
797  def @G=gmsring(t,"s");
798  setring(@G);
799
800  int mu=ncols(gmsbasis);
801  matrix A;
802  ideal e;
803  intvec m;
804
805  def A0,r,H,H0,K0=saturate();
806  e,m,A0,r=eigenval(A0,r,H,K0);
807  A,A0,r,H,H0,e,m,K0=transform(A,A0,r,H,H0,e,m,K0,0,1);
808
809  dbprint(printlevel-voice+2,"// compute weight filtration basis");
810  list l=jordanbasis(A,e,m);
811  def U,v=l[1..2];
812  kill l;
813  vector u0;
814  int v0;
815  int i,j,k,l;
816  for(k,l=1,1;l<=ncols(e);k,l=k+m[l],l+1)
817  {
818    for(i=k+m[l]-1;i>=k+1;i--)
819    {
820      for(j=i-1;j>=k;j--)
821      {
822        if(v[i]>v[j])
823        {
824          v0=v[i];v[i]=v[j];v[j]=v0;
825          u0=U[i];U[i]=U[j];U[j]=u0;
826        }
827      }
828    }
829  }
830
831  dbprint(printlevel-voice+2,"// transform to weight filtration basis");
832  matrix V=inverse(U);
833  A=V*A*U;
834  dbprint(printlevel-voice+2,"// compute standard basis of H''");
835  H=H*U;
836  H0=std(V*H0);
837
838  dbprint(printlevel-voice+2,"// compute spectral pairs");
839  ideal a;
840  intvec w;
841  for(i=1;i<=mu;i++)
842  {
843    j=leadexp(H0[i])[nvars(basering)+1];
844    a[i]=A[j,j]+ord(H0[i])/deg(var(1))-1;
845    w[i]=v[j]+n;
846  }
847  H=H*H0;
848  H=simplify(jet(H/var(1)^(mindegree(H)/deg(var(1))),0),1);
849
850  kill l;
851  list l=sppnf(list(a,w,H))+list(gmsbasis,gmsstd);
852  setring(@R);
853  list l=imap(@G,l);
854  kill @G,gmsmaxdeg;
855  attrib(l[5],"isSB",1);
856
857  return(l);
858}
859example
860{ "EXAMPLE:"; echo=2;
861  ring R=0,(x,y),ds;
862  poly t=x5+x2y2+y5;
863  vwfilt(t);
864}
865///////////////////////////////////////////////////////////////////////////////
866
867static proc fsplit(ideal e0,intvec m0,matrix A,module H,module H0)
868{
869  int mu=ncols(gmsbasis);
870
871  dbprint(printlevel-voice+2,"// compute standard basis of H''");
872  H0=std(H0);
873  H0=simplify(H0,1);
874
875  dbprint(printlevel-voice+2,"// compute Hodge filtration");
876  int i,j,k;
877  ideal e;
878  intvec m;
879  e[mu]=0;
880  for(i=1;i<=ncols(e0);i++)
881  {
882    for(j=m0[i];j>=1;j--)
883    {
884      k++;
885      e[k]=e0[i];
886      m[k]=i;
887    }
888  }
889
890  number n,n0;
891  vector v,v0;
892  list F;
893  for(i=ncols(e0);i>=1;i--)
894  {
895    F[i]=module(matrix(0,mu,1));
896  }
897  for(i=mu;i>=1;i--)
898  {
899    v=H0[i];
900    v0=lead(v);
901    n0=leadcoef(e[leadexp(v0)[nvars(basering)+1]])+leadexp(v0)[1];
902    v=v-lead(v);
903    while(v!=0)
904    {
905      n=leadcoef(e[leadexp(v)[nvars(basering)+1]])+leadexp(v)[1];
906      if(n==n0)
907      {
908        v0=v0+lead(v);
909        v=v-lead(v);
910      }
911      else
912      {
913        v=0;
914      }
915    }
916    j=m[leadexp(v0)[nvars(basering)+1]];
917    F[j]=F[j]+v0;
918  }
919
920  dbprint(printlevel-voice+2,"// compute splitting of Hodge filtration");
921  matrix A0=jet(A,0);
922  module U,U0,U1,U2;
923  matrix N;
924  s=0;
925  for(i=size(F);i>=1;i--)
926  {
927    N=A0-e0[i];
928    U0=0;
929    while(size(F[i])>0)
930    {
931      U1=jet(F[i],0);
932      k=0;
933      while(size(U1)>0)
934      {
935        for(j=ncols(U1);j>=1;j--)
936        {
937          if(size(reduce(U1[j],std(U0)))>0)
938          {
939            U0=U0+U1[j];
940          }
941        }
942        U1=N*U1;
943        k++;
944      }
945      F[i]=module(F[i]/var(1));
946    }
947    U=U0+U;
948  }
949
950  dbprint(printlevel-voice+2,"// transform to Hodge splitting basis");
951  H=H*U;
952  H0=lift(U,H0);
953  A=lift(U,A*U);
954
955  return(e,A,H,H0);
956}
957///////////////////////////////////////////////////////////////////////////////
958
959static proc glift(ideal e,matrix A,module H,module H0,int K)
960{
961  poly s=var(1);
962  int mu=ncols(gmsbasis);
963
964  dbprint(printlevel-voice+2,"// compute standard basis of H''");
965  H0=std(H0);
966  H0=simplify(H0,1);
967
968  int i,j,k;
969  ideal v;
970  for(i=mu;i>=1;i--)
971  {
972    v[i]=e[leadexp(H0[i])[nvars(basering)+1]]+leadexp(H0[i])[1];
973  }
974
975  dbprint(printlevel-voice+2,
976    "// compute matrix A0 of t w.r.t. good basis H0 of H''");
977  number c;
978  matrix h0[mu][1];
979  matrix m[mu][1];
980  matrix a0[mu][1];
981  matrix A0[mu][mu];
982  module M=H0;
983  module N=jet(s*A*matrix(H0)+s^2*diff(matrix(H0),s),K+1);
984  while(size(N)>0)
985  {
986    j=mu;
987    for(k=mu-1;k>=1;k--)
988    {
989      if(N[k]>N[j])
990      {
991        j=k;
992      }
993    }
994    i=mu;
995    while(leadexp(M[i])[nvars(basering)+1]!=leadexp(N[j])[nvars(basering)+1])
996    {
997      i--;
998    }
999    k=leadexp(N[j])[1]-leadexp(M[i])[1];
1000    if(k==0||i==j)
1001    {
1002      dbprint(printlevel-voice+3,"// compute A0["+string(i)+","+string(j)+"]");
1003      c=leadcoef(N[j])/leadcoef(M[i]);
1004      A0[i,j]=A0[i,j]+c*s^k;
1005      N[j]=jet(N[j]-c*s^k*M[i],K+1);
1006    }
1007    else
1008    {
1009      dbprint(printlevel-voice+3,
1010        "// reduce H0["+string(j)+"] with H0["+string(i)+"]");
1011      c=leadcoef(N[j])/leadcoef(M[i])/(1-k-leadcoef(v[i])+leadcoef(v[j]));
1012      H0[j]=H0[j]+c*s^(k-1)*H0[i];
1013      M[j]=M[j]+c*s^(k-1)*M[i];
1014      h0=c*s^(k-1)*H0[i];
1015      N[j]=N[j]+jet(s*A*h0+s^2*diff(h0,s),K+1)[1];
1016      m=M[i];
1017      a0=transpose(A0)[j];
1018      N=N-jet(c*s^(k-1)*m*transpose(a0),K+1);
1019    }
1020  }
1021
1022  H0=H*H0;
1023  H0=H0/var(1)^(mindegree(H0)/deg(var(1)));
1024
1025  return(A0,H0);
1026}
1027///////////////////////////////////////////////////////////////////////////////
1028
1029proc tmatrix(poly t)
1030"USAGE:    tmatrix(t); poly t
1031ASSUME:   characteristic 0; local degree ordering;
1032          isolated critical point 0 of t
1033RETURN:
1034@format
1035list l=A0,A1,T,M;
1036  matrix A0,A1;  t=A0+s*A1+s^2*(d/ds) on H'' w.r.t. C[[s]]-basis M*T
1037  module T;  C-basis of C^mu
1038  ideal M;  monomial C-basis of H''/sH''
1039@end format
1040KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
1041          mixed Hodge structure; V-filtration; weight filtration;
1042          monodromy; spectrum; spectral pairs; good basis
1043EXAMPLE:  example tmatrix; shows examples
1044"
1045{
1046  def @R=basering;
1047  int n=nvars(@R)-1;
1048  def @G=gmsring(t,"s");
1049  setring(@G);
1050
1051  int mu=ncols(gmsbasis);
1052  matrix A;
1053  module U0;
1054  ideal e;
1055  intvec m;
1056
1057  def A0,r,H,H0,K0=saturate();
1058  e,m,A0,r=eigenval(A0,r,H,K0);
1059  A,A0,r,H,H0,e,m,K0=transform(A,A0,r,H,H0,e,m,K0,K0+int(e[ncols(e)]-e[1]),1);
1060  A,H0=glift(fsplit(e,m,A,H,H0),K0);
1061
1062  A0=jet(A,0);
1063  A=jet(A/var(1),0);
1064
1065  list l=A0,A,H0,gmsbasis;
1066  setring(@R);
1067  list l=imap(@G,l);
1068  kill @G,gmsmaxdeg;
1069
1070  return(l);
1071}
1072example
1073{ "EXAMPLE:"; echo=2;
1074  ring R=0,(x,y),ds;
1075  poly t=x5+x2y2+y5;
1076  list l=tmatrix(t);
1077  print(l[1]);
1078  print(l[2]);
1079  print(l[3]);
1080  print(l[4]);
1081}
1082///////////////////////////////////////////////////////////////////////////////
1083
1084proc endvfilt(list v)
1085"USAGE:   endvfilt(v); list v
1086ASSUME:  v returned by vfilt
1087RETURN:
1088@format
1089list ev;  V-filtration on Jacobian algebra
1090  ideal ev[1];
1091    number ev[1][i];  i-th V-filtration index
1092  intvec ev[2];
1093    int ev[2][i];  i-th multiplicity
1094  list ev[3];
1095    module ev[3][i];  vector space of i-th graded part in terms of ev[4]
1096  ideal ev[4];  monomial vector space basis of Jacobian algebra
1097  ideal ev[5];  standard basis of Jacobian ideal
1098@end format
1099KEYWORDS: singularities; Gauss-Manin system; Brieskorn lattice;
1100          mixed Hodge structure; V-filtration; endomorphism filtration
1101EXAMPLE:  example endvfilt; shows examples
1102"
1103{
1104  def a,d,V,m,g=v[1..5];
1105  attrib(g,"isSB",1);
1106  int mu=ncols(m);
1107
1108  module V0=V[1];
1109  for(int i=2;i<=size(V);i++)
1110  {
1111    V0=V0,V[i];
1112  }
1113
1114  dbprint(printlevel-voice+2,"// compute multiplication in Jacobian algebra");
1115  list M;
1116  module U=freemodule(ncols(m));
1117  for(i=ncols(m);i>=1;i--)
1118  {
1119    M[i]=division(coeffs(reduce(m[i]*m,g,U),m)*V0,V0)[1];
1120  }
1121
1122  int j,k,i0,j0,i1,j1;
1123  number b0=number(a[1]-a[ncols(a)]);
1124  number b1,b2;
1125  matrix M0;
1126  module L;
1127  list v0=freemodule(ncols(m));
1128  ideal a0=b0;
1129  list l;
1130
1131  while(b0<number(a[ncols(a)]-a[1]))
1132  {
1133    dbprint(printlevel-voice+2,"// find next possible index");
1134    b1=number(a[ncols(a)]-a[1]);
1135    for(j=ncols(a);j>=1;j--)
1136    {
1137      for(i=ncols(a);i>=1;i--)
1138      {
1139        b2=number(a[i]-a[j]);
1140        if(b2>b0&&b2<b1)
1141        {
1142          b1=b2;
1143        }
1144        else
1145        {
1146          if(b2<=b0)
1147          {
1148            i=0;
1149          }
1150        }
1151      }
1152    }
1153    b0=b1;
1154
1155    l=ideal();
1156    for(k=ncols(m);k>=2;k--)
1157    {
1158      l=l+list(ideal());
1159    }
1160
1161    dbprint(printlevel-voice+2,"// collect conditions for EV["+string(b0)+"]");
1162    j=ncols(a);
1163    j0=mu;
1164    while(j>=1)
1165    {
1166      i0=1;
1167      i=1;
1168      while(i<ncols(a)&&a[i]<a[j]+b0)
1169      {
1170        i0=i0+d[i];
1171        i++;
1172      }
1173      if(a[i]<a[j]+b0)
1174      {
1175        i0=i0+d[i];
1176        i++;
1177      }
1178      for(k=1;k<=ncols(m);k++)
1179      {
1180        M0=M[k];
1181        if(i0>1)
1182        {
1183          l[k]=l[k],M0[1..i0-1,j0-d[j]+1..j0];
1184        }
1185      }
1186      j0=j0-d[j];
1187      j--;
1188    }
1189
1190    dbprint(printlevel-voice+2,"// compose condition matrix");
1191    L=transpose(module(l[1]));
1192    for(k=2;k<=ncols(m);k++)
1193    {
1194      L=L,transpose(module(l[k]));
1195    }
1196
1197    dbprint(printlevel-voice+2,"// compute kernel of condition matrix");
1198    v0=v0+list(syz(L));
1199    a0=a0,b0;
1200  }
1201
1202  dbprint(printlevel-voice+2,"// compute graded parts");
1203  option(redSB);
1204  for(i=1;i<size(v0);i++)
1205  {
1206    v0[i+1]=std(v0[i+1]);
1207    v0[i]=std(reduce(v0[i],v0[i+1]));
1208  }
1209  option(noredSB);
1210
1211  dbprint(printlevel-voice+2,"// remove trivial graded parts");
1212  i=1;
1213  while(size(v0[i])==0)
1214  {
1215    i++;
1216  }
1217  list v1=v0[i];
1218  intvec d1=size(v0[i]);
1219  ideal a1=a0[i];
1220  i++;
1221  while(i<=size(v0))
1222  {
1223    if(size(v0[i])>0)
1224    {
1225      v1=v1+list(v0[i]);
1226      d1=d1,size(v0[i]);
1227      a1=a1,a0[i];
1228    }
1229    i++;
1230  }
1231  return(list(a1,d1,v1,m,g));
1232}
1233example
1234{ "EXAMPLE:"; echo=2;
1235  ring R=0,(x,y),ds;
1236  poly t=x5+x2y2+y5;
1237  endvfilt(vfilt(t));
1238}
1239///////////////////////////////////////////////////////////////////////////////
1240
1241proc sppnf(list sp)
1242"USAGE:   sppnf(list(a,w[,m])); ideal a, intvec w, intvec m
1243ASSUME:  ncols(e)==size(w)==size(m)
1244RETURN:  order (a[i][,w[i]]) with multiplicity m[i] lexicographically
1245EXAMPLE: example sppnf; shows examples
1246"
1247{
1248  ideal a=sp[1];
1249  intvec w=sp[2];
1250  int n=ncols(a);
1251  intvec m;
1252  list V;
1253  module v;
1254  int i,j;
1255  for(i=3;i<=size(sp);i++)
1256  {
1257    if(typeof(sp[i])=="intvec")
1258    {
1259      m=sp[i];
1260    }
1261    if(typeof(sp[i])=="module")
1262    {
1263      v=sp[i];
1264      for(j=n;j>=1;j--)
1265      {
1266        V[j]=module(v[j]);
1267      }
1268    }
1269    if(typeof(sp[i])=="list")
1270    {
1271      V=sp[i];
1272    }
1273  }
1274  if(m==0)
1275  {
1276    for(i=n;i>=1;i--)
1277    {
1278      m[i]=1;
1279    }
1280  }
1281
1282  int k;
1283  ideal a0;
1284  intvec w0,m0;
1285  list V0;
1286  number a1;
1287  int w1,m1;
1288  for(i=n;i>=1;i--)
1289  {
1290    if(m[i]!=0)
1291    {
1292      for(j=i-1;j>=1;j--)
1293      {
1294        if(m[j]!=0)
1295        {
1296          if(number(a[i])>number(a[j])||
1297            (number(a[i])==number(a[j])&&w[i]<w[j]))
1298          {
1299            a1=number(a[i]);
1300            a[i]=a[j];
1301            a[j]=a1;
1302            w1=w[i];
1303            w[i]=w[j];
1304            w[j]=w1;
1305            m1=m[i];
1306            m[i]=m[j];
1307            m[j]=m1;
1308            if(size(V)>0)
1309            {
1310              v=V[i];
1311              V[i]=V[j];
1312              V[j]=v;
1313            }
1314          }
1315          if(number(a[i])==number(a[j])&&w[i]==w[j])
1316          {
1317            m[i]=m[i]+m[j];
1318            m[j]=0;
1319            if(size(V)>0)
1320            {
1321              V[i]=V[i]+V[j];
1322            }
1323          }
1324        }
1325      }
1326      k++;
1327      a0[k]=a[i];
1328      w0[k]=w[i];
1329      m0[k]=m[i];
1330      if(size(V)>0)
1331      {
1332        V0[k]=V[i];
1333      }
1334    }
1335  }
1336
1337  if(size(V0)>0)
1338  {
1339    n=size(V0);
1340    module U=std(V0[n]);
1341    for(i=n-1;i>=1;i--)
1342    {
1343      V0[i]=simplify(reduce(V0[i],U),1);
1344      if(i>=2)
1345      {
1346        U=std(U+V0[i]);
1347      }
1348    }
1349  }
1350
1351  if(k>0)
1352  {
1353    sp=a0,w0,m0;
1354    if(size(V0)>0)
1355    {
1356      sp[4]=V0;
1357    }
1358  }
1359  return(sp);
1360}
1361example
1362{ "EXAMPLE:"; echo=2;
1363  ring R=0,(x,y),ds;
1364  list sp=list(ideal(-1/2,-3/10,-3/10,-1/10,-1/10,0,1/10,1/10,3/10,3/10,1/2),
1365    intvec(2,1,1,1,1,1,1,1,1,1,0));
1366  sppprint(sppnf(sp));
1367}
1368///////////////////////////////////////////////////////////////////////////////
1369
1370proc sppprint(list spp)
1371"USAGE:   sppprint(spp); list spp
1372RETURN:  string s;  spectral pairs spp
1373EXAMPLE: example sppprint; shows examples
1374"
1375{
1376  string s;
1377  for(int i=1;i<size(spp[3]);i++)
1378  {
1379    s=s+"(("+string(spp[1][i])+","+string(spp[2][i])+"),"
1380      +string(spp[3][i])+"),";
1381  }
1382  s=s+"(("+string(spp[1][i])+","+string(spp[2][i])+"),"+string(spp[3][i])+")";
1383  return(s);
1384}
1385example
1386{ "EXAMPLE:"; echo=2;
1387  ring R=0,(x,y),ds;
1388  list spp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(2,1,1,1,1,1,0),
1389    intvec(1,2,2,1,2,2,1));
1390  sppprint(spp);
1391}
1392///////////////////////////////////////////////////////////////////////////////
1393
1394proc spadd(list sp1,list sp2)
1395"USAGE:   spadd(sp1,sp2); list sp1, list sp2
1396RETURN:  list sp;  sum of spectra sp1 and sp2
1397EXAMPLE: example spadd; shows examples
1398"
1399{
1400  ideal s;
1401  intvec m;
1402  int i,i1,i2=1,1,1;
1403  while(i1<=size(sp1[2])||i2<=size(sp2[2]))
1404  {
1405    if(i1<=size(sp1[2]))
1406    {
1407      if(i2<=size(sp2[2]))
1408      {
1409        if(number(sp1[1][i1])<number(sp2[1][i2]))
1410        {
1411          s[i]=sp1[1][i1];
1412          m[i]=sp1[2][i1];
1413          i++;
1414          i1++;
1415        }
1416        else
1417        {
1418          if(number(sp1[1][i1])>number(sp2[1][i2]))
1419          {
1420            s[i]=sp2[1][i2];
1421            m[i]=sp2[2][i2];
1422            i++;
1423            i2++;
1424          }
1425          else
1426          {
1427            if(sp1[2][i1]+sp2[2][i2]!=0)
1428            {
1429              s[i]=sp1[1][i1];
1430              m[i]=sp1[2][i1]+sp2[2][i2];
1431              i++;
1432            }
1433            i1++;
1434            i2++;
1435          }
1436        }
1437      }
1438      else
1439      {
1440        s[i]=sp1[1][i1];
1441        m[i]=sp1[2][i1];
1442        i++;
1443        i1++;
1444      }
1445    }
1446    else
1447    {
1448      s[i]=sp2[1][i2];
1449      m[i]=sp2[2][i2];
1450      i++;
1451      i2++;
1452    }
1453  }
1454  return(list(s,m));
1455}
1456example
1457{ "EXAMPLE:"; echo=2;
1458  ring R=0,(x,y),ds;
1459  list sp1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1460  spprint(sp1);
1461  list sp2=list(ideal(-1/6,1/6),intvec(1,1));
1462  spprint(sp2);
1463  spprint(spadd(sp1,sp2));
1464}
1465///////////////////////////////////////////////////////////////////////////////
1466
1467proc spsub(list sp1,list sp2)
1468"USAGE:   spsub(sp1,sp2); list sp1, list sp2
1469RETURN:  list sp;  difference of spectra sp1 and sp2
1470EXAMPLE: example spsub; shows examples
1471"
1472{
1473  return(spadd(sp1,spmul(sp2,-1)));
1474}
1475example
1476{ "EXAMPLE:"; echo=2;
1477  ring R=0,(x,y),ds;
1478  list sp1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1479  spprint(sp1);
1480  list sp2=list(ideal(-1/6,1/6),intvec(1,1));
1481  spprint(sp2);
1482  spprint(spsub(sp1,sp2));
1483}
1484///////////////////////////////////////////////////////////////////////////////
1485
1486proc spmul(list #)
1487"USAGE:   spmul(sp0,k); list sp0, int[vec] k
1488RETURN:  list sp;  linear combination of spectra sp0 with coefficients k
1489EXAMPLE: example spmul; shows examples
1490"
1491{
1492  if(size(#)==2)
1493  {
1494    if(typeof(#[1])=="list")
1495    {
1496      if(typeof(#[2])=="int")
1497      {
1498        return(list(#[1][1],#[1][2]*#[2]));
1499      }
1500      if(typeof(#[2])=="intvec")
1501      {
1502        list sp0=list(ideal(),intvec(0));
1503        for(int i=size(#[2]);i>=1;i--)
1504        {
1505          sp0=spadd(sp0,spmul(#[1][i],#[2][i]));
1506        }
1507        return(sp0);
1508      }
1509    }
1510  }
1511  return(list(ideal(),intvec(0)));
1512}
1513example
1514{ "EXAMPLE:"; echo=2;
1515  ring R=0,(x,y),ds;
1516  list sp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1517  spprint(sp);
1518  spprint(spmul(sp,2));
1519  list sp1=list(ideal(-1/6,1/6),intvec(1,1));
1520  spprint(sp1);
1521  list sp2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
1522  spprint(sp2);
1523  spprint(spmul(list(sp1,sp2),intvec(1,2)));
1524}
1525///////////////////////////////////////////////////////////////////////////////
1526
1527proc spissemicont(list sp,list #)
1528"USAGE:   spissemicont(sp[,1]); list sp, int opt
1529RETURN:
1530@format
1531int k=
1532  1;  if sum of sp is positive on all intervals [a,a+1) [and (a,a+1)]
1533  0;  if sum of sp is negative on some interval [a,a+1) [or (a,a+1)]
1534@end format
1535EXAMPLE: example spissemicont; shows examples
1536"
1537{
1538  int opt=0;
1539  if(size(#)>0)
1540  {
1541    if(typeof(#[1])=="int")
1542    {
1543      opt=1;
1544    }
1545  }
1546  int i,j,k;
1547  i=1;
1548  while(i<=size(sp[2])-1)
1549  {
1550    j=i+1;
1551    k=0;
1552    while(j+1<=size(sp[2])&&number(sp[1][j])<=number(sp[1][i])+1)
1553    {
1554      if(opt==0||number(sp[1][j])<number(sp[1][i])+1)
1555      {
1556        k=k+sp[2][j];
1557      }
1558      j++;
1559    }
1560    if(j==size(sp[2])&&number(sp[1][j])<=number(sp[1][i])+1)
1561    {
1562      if(opt==0||number(sp[1][j])<number(sp[1][i])+1)
1563      {
1564        k=k+sp[2][j];
1565      }
1566    }
1567    if(k<0)
1568    {
1569      return(0);
1570    }
1571    i++;
1572  }
1573  return(1);
1574}
1575example
1576{ "EXAMPLE:"; echo=2;
1577  ring R=0,(x,y),ds;
1578  list sp1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1579  spprint(sp1);
1580  list sp2=list(ideal(-1/6,1/6),intvec(1,1));
1581  spprint(sp2);
1582  spissemicont(spsub(sp1,spmul(sp2,3)));
1583  spissemicont(spsub(sp1,spmul(sp2,4)));
1584}
1585///////////////////////////////////////////////////////////////////////////////
1586
1587proc spsemicont(list sp0,list sp,list #)
1588"USAGE:   spsemicont(sp0,sp,k[,1]); list sp0, list sp
1589RETURN:
1590@format
1591list l;
1592  intvec l[i];  if the spectra sp0 occur with multiplicities k
1593                in a deformation of a [quasihomogeneous] singularity
1594                with spectrum sp then k<=l[i]
1595@end format
1596EXAMPLE: example spsemicont; shows examples
1597"
1598{
1599  list l,l0;
1600  int i,j,k;
1601  while(spissemicont(sp0,#))
1602  {
1603    if(size(sp)>1)
1604    {
1605      l0=spsemicont(sp0,list(sp[1..size(sp)-1]));
1606      for(i=1;i<=size(l0);i++)
1607      {
1608        if(size(l)>0)
1609        {
1610          j=1;
1611          while(j<size(l)&&l[j]!=l0[i])
1612          {
1613            j++;
1614          }
1615          if(l[j]==l0[i])
1616          {
1617            l[j][size(sp)]=k;
1618          }
1619          else
1620          {
1621            l0[i][size(sp)]=k;
1622            l=l+list(l0[i]);
1623          }
1624        }
1625        else
1626        {
1627          l=l0;
1628        }
1629      }
1630    }
1631    sp0=spsub(sp0,sp[size(sp)]);
1632    k++;
1633  }
1634  if(size(sp)>1)
1635  {
1636    return(l);
1637  }
1638  else
1639  {
1640    return(list(intvec(k-1)));
1641  }
1642}
1643example
1644{ "EXAMPLE:"; echo=2;
1645  ring R=0,(x,y),ds;
1646  list sp0=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1647  spprint(sp0);
1648  list sp1=list(ideal(-1/6,1/6),intvec(1,1));
1649  spprint(sp1);
1650  list sp2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
1651  spprint(sp2);
1652  list sp=sp1,sp2;
1653  list l=spsemicont(sp0,sp);
1654  l;
1655  spissemicont(spsub(sp0,spmul(sp,l[1])));
1656  spissemicont(spsub(sp0,spmul(sp,l[1]-1)));
1657  spissemicont(spsub(sp0,spmul(sp,l[1]+1)));
1658}
1659///////////////////////////////////////////////////////////////////////////////
1660
1661proc spmilnor(list sp)
1662"USAGE:   spmilnor(sp); list sp
1663RETURN:  int mu;  Milnor number of spectrum sp
1664EXAMPLE: example spmilnor; shows examples
1665"
1666{
1667  return(sum(sp[2]));
1668}
1669example
1670{ "EXAMPLE:"; echo=2;
1671  ring R=0,(x,y),ds;
1672  list sp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1673  spprint(sp);
1674  spmilnor(sp);
1675}
1676///////////////////////////////////////////////////////////////////////////////
1677
1678proc spgeomgenus(list sp)
1679"USAGE:   spgeomgenus(sp); list sp
1680RETURN:  int g;  geometrical genus of spectrum sp
1681EXAMPLE: example spgeomgenus; shows examples
1682"
1683{
1684  int g=0;
1685  int i=1;
1686  while(i+1<=size(sp[2])&&number(sp[1][i])<=number(0))
1687  {
1688    g=g+sp[2][i];
1689    i++;
1690  }
1691  if(i==size(sp[2])&&number(sp[1][i])<=number(0))
1692  {
1693    g=g+sp[2][i];
1694  }
1695  return(g);
1696}
1697example
1698{ "EXAMPLE:"; echo=2;
1699  ring R=0,(x,y),ds;
1700  list sp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1701  spprint(sp);
1702  spgeomgenus(sp);
1703}
1704///////////////////////////////////////////////////////////////////////////////
1705
1706proc spgamma(list sp)
1707"USAGE:   spgamma(sp); list sp
1708RETURN:  number gamma;  gamma invariant of spectrum sp
1709EXAMPLE: example spgamma; shows examples
1710"
1711{
1712  int i,j;
1713  number g=0;
1714  for(i=1;i<=ncols(sp[1]);i++)
1715  {
1716    for(j=1;j<=sp[2][i];j++)
1717    {
1718      g=g+(number(sp[1][i])-number(nvars(basering)-2)/2)^2;
1719    }
1720  }
1721  g=-g/4+sum(sp[2])*number(sp[1][ncols(sp[1])]-sp[1][1])/48;
1722  return(g);
1723}
1724example
1725{ "EXAMPLE:"; echo=2;
1726  ring R=0,(x,y),ds;
1727  list sp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1728  spprint(sp);
1729  spgamma(sp);
1730}
1731///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.