source: git/Singular/LIB/gaussman.lib @ 96967e

spielwiese
Last change on this file since 96967e was 96967e, checked in by Gert-Martin Greuel <greuel@…>, 23 years ago
* GMG: Kosmetik fuer html-Hilfe git-svn-id: file:///usr/local/Singular/svn/trunk@4989 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 21.0 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2
3version="$Id: gaussman.lib,v 1.19 2000-12-23 17:11:30 greuel Exp $";
4category="Singularities";
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[,...]);    monodromy matrix of f, spectrum of monodromy of f
15 vfiltration(f[,...]);  V-filtration of f on H''/H', singularity spectrum of f
16 vfiltjacalg(...);      V-filtration on Jacobian algebra
17 gamma(...);            Hertling's gamma invariant
18 gamma4(...);           Hertling's gamma4 invariant
19
20SEE ALSO: mondromy.lib, spectrum.lib, jordan.lib
21
22KEYWORDS: singularities; Gauss-Manin connection; monodromy; spectrum;
23          Brieskorn lattice; Hodge filtration; V-filtration;
24";
25
26LIB "jordan.lib";
27
28///////////////////////////////////////////////////////////////////////////////
29
30static proc maxintdif(ideal e)
31{
32  dbprint(printlevel-voice+2,"//gaussman::maxintdif");
33  int i,j,id;
34  int mid=0;
35  for(i=ncols(e);i>=1;i--)
36  {
37    for(j=i-1;j>=1;j--)
38    {
39      id=int(e[i]-e[j]);
40      if(id<0)
41      {
42        id=-id;
43      }
44      if(id>mid)
45      {
46        mid=id;
47      }
48    }
49  }
50  return(mid);
51}
52///////////////////////////////////////////////////////////////////////////////
53
54static proc maxorddif(matrix H)
55{
56  dbprint(printlevel-voice+2,"//gaussman::maxorddif");
57  int i,j,d;
58  int d0,d1=-1,-1;
59  for(i=nrows(H);i>=1;i--)
60  {
61    for(j=ncols(H);j>=1;j--)
62    {
63      d=ord(H[i,j]);
64      if(d>=0)
65      {
66        if(d0<0||d<d0)
67        {
68          d0=d;
69        }
70        if(d1<0||d>d1)
71        {
72          d1=d;
73        }
74      }
75    }
76  }
77  return(d1-d0);
78}
79///////////////////////////////////////////////////////////////////////////////
80
81static proc invunit(poly u,int n)
82{
83  dbprint(printlevel-voice+2,"//gaussman::invunit");
84  if(ord(u)>0)
85  {
86    ERROR("no unit");
87  }
88  poly u0=jet(u,0);
89  u=jet(1-u/u0,n);
90  poly w=u;
91  poly v=1+u;
92  for(int i=n div ord(u);i>1;i--)
93  {
94    w=jet(w*u,n);
95    v=v+w;
96  }
97  v=jet(v,n)/u0;
98  return(v);
99}
100///////////////////////////////////////////////////////////////////////////////
101
102static proc expand(matrix M,matrix U,int n)
103{
104  dbprint(printlevel-voice+2,"//gaussman::expand");
105  for(int i=ncols(U);i>=1;i--)
106  {
107    U[i,i]=invunit(U[i,i],n-ord(M[i]));
108  }
109  return(jet(M*U,n));
110}
111///////////////////////////////////////////////////////////////////////////////
112
113static proc redNF(module M,module N,list #)
114{
115  matrix U=freemodule(ncols(M));
116  if(size(#)>0)
117  {
118    if(typeof(#[1])=="matrix")
119    {
120      U=#[1];
121    }
122  }
123
124  for(int i=ncols(M);i>=1;i--)
125  {
126    M[i]=M[i]/lead(U[i,i]);
127    U[i,i]=U[i,i]/lead(U[i,i]);
128  }
129  N=std(N);
130  module redNFM=matrix(0,nrows(M),ncols(M));
131  dbprint(printlevel-voice+2,"//gaussman::redNF: reduce");
132  module weakNFM=reduce(M,N);
133  while(size(weakNFM)>0)
134  {
135    redNFM=matrix(redNFM)+matrix(lead(weakNFM));
136    M=matrix(M)-matrix(lead(weakNFM))*U;
137    dbprint(printlevel-voice+2,"//gaussman::redNF: reduce");
138    weakNFM=reduce(M,N);
139  }
140  return(redNFM);
141}
142///////////////////////////////////////////////////////////////////////////////
143
144proc monodromy(poly f,list #)
145"USAGE:    monodromy(f[,mode]); poly f, int mode[=1]
146ASSUME:   local ordering, f isolated singularity at 0
147RETURN:   if mode=0 :
148            matrix M : exp(-2*pi*i*M) monodromy matrix of f
149          if mode=1 :
150            ideal e : exp(-2*pi*i*e) spectrum of monodromy of f
151SEE ALSO: monodromy.lib, jordan.lib
152KEYWORDS: singularities; Gauss-Manin connection; monodromy;
153          Brieskorn lattice
154EXAMPLE:  example monodromy; shows an example
155"
156{
157  int mode=1;
158  if(size(#)>0)
159  {
160    if(typeof(#[1])=="int")
161    {
162      mode=#[1];
163    }
164  }
165
166  int i,j;
167  int n=nvars(basering)-1;
168  for(i=n+1;i>=1;i--)
169  {
170    if(var(i)>1)
171    {
172      ERROR("no local ordering");
173    }
174  }
175  ideal J=jacob(f);
176  ideal sJ=std(J);
177  if(vdim(sJ)<=0)
178  {
179    if(vdim(sJ)==0)
180    {
181      ERROR("no singularity at 0");
182    }
183    else
184    {
185      ERROR("no isolated singularity at 0");
186    }
187  }
188  ideal m=kbase(sJ);
189  int mu,modm=ncols(m),maxorddif(m);
190
191  ideal w=f*m;
192  matrix U=freemodule(mu);
193  matrix A0[mu][mu],A,C,D;
194  list l;
195  module H,dH=freemodule(mu),freemodule(mu);
196  module H0;
197  int sdH=1;
198  int k=-1;
199  int K,N,mide;
200
201  while(k<K||sdH>0)
202  {
203    k++;
204    dbprint(printlevel-voice+2,"//gaussman::monodromy: k="+string(k));
205
206    dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
207    C=coeffs(redNF(w,sJ,U),m);
208    A0=A0+C*var(1)^k;
209
210    if(sdH>0)
211    {
212      H0=H;
213      dbprint(printlevel-voice+2,"//gaussman::monodromy: compute dH");
214      dH=jet(module(A0*dH+var(1)^2*diff(matrix(dH),var(1))),k);
215      H=H*var(1)+dH;
216
217      dbprint(printlevel-voice+2,"//gaussman::monodromy: test dH==0");
218      sdH=size(reduce(H,std(H0*var(1))));
219      if(sdH>0)
220      {
221        A0=A0-var(1);
222      }
223      else
224      {
225        dbprint(printlevel-voice+2,
226          "//gaussman::monodromy: compute basis of saturation");
227        H=minbase(H0);
228        int modH=maxorddif(H);
229        K=modH+1;
230      }
231    }
232
233    if(k==K&&sdH==0)
234    {
235      N=k-modH;
236      dbprint(printlevel-voice+2,
237        "//gaussman::monodromy: compute A on saturation");
238      l=division(H*var(1),A0*H+var(1)^2*diff(matrix(H),var(1)));
239      A=expand(l[1],l[2],N-1);
240      dbprint(printlevel-voice+2,
241        "//gaussman::monodromy: compute eigenvalues e of A");
242      ideal e=jordan(A,-1)[1];
243      dbprint(printlevel-voice+2,"//gaussman::monodromy: e="+string(e));
244      if(mode==1)
245      {
246        mide=maxintdif(e);
247        K=K+mide;
248      }
249    }
250
251    if(k<K||sdH>0)
252    {
253      dbprint(printlevel-voice+2,"//gaussman::monodromy: divide by J");
254      l=division(J,ideal(matrix(w)-matrix(m)*C*U));
255      D=l[1];
256
257      dbprint(printlevel-voice+2,"//gaussman::monodromy: compute w/U");
258      for(j=mu;j>=1;j--)
259      {
260        if(l[2][j,j]!=0)
261        {
262          dbprint(printlevel-voice+2,
263            "//gaussman::monodromy: compute U["+string(j)+"]");
264          U[j,j]=U[j,j]*l[2][j,j];
265        }
266        dbprint(printlevel-voice+2,
267          "//gaussman::monodromy: compute w["+string(j)+"]");
268        w[j]=0;
269        for(i=n+1;i>=1;i--)
270        {
271          w[j]=w[j]+U[j,j]*diff(D[i,j],var(i))-diff(U[j,j],var(i))*D[i,j];
272        }
273      }
274      U=U*U;
275    }
276  }
277
278  while(mide>0)
279  {
280    dbprint(printlevel-voice+2,"//gaussman::monodromy: mide="+string(mide));
281
282    intvec b;
283    U=0;
284    module dU;
285    A0=jet(A,0);
286    matrix A0e;
287    for(i=ncols(e);i>=1;i--)
288    {
289      A0e=freemodule(mu);
290      for(j=n;j>=0;j--) // Potenzen von Matrizen?
291      {
292        A0e=A0e*(A0-e[i]);
293      }
294      dU=syz(A0e);
295      U=dU+U;
296      b[i]=size(dU);
297    }
298    A=division(U,A*U)[1];
299
300    intvec ide;
301    ide[mu]=0;
302    for(i=ncols(e);i>=1;i--)
303    {
304      for(j=i-1;j>=1;j--)
305      {
306        k=int(e[j]-e[i]);
307        if(k>ide[i])
308        {
309          ide[i]=k;
310        }
311        if(-k>ide[j])
312        {
313          ide[j]=-k;
314        }
315      }
316    }
317    for(j,k=ncols(e),mu;j>=1;j--)
318    {
319      for(i=b[j];i>=1;i--)
320      {
321        ide[k]=ide[j];
322        k--;
323      }
324    }
325
326    for(i=mu;i>=1;i--)
327    {
328      if(ide[i]>0)
329      {
330        A[i,i]=A[i,i]+1;
331        e[i]=e[i]+1;
332      }
333      for(j=mu;j>=1;j--)
334      {
335        if(ide[i]==0&&ide[j]>0)
336        {
337          A[i,j]=A[i,j]*var(1);
338        }
339        else
340        {
341          if(ide[i]>0&&ide[j]==0)
342          {
343            A[i,j]=A[i,j]/var(1);
344          }
345        }
346      }
347    }
348    mide--;
349  }
350
351  if(mode==1)
352  {
353    return(jet(A,0));
354  }
355  else
356  {
357    return(e);
358  }
359}
360example
361{ "EXAMPLE:"; echo=2;
362  ring R=0,(x,y),ds;
363  poly f=x5+x2y2+y5;
364  matrix M=monodromy(f);
365  print(M);
366}
367///////////////////////////////////////////////////////////////////////////////
368
369proc vfiltration(poly f,list #)
370"USAGE:    vfiltration(f[,mode]); poly f, int mode[default=1]
371ASSUME:   local ordering, f isolated singularity at 0
372RETURN:   list l:
373          @format
374          if mode=0 or mode=1:
375            l[1]: ideal, spectral numbers in increasing order
376            l[2]: intvec
377                  l[2][i]: int, multiplicity of spectral number l[1][i]
378          if mode=1 :
379            l[3]: list
380                  l[3][i]: module, vector space basis of l[1][i]-th graded
381                           part of the V-filtration on H''/H' in terms of l[4]
382            l[4]: ideal, monomial vector space basis of H''/H'
383            l[5]: ideal, standard basis of Jacobian ideal
384          @end format
385NOTE:     H' and H'' denote Brieskorn lattices
386SEE ALSO: spectrum.lib
387KEYWORDS: singularities; Gauss-Manin connection; spectrum;
388          Brieskorn lattice; Hodge filtration; V-filtration;
389EXAMPLE:  example vfiltration; shows an example
390"
391{
392  int mode=1;
393  if(size(#)>0)
394  {
395    if(typeof(#[1])=="int")
396    {
397      mode=#[1];
398    }
399  }
400
401  int i,j;
402  int n=nvars(basering)-1;
403  for(i=n+1;i>=1;i--)
404  {
405    if(var(i)>1)
406    {
407      ERROR("no local ordering");
408    }
409  }
410  ideal J=jacob(f);
411  ideal sJ=std(J);
412  if(vdim(sJ)<=0)
413  {
414    if(vdim(sJ)==0)
415    {
416      ERROR("no singularity at 0");
417    }
418    else
419    {
420      ERROR("no isolated singularity at 0");
421    }
422  }
423  ideal m=kbase(sJ);
424  int mu,modm=ncols(m),maxorddif(m);
425
426  ideal w=f*m;
427  matrix U=freemodule(mu);
428  matrix A[mu][mu],C,D;
429  list l;
430  module H,dH=freemodule(mu),freemodule(mu);
431  module H0;
432  int sdH=1;
433  int k=-1;
434  int K;
435
436  while(k<K||sdH>0)
437  {
438    k++;
439    dbprint(printlevel-voice+2,"//gaussman::vfiltration: k="+string(k));
440
441    dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute C");
442    C=coeffs(redNF(w,sJ,U),m);
443    A=A+C*var(1)^k;
444
445    if(sdH>0)
446    {
447      H0=H;
448      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute dH");
449      dH=jet(module(A*dH+var(1)^2*diff(matrix(dH),var(1))),k);
450      H=H*var(1)+dH;
451
452      dbprint(printlevel-voice+2,"//gaussman::vfiltration: test dH==0");
453      sdH=size(reduce(H,std(H0*var(1))));
454      if(sdH>0)
455      {
456        A=A-var(1);
457      }
458      else
459      {
460        dbprint(printlevel-voice+2,
461          "//gaussman::vfiltration: compute basis of saturation");
462        H=minbase(H0);
463        int modH=maxorddif(H);
464        if(k<n)
465        {
466          K=modH+n+1;
467        }
468        else
469        {
470          K=modH+k+1;
471        }
472        H0=freemodule(mu)*var(1)^k;
473      }
474    }
475
476    if(k<K||sdH>0)
477    {
478      dbprint(printlevel-voice+2,"//gaussman::vfiltration: divide by J");
479      l=division(J,ideal(matrix(w)-matrix(m)*C*U));
480      D=l[1];
481
482      dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute w/U");
483      for(j=mu;j>=1;j--)
484      {
485        if(l[2][j,j]!=0)
486        {
487          dbprint(printlevel-voice+2,
488            "//gaussman::vfiltration: compute U["+string(j)+"]");
489          U[j,j]=U[j,j]*l[2][j,j];
490        }
491        dbprint(printlevel-voice+2,
492          "//gaussman::vfiltration: compute w["+string(j)+"]");
493        w[j]=0;
494        for(i=n+1;i>=1;i--)
495        {
496          w[j]=w[j]+U[j,j]*diff(D[i,j],var(i))-diff(U[j,j],var(i))*D[i,j];
497        }
498      }
499      U=U*U;
500    }
501  }
502  int N=k-modH;
503
504  dbprint(printlevel-voice+2,
505    "//gaussman::vfiltration: transform H0 to saturation");
506  l=division(H,H0);
507  H0=expand(l[1],l[2],N-1);
508
509  dbprint(printlevel-voice+2,
510    "//gaussman::vfiltration: compute H0 as vector space V0");
511  dbprint(printlevel-voice+2,
512    "//gaussman::vfiltration: compute H1 as vector space V1");
513  poly p;
514  int i0,j0,i1,j1;
515  matrix V0[mu*N][mu*N];
516  matrix V1[mu*N][mu*(N-1)];
517  for(i0=mu;i0>=1;i0--)
518  {
519    for(i1=mu;i1>=1;i1--)
520    {
521      p=H0[i1,i0];
522      while(p!=0)
523      {
524        j1=leadexp(p)[1];
525        for(j0=N-j1-1;j0>=0;j0--)
526        {
527          V0[i1+(j1+j0)*mu,i0+j0*mu]=V0[i1+(j1+j0)*mu,i0+j0*mu]+leadcoef(p);
528          if(j1+j0+1<N)
529          {
530            V1[i1+(j1+j0+1)*mu,i0+j0*mu]=
531            V1[i1+(j1+j0+1)*mu,i0+j0*mu]+leadcoef(p);
532          }
533        }
534        p=p-lead(p);
535      }
536    }
537  }
538
539  dbprint(printlevel-voice+2,
540    "//gaussman::vfiltration: compute A on saturation");
541  l=division(H*var(1),A*H+var(1)^2*diff(matrix(H),var(1)));
542  A=expand(l[1],l[2],N-1);
543
544  dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute matrix M of A");
545  matrix M[mu*N][mu*N];
546  for(i0=mu;i0>=1;i0--)
547  {
548    for(i1=mu;i1>=1;i1--)
549    {
550      p=A[i1,i0];
551      while(p!=0)
552      {
553        j1=leadexp(p)[1];
554        for(j0=N-j1-1;j0>=0;j0--)
555        {
556          M[i1+(j0+j1)*mu,i0+j0*mu]=leadcoef(p);
557        }
558        p=p-lead(p);
559      }
560    }
561  }
562  for(i0=mu;i0>=1;i0--)
563  {
564    for(j0=N-1;j0>=0;j0--)
565    {
566      M[i0+j0*mu,i0+j0*mu]=M[i0+j0*mu,i0+j0*mu]+j0;
567    }
568  }
569
570  dbprint(printlevel-voice+2,
571    "//gaussman::vfiltration: compute eigenvalues eA of A");
572  ideal eA=jordan(A,-1)[1];
573  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eA="+string(eA));
574
575  dbprint(printlevel-voice+2,
576    "//gaussman::vfiltration: compute eigenvalues eM of M");
577  ideal eM;
578  k=0;
579  intvec u;
580  for(i=N;i>=1;i--)
581  {
582    u[i]=1;
583  }
584  i0=1;
585  while(u[N]<=ncols(eA))
586  {
587    for(i,i1=i0+1,i0;i<=N;i++)
588    {
589      if(eA[u[i]]+i<eA[u[i1]]+i1)
590      {
591        i1=i;
592      }
593    }
594    k++;
595    eM[k]=eA[u[i1]]+i1-1;
596    u[i1]=u[i1]+1;
597    if(u[i1]>ncols(eA))
598    {
599      i0=i1+1;
600    }
601  }
602  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eM="+string(eM));
603
604  dbprint(printlevel-voice+2,
605    "//gaussman::vfiltration: compute V-filtration on H0/H1");
606  ideal s;
607  k=0;
608  list V;
609  matrix Me;
610  V[ncols(eM)+1]=std(V1);
611  intvec d;
612  if(mode==0)
613  {
614    for(i=ncols(eM);number(eM[i])-1>number(n-1)/2;i--)
615    {
616      Me=freemodule(mu*N);
617      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
618      {
619        Me=Me*(M-eM[i]);
620      }
621
622      dbprint(printlevel-voice+2,
623        "//gaussman::vfiltration: compute V["+string(i)+"]");
624      V1=module(V1)+syz(Me);
625      V[i]=std(intersect(V1,V0));
626
627      if(size(V[i])>size(V[i+1]))
628      {
629        k++;
630        s[k]=eM[i]-1;
631        d[k]=size(V[i])-size(V[i+1]);
632      }
633    }
634
635    dbprint(printlevel-voice+2,
636      "//gaussman::vfiltration: symmetry index found");
637    j=k;
638    if(number(eM[i])-1==number(n-1)/2)
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        k++;
654        s[k]=eM[i]-1;
655        d[k]=size(V[i])-size(V[i+1]);
656      }
657    }
658
659    dbprint(printlevel-voice+2,"//gaussman::vfiltration: apply symmetry");
660    while(j>=1)
661    {
662      k++;
663      s[k]=s[j];
664      s[j]=n-1-s[k];
665      d[k]=d[j];
666      j--;
667    }
668
669    return(list(s,d));
670  }
671  else
672  {
673    list v;
674    j=-1;
675    for(i=ncols(eM);i>=1;i--)
676    {
677      Me=freemodule(mu*N);
678      for(i0=n;i0>=0;i0--) // Potenzen von Matrizen?
679      {
680        Me=Me*(M-eM[i]);
681      }
682
683      dbprint(printlevel-voice+2,
684        "//gaussman::vfiltration: compute V["+string(i)+"]");
685      V1=module(V1)+syz(Me);
686      V[i]=std(intersect(V1,V0));
687
688      if(size(V[i])>size(V[i+1]))
689      {
690        if(number(eM[i]-1)>=number(n-1)/2)
691        {
692          k++;
693          s[k]=eM[i]-1;
694          dbprint(printlevel-voice+2,
695            "//gaussman::vfiltration: transform to V0");
696          v[k]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
697        }
698        else
699        {
700          if(j<0)
701          {
702            if(s[k]==number(n-1)/2)
703            {
704              j=k-1;
705            }
706            else
707            {
708              j=k;
709            }
710          }
711          k++;
712          s[k]=s[j];
713          s[j]=eM[i]-1;
714          v[k]=v[j];
715          dbprint(printlevel-voice+2,
716            "//gaussman::vfiltration: transform to V0");
717          v[j]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1];
718          j--;
719        }
720      }
721    }
722
723    dbprint(printlevel-voice+2,
724      "//gaussman::vfiltration: compute graded parts");
725    option(redSB);
726    for(k=1;k<size(v);k++)
727    {
728      v[k]=std(reduce(v[k],std(v[k+1])));
729      d[k]=size(v[k]);
730    }
731    v[k]=std(v[k]);
732    d[k]=size(v[k]);
733
734    return(list(s,d,v,m,sJ));
735  }
736}
737example
738{ "EXAMPLE:"; echo=2;
739  ring R=0,(x,y),ds;
740  poly f=x5+x2y2+y5;
741  list l=vfiltration(f);
742  print(l);
743}
744///////////////////////////////////////////////////////////////////////////////
745
746proc vfiltjacalg(list l)
747"USAGE:   vfiltjacalg(vfiltration(f));
748ASSUME:  local ordering, f isolated singularity at 0
749RETURN:  list l:
750         @format
751           l[1]: ideal, spectral numbers of the V-filtration on the
752                 Jacobian algebra in increasing order
753           l[2]: intvec
754               l[2][i]: int, multiplicity of spectral number l[1][i]
755           l[3]: list
756               l[3][i]: module, vector space basis of l[1][i]-th graded part
757                        of the V-filtration on the Jacobian algebra in terms
758                        of l[4]
759           l[4]: ideal, monomial vector space basis of the Jacobian algebra
760           l[5]: ideal, standard basis of Jacobian ideal
761           @end format
762EXAMPLE: example vfiltjacalg; shows an example
763"
764{
765  def s,d,v,m,sJ=l[1..5];
766  int mu=ncols(m);
767
768  int i,j,k;
769  module V=v[1];
770  for(i=2;i<=size(v);i++)
771  {
772    V=V,v[i];
773  }
774
775  dbprint(printlevel-voice+2,
776    "//gaussman::vfiltjacalg: compute multiplication in Jacobian algebra");
777  list M;
778  for(i=ncols(m);i>=1;i--)
779  {
780    M[i]=lift(V,coeffs(redNF(m[i]*m,sJ),m)*V);
781  }
782
783  int i0,j0,i1,j1;
784  number r0=number(s[1]-s[ncols(s)]);
785  number r1,r2;
786  matrix M0;
787  module L;
788  list v0=freemodule(ncols(m));
789  ideal s0=r0;
790
791  while(r0<number(s[ncols(s)]-s[1]))
792  {
793    dbprint(printlevel-voice+2,
794      "//gaussman::vfiltjacalg: find next possible index");
795    r1=number(s[ncols(s)]-s[1]);
796    for(j=ncols(s);j>=1;j--)
797    {
798      for(i=ncols(s);i>=1;i--)
799      {
800        r2=number(s[i]-s[j]);
801        if(r2>r0&&r2<r1)
802        {
803          r1=r2;
804        }
805        else
806        {
807          if(r2<=r0)
808          {
809            i=0;
810          }
811        }
812      }
813    }
814    r0=r1;
815
816    l=ideal();
817    for(k=ncols(m);k>=2;k--)
818    {
819      l=l+list(ideal());
820    }
821
822    dbprint(printlevel-voice+2,
823      "//gaussman::vfiltjacalg: collect conditions for V["+string(r0)+"]");
824    j=ncols(s);
825    j0=mu;
826    while(j>=1)
827    {
828      i0=1;
829      i=1;
830      while(i<ncols(s)&&s[i]<s[j]+r0)
831      {
832        i0=i0+d[i];
833        i++;
834      }
835      if(s[i]<s[j]+r0)
836      {
837        i0=i0+d[i];
838        i++;
839      }
840      for(k=1;k<=ncols(m);k++)
841      {
842        M0=M[k];
843        if(i0>1)
844        {
845          l[k]=l[k],M0[1..i0-1,j0-d[j]+1..j0];
846        }
847      }
848      j0=j0-d[j];
849      j--;
850    }
851
852    dbprint(printlevel-voice+2,
853      "//gaussman::vfiltjacalg: compose condition matrix");
854    L=transpose(module(l[1]));
855    for(k=2;k<=ncols(m);k++)
856    {
857      L=L,transpose(module(l[k]));
858    }
859
860    dbprint(printlevel-voice+2,
861      "//gaussman::vfiltjacalg: compute kernel of condition matrix");
862    v0=v0+list(syz(L));
863    s0=s0,r0;
864  }
865
866  dbprint(printlevel-voice+2,"//gaussman::vfiltjacalg: compute graded parts");
867  option(redSB);
868  for(i=1;i<size(v0);i++)
869  {
870    v0[i+1]=std(v0[i+1]);
871    v0[i]=std(reduce(v0[i],v0[i+1]));
872  }
873
874  dbprint(printlevel-voice+2,
875    "//gaussman::vfiltjacalg: remove trivial graded parts");
876  i=1;
877  while(size(v0[i])==0)
878  {
879    i++;
880  }
881  list v1=v0[i];
882  intvec d1=size(v0[i]);
883  ideal s1=s0[i];
884  i++;
885  while(i<=size(v0))
886  {
887    if(size(v0[i])>0)
888    {
889      v1=v1+list(v0[i]);
890      d1=d1,size(v0[i]);
891      s1=s1,s0[i];
892    }
893    i++;
894  }
895  return(list(s1,d1,v1,m));
896}
897example
898{ "EXAMPLE:"; echo=2;
899  ring R=0,(x,y),ds;
900  poly f=x5+x2y2+y5;
901  vfiltjacalg(vfiltration(f));
902}
903///////////////////////////////////////////////////////////////////////////////
904
905proc gamma(list l)
906"USAGE:   gamma(vfiltration(f,0)); poly f
907ASSUME:  local ordering, f isolated singularity at 0
908RETURN:  number g : Hertling's gamma invariant
909EXAMPLE: example gamma; shows an example
910"
911{
912  ideal s=l[1];
913  intvec d=l[2];
914  int n=nvars(basering)-1;
915  number g=0;
916  int i,j;
917  for(i=1;i<=ncols(s);i++)
918  {
919    for(j=1;j<=d[i];j++)
920    {
921      g=g+(number(s[i])-number(n-1)/2)^2;
922    }
923  }
924  g=-g/4+sum(d)*number(s[ncols(s)]-s[1])/48;
925  return(g);
926}
927example
928{ "EXAMPLE:"; echo=2;
929  ring R=0,(x,y),ds;
930  poly f=x5+x2y2+y5;
931  gamma(vfiltration(f,0));
932}
933///////////////////////////////////////////////////////////////////////////////
934
935proc gamma4(list l)
936"USAGE:   gamma4(vfiltration(f,0)); poly f
937ASSUME:  local ordering, f isolated singularity at 0
938RETURN:  number g4 : Hertling's gamma4 invariant
939EXAMPLE: example gamma4; shows an example
940"
941{
942  ideal s=l[1];
943  intvec d=l[2];
944  int n=nvars(basering)-1;
945  number g4=0;
946  int i,j;
947  for(i=1;i<=ncols(s);i++)
948  {
949    for(j=1;j<=d[i];j++)
950    {
951      g4=g4+(number(s[i])-number(n-1)/2)^4;
952    }
953  }
954  g4=g4-(number(s[ncols(s)]-s[1])/12-1/30)*
955    (sum(d)*number(s[ncols(s)]-s[1])/4-24*gamma(l));
956 return(g4);
957}
958example
959{ "EXAMPLE:"; echo=2;
960  ring R=0,(x,y),ds;
961  poly f=x5+x2y2+y5;
962  gamma4(vfiltration(f,0));
963}
964///////////////////////////////////////////////////////////////////////////////
965
966proc tst_gaussm(poly f)
967{
968  echo=2;
969  basering;
970  f;
971  print(monodromy(f));
972  list l=vfiltration(f);
973  l;
974  vfiltjacalg(l);
975  gamma(l);
976  gamma4(l);
977}
978///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.