source: git/Singular/LIB/gmspoly.lib @ 91c7192

spielwiese
Last change on this file since 91c7192 was 91c7192, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: force matrix format git-svn-id: file:///usr/local/Singular/svn/trunk@10628 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.5 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: gmspoly.lib,v 1.16 2008-03-18 15:50:14 Singular Exp $";
3category="Singularities";
4
5info="
6LIBRARY:  gmspoly.lib  Gauss-Manin System of Tame Polynomials
7
8AUTHOR:   Mathias Schulze, email: mschulze@mathematik.uni-kl.de
9
10OVERVIEW: A library to compute invariants related to the Gauss-Manin system
11          of a cohomologically tame polynomial
12
13PROCEDURES:
14  isTame(f);  test if the polynomial f is tame
15  goodBasis(f);  a good basis of the Brieskorn lattice of a cohomologically tame f
16
17SEE ALSO: gmssing_lib
18
19KEYWORDS: tame polynomial; Gauss-Manin system; Brieskorn lattice;
20          mixed Hodge structure; V-filtration; weight filtration;
21          monodromy; spectrum; spectral pairs; good basis
22";
23
24LIB "linalg.lib";
25LIB "ring.lib";
26
27///////////////////////////////////////////////////////////////////////////////
28
29static proc mindegree(matrix A)
30{
31  int d=0;
32
33  while(A/var(1)^(d+1)*var(1)^(d+1)==A)
34  {
35    d++;
36  }
37
38  return(d);
39}
40///////////////////////////////////////////////////////////////////////////////
41
42static proc maxdegree(matrix A)
43{
44  int d=0;
45  matrix N[nrows(A)][ncols(A)];
46
47  while(A/var(1)^(d+1)!=N)
48  {
49    d++;
50  }
51
52  return(d);
53}
54///////////////////////////////////////////////////////////////////////////////
55
56proc isTame(poly f)
57"USAGE:    isTame(f); poly f
58ASSUME:   basering has no variables named w(1),w(2),...
59RETURN:
60@format
61int k=
62  0;  if f is not tame
63  1;  if f is tame
64@end format
65KEYWORDS: tame polynomial
66EXAMPLE:  example isTame; shows examples
67"
68{
69  int d=vdim(std(jacob(f)));
70  def @X=basering;
71  int n=nvars(@X);
72  def @WX=changechar("(0,w(1.."+string(n)+"))");
73  setring @WX;
74  ideal J=jacob(imap(@X,f));
75  int i;
76  for(i=1;i<=n;i++)
77  {
78    J[i]=J[i]+w(i);
79  }
80  int D=vdim(std(J));
81
82  setring(@X);
83  kill @WX;
84
85  return(d>0&&d==D);
86}
87example
88{ "EXAMPLE:"; echo=2;
89  ring R=0,(x,y),dp;
90  isTame(x2y+x);
91  isTame(x3+y3+xy);
92}
93///////////////////////////////////////////////////////////////////////////////
94
95static proc chart(matrix A)
96{
97  A=ideal(homog(transpose(ideal(A)),var(2)));
98  def r=basering;
99  map h=r,1,var(1);
100  return(h(A));
101}
102///////////////////////////////////////////////////////////////////////////////
103
104static proc pidbasis(module M0,module M)
105{
106  int m=nrows(M);
107  int n=ncols(M);
108
109  module L,N;
110  module T=freemodule(m);
111  while(matrix(L)!=matrix(M))
112  {
113    L=M;
114
115    M=T,M;
116    N=transpose(std(transpose(M)));
117    T=N[1..m];
118    M=N[m+1..ncols(N)];
119
120    M=freemodule(n),transpose(M);
121    N=std(transpose(M));
122    N=transpose(simplify(N,1));
123    M=N[n+1..ncols(N)];
124    M=transpose(M);
125  }
126
127  if(maxdegree(M)>0)
128  {
129    print("   ? module not free");
130    return(module());
131  }
132
133  attrib(M,"isSB",1);
134  N=lift(T,simplify(reduce(M0,M),2));
135
136  return(N);
137}
138///////////////////////////////////////////////////////////////////////////////
139
140static proc vfilt(matrix B,int d)
141{
142  int mu=ncols(B);
143
144  module V=freemodule(mu);
145  module V0=var(1)^(d-1)*freemodule(mu);
146  attrib(V0,"isSB",1);
147  module V1=B;
148  option(redSB);
149  while(size(reduce(V1,V0))>0)
150  {
151    V=std(V0+V1);
152    V0=var(1)^(d-1)*V;
153    attrib(V0,"isSB",1);
154    V1=B*matrix(V1)-var(1)^d*diff(matrix(V1),var(1));
155  }
156  option("noredSB");
157
158  B=lift(V0,B*matrix(V)-var(1)^d*diff(matrix(V),var(1)));
159  list l=eigenvals(B);
160  def e0,s0=l[1..2];
161
162  module U;
163  int i,j,i0,j0,i1,j1,k;
164  for(k=int(e0[ncols(e0)]-e0[1]);k>=1;k--)
165  {
166    U=0;
167    for(i=1;i<=ncols(e0);i++)
168    {
169      U=U+syz(power(jet(B,0)-e0[i],s0[i]));
170    }
171    B=lift(U,B*U);
172    V=V*U;
173
174    for(i0,i=1,1;i0<=ncols(e0);i0++)
175    {
176      for(i1=1;i1<=s0[i0];i1,i=i1+1,i+1)
177      {
178        for(j0,j=1,1;j0<=ncols(e0);j0++)
179        {
180          for(j1=1;j1<=s0[j0];j1,j=j1+1,j+1)
181          {
182            if(leadcoef(e0[i0]-e0[1])>=1&&leadcoef(e0[j0]-e0[1])<1)
183            {
184              B[i,j]=B[i,j]/var(1);
185            }
186            if(leadcoef(e0[i0]-e0[1])<1&&leadcoef(e0[j0]-e0[1])>=1)
187            {
188              B[i,j]=B[i,j]*var(1);
189            }
190          }
191        }
192      }
193    }
194
195    for(i0,i=1,1;i0<=ncols(e0);i0++)
196    {
197      if(leadcoef(e0[i0]-e0[1])>=1)
198      {
199        for(i1=1;i1<=s0[i0];i1,i=i1+1,i+1)
200        {
201          B[i,i]=B[i,i]-1;
202          V[i]=V[i]*var(1);
203        }
204        e0[i0]=e0[i0]-1;
205      }
206      else
207      {
208        i=i+s0[i0];
209      }
210    }
211
212    l=spnf(list(e0,s0));
213    e0,s0=l[1..2];
214  }
215
216  U=0;
217  for(i=1;i<=ncols(e0);i++)
218  {
219    U=U+syz(power(jet(B,0)-e0[i],s0[i]));
220  }
221  B=lift(U,B*U);
222  V=V*U;
223
224  d=mindegree(V);
225  V=V/var(1)^d;
226  B=B+d*matrix(freemodule(mu));
227  for(i=ncols(e0);i>=1;i--)
228  {
229    e0[i]=e0[i]+d;
230  }
231
232  return(e0,s0,V,B);
233}
234///////////////////////////////////////////////////////////////////////////////
235
236static proc spec(ideal e0,intvec s0,module V,matrix B)
237{
238  int mu=ncols(B);
239
240  int i,j,k;
241
242  int d=maxdegree(V);
243  int d0=d;
244  V=chart(V);
245  module U=std(V);
246  while(size(reduce(var(1)^d*freemodule(mu),U))>0)
247  {
248    d++;
249  }
250  if(d>d0)
251  {
252    k=d-d0;
253    B=B-k*freemodule(mu);
254    for(i=1;i<=ncols(e0);i++)
255    {
256      e0[i]=e0[i]-k;
257    }
258  }
259  module G=lift(V,var(1)^d*freemodule(mu));
260  G=std(G);
261  G=simplify(G,1);
262
263  ideal e;
264  intvec s;
265  e[mu]=0;
266  for(j,k=1,1;j<=ncols(e0);j++)
267  {
268    for(i=s0[j];i>=1;i,k=i-1,k+1)
269    {
270      e[k]=e0[j];
271      s[k]=j;
272    }
273  }
274
275  ideal a;
276  a[mu]=0;
277  for(i=1;i<=mu;i++)
278  {
279    a[i]=leadcoef(e[leadexp(G[i])[nvars(basering)+1]])+leadexp(G[i])[1];
280  }
281
282  return(a,e0,e,s,V,B,G);
283}
284///////////////////////////////////////////////////////////////////////////////
285
286static proc fsplit(ideal e0,ideal e,intvec s,module V,matrix B,module G)
287{
288  int mu=ncols(e);
289
290  int i,j,k;
291
292  number n,n0;
293  vector v,v0;
294  list F;
295  for(i=ncols(e0);i>=1;i--)
296  {
297    F[i]=module(matrix(0,mu,1));
298  }
299  for(i=mu;i>=1;i--)
300  {
301    v=G[i];
302    v0=lead(v);
303    n0=leadcoef(e[leadexp(v0)[nvars(basering)+1]])+leadexp(v0)[1];
304    v=v-lead(v);
305    while(v!=0)
306    {
307      n=leadcoef(e[leadexp(v)[nvars(basering)+1]])+leadexp(v)[1];
308      if(n==n0)
309      {
310        v0=v0+lead(v);
311        v=v-lead(v);
312      }
313      else
314      {
315        v=0;
316      }
317    }
318    j=s[leadexp(v0)[nvars(basering)+1]];
319    F[j]=F[j]+v0;
320  }
321
322  matrix B0=jet(B,0);
323  module U,U0,U1,U2;
324  matrix N;
325  for(i=size(F);i>=1;i--)
326  {
327    N=B0-e0[i];
328    U0=0;
329    while(size(F[i])>0)
330    {
331      k=0;
332      U1=jet(F[i],0);
333      while(size(U1)>0)
334      {
335        for(j=ncols(U1);j>=1;j--)
336        {
337          if(size(reduce(U1[j],std(U0)))>0)
338          {
339            U0=U1[j]+U0;
340          }
341        }
342        U1=N*U1;
343        k++;
344      }
345      F[i]=module(F[i]/var(1));
346    }
347    U=U0+U;
348  }
349
350  V=V*U;
351  G=lift(U,G);
352  B=lift(U,B*U);
353
354  return(e,V,B,G);
355}
356///////////////////////////////////////////////////////////////////////////////
357
358static proc glift(ideal e,module V,matrix B,module G)
359{
360  int mu=ncols(e);
361
362  int d=maxdegree(B);
363  B=chart(B);
364  G=std(G);
365  G=simplify(G,1);
366
367  int i,j,k;
368
369  ideal v;
370  for(i=mu;i>=1;i--)
371  {
372    v[i]=e[leadexp(G[i])[nvars(basering)+1]]+leadexp(G[i])[1];
373  }
374
375  number c;
376  matrix g[mu][1];
377  matrix m[mu][1];
378  matrix a[mu][1];
379  matrix A[mu][mu];
380  module M=var(1)^d*G;
381  module N=var(1)*B*matrix(G)+var(1)^(d+2)*diff(matrix(G),var(1));
382  while(size(N)>0)
383  {
384    j=mu;
385    for(k=mu-1;k>=1;k--)
386    {
387      if(N[k]>N[j])
388      {
389        j=k;
390      }
391    }
392
393    i=mu;
394    while(leadexp(M[i])[nvars(basering)+1]!=leadexp(N[j])[nvars(basering)+1])
395    {
396      i--;
397    }
398
399    k=leadexp(N[j])[1]-leadexp(M[i])[1];
400    if(k==0||i==j)
401    {
402      c=leadcoef(N[j])/leadcoef(M[i]);
403      A[i,j]=A[i,j]+c*var(1)^k;
404      N[j]=N[j]-c*var(1)^k*M[i];
405    }
406    else
407    {
408      c=leadcoef(N[j])/leadcoef(M[i])/(1-k-leadcoef(v[i])+leadcoef(v[j]));
409      G[j]=G[j]+c*var(1)^(k-1)*G[i];
410      M[j]=M[j]+c*var(1)^(k-1)*M[i];
411      g=c*var(1)^(k-1)*G[i];
412      N[j]=N[j]+(var(1)*B*g+var(1)^(d+2)*diff(g,var(1)))[1];
413      m=M[i];
414      a=transpose(A)[j];
415      N=N-c*var(1)^(k-1)*m*transpose(a);
416    }
417  }
418
419  G=V*G;
420  G=G/var(1)^mindegree(G);
421
422  return(G,A);
423}
424///////////////////////////////////////////////////////////////////////////////
425
426proc goodBasis(poly f)
427"USAGE:    goodBasis(f); poly f
428ASSUME:   f is cohomologically tame
429RETURN:
430@format
431ring R;  basering with new variable s
432  ideal b;  [matrix(b)] is a good basis of the Brieskorn lattice
433  matrix A;  A(s)=A0+s*A1 and t[matrix(b)]=[matrix(b)](A(s)+s^2*(d/ds))
434@end format
435KEYWORDS: tame polynomial; Gauss-Manin system; Brieskorn lattice;
436          mixed Hodge structure; V-filtration; weight filtration;
437          monodromy; spectrum; spectral pairs; good basis
438EXAMPLE:  example goodBasis; shows examples
439"
440{
441  def @X=basering;
442  int n=nvars(@X);
443  ideal J=jacob(f);
444
445  if(vdim(std(J))<=0)
446  {
447    ERROR("input is not cohomologically tame");
448  }
449
450  int i,j,k,l;
451
452  ideal X=maxideal(1);
453  string c="ring @XS=0,(s,"+varstr(@X)+"),(C,dp(1),dp(n));";
454  execute(c);
455  poly f=imap(@X,f);
456  ideal J=imap(@X,J);
457  ideal JS=std(J+var(1));
458  ideal b0=kbase(JS);
459  int mu=ncols(b0);
460  ideal X=imap(@X,X);
461
462  ideal b;
463  matrix A;
464  module B,B0;
465  ideal K,L,M=1,J,1;
466  ideal K0,L0,M0=X,X,X;
467  module K1,L1,K2,L2;
468  module LL1;
469  for(i=1;i<=deg(f)-1;i++)
470  {
471    M=M,M0;
472    M0=M0*X;
473  }
474
475  ring @S=0,(s,t),(dp,C);
476  number a0;
477  ideal a;
478  int d;
479  ideal e,e0;
480  intvec s,s0;
481  matrix A,B;
482  module V,G;
483
484  while(2*a0!=mu*n)
485  {
486    setring(@XS);
487
488    B=0;
489    while(size(B)<mu||size(B0)<mu||maxdegree(b)+deg(f)>k)
490    {
491      k++;
492      K=K,K0;
493      K0=K0*X;
494
495      B=0;
496      while(size(B)==0||size(B)>mu)
497      {
498        l++;
499        for(i=1;i<=size(L0);i++)
500        {
501          for(j=1;j<=n;j++)
502          {
503            L=L,J[j]*L0[i]-var(1)*diff(L0[i],var(j+1));
504          }
505        }
506        L0=L0*X;
507        M=M,M0;
508        M0=M0*X;
509
510        K1=coeffs(K,K,product(X));
511        L1=std(coeffs(L,M,product(X)));
512        LL1=jet(lead(L1),0);
513        attrib(LL1,"isSB",1);
514        K2=simplify(reduce(K1,LL1),2);
515        L2=intersect(K2,L1);
516
517        B=pidbasis(K2,L2);
518      }
519      B0=std(coeffs(reduce(matrix(K,nrows(K),nrows(B))*B,JS),b0));
520      b=matrix(K,nrows(K),nrows(B))*B;
521    }
522
523    A=lift(B,reduce(coeffs(f*b+var(1)^2*diff(b,var(1)),M,product(X)),L1));
524    d=maxdegree(A);
525    A=chart(A);
526
527    setring(@S);
528
529    e0,s0,V,B=vfilt(imap(@XS,A),d);
530    a,e0,e,s,V,B,G=spec(e0,s0,V,B);
531
532    a0=leadcoef(a[1]);
533    for(i=2;i<=mu;i++)
534    {
535      a0=a0+leadcoef(a[i]);
536    }
537  }
538
539  G,A=glift(fsplit(e0,e,s,V,B,G));
540
541  setring(@XS);
542  b=matrix(b)*imap(@S,G);
543  A=imap(@S,A);
544  export(b,A);
545  kill @S;
546
547  setring(@X);
548  return(@XS);
549}
550example
551{ "EXAMPLE:"; echo=2;
552  ring R=0,(x,y,z),dp;
553  poly f=x+y+z+x2y2z2;
554  def Rs=goodBasis(f);
555  setring(Rs);
556  b;
557  print(jet(A,0));
558  print(jet(A/var(1),0));
559}
560///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.