source: git/Singular/LIB/gaussman.lib @ 8960ec

spielwiese
Last change on this file since 8960ec was 8960ec, checked in by Mathias Schulze <mschulze@…>, 23 years ago
*mschulze: added spmilnor git-svn-id: file:///usr/local/Singular/svn/trunk@5327 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.5 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: gaussman.lib,v 1.36 2001-03-20 18:19:48 mschulze Exp $";
3category="Singularities";
4
5info="
6LIBRARY:  gaussman.lib  Gauss-Manin Connection of a Singularity
7
8AUTHOR:   Mathias Schulze, email: mschulze@mathematik.uni-kl.de
9
10OVERVIEW: A library to compute invariants related to the Gauss-Manin connection
11          of a an isolated hypersurface singularity
12
13PROCEDURES:
14 monodromy(f[,opt]);        monodromy matrix, spectrum of monodromy of f
15 vfiltration(f[,opt]);      V-filtration on H''/H', singularity spectrum of f
16 spectrum(f);               singularity spectrum of f
17 endfilt(poly f,list V);    endomorphism filtration of filtration V
18 spprint(list S);           print spectrum S
19 spadd(list S1,list S2);    sum of spectra S1 and S2
20 spsub(list S1,list S2);    difference of spectra S1 and S2
21 spmul(list S,int k);       product of spectrum S and integer k
22 spmul(list S,intvec k);    linear combination of spectra S with coefficients k
23 spissemicont(list S[,opt]);        test spectrum S for semicontinuity
24 spsemicont(list S0,list S[,opt]);  relative semicontinuity of spectra S0 and S
25 spmilnor(list S);          milnor number of spectrum S
26 spgeomgenus(list S);       geometrical genus of spectrum S
27 spgamma(list S);           gamma invariant of spectrum S
28
29SEE ALSO: mondromy_lib, spectrum_lib
30
31KEYWORDS: singularities; Gauss-Manin connection; monodromy; spectrum;
32          Brieskorn lattice; Hodge filtration; V-filtration
33";
34
35LIB "linalg.lib";
36
37///////////////////////////////////////////////////////////////////////////////
38
39static proc maxintdif(ideal e)
40{
41  dbprint(printlevel-voice+2,"//gaussman::maxintdif");
42  int i,j,id;
43  int mid=0;
44  for(i=ncols(e);i>=1;i--)
45  {
46    for(j=i-1;j>=1;j--)
47    {
48      id=int(e[i]-e[j]);
49      if(id<0)
50      {
51        id=-id;
52      }
53      if(id>mid)
54      {
55        mid=id;
56      }
57    }
58  }
59  return(mid);
60}
61///////////////////////////////////////////////////////////////////////////////
62
63static proc maxorddif(matrix H)
64{
65  dbprint(printlevel-voice+2,"//gaussman::maxorddif");
66  int i,j,d;
67  int d0,d1=-1,-1;
68  for(i=nrows(H);i>=1;i--)
69  {
70    for(j=ncols(H);j>=1;j--)
71    {
72      d=ord(H[i,j]);
73      if(d>=0)
74      {
75        if(d0<0||d<d0)
76        {
77          d0=d;
78        }
79        if(d1<0||d>d1)
80        {
81          d1=d;
82        }
83      }
84    }
85  }
86  return(d1-d0);
87}
88///////////////////////////////////////////////////////////////////////////////
89
90proc monodromy(poly f,list #)
91"USAGE:    monodromy(f[,opt]); poly f, int opt
92ASSUME:   basering has characteristic 0 and local ordering,
93          f has isolated singularity at 0
94RETURN:
95@format
96if opt==0:
97  matrix M: exp(-2*pi*i*M) is a monodromy matrix of f
98if opt==1:
99  ideal e: exp(-2*pi*i*e) are the eigenvalues of the monodromy of f
100default: opt=1
101@end format
102SEE ALSO: mondromy_lib
103KEYWORDS: singularities; Gauss-Manin connection; monodromy
104EXAMPLE:  example monodromy; shows an example
105"
106{
107  if(charstr(basering)!="0")
108  {
109    ERROR("characteristic 0 expected");
110  }
111  int n=nvars(basering)-1;
112  for(int i=n+1;i>=1;i--)
113  {
114    if(var(i)>1)
115    {
116      ERROR("local ordering expected");
117    }
118  }
119  ideal J=jacob(f);
120  ideal sJ=std(J);
121  if(vdim(sJ)<=0)
122  {
123    if(vdim(sJ)==0)
124    {
125      ERROR("singularity at 0 expected");
126    }
127    else
128    {
129      ERROR("isolated singularity at 0 expected");
130    }
131  }
132  int opt=1;
133  if(size(#)>0)
134  {
135    if(typeof(#[1])=="int")
136    {
137      opt=#[1];
138    }
139  }
140
141  ideal m=kbase(sJ);
142  int mu,modm=ncols(m),maxorddif(m);
143
144  ideal w=f*m;
145  matrix U=freemodule(mu);
146  matrix A0[mu][mu],A,C,D;
147  list l;
148  module H,dH=freemodule(mu),freemodule(mu);
149  module H0;
150  int sdH=1;
151  int k=-1;
152  int j,K,N,mide;
153
154  while(k<K||sdH>0)
155  {
156    k++;
157    dbprint(printlevel-voice+2,"//gaussman::monodromy: k="+string(k));
158
159    dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
160    C=coeffs(reduce(w,U,sJ),m);
161    A0=A0+C*var(1)^k;
162
163    if(sdH>0)
164    {
165      H0=H;
166      dbprint(printlevel-voice+2,"//gaussman::monodromy: compute dH");
167      dH=jet(module(A0*dH+var(1)^2*diff(matrix(dH),var(1))),k);
168      H=H*var(1)+dH;
169
170      dbprint(printlevel-voice+2,"//gaussman::monodromy: test dH==0");
171      sdH=size(reduce(H,std(H0*var(1))));
172      if(sdH>0)
173      {
174        A0=A0-var(1);
175      }
176      else
177      {
178        dbprint(printlevel-voice+2,
179          "//gaussman::monodromy: compute basis of saturation");
180        H=minbase(H0);
181        int modH=maxorddif(H);
182        K=modH+1;
183      }
184    }
185
186    if(k==K&&sdH==0)
187    {
188      N=k-modH;
189      dbprint(printlevel-voice+2,
190        "//gaussman::monodromy: compute A on saturation");
191      l=division(H*var(1),A0*H+var(1)^2*diff(matrix(H),var(1)));
192      A=jet(l[1],l[2],N-1);
193      if(mide==0)
194      {
195        dbprint(printlevel-voice+2,
196          "//gaussman::monodromy: compute eigenvalues e and"+
197          "multiplicities b of A");
198        l=system("eigenval",jet(A,0));
199        ideal e=l[1];
200        intvec b=l[2];
201        dbprint(printlevel-voice+2,"//gaussman::monodromy: e="+string(e));
202        dbprint(printlevel-voice+2,"//gaussman::monodromy: b="+string(b));
203        if(opt==1)
204        {
205          mide=maxintdif(e);
206          K=K+mide;
207        }
208      }
209    }
210
211    if(k<K||sdH>0)
212    {
213      dbprint(printlevel-voice+2,"//gaussman::monodromy: divide by J");
214      l=division(J,ideal(matrix(w)-matrix(m)*C*U));
215      D=l[1];
216
217      dbprint(printlevel-voice+2,"//gaussman::monodromy: compute w/U");
218      for(j=mu;j>=1;j--)
219      {
220        if(l[2][j,j]!=0)
221        {
222          dbprint(printlevel-voice+2,
223            "//gaussman::monodromy: compute U["+string(j)+"]");
224          U[j,j]=U[j,j]*l[2][j,j];
225        }
226        dbprint(printlevel-voice+2,
227          "//gaussman::monodromy: compute w["+string(j)+"]");
228        w[j]=0;
229        for(i=n+1;i>=1;i--)
230        {
231          w[j]=w[j]+U[j,j]*diff(D[i,j],var(i))-diff(U[j,j],var(i))*D[i,j];
232        }
233      }
234      U=U*U;
235    }
236  }
237
238  if(mide>0)
239  {
240    intvec ide;
241    ide[mu]=0;
242    module dU;
243    matrix A0e;
244    for(i=ncols(e);i>=1;i--)
245    {
246      for(j=i-1;j>=1;j--)
247      {
248        k=int(e[j]-e[i]);
249        if(k>ide[i])
250        {
251          ide[i]=k;
252        }
253        if(-k>ide[j])
254        {
255          ide[j]=-k;
256        }
257      }
258    }
259    for(j,k=ncols(e),mu;j>=1;j--)
260    {
261      for(i=b[j];i>=1;i--)
262      {
263        ide[k]=ide[j];
264        k--;
265      }
266    }
267  }
268  while(mide>0)
269  {
270    dbprint(printlevel-voice+2,"//gaussman::monodromy: mide="+string(mide));
271
272    U=0;
273    A0=jet(A,0);
274    for(i=ncols(e);i>=1;i--)
275    {
276      A0e=freemodule(mu);
277      for(j=n;j>=0;j--) // Potenzen von Matrizen?
278      {
279        A0e=A0e*(A0-e[i]);
280      }
281      dU=syz(A0e);
282      U=dU+U;
283    }
284    A=division(U,A*U)[1];
285
286    for(i=mu;i>=1;i--)
287    {
288      for(j=mu;j>=1;j--)
289      {
290        if(ide[i]==0&&ide[j]>0)
291        {
292          A[i,j]=A[i,j]*var(1);
293        }
294        else
295        {
296          if(ide[i]>0&&ide[j]==0)
297          {
298            A[i,j]=A[i,j]/var(1);
299          }
300        }
301      }
302    }
303    for(i=mu;i>=1;i--)
304    {
305      if(ide[i]>0)
306      {
307        A[i,i]=A[i,i]+1;
308        e[i]=e[i]+1;
309        ide[i]=ide[i]-1;
310      }
311    }
312    mide--;
313  }
314
315  if(opt==1)
316  {
317    return(jet(A,0));
318  }
319  else
320  {
321    return(e);
322  }
323}
324example
325{ "EXAMPLE:"; echo=2;
326  ring R=0,(x,y),ds;
327  poly f=x5+x2y2+y5;
328  print(monodromy(f));
329}
330///////////////////////////////////////////////////////////////////////////////
331
332proc vfiltration(poly f,list #)
333"USAGE:    vfiltration(f[,opt]); poly f, int opt
334ASSUME:   basering has characteristic 0 and local ordering,
335          f has isolated singularity at 0
336RETURN:
337@format
338list V: V-filtration of f on H''/H'
339if opt==0 or opt==1:
340  ideal V[1]: spectral numbers in increasing order
341  intvec V[2]:
342    int V[2][i]: multiplicity of spectral number V[1][i]
343if opt==1:
344  list V[3]:
345    module V[3][i]: vector space basis of V[1][i]-th graded part
346                    in terms of V[4]
347  ideal V[4]: monomial vector space basis
348default: opt=1
349@end format
350NOTE:     H' and H'' denote the Brieskorn lattices
351SEE ALSO: spectrum_lib
352KEYWORDS: singularities; Gauss-Manin connection;
353          Brieskorn lattice; Hodge filtration; V-filtration; spectrum
354EXAMPLE:  example vfiltration; shows an example
355"
356{
357  if(charstr(basering)!="0")
358  {
359    ERROR("characteristic 0 expected");
360  }
361  int n=nvars(basering)-1;
362  for(int i=n+1;i>=1;i--)
363  {
364    if(var(i)>1)
365    {
366      ERROR("local ordering expected");
367    }
368  }
369  ideal J=jacob(f);
370  ideal sJ=std(J);
371  if(vdim(sJ)<=0)
372  {
373    if(vdim(sJ)==0)
374    {
375      ERROR("singularity at 0 expected");
376    }
377    else
378    {
379      ERROR("isolated singularity at 0 expected");
380    }
381  }
382  int opt=1;
383  if(size(#)>0)
384  {
385    if(typeof(#[1])=="int")
386    {
387      opt=#[1];
388    }
389  }
390
391  ideal m=kbase(sJ);
392  int mu,modm=ncols(m),maxorddif(m);
393
394  ideal w=f*m;
395  matrix U=freemodule(mu);
396  matrix A[mu][mu],C,D;
397  list l;
398  module H,dH=freemodule(mu),freemodule(mu);
399  module H0;
400  int sdH=1;
401  int k=-1;
402  int j,K;
403
404  while(k<K||sdH>0)
405  {
406    k++;
407    dbprint(printlevel-voice+2,"//gaussman::vfiltration: k="+string(k));
408
409    dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute C");
410    C=coeffs(reduce(w,U,sJ),m);
411    A=A+C*var(1)^k;
412
413    if(sdH>0)
414    {
415      H0=H;
416      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute dH");
417      dH=jet(module(A*dH+var(1)^2*diff(matrix(dH),var(1))),k);
418      H=H*var(1)+dH;
419
420      dbprint(printlevel-voice+2,"//gaussman::vfiltration: test dH==0");
421      sdH=size(reduce(H,std(H0*var(1))));
422      if(sdH>0)
423      {
424        A=A-var(1);
425      }
426      else
427      {
428        dbprint(printlevel-voice+2,
429          "//gaussman::vfiltration: compute basis of saturation");
430        H=minbase(H0);
431        int modH=maxorddif(H);
432        K=modH+n+1;
433        H0=freemodule(mu)*var(1)^k;
434      }
435    }
436
437    if(k<K||sdH>0)
438    {
439      dbprint(printlevel-voice+2,"//gaussman::vfiltration: divide by J");
440      l=division(J,ideal(matrix(w)-matrix(m)*C*U));
441      D=l[1];
442
443      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute w/U");
444      for(j=mu;j>=1;j--)
445      {
446        if(l[2][j,j]!=0)
447        {
448          dbprint(printlevel-voice+2,
449            "//gaussman::vfiltration: compute U["+string(j)+"]");
450          U[j,j]=U[j,j]*l[2][j,j];
451        }
452        dbprint(printlevel-voice+2,
453          "//gaussman::vfiltration: compute w["+string(j)+"]");
454        w[j]=0;
455        for(i=n+1;i>=1;i--)
456        {
457          w[j]=w[j]+U[j,j]*diff(D[i,j],var(i))-diff(U[j,j],var(i))*D[i,j];
458        }
459      }
460      U=U*U;
461    }
462  }
463  int N=k-modH;
464
465  dbprint(printlevel-voice+2,
466    "//gaussman::vfiltration: transform H0 to saturation");
467  l=division(H,H0);
468  H0=jet(l[1],l[2],N-1);
469
470  dbprint(printlevel-voice+2,
471    "//gaussman::vfiltration: compute H0 as vector space V0");
472  dbprint(printlevel-voice+2,
473    "//gaussman::vfiltration: compute H1 as vector space V1");
474  poly p;
475  int i0,j0,i1,j1;
476  matrix V0[mu*N][mu*N];
477  matrix V1[mu*N][mu*(N-1)];
478  for(i0=mu;i0>=1;i0--)
479  {
480    for(i1=mu;i1>=1;i1--)
481    {
482      p=H0[i1,i0];
483      while(p!=0)
484      {
485        j1=leadexp(p)[1];
486        for(j0=N-j1-1;j0>=0;j0--)
487        {
488          V0[i1+(j1+j0)*mu,i0+j0*mu]=V0[i1+(j1+j0)*mu,i0+j0*mu]+leadcoef(p);
489          if(j1+j0+1<N)
490          {
491            V1[i1+(j1+j0+1)*mu,i0+j0*mu]=
492            V1[i1+(j1+j0+1)*mu,i0+j0*mu]+leadcoef(p);
493          }
494        }
495        p=p-lead(p);
496      }
497    }
498  }
499
500  dbprint(printlevel-voice+2,
501    "//gaussman::vfiltration: compute A on saturation");
502  l=division(H*var(1),A*H+var(1)^2*diff(matrix(H),var(1)));
503  A=jet(l[1],l[2],N-1);
504
505  dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute matrix M of A");
506  matrix M[mu*N][mu*N];
507  for(i0=mu;i0>=1;i0--)
508  {
509    for(i1=mu;i1>=1;i1--)
510    {
511      p=A[i1,i0];
512      while(p!=0)
513      {
514        j1=leadexp(p)[1];
515        for(j0=N-j1-1;j0>=0;j0--)
516        {
517          M[i1+(j0+j1)*mu,i0+j0*mu]=leadcoef(p);
518        }
519        p=p-lead(p);
520      }
521    }
522  }
523  for(i0=mu;i0>=1;i0--)
524  {
525    for(j0=N-1;j0>=0;j0--)
526    {
527      M[i0+j0*mu,i0+j0*mu]=M[i0+j0*mu,i0+j0*mu]+j0;
528    }
529  }
530
531  dbprint(printlevel-voice+2,
532    "//gaussman::vfiltration: compute eigenvalues eA of A");
533  ideal eA=system("eigenval",jet(A,0))[1];
534  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eA="+string(eA));
535
536  dbprint(printlevel-voice+2,
537    "//gaussman::vfiltration: compute eigenvalues eM of M");
538  ideal eM;
539  k=0;
540  intvec u;
541  for(i=N;i>=1;i--)
542  {
543    u[i]=1;
544  }
545  i0=1;
546  while(u[N]<=ncols(eA))
547  {
548    for(i,i1=i0+1,i0;i<=N;i++)
549    {
550      if(eA[u[i]]+i<eA[u[i1]]+i1)
551      {
552        i1=i;
553      }
554    }
555    k++;
556    eM[k]=eA[u[i1]]+i1-1;
557    u[i1]=u[i1]+1;
558    if(u[i1]>ncols(eA))
559    {
560      i0=i1+1;
561    }
562  }
563  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eM="+string(eM));
564
565  dbprint(printlevel-voice+2,
566    "//gaussman::vfiltration: compute V-filtration on H0/H1");
567  ideal a;
568  k=0;
569  list V;
570  matrix Me;
571  V[ncols(eM)+1]=std(V1);
572  intvec d;
573  if(opt==0)
574  {
575    for(i=ncols(eM);number(eM[i])-1>number(n-1)/2;i--)
576    {
577      Me=freemodule(mu*N);
578      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
579      {
580        Me=Me*(M-eM[i]);
581      }
582
583      dbprint(printlevel-voice+2,
584        "//gaussman::vfiltration: compute V["+string(i)+"]");
585      V1=module(V1)+syz(Me);
586      V[i]=std(intersect(V1,V0));
587
588      if(size(V[i])>size(V[i+1]))
589      {
590        k++;
591        a[k]=eM[i]-1;
592        d[k]=size(V[i])-size(V[i+1]);
593      }
594    }
595
596    dbprint(printlevel-voice+2,
597      "//gaussman::vfiltration: symmetry index found");
598    j=k;
599
600    if(number(eM[i])-1==number(n-1)/2)
601    {
602      Me=freemodule(mu*N);
603      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
604      {
605        Me=Me*(M-eM[i]);
606      }
607
608      dbprint(printlevel-voice+2,
609        "//gaussman::vfiltration: compute V["+string(i)+"]");
610      V1=module(V1)+syz(Me);
611      V[i]=std(intersect(V1,V0));
612
613      if(size(V[i])>size(V[i+1]))
614      {
615        k++;
616        a[k]=eM[i]-1;
617        d[k]=size(V[i])-size(V[i+1]);
618      }
619    }
620
621    dbprint(printlevel-voice+2,"//gaussman::vfiltration: apply symmetry");
622    while(j>=1)
623    {
624      k++;
625      a[k]=a[j];
626      a[j]=n-1-a[k];
627      d[k]=d[j];
628      j--;
629    }
630
631    return(list(a,d));
632  }
633  else
634  {
635    list v;
636    j=-1;
637    for(i=ncols(eM);i>=1;i--)
638    {
639      Me=freemodule(mu*N);
640      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
641      {
642        Me=Me*(M-eM[i]);
643      }
644
645      dbprint(printlevel-voice+2,
646        "//gaussman::vfiltration: compute V["+string(i)+"]");
647      V1=module(V1)+syz(Me);
648      V[i]=std(intersect(V1,V0));
649
650      if(size(V[i])>size(V[i+1]))
651      {
652        if(number(eM[i]-1)>=number(n-1)/2)
653        {
654          k++;
655          a[k]=eM[i]-1;
656          dbprint(printlevel-voice+2,
657            "//gaussman::vfiltration: transform to V0");
658          v[k]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
659        }
660        else
661        {
662          if(j<0)
663          {
664            if(a[k]==number(n-1)/2)
665            {
666              j=k-1;
667            }
668            else
669            {
670              j=k;
671            }
672          }
673          k++;
674          a[k]=a[j];
675          a[j]=eM[i]-1;
676          v[k]=v[j];
677          dbprint(printlevel-voice+2,
678            "//gaussman::vfiltration: transform to V0");
679          v[j]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
680          j--;
681        }
682      }
683    }
684
685    dbprint(printlevel-voice+2,
686      "//gaussman::vfiltration: compute graded parts");
687    option(redSB);
688    for(k=1;k<size(v);k++)
689    {
690      v[k]=std(reduce(v[k],std(v[k+1])));
691      d[k]=size(v[k]);
692    }
693    v[k]=std(v[k]);
694    d[k]=size(v[k]);
695
696    return(list(a,d,v,m));
697  }
698}
699example
700{ "EXAMPLE:"; echo=2;
701  ring R=0,(x,y),ds;
702  poly f=x5+x2y2+y5;
703  vfiltration(f);
704}
705///////////////////////////////////////////////////////////////////////////////
706
707proc spectrum(poly f)
708"USAGE:    spectrum(f); poly f
709ASSUME:   basering has characteristic 0 and local ordering,
710          f has isolated singularity at 0
711RETURN:
712@format
713list S: singularity spectrum of f
714  ideal S[1]: spectral numbers in increasing order
715  intvec S[2]:
716    int S[2][i]: multiplicity of spectral number S[1][i]
717@end format
718SEE ALSO: spectrum_lib
719KEYWORDS: singularities; Gauss-Manin connection; spectrum
720EXAMPLE:  example spectrum; shows an example
721"
722{
723  return(vfiltration(f,0));
724}
725example
726{ "EXAMPLE:"; echo=2;
727  ring R=0,(x,y),ds;
728  poly f=x5+x2y2+y5;
729  spprint(spectrum(f));
730}
731///////////////////////////////////////////////////////////////////////////////
732
733proc endfilt(poly f,list V)
734"USAGE:   endfilt(f,V); poly f, list V
735ASSUME:  basering has characteristic 0 and local ordering,
736         f has isolated singularity at 0
737RETURN:
738@format
739list V1: endomorphim filtration of V on the Jacobian algebra of f
740  ideal V1[1]: spectral numbers in increasing order
741  intvec V1[2]:
742    int V1[2][i]: multiplicity of spectral number V1[1][i]
743  list V1[3]:
744    module V1[3][i]: vector space basis of the V1[1][i]-th graded part
745                     in terms of V1[4]
746  ideal V1[4]: monomial vector space basis
747@end format
748SEE ALSO: spectrum_lib
749KEYWORDS: singularities; Gauss-Manin connection; spectrum;
750          Brieskorn lattice; Hodge filtration; V-filtration
751EXAMPLE: example endfilt; shows an example
752"
753{
754  if(charstr(basering)!="0")
755  {
756    ERROR("characteristic 0 expected");
757  }
758  int n=nvars(basering)-1;
759  for(int i=n+1;i>=1;i--)
760  {
761    if(var(i)>1)
762    {
763      ERROR("local ordering expected");
764    }
765  }
766  ideal sJ=std(jacob(f));
767  if(vdim(sJ)<=0)
768  {
769    if(vdim(sJ)==0)
770    {
771      ERROR("singularity at 0 expected");
772    }
773    else
774    {
775      ERROR("isolated singularity at 0 expected");
776    }
777  }
778
779  def a,d,v,m=V[1..4];
780  int mu=ncols(m);
781
782  module V0=v[1];
783  for(i=2;i<=size(v);i++)
784  {
785    V0=V0,v[i];
786  }
787
788  dbprint(printlevel-voice+2,
789    "//gaussman::endfilt: compute multiplication in Jacobian algebra");
790  list M;
791  matrix U=freemodule(ncols(m));
792  for(i=ncols(m);i>=1;i--)
793  {
794    M[i]=lift(V0,coeffs(reduce(m[i]*m,U,sJ),m)*V0);
795  }
796
797  int j,k,i0,j0,i1,j1;
798  number b0=number(a[1]-a[ncols(a)]);
799  number b1,b2;
800  matrix M0;
801  module L;
802  list v0=freemodule(ncols(m));
803  ideal a0=b0;
804
805  while(b0<number(a[ncols(a)]-a[1]))
806  {
807    dbprint(printlevel-voice+2,
808      "//gaussman::endfilt: find next possible index");
809    b1=number(a[ncols(a)]-a[1]);
810    for(j=ncols(a);j>=1;j--)
811    {
812      for(i=ncols(a);i>=1;i--)
813      {
814        b2=number(a[i]-a[j]);
815        if(b2>b0&&b2<b1)
816        {
817          b1=b2;
818        }
819        else
820        {
821          if(b2<=b0)
822          {
823            i=0;
824          }
825        }
826      }
827    }
828    b0=b1;
829
830    list l=ideal();
831    for(k=ncols(m);k>=2;k--)
832    {
833      l=l+list(ideal());
834    }
835
836    dbprint(printlevel-voice+2,
837      "//gaussman::endfilt: collect conditions for V1["+string(b0)+"]");
838    j=ncols(a);
839    j0=mu;
840    while(j>=1)
841    {
842      i0=1;
843      i=1;
844      while(i<ncols(a)&&a[i]<a[j]+b0)
845      {
846        i0=i0+d[i];
847        i++;
848      }
849      if(a[i]<a[j]+b0)
850      {
851        i0=i0+d[i];
852        i++;
853      }
854      for(k=1;k<=ncols(m);k++)
855      {
856        M0=M[k];
857        if(i0>1)
858        {
859          l[k]=l[k],M0[1..i0-1,j0-d[j]+1..j0];
860        }
861      }
862      j0=j0-d[j];
863      j--;
864    }
865
866    dbprint(printlevel-voice+2,
867      "//gaussman::endfilt: compose condition matrix");
868    L=transpose(module(l[1]));
869    for(k=2;k<=ncols(m);k++)
870    {
871      L=L,transpose(module(l[k]));
872    }
873
874    dbprint(printlevel-voice+2,
875      "//gaussman::endfilt: compute kernel of condition matrix");
876    v0=v0+list(syz(L));
877    a0=a0,b0;
878  }
879
880  dbprint(printlevel-voice+2,"//gaussman::endfilt: compute graded parts");
881  option(redSB);
882  for(i=1;i<size(v0);i++)
883  {
884    v0[i+1]=std(v0[i+1]);
885    v0[i]=std(reduce(v0[i],v0[i+1]));
886  }
887
888  dbprint(printlevel-voice+2,
889    "//gaussman::endfilt: remove trivial graded parts");
890  i=1;
891  while(size(v0[i])==0)
892  {
893    i++;
894  }
895  list v1=v0[i];
896  intvec d1=size(v0[i]);
897  ideal a1=a0[i];
898  i++;
899  while(i<=size(v0))
900  {
901    if(size(v0[i])>0)
902    {
903      v1=v1+list(v0[i]);
904      d1=d1,size(v0[i]);
905      a1=a1,a0[i];
906    }
907    i++;
908  }
909  return(list(a1,d1,v1,m));
910}
911example
912{ "EXAMPLE:"; echo=2;
913  ring R=0,(x,y),ds;
914  poly f=x5+x2y2+y5;
915  endfilt(f,vfiltration(f));
916}
917///////////////////////////////////////////////////////////////////////////////
918
919proc spprint(list S)
920"USAGE:   spprint(S); list S
921RETURN:  string: spectrum S
922EXAMPLE: example spprint; shows an example
923"
924{
925  string s;
926  for(int i=1;i<size(S[2]);i++)
927  {
928    s=s+"("+string(S[1][i])+","+string(S[2][i])+"),";
929  }
930  s=s+"("+string(S[1][i])+","+string(S[2][i])+")";
931  return(s);
932}
933example
934{ "EXAMPLE:"; echo=2;
935  ring R=0,(x,y),ds;
936  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
937  spprint(S);
938}
939///////////////////////////////////////////////////////////////////////////////
940
941proc spadd(list S1,list S2)
942"USAGE:   spadd(S1,S2); list S1,S2
943RETURN:  list: sum of spectra S1 and S2
944EXAMPLE: example spadd; shows an example
945"
946{
947  ideal s;
948  intvec m;
949  int i,i1,i2=1,1,1;
950  while(i1<=size(S1[2])||i2<=size(S2[2]))
951  {
952    if(i1<=size(S1[2]))
953    {
954      if(i2<=size(S2[2]))
955      {
956        if(number(S1[1][i1])<number(S2[1][i2]))
957        {
958          s[i]=S1[1][i1];
959          m[i]=S1[2][i1];
960          i++;
961          i1++;
962        }
963        else
964        {
965          if(number(S1[1][i1])>number(S2[1][i2]))
966          {
967            s[i]=S2[1][i2];
968            m[i]=S2[2][i2];
969            i++;
970            i2++;
971          }
972          else
973          {
974            if(S1[2][i1]+S2[2][i2]!=0)
975            {
976              s[i]=S1[1][i1];
977              m[i]=S1[2][i1]+S2[2][i2];
978              i++;
979            }
980            i1++;
981            i2++;
982          }
983        }
984      }
985      else
986      {
987        s[i]=S1[1][i1];
988        m[i]=S1[2][i1];
989        i++;
990        i1++;
991      }
992    }
993    else
994    {
995      s[i]=S2[1][i2];
996      m[i]=S2[2][i2];
997      i++;
998      i2++;
999    }
1000  }
1001  return(list(s,m));
1002}
1003example
1004{ "EXAMPLE:"; echo=2;
1005  ring R=0,(x,y),ds;
1006  list S1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1007  spprint(S1);
1008  list S2=list(ideal(-1/6,1/6),intvec(1,1));
1009  spprint(S2);
1010  spprint(spadd(S1,S2));
1011}
1012///////////////////////////////////////////////////////////////////////////////
1013
1014proc spsub(list S1,list S2)
1015"USAGE:   spsub(S1,S2); list S1,S2
1016RETURN:  list: difference of spectra S1 and S2
1017EXAMPLE: example spsub; shows an example
1018"
1019{
1020  return(spadd(S1,spmul(S2,-1)));
1021}
1022example
1023{ "EXAMPLE:"; echo=2;
1024  ring R=0,(x,y),ds;
1025  list S1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1026  spprint(S1);
1027  list S2=list(ideal(-1/6,1/6),intvec(1,1));
1028  spprint(S2);
1029  spprint(spsub(S1,S2));
1030}
1031///////////////////////////////////////////////////////////////////////////////
1032
1033proc spmul(list #)
1034"USAGE:
1035@format
10361) spmul(S,k); list S, int k
10372) spmul(S,k); list S, intvec k
1038@end format
1039RETURN:
1040@format
10411) list: product of spectrum S and integer k
10422) list: linear combination of spectra S with coefficients k
1043@end format
1044EXAMPLE: example spmul; shows an example
1045"
1046{
1047  if(size(#)==2)
1048  {
1049    if(typeof(#[1])=="list")
1050    {
1051      if(typeof(#[2])=="int")
1052      {
1053        return(list(#[1][1],#[1][2]*#[2]));
1054      }
1055      if(typeof(#[2])=="intvec")
1056      {
1057        list S0=list(ideal(),intvec(0));
1058        for(int i=size(#[2]);i>=1;i--)
1059        {
1060          S0=spadd(S0,spmul(#[1][i],#[2][i]));
1061        }
1062        return(S0);
1063      }
1064    }
1065  }
1066  return(list(ideal(),intvec(0)));
1067}
1068example
1069{ "EXAMPLE:"; echo=2;
1070  ring R=0,(x,y),ds;
1071  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1072  spprint(S);
1073  spprint(spmul(S,2));
1074  list S1=list(ideal(-1/6,1/6),intvec(1,1));
1075  spprint(S1);
1076  list S2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
1077  spprint(S2);
1078  spprint(spmul(list(S1,S2),intvec(1,2)));
1079}
1080///////////////////////////////////////////////////////////////////////////////
1081
1082proc spissemicont(list S,list #)
1083"USAGE:   spissemicont(S[,opt]); list S, int opt
1084RETURN:
1085@format
1086int k=
1087if opt==0:
1088  1, if sum of spectrum S over all intervals [a,a+1) is positive
1089  0, if sum of spectrum S over some interval [a,a+1) is negative
1090if opt==1:
1091  1, if sum of spectrum S over all intervals [a,a+1) and (a,a+1) is positive
1092  0, if sum of spectrum S over some interval [a,a+1) or (a,a+1) is negative
1093default: opt=0
1094@end format
1095EXAMPLE: example spissemicont; shows an example
1096"
1097{
1098  int opt=0;
1099  if(size(#)>0)
1100  {
1101    if(typeof(#[1])=="int")
1102    {
1103      opt=1;
1104    }
1105  }
1106  int i,j,k=1,1,0;
1107  while(j<=size(S[2]))
1108  {
1109    while(j+1<=size(S[2])&&S[1][j]<S[1][i]+1)
1110    {
1111      k=k+S[2][j];
1112      j++;
1113    }
1114    if(j==size(S[2])&&S[1][j]<S[1][i]+1)
1115    {
1116      k=k+S[2][j];
1117      j++;
1118    }
1119    if(k<0)
1120    {
1121      return(0);
1122    }
1123    k=k-S[2][i];
1124    if(k<0&&opt==1)
1125    {
1126      return(0);
1127    }
1128    i++;
1129  }
1130  return(1);
1131}
1132example
1133{ "EXAMPLE:"; echo=2;
1134  ring R=0,(x,y),ds;
1135  list S1=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1136  spprint(S1);
1137  list S2=list(ideal(-1/6,1/6),intvec(1,1));
1138  spprint(S2);
1139  spissemicont(spsub(S1,spmul(S2,5)));
1140  spissemicont(spsub(S1,spmul(S2,5)),1);
1141  spissemicont(spsub(S1,spmul(S2,6)));
1142}
1143///////////////////////////////////////////////////////////////////////////////
1144
1145proc spsemicont(list S0,list S,list #)
1146"USAGE:   spsemicont(S,k[,opt]); list S0, list S, int opt
1147RETURN:  list of intvecs l:
1148         spissemicont(sub(S0,spmul(S,k)),opt)==1 iff k<=l[i] for some i
1149NOTE:    if the spectra S occur with multiplicities k in a deformation
1150         of the [quasihomogeneous] spectrum S0 then
1151         spissemicont(sub(S0,spmul(S,k))[,1])==1
1152EXAMPLE: example spsemicont; shows an example
1153"
1154{
1155  list l,l0;
1156  int i,k;
1157  while(spissemicont(S0,#))
1158  {
1159    if(size(S)>1)
1160    {
1161      l0=spsemicont(S0,list(S[1..size(S)-1]));
1162      for(i=size(l0);i>=1;i--)
1163      {
1164        l0[i][size(S)]=k;
1165      }
1166      l=l+l0;
1167    }
1168    S0=spsub(S0,S[size(S)]);
1169    k++;
1170  }
1171  if(size(S)>1)
1172  {
1173    return(l);
1174  }
1175  else
1176  {
1177    return(list(intvec(k-1)));
1178  }
1179}
1180example
1181{ "EXAMPLE:"; echo=2;
1182  ring R=0,(x,y),ds;
1183  list S0=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1184  spprint(S0);
1185  list S1=list(ideal(-1/6,1/6),intvec(1,1));
1186  spprint(S1);
1187  list S2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
1188  spprint(S2);
1189  list S=S1,S2;
1190  list l=spsemicont(S0,S);
1191  l;
1192  spissemicont(spsub(S0,spmul(S,l[1])));
1193  spissemicont(spsub(S0,spmul(S,l[1]-1)));
1194  spissemicont(spsub(S0,spmul(S,l[1]+1)));
1195}
1196///////////////////////////////////////////////////////////////////////////////
1197
1198proc spmilnor(list S)
1199"USAGE:   spmilnor(S); list S
1200RETURN:  int: Milnor number of spectrum S
1201EXAMPLE: example spmilnor; shows an example
1202"
1203{
1204  return(sum(S[2]));
1205}
1206example
1207{ "EXAMPLE:"; echo=2;
1208  ring R=0,(x,y),ds;
1209  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1210  spprint(S);
1211  spmilnor(S);
1212}
1213///////////////////////////////////////////////////////////////////////////////
1214
1215proc spgeomgenus(list S)
1216"USAGE:   spgeomgenus(S); list S
1217RETURN:  int: geometrical genus of spectrum S
1218EXAMPLE: example spgeomgenus; shows an example
1219"
1220{
1221  int g=0;
1222  int i=1;
1223  while(i+1<=size(S[2])&&number(S[1][i])<=number(0))
1224  {
1225    g=g+S[2][i];
1226    i++;
1227  }
1228  if(i==size(S[2])&&number(S[1][i])<=number(0))
1229  {
1230    g=g+S[2][i];
1231  }
1232  return(g);
1233}
1234example
1235{ "EXAMPLE:"; echo=2;
1236  ring R=0,(x,y),ds;
1237  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1238  spprint(S);
1239  spgeomgenus(S);
1240}
1241///////////////////////////////////////////////////////////////////////////////
1242
1243proc spgamma(list S)
1244"USAGE:   spgamma(S); list S
1245RETURN:  number: gamma invariant of spectrum S
1246EXAMPLE: example spgamma; shows an example
1247"
1248{
1249  int i,j;
1250  number g=0;
1251  for(i=1;i<=ncols(S[1]);i++)
1252  {
1253    for(j=1;j<=S[2][i];j++)
1254    {
1255      g=g+(number(S[1][i])-number(nvars(basering)-2)/2)^2;
1256    }
1257  }
1258  g=-g/4+sum(S[2])*number(S[1][ncols(S[1])]-S[1][1])/48;
1259  return(g);
1260}
1261example
1262{ "EXAMPLE:"; echo=2;
1263  ring R=0,(x,y),ds;
1264  list S=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1265  spprint(S);
1266  spgamma(S);
1267}
1268///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.