source: git/Singular/LIB/gaussman.lib @ 38f6b33

spielwiese
Last change on this file since 38f6b33 was 38f6b33, checked in by Mathias Schulze <mschulze@…>, 23 years ago
*mschulze: hopefully no more redundant elements in return list of spsemicont git-svn-id: file:///usr/local/Singular/svn/trunk@5401 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.8 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: gaussman.lib,v 1.37 2001-04-26 12:05:15 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  intvec V[1]: numerators of spectral numbers
341  intvec V[2]: denominators of spectral numbers
342  intvec V[3]:
343    int V[3][i]: multiplicity of spectral number V[1][i]/V[2][i]
344if opt==1:
345  list V[4]:
346    module V[4][i]: vector space basis of V[1][i]/V[2][i]-th graded part
347                    in terms of V[5]
348  ideal V[5]: monomial vector space basis
349default: opt=1
350@end format
351NOTE:     H' and H'' denote the Brieskorn lattices
352SEE ALSO: spectrum_lib
353KEYWORDS: singularities; Gauss-Manin connection;
354          Brieskorn lattice; Hodge filtration; V-filtration; spectrum
355EXAMPLE:  example vfiltration; shows an example
356"
357{
358  if(charstr(basering)!="0")
359  {
360    ERROR("characteristic 0 expected");
361  }
362  int n=nvars(basering)-1;
363  for(int i=n+1;i>=1;i--)
364  {
365    if(var(i)>1)
366    {
367      ERROR("local ordering expected");
368    }
369  }
370  ideal J=jacob(f);
371  ideal sJ=std(J);
372  if(vdim(sJ)<=0)
373  {
374    if(vdim(sJ)==0)
375    {
376      ERROR("singularity at 0 expected");
377    }
378    else
379    {
380      ERROR("isolated singularity at 0 expected");
381    }
382  }
383  int opt=1;
384  if(size(#)>0)
385  {
386    if(typeof(#[1])=="int")
387    {
388      opt=#[1];
389    }
390  }
391
392  ideal m=kbase(sJ);
393  int mu,modm=ncols(m),maxorddif(m);
394
395  ideal w=f*m;
396  matrix U=freemodule(mu);
397  matrix A[mu][mu],C,D;
398  list l;
399  module H,dH=freemodule(mu),freemodule(mu);
400  module H0;
401  int sdH=1;
402  int k=-1;
403  int j,K;
404
405  while(k<K||sdH>0)
406  {
407    k++;
408    dbprint(printlevel-voice+2,"//gaussman::vfiltration: k="+string(k));
409
410    dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute C");
411    C=coeffs(reduce(w,U,sJ),m);
412    A=A+C*var(1)^k;
413
414    if(sdH>0)
415    {
416      H0=H;
417      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute dH");
418      dH=jet(module(A*dH+var(1)^2*diff(matrix(dH),var(1))),k);
419      H=H*var(1)+dH;
420
421      dbprint(printlevel-voice+2,"//gaussman::vfiltration: test dH==0");
422      sdH=size(reduce(H,std(H0*var(1))));
423      if(sdH>0)
424      {
425        A=A-var(1);
426      }
427      else
428      {
429        dbprint(printlevel-voice+2,
430          "//gaussman::vfiltration: compute basis of saturation");
431        H=minbase(H0);
432        int modH=maxorddif(H);
433        K=modH+n+1;
434        H0=freemodule(mu)*var(1)^k;
435      }
436    }
437
438    if(k<K||sdH>0)
439    {
440      dbprint(printlevel-voice+2,"//gaussman::vfiltration: divide by J");
441      l=division(J,ideal(matrix(w)-matrix(m)*C*U));
442      D=l[1];
443
444      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute w/U");
445      for(j=mu;j>=1;j--)
446      {
447        if(l[2][j,j]!=0)
448        {
449          dbprint(printlevel-voice+2,
450            "//gaussman::vfiltration: compute U["+string(j)+"]");
451          U[j,j]=U[j,j]*l[2][j,j];
452        }
453        dbprint(printlevel-voice+2,
454          "//gaussman::vfiltration: compute w["+string(j)+"]");
455        w[j]=0;
456        for(i=n+1;i>=1;i--)
457        {
458          w[j]=w[j]+U[j,j]*diff(D[i,j],var(i))-diff(U[j,j],var(i))*D[i,j];
459        }
460      }
461      U=U*U;
462    }
463  }
464  int N=k-modH;
465
466  dbprint(printlevel-voice+2,
467    "//gaussman::vfiltration: transform H0 to saturation");
468  l=division(H,H0);
469  H0=jet(l[1],l[2],N-1);
470
471  dbprint(printlevel-voice+2,
472    "//gaussman::vfiltration: compute H0 as vector space V0");
473  dbprint(printlevel-voice+2,
474    "//gaussman::vfiltration: compute H1 as vector space V1");
475  poly p;
476  int i0,j0,i1,j1;
477  matrix V0[mu*N][mu*N];
478  matrix V1[mu*N][mu*(N-1)];
479  for(i0=mu;i0>=1;i0--)
480  {
481    for(i1=mu;i1>=1;i1--)
482    {
483      p=H0[i1,i0];
484      while(p!=0)
485      {
486        j1=leadexp(p)[1];
487        for(j0=N-j1-1;j0>=0;j0--)
488        {
489          V0[i1+(j1+j0)*mu,i0+j0*mu]=V0[i1+(j1+j0)*mu,i0+j0*mu]+leadcoef(p);
490          if(j1+j0+1<N)
491          {
492            V1[i1+(j1+j0+1)*mu,i0+j0*mu]=
493            V1[i1+(j1+j0+1)*mu,i0+j0*mu]+leadcoef(p);
494          }
495        }
496        p=p-lead(p);
497      }
498    }
499  }
500
501  dbprint(printlevel-voice+2,
502    "//gaussman::vfiltration: compute A on saturation");
503  l=division(H*var(1),A*H+var(1)^2*diff(matrix(H),var(1)));
504  A=jet(l[1],l[2],N-1);
505
506  dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute matrix M of A");
507  matrix M[mu*N][mu*N];
508  for(i0=mu;i0>=1;i0--)
509  {
510    for(i1=mu;i1>=1;i1--)
511    {
512      p=A[i1,i0];
513      while(p!=0)
514      {
515        j1=leadexp(p)[1];
516        for(j0=N-j1-1;j0>=0;j0--)
517        {
518          M[i1+(j0+j1)*mu,i0+j0*mu]=leadcoef(p);
519        }
520        p=p-lead(p);
521      }
522    }
523  }
524  for(i0=mu;i0>=1;i0--)
525  {
526    for(j0=N-1;j0>=0;j0--)
527    {
528      M[i0+j0*mu,i0+j0*mu]=M[i0+j0*mu,i0+j0*mu]+j0;
529    }
530  }
531
532  dbprint(printlevel-voice+2,
533    "//gaussman::vfiltration: compute eigenvalues eA of A");
534  ideal eA=system("eigenval",jet(A,0))[1];
535  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eA="+string(eA));
536
537  dbprint(printlevel-voice+2,
538    "//gaussman::vfiltration: compute eigenvalues eM of M");
539  ideal eM;
540  k=0;
541  intvec u;
542  for(i=N;i>=1;i--)
543  {
544    u[i]=1;
545  }
546  i0=1;
547  while(u[N]<=ncols(eA))
548  {
549    for(i,i1=i0+1,i0;i<=N;i++)
550    {
551      if(eA[u[i]]+i<eA[u[i1]]+i1)
552      {
553        i1=i;
554      }
555    }
556    k++;
557    eM[k]=eA[u[i1]]+i1-1;
558    u[i1]=u[i1]+1;
559    if(u[i1]>ncols(eA))
560    {
561      i0=i1+1;
562    }
563  }
564  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eM="+string(eM));
565
566  dbprint(printlevel-voice+2,
567    "//gaussman::vfiltration: compute V-filtration on H0/H1");
568  ideal a;
569  k=0;
570  list V;
571  matrix Me;
572  V[ncols(eM)+1]=std(V1);
573  intvec d;
574  if(opt==0)
575  {
576    for(i=ncols(eM);number(eM[i])-1>number(n-1)/2;i--)
577    {
578      Me=freemodule(mu*N);
579      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
580      {
581        Me=Me*(M-eM[i]);
582      }
583
584      dbprint(printlevel-voice+2,
585        "//gaussman::vfiltration: compute V["+string(i)+"]");
586      V1=module(V1)+syz(Me);
587      V[i]=std(intersect(V1,V0));
588
589      if(size(V[i])>size(V[i+1]))
590      {
591        k++;
592        a[k]=eM[i]-1;
593        d[k]=size(V[i])-size(V[i+1]);
594      }
595    }
596
597    dbprint(printlevel-voice+2,
598      "//gaussman::vfiltration: symmetry index found");
599    j=k;
600
601    if(number(eM[i])-1==number(n-1)/2)
602    {
603      Me=freemodule(mu*N);
604      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
605      {
606        Me=Me*(M-eM[i]);
607      }
608
609      dbprint(printlevel-voice+2,
610        "//gaussman::vfiltration: compute V["+string(i)+"]");
611      V1=module(V1)+syz(Me);
612      V[i]=std(intersect(V1,V0));
613
614      if(size(V[i])>size(V[i+1]))
615      {
616        k++;
617        a[k]=eM[i]-1;
618        d[k]=size(V[i])-size(V[i+1]);
619      }
620    }
621
622    dbprint(printlevel-voice+2,"//gaussman::vfiltration: apply symmetry");
623    while(j>=1)
624    {
625      k++;
626      a[k]=a[j];
627      a[j]=n-1-a[k];
628      d[k]=d[j];
629      j--;
630    }
631
632    return(list(a,d));
633  }
634  else
635  {
636    list v;
637    j=-1;
638    for(i=ncols(eM);i>=1;i--)
639    {
640      Me=freemodule(mu*N);
641      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
642      {
643        Me=Me*(M-eM[i]);
644      }
645
646      dbprint(printlevel-voice+2,
647        "//gaussman::vfiltration: compute V["+string(i)+"]");
648      V1=module(V1)+syz(Me);
649      V[i]=std(intersect(V1,V0));
650
651      if(size(V[i])>size(V[i+1]))
652      {
653        if(number(eM[i]-1)>=number(n-1)/2)
654        {
655          k++;
656          a[k]=eM[i]-1;
657          dbprint(printlevel-voice+2,
658            "//gaussman::vfiltration: transform to V0");
659          v[k]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
660        }
661        else
662        {
663          if(j<0)
664          {
665            if(a[k]==number(n-1)/2)
666            {
667              j=k-1;
668            }
669            else
670            {
671              j=k;
672            }
673          }
674          k++;
675          a[k]=a[j];
676          a[j]=eM[i]-1;
677          v[k]=v[j];
678          dbprint(printlevel-voice+2,
679            "//gaussman::vfiltration: transform to V0");
680          v[j]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
681          j--;
682        }
683      }
684    }
685
686    dbprint(printlevel-voice+2,
687      "//gaussman::vfiltration: compute graded parts");
688    option(redSB);
689    for(k=1;k<size(v);k++)
690    {
691      v[k]=std(reduce(v[k],std(v[k+1])));
692      d[k]=size(v[k]);
693    }
694    v[k]=std(v[k]);
695    d[k]=size(v[k]);
696
697    return(list(a,d,v,m));
698  }
699}
700example
701{ "EXAMPLE:"; echo=2;
702  ring R=0,(x,y),ds;
703  poly f=x5+x2y2+y5;
704  vfiltration(f);
705}
706///////////////////////////////////////////////////////////////////////////////
707
708proc spectrum(poly f)
709"USAGE:    spectrum(f); poly f
710ASSUME:   basering has characteristic 0 and local ordering,
711          f has isolated singularity at 0
712RETURN:
713@format
714list S: singularity spectrum of f
715  ideal S[1]: spectral numbers in increasing order
716  intvec S[2]:
717    int S[2][i]: multiplicity of spectral number S[1][i]
718@end format
719SEE ALSO: spectrum_lib
720KEYWORDS: singularities; Gauss-Manin connection; spectrum
721EXAMPLE:  example spectrum; shows an example
722"
723{
724  return(vfiltration(f,0));
725}
726example
727{ "EXAMPLE:"; echo=2;
728  ring R=0,(x,y),ds;
729  poly f=x5+x2y2+y5;
730  spprint(spectrum(f));
731}
732///////////////////////////////////////////////////////////////////////////////
733
734proc endfilt(poly f,list V)
735"USAGE:   endfilt(f,V); poly f, list V
736ASSUME:  basering has characteristic 0 and local ordering,
737         f has isolated singularity at 0
738RETURN:
739@format
740list V1: endomorphim filtration of V on the Jacobian algebra of f
741  ideal V1[1]: spectral numbers in increasing order
742  intvec V1[2]:
743    int V1[2][i]: multiplicity of spectral number V1[1][i]
744  list V1[3]:
745    module V1[3][i]: vector space basis of the V1[1][i]-th graded part
746                     in terms of V1[4]
747  ideal V1[4]: monomial vector space basis
748@end format
749SEE ALSO: spectrum_lib
750KEYWORDS: singularities; Gauss-Manin connection; spectrum;
751          Brieskorn lattice; Hodge filtration; V-filtration
752EXAMPLE: example endfilt; shows an example
753"
754{
755  if(charstr(basering)!="0")
756  {
757    ERROR("characteristic 0 expected");
758  }
759  int n=nvars(basering)-1;
760  for(int i=n+1;i>=1;i--)
761  {
762    if(var(i)>1)
763    {
764      ERROR("local ordering expected");
765    }
766  }
767  ideal sJ=std(jacob(f));
768  if(vdim(sJ)<=0)
769  {
770    if(vdim(sJ)==0)
771    {
772      ERROR("singularity at 0 expected");
773    }
774    else
775    {
776      ERROR("isolated singularity at 0 expected");
777    }
778  }
779
780  def a,d,v,m=V[1..4];
781  int mu=ncols(m);
782
783  module V0=v[1];
784  for(i=2;i<=size(v);i++)
785  {
786    V0=V0,v[i];
787  }
788
789  dbprint(printlevel-voice+2,
790    "//gaussman::endfilt: compute multiplication in Jacobian algebra");
791  list M;
792  matrix U=freemodule(ncols(m));
793  for(i=ncols(m);i>=1;i--)
794  {
795    M[i]=lift(V0,coeffs(reduce(m[i]*m,U,sJ),m)*V0);
796  }
797
798  int j,k,i0,j0,i1,j1;
799  number b0=number(a[1]-a[ncols(a)]);
800  number b1,b2;
801  matrix M0;
802  module L;
803  list v0=freemodule(ncols(m));
804  ideal a0=b0;
805
806  while(b0<number(a[ncols(a)]-a[1]))
807  {
808    dbprint(printlevel-voice+2,
809      "//gaussman::endfilt: find next possible index");
810    b1=number(a[ncols(a)]-a[1]);
811    for(j=ncols(a);j>=1;j--)
812    {
813      for(i=ncols(a);i>=1;i--)
814      {
815        b2=number(a[i]-a[j]);
816        if(b2>b0&&b2<b1)
817        {
818          b1=b2;
819        }
820        else
821        {
822          if(b2<=b0)
823          {
824            i=0;
825          }
826        }
827      }
828    }
829    b0=b1;
830
831    list l=ideal();
832    for(k=ncols(m);k>=2;k--)
833    {
834      l=l+list(ideal());
835    }
836
837    dbprint(printlevel-voice+2,
838      "//gaussman::endfilt: collect conditions for V1["+string(b0)+"]");
839    j=ncols(a);
840    j0=mu;
841    while(j>=1)
842    {
843      i0=1;
844      i=1;
845      while(i<ncols(a)&&a[i]<a[j]+b0)
846      {
847        i0=i0+d[i];
848        i++;
849      }
850      if(a[i]<a[j]+b0)
851      {
852        i0=i0+d[i];
853        i++;
854      }
855      for(k=1;k<=ncols(m);k++)
856      {
857        M0=M[k];
858        if(i0>1)
859        {
860          l[k]=l[k],M0[1..i0-1,j0-d[j]+1..j0];
861        }
862      }
863      j0=j0-d[j];
864      j--;
865    }
866
867    dbprint(printlevel-voice+2,
868      "//gaussman::endfilt: compose condition matrix");
869    L=transpose(module(l[1]));
870    for(k=2;k<=ncols(m);k++)
871    {
872      L=L,transpose(module(l[k]));
873    }
874
875    dbprint(printlevel-voice+2,
876      "//gaussman::endfilt: compute kernel of condition matrix");
877    v0=v0+list(syz(L));
878    a0=a0,b0;
879  }
880
881  dbprint(printlevel-voice+2,"//gaussman::endfilt: compute graded parts");
882  option(redSB);
883  for(i=1;i<size(v0);i++)
884  {
885    v0[i+1]=std(v0[i+1]);
886    v0[i]=std(reduce(v0[i],v0[i+1]));
887  }
888
889  dbprint(printlevel-voice+2,
890    "//gaussman::endfilt: remove trivial graded parts");
891  i=1;
892  while(size(v0[i])==0)
893  {
894    i++;
895  }
896  list v1=v0[i];
897  intvec d1=size(v0[i]);
898  ideal a1=a0[i];
899  i++;
900  while(i<=size(v0))
901  {
902    if(size(v0[i])>0)
903    {
904      v1=v1+list(v0[i]);
905      d1=d1,size(v0[i]);
906      a1=a1,a0[i];
907    }
908    i++;
909  }
910  return(list(a1,d1,v1,m));
911}
912example
913{ "EXAMPLE:"; echo=2;
914  ring R=0,(x,y),ds;
915  poly f=x5+x2y2+y5;
916  endfilt(f,vfiltration(f));
917}
918///////////////////////////////////////////////////////////////////////////////
919
920proc spprint(list S)
921"USAGE:   spprint(S); list S
922RETURN:  string: spectrum S
923EXAMPLE: example spprint; shows an example
924"
925{
926  string s;
927  for(int i=1;i<size(S[2]);i++)
928  {
929    s=s+"("+string(S[1][i])+","+string(S[2][i])+"),";
930  }
931  s=s+"("+string(S[1][i])+","+string(S[2][i])+")";
932  return(s);
933}
934example
935{ "EXAMPLE:"; echo=2;
936  ring R=0,(x,y),ds;
937  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));
938  spprint(S);
939}
940///////////////////////////////////////////////////////////////////////////////
941
942proc spadd(list S1,list S2)
943"USAGE:   spadd(S1,S2); list S1,S2
944RETURN:  list: sum of spectra S1 and S2
945EXAMPLE: example spadd; shows an example
946"
947{
948  ideal s;
949  intvec m;
950  int i,i1,i2=1,1,1;
951  while(i1<=size(S1[2])||i2<=size(S2[2]))
952  {
953    if(i1<=size(S1[2]))
954    {
955      if(i2<=size(S2[2]))
956      {
957        if(number(S1[1][i1])<number(S2[1][i2]))
958        {
959          s[i]=S1[1][i1];
960          m[i]=S1[2][i1];
961          i++;
962          i1++;
963        }
964        else
965        {
966          if(number(S1[1][i1])>number(S2[1][i2]))
967          {
968            s[i]=S2[1][i2];
969            m[i]=S2[2][i2];
970            i++;
971            i2++;
972          }
973          else
974          {
975            if(S1[2][i1]+S2[2][i2]!=0)
976            {
977              s[i]=S1[1][i1];
978              m[i]=S1[2][i1]+S2[2][i2];
979              i++;
980            }
981            i1++;
982            i2++;
983          }
984        }
985      }
986      else
987      {
988        s[i]=S1[1][i1];
989        m[i]=S1[2][i1];
990        i++;
991        i1++;
992      }
993    }
994    else
995    {
996      s[i]=S2[1][i2];
997      m[i]=S2[2][i2];
998      i++;
999      i2++;
1000    }
1001  }
1002  return(list(s,m));
1003}
1004example
1005{ "EXAMPLE:"; echo=2;
1006  ring R=0,(x,y),ds;
1007  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));
1008  spprint(S1);
1009  list S2=list(ideal(-1/6,1/6),intvec(1,1));
1010  spprint(S2);
1011  spprint(spadd(S1,S2));
1012}
1013///////////////////////////////////////////////////////////////////////////////
1014
1015proc spsub(list S1,list S2)
1016"USAGE:   spsub(S1,S2); list S1,S2
1017RETURN:  list: difference of spectra S1 and S2
1018EXAMPLE: example spsub; shows an example
1019"
1020{
1021  return(spadd(S1,spmul(S2,-1)));
1022}
1023example
1024{ "EXAMPLE:"; echo=2;
1025  ring R=0,(x,y),ds;
1026  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));
1027  spprint(S1);
1028  list S2=list(ideal(-1/6,1/6),intvec(1,1));
1029  spprint(S2);
1030  spprint(spsub(S1,S2));
1031}
1032///////////////////////////////////////////////////////////////////////////////
1033
1034proc spmul(list #)
1035"USAGE:
1036@format
10371) spmul(S,k); list S, int k
10382) spmul(S,k); list S, intvec k
1039@end format
1040RETURN:
1041@format
10421) list: product of spectrum S and integer k
10432) list: linear combination of spectra S with coefficients k
1044@end format
1045EXAMPLE: example spmul; shows an example
1046"
1047{
1048  if(size(#)==2)
1049  {
1050    if(typeof(#[1])=="list")
1051    {
1052      if(typeof(#[2])=="int")
1053      {
1054        return(list(#[1][1],#[1][2]*#[2]));
1055      }
1056      if(typeof(#[2])=="intvec")
1057      {
1058        list S0=list(ideal(),intvec(0));
1059        for(int i=size(#[2]);i>=1;i--)
1060        {
1061          S0=spadd(S0,spmul(#[1][i],#[2][i]));
1062        }
1063        return(S0);
1064      }
1065    }
1066  }
1067  return(list(ideal(),intvec(0)));
1068}
1069example
1070{ "EXAMPLE:"; echo=2;
1071  ring R=0,(x,y),ds;
1072  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));
1073  spprint(S);
1074  spprint(spmul(S,2));
1075  list S1=list(ideal(-1/6,1/6),intvec(1,1));
1076  spprint(S1);
1077  list S2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
1078  spprint(S2);
1079  spprint(spmul(list(S1,S2),intvec(1,2)));
1080}
1081///////////////////////////////////////////////////////////////////////////////
1082
1083proc spissemicont(list S,list #)
1084"USAGE:   spissemicont(S[,opt]); list S, int opt
1085RETURN:
1086@format
1087int k=
1088if opt==0:
1089  1, if sum of spectrum S over all intervals [a,a+1) is positive
1090  0, if sum of spectrum S over some interval [a,a+1) is negative
1091if opt==1:
1092  1, if sum of spectrum S over all intervals [a,a+1) and (a,a+1) is positive
1093  0, if sum of spectrum S over some interval [a,a+1) or (a,a+1) is negative
1094default: opt=0
1095@end format
1096EXAMPLE: example spissemicont; shows an example
1097"
1098{
1099  int opt=0;
1100  if(size(#)>0)
1101  {
1102    if(typeof(#[1])=="int")
1103    {
1104      opt=1;
1105    }
1106  }
1107  int i,j,k=1,1,0;
1108  while(j<=size(S[2]))
1109  {
1110    while(j+1<=size(S[2])&&S[1][j]<S[1][i]+1)
1111    {
1112      k=k+S[2][j];
1113      j++;
1114    }
1115    if(j==size(S[2])&&S[1][j]<S[1][i]+1)
1116    {
1117      k=k+S[2][j];
1118      j++;
1119    }
1120    if(k<0)
1121    {
1122      return(0);
1123    }
1124    k=k-S[2][i];
1125    if(k<0&&opt==1)
1126    {
1127      return(0);
1128    }
1129    i++;
1130  }
1131  return(1);
1132}
1133example
1134{ "EXAMPLE:"; echo=2;
1135  ring R=0,(x,y),ds;
1136  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));
1137  spprint(S1);
1138  list S2=list(ideal(-1/6,1/6),intvec(1,1));
1139  spprint(S2);
1140  spissemicont(spsub(S1,spmul(S2,5)));
1141  spissemicont(spsub(S1,spmul(S2,5)),1);
1142  spissemicont(spsub(S1,spmul(S2,6)));
1143}
1144///////////////////////////////////////////////////////////////////////////////
1145
1146proc spsemicont(list S0,list S,list #)
1147"USAGE:   spsemicont(S,k[,opt]); list S0, list S, int opt
1148RETURN:  list of intvecs l:
1149         spissemicont(sub(S0,spmul(S,k)),opt)==1 iff k<=l[i] for some i
1150NOTE:    if the spectra S occur with multiplicities k in a deformation
1151         of the [quasihomogeneous] spectrum S0 then
1152         spissemicont(sub(S0,spmul(S,k))[,1])==1
1153EXAMPLE: example spsemicont; shows an example
1154"
1155{
1156  list l,l0;
1157  int i,j,k;
1158  while(spissemicont(S0,#))
1159  {
1160    if(size(S)>1)
1161    {
1162      l0=spsemicont(S0,list(S[1..size(S)-1]));
1163      for(i=1;i<=size(l0);i++)
1164      {
1165        if(size(l)>0)
1166        {
1167          j=1;
1168          while(j<size(l)&&l[j]!=l0[i])
1169          {
1170            j++;
1171          }
1172          if(l[j]==l0[i])
1173          {
1174            l[j][size(S)]=k;
1175          }
1176          else
1177          {
1178            l0[i][size(S)]=k;
1179            l=l+list(l0[i]);
1180          }
1181        }
1182        else
1183        {
1184          l=l0;
1185        }
1186      }
1187    }
1188    S0=spsub(S0,S[size(S)]);
1189    k++;
1190  }
1191  if(size(S)>1)
1192  {
1193    return(l);
1194  }
1195  else
1196  {
1197    return(list(intvec(k-1)));
1198  }
1199}
1200example
1201{ "EXAMPLE:"; echo=2;
1202  ring R=0,(x,y),ds;
1203  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));
1204  spprint(S0);
1205  list S1=list(ideal(-1/6,1/6),intvec(1,1));
1206  spprint(S1);
1207  list S2=list(ideal(-1/3,0,1/3),intvec(1,2,1));
1208  spprint(S2);
1209  list S=S1,S2;
1210  list l=spsemicont(S0,S);
1211  l;
1212  spissemicont(spsub(S0,spmul(S,l[1])));
1213  spissemicont(spsub(S0,spmul(S,l[1]-1)));
1214  spissemicont(spsub(S0,spmul(S,l[1]+1)));
1215}
1216///////////////////////////////////////////////////////////////////////////////
1217
1218proc spmilnor(list S)
1219"USAGE:   spmilnor(S); list S
1220RETURN:  int: Milnor number of spectrum S
1221EXAMPLE: example spmilnor; shows an example
1222"
1223{
1224  return(sum(S[2]));
1225}
1226example
1227{ "EXAMPLE:"; echo=2;
1228  ring R=0,(x,y),ds;
1229  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));
1230  spprint(S);
1231  spmilnor(S);
1232}
1233///////////////////////////////////////////////////////////////////////////////
1234
1235proc spgeomgenus(list S)
1236"USAGE:   spgeomgenus(S); list S
1237RETURN:  int: geometrical genus of spectrum S
1238EXAMPLE: example spgeomgenus; shows an example
1239"
1240{
1241  int g=0;
1242  int i=1;
1243  while(i+1<=size(S[2])&&number(S[1][i])<=number(0))
1244  {
1245    g=g+S[2][i];
1246    i++;
1247  }
1248  if(i==size(S[2])&&number(S[1][i])<=number(0))
1249  {
1250    g=g+S[2][i];
1251  }
1252  return(g);
1253}
1254example
1255{ "EXAMPLE:"; echo=2;
1256  ring R=0,(x,y),ds;
1257  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));
1258  spprint(S);
1259  spgeomgenus(S);
1260}
1261///////////////////////////////////////////////////////////////////////////////
1262
1263proc spgamma(list S)
1264"USAGE:   spgamma(S); list S
1265RETURN:  number: gamma invariant of spectrum S
1266EXAMPLE: example spgamma; shows an example
1267"
1268{
1269  int i,j;
1270  number g=0;
1271  for(i=1;i<=ncols(S[1]);i++)
1272  {
1273    for(j=1;j<=S[2][i];j++)
1274    {
1275      g=g+(number(S[1][i])-number(nvars(basering)-2)/2)^2;
1276    }
1277  }
1278  g=-g/4+sum(S[2])*number(S[1][ncols(S[1])]-S[1][1])/48;
1279  return(g);
1280}
1281example
1282{ "EXAMPLE:"; echo=2;
1283  ring R=0,(x,y),ds;
1284  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));
1285  spprint(S);
1286  spgamma(S);
1287}
1288///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.