source: git/Singular/LIB/gmssing.lib @ 4c20ee

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