source: git/Singular/LIB/mondromy.lib @ 4ad2bb

spielwiese
Last change on this file since 4ad2bb was 4ad2bb, checked in by Mathias Schulze <mschulze@…>, 25 years ago
*** empty log message *** git-svn-id: file:///usr/local/Singular/svn/trunk@3161 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.2 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2
3version="$Id: mondromy.lib,v 1.3 1999-06-22 15:38:54 mschulze Exp $";
4info="
5LIBRARY: mondromy.lib  PROCEDURES TO COMPUTE THE MONODROMY OF A SINGULARITY
6AUTHOR:  Mathias Schulze, email: mschulze@mathematik.uni-kl.de
7
8 invunit(u,n);         series inverse of polynomial u up to order n
9 detadj(U);            determinant and adjoint matrix of square matrix U
10 jacoblift(f);         lifts f^kappa in jacob(f) with minimal kappa
11 monodromy(f[,opt]);   monodromy of isolated hypersurface singularity f
12 H''basis(f);          basis of Brieskorn lattice H''
13";
14
15LIB "ring.lib";
16LIB "sing.lib";
17LIB "jordan.lib";
18///////////////////////////////////////////////////////////////////////////////
19
20static proc pcvladdl(list l1,list l2)
21{
22  return(system("pcvLAddL",l1,l2));
23}
24
25static proc pcvpmull(poly p,list l)
26{
27  return(system("pcvPMulL",p,l));
28}
29
30static proc pcvmindeg(list #)
31{
32  return(system("pcvMinDeg",#[1]));
33}
34
35static proc pcvp2cv(list l,int i0,int i1)
36{
37  return(system("pcvP2CV",l,i0,i1));
38}
39
40static proc pcvcv2p(list l,int i0,int i1)
41{
42  return(system("pcvCV2P",l,i0,i1));
43}
44
45static proc pcvdim(int i0,int i1)
46{
47  return(system("pcvDim",i0,i1));
48}
49
50static proc pcvbasis(int i0,int i1)
51{
52  return(system("pcvBasis",i0,i1));
53}
54
55static proc pcvinit()
56{
57  if(system("with","DynamicLoading"))
58  {
59    pcvlal=Pcv::LAddL;
60    pcvpmull=Pcv::PMulL;
61    pcvmindeg=Pcv::MinDeg;
62    pcvp2cv=Pcv::P2CV;
63    pcvcv2p=Pcv::CV2P;
64    pcvdim=Pcv::Dim;
65    pcvbasis=Pcv::Basis;
66  }
67}
68
69///////////////////////////////////////////////////////////////////////////////
70
71static proc min(intvec v)
72{
73  int m=v[1];
74  int i;
75  for(i=2;i<=size(v);i++)
76  {
77    if(m>v[i])
78    {
79      m=v[i];
80    }
81  }
82  return(m);
83}
84///////////////////////////////////////////////////////////////////////////////
85
86static proc max(intvec v)
87{
88  int m=v[1];
89  int i;
90  for(i=2;i<=size(v);i++)
91  {
92    if(m<v[i])
93    {
94      m=v[i];
95    }
96  }
97  return(m);
98}
99///////////////////////////////////////////////////////////////////////////////
100
101static proc mdivp(matrix m,poly p)
102{
103  int i,j;
104  for(i=nrows(m);i>=1;i--)
105  {
106    for(j=ncols(m);j>=1;j--)
107    {
108      m[i,j]=m[i,j]/p;
109    }
110  }
111  return(m); 
112}
113///////////////////////////////////////////////////////////////////////////////
114
115proc codimV(list V,int N)
116{
117  int codim=pcvdim(0,N);
118  if(size(V)>0)
119  {
120    dbprint(printlevel-voice+2,"//vector space dimension: "+string(codim));
121    dbprint(printlevel-voice+2,
122      "//number of subspace generators: "+string(size(V)));
123    int t=timer;
124    codim=codim-ncols(interred(module(V[1..size(V)])));
125    dbprint(printlevel-voice+2,"//codimension: "+string(codim));
126  }
127  return(codim);
128}
129///////////////////////////////////////////////////////////////////////////////
130
131proc quotV(list V,int N)
132{
133  module Q=freemodule(pcvdim(0,N));
134  if(size(V)>0)
135  {
136    dbprint(printlevel-voice+2,"//vector space dimension: "+string(nrows(Q)));
137    dbprint(printlevel-voice+2,
138      "//number of subspace generators: "+string(size(V)));
139    int t=timer;
140    Q=interred(reduce(std(Q),std(module(V[1..size(V)]))));
141  }
142  return(list(Q[1..size(Q)]));
143}
144///////////////////////////////////////////////////////////////////////////////
145
146proc invunit(poly u,int n)
147"USAGE:   invunit(u,n); u poly, n int
148ASSUME:  The polynomial u is a series unit.
149RETURN:  The procedure returns the series inverse of u up to order n
150         or a zero polynomial if u is no series unit.
151DISPLAY: The procedure displays comments if printlevel>=1.
152EXAMPLE: example invunit; shows an example.
153"
154{
155  if(pcvmindeg(u)==0)
156  {
157    def br=basering;
158    changeord("pr","dp");
159    qring qr=std(maxideal(n+1));
160    kill pr;
161
162    dbprint(printlevel-voice+2,"//computing inverse...");
163    int t=timer;
164    poly v=lift(fetch(br,u),1)[1,1];
165    dbprint(printlevel-voice+2,"//...inverse computed ["+string(timer-t)+
166      " secs, "+string(memory(1))+" bytes]");
167
168    setring br;
169    return(fetch(qr,v));
170  }
171  else
172  {
173    print("//no series unit");
174    return(poly(0));
175  }
176}
177example
178{ "EXAMPLE:"; echo=2;
179  ring R=0,x,dp;
180  invunit(1+x,10);
181}
182///////////////////////////////////////////////////////////////////////////////
183
184proc detadj(module U)
185"USAGE:   detadj(U); U matrix
186ASSUME:  U is a square matrix with non zero determinant.
187RETURN:  The procedure returns a list with at most 2 entries.
188         If U is not a sqaure matrix, the list is empty.
189         If U is a sqaure matrix, then the first entry is the determinant of U.
190         If U is a square matrix and the determinant of U not zero,
191         then the second entry is the adjoint matrix of U.
192DISPLAY: The procedure displays comments if printlevel>=1.
193EXAMPLE: example detadj; shows an example.
194"
195{
196  if(nrows(U)==ncols(U))
197  {
198    dbprint(printlevel-voice+2,"//computing determinant...");
199    int t=timer;
200    poly detU=det(U);
201    dbprint(printlevel-voice+2,"//...determinant computed ["+string(timer-t)+
202      " secs, "+string(memory(1))+" bytes]");
203
204    if(detU==0)
205    {
206      print("//determinant zero");
207      return(list(detU));
208    }
209    else
210    {
211      def br=basering;
212      changeord("pr","dp");
213      matrix U=fetch(br,U);
214      poly detU=fetch(br,detU);
215
216      dbprint(printlevel-voice+2,"//computing adjoint matrix...");
217      t=timer;
218      matrix adjU=lift(U,detU*freemodule(nrows(U)));
219      dbprint(printlevel-voice+2,"//...adjoint matrix computed ["
220        +string(timer-t)+" secs, "+string(memory(1))+" bytes]");
221
222      setring br;
223      matrix adjU=fetch(pr,adjU);
224      kill pr;
225    }
226  }
227  else
228  {
229    print("//no square matrix");
230    return(list());
231  }
232
233  return(list(detU,adjU));
234}
235example
236{ "EXAMPLE:"; echo=2;
237  ring R=0,x,dp;
238  matrix U[2][2]=1,1+x,1+x2,1+x3;
239  list daU=detadj(U);
240  daU[1];
241  print(daU[2]);
242}
243///////////////////////////////////////////////////////////////////////////////
244
245proc jacoblift(poly f)
246"USAGE:   jacoblift(f); f poly
247ASSUME:  The polynomial f in a series ring (local ordering) defines
248         an isolated hypersurface singularity.
249RETURN:  The procedure returns a list with entries kappa, xi, u of type
250         int, vector, poly such that kappa is minimal with f^kappa in jacob(f),
251         u is a unit, and u*f^kappa=(matrix(jacob(f))*xi)[1,1].
252DISPLAY: The procedure displays comments if printlevel>=1.
253EXAMPLE: example jacoblift; shows an example.
254"
255{
256  dbprint(printlevel-voice+2,"//computing kappa...");
257  int t=timer;
258  ideal jf=jacob(f);
259  ideal sjf=std(jf);
260  int kappa=1;
261  poly fkappa=f;
262  while(reduce(fkappa,sjf)!=0)
263  {
264    dbprint(printlevel-voice+2,"//kappa="+string(kappa));
265    kappa++;
266    fkappa=fkappa*f;
267  }
268  dbprint(printlevel-voice+2,"//kappa="+string(kappa));
269  dbprint(printlevel-voice+2,"//...kappa computed ["+string(timer-t)+" secs, "
270    +string(memory(1))+" bytes]");
271
272  dbprint(printlevel-voice+2,"//computing xi...");
273  t=timer;
274  vector xi=lift(jf,fkappa)[1];
275  dbprint(printlevel-voice+2,"//...xi computed ["+string(timer-t)+" secs, "
276    +string(memory(1))+" bytes]");
277
278  dbprint(printlevel-voice+2,"//computing u...");
279  t=timer;
280  poly u=(matrix(jf)*xi)[1,1]/fkappa;
281  dbprint(printlevel-voice+2,"//...u computed ["+string(timer-t)+" secs, "
282    +string(memory(1))+" bytes]");
283
284  return(list(kappa,xi,u));
285}
286example
287{ "EXAMPLE:"; echo=2;
288  ring R=0,(x,y),ds;
289  poly f=x2y2+x6+y6;
290  jacoblift(f);
291}
292///////////////////////////////////////////////////////////////////////////////
293
294static proc getdeltaP1(poly f,int K,int N,int dN)
295{
296  return(pcvpmull(f^K,pcvbasis(0,N+dN-K*pcvmindeg(f))));
297}
298///////////////////////////////////////////////////////////////////////////////
299
300static proc getdeltaP2(poly f,int N,int dN)
301{
302  def of,jf=pcvmindeg(f),jacob(f);
303  list b=pcvbasis(N-of+2,N+dN-of+2);
304  list P2;
305  P2[size(b)*((nvars(basering)-1)*nvars(basering))/2]=0;
306  int i,j,k,l;
307  intvec alpha;
308  for(k,l=1,1;k<=size(b);k++)
309  {
310    alpha=leadexp(b[k]);
311    for(i=nvars(basering)-1;i>=1;i--)
312    {
313      for(j=nvars(basering);j>i;j--)
314      {
315        P2[l]=alpha[i]*jf[j]*(b[k]/var(i))-alpha[j]*jf[i]*(b[k]/var(j));
316        l++;
317      }
318    }
319  }
320  return(P2);
321}
322///////////////////////////////////////////////////////////////////////////////
323
324static proc getdeltaPe(poly f,list e,int K,int dK)
325{
326  int k;
327  list Pe,fke;
328  for(k,fke=K,pcvpmull(f^K,e);k<K+dK;k,fke=k+1,pcvpmull(f,fke))
329  {
330    Pe=Pe+fke;
331  }
332  return(Pe);
333}
334///////////////////////////////////////////////////////////////////////////////
335
336static proc incK(poly f,int mu,int K,int deltaK,int N,
337  list e,list P1,list P2,list Pe,list V1,list V2,list Ve)
338{
339  int deltaN=deltaK*pcvmindeg(f);
340
341  list deltaP1;
342  P1=pcvpmull(f^deltaK,P1);
343  V1=pcvp2cv(P1,0,N+deltaN);
344
345  list deltaP2=getdeltaP2(f,N,deltaN);
346  V2=pcvladdl(V2,pcvp2cv(P2,N,N+deltaN))+pcvp2cv(deltaP2,0,N+deltaN);
347  P2=P2+deltaP2;
348
349  list deltaPe=getdeltaPe(f,e,K,deltaK);
350  Ve=pcvladdl(Ve,pcvp2cv(Pe,N,N+deltaN))+pcvp2cv(deltaPe,0,N+deltaN);
351  Pe=Pe+deltaPe;
352
353  K=K+deltaK;
354  dbprint(printlevel-voice+2,"//K="+string(K));
355
356  N=N+deltaN;
357  dbprint(printlevel-voice+2,"//N="+string(N));
358
359  deltaN=1;
360  dbprint(printlevel-voice+2,"//computing codimension of");
361  dbprint(printlevel-voice+2,"//df^dOmega^(n-1)+f^K*Omega^(n+1) in "
362    +"Omega^(n+1) mod m^N*Omega^(n+1)...");
363  int t=timer;
364  while(codimV(V1+V2,N)<K*mu)
365  {
366    dbprint(printlevel-voice+2,"//...codimension computed ["+string(timer-t)
367      +" secs, "+string(memory(1))+" bytes]");
368
369    deltaP1=getdeltaP1(f,K,N,deltaN);
370    V1=pcvladdl(V1,pcvp2cv(P1,N,N+deltaN))+pcvp2cv(deltaP1,0,N+deltaN);
371    P1=P1+deltaP1;
372
373    deltaP2=getdeltaP2(f,N,deltaN);
374    V2=pcvladdl(V2,pcvp2cv(P2,N,N+deltaN))+pcvp2cv(deltaP2,0,N+deltaN);
375    P2=P2+deltaP2;
376
377    Ve=pcvladdl(Ve,pcvp2cv(Pe,N,N+deltaN));
378
379    N=N+deltaN;
380    dbprint(printlevel-voice+2,"//N="+string(N));
381
382    dbprint(printlevel-voice+2,"//computing codimension of");
383    dbprint(printlevel-voice+2,"//df^dOmega^(n-1)+f^K*Omega^(n+1) in "
384      +"Omega^(n+1) mod m^N*Omega^(n+1)...");
385    t=timer;
386  }
387  dbprint(printlevel-voice+2,"//...codimension computed ["+string(timer-t)
388    +" secs, "+string(memory(1))+" bytes]");
389
390  return(K,N,P1,P2,Pe,V1,V2,Ve);
391}
392///////////////////////////////////////////////////////////////////////////////
393
394static proc nablaK(poly f,int kappa,vector xi,poly u,int N,int prevN,
395  list Vnablae,list e)
396{
397  xi=jet(xi,N);
398  u=invunit(u,N);
399  poly fkappa=kappa*f^(kappa-1);
400
401  poly p,q;
402  list nablae;
403  int i,j;
404  for(i=1;i<=size(e);i++)
405  {
406    for(j,p=nvars(basering),0;j>=1;j--)
407    {
408      q=jet(e[i]*xi[j],N);
409      if(q!=0)
410      {
411        p=p+diff(q*jet(u,N-pcvmindeg(q)),var(j));
412      }
413    }
414    nablae=nablae+list(p-jet(fkappa*e[i],N-1));
415  }
416
417  return(pcvladdl(Vnablae,pcvp2cv(nablae,prevN,N-prevN)));
418}
419///////////////////////////////////////////////////////////////////////////////
420
421static proc MK(poly f,int mu,int kappa,vector xi,poly u,
422  int K,int N,int prevN,list e,list V1,list V2,list Ve,list Vnablae)
423{
424  dbprint(printlevel-voice+2,"//computing nabla(e)...");
425  int t=timer;
426  Vnablae=nablaK(f,kappa,xi,u,N,prevN,Vnablae,e);
427  dbprint(printlevel-voice+2,"//...nabla(e) computed ["+string(timer-t)
428    +" secs, "+string(memory(1))+" bytes]");
429
430  dbprint(printlevel-voice+2,
431    "//lifting nabla(e) to C-basis of H''/t^KH''...");
432  list V=Ve+V1+V2;
433  module W=module(V[1..size(V)]);
434  dbprint(printlevel-voice+2,"//vector space dimension: "+string(nrows(W)));
435  dbprint(printlevel-voice+2,"//number of generators: "+string(ncols(W)));
436  t=timer;
437  matrix C=lift(W,module(Vnablae[1..size(Vnablae)]));
438  dbprint(printlevel-voice+2,"//...nabla(e) lifted ["+string(timer-t)
439    +" secs, "+string(memory(1))+" bytes]");
440
441  dbprint(printlevel-voice+2,"//computing e-lift of nabla(e)...");
442  t=timer;
443  int i1,i2,j,k;
444  matrix M[mu][mu];
445  for(j=1;j<=mu;j++)
446  {
447    for(k,i2=0,1;k<K;k++)
448    {
449      for(i1=1;i1<=mu;i1,i2=i1+1,i2+1)
450      {
451        M[i1,j]=M[i1,j]+C[i2,j]*var(1)^k;
452      }
453    }
454  }
455  dbprint(printlevel-voice+2,"//...e-lift of nabla(e) computed ["
456    +string(timer-t)+" secs, "+string(memory(1))+" bytes]");
457
458  return(M,N,Vnablae);
459}
460///////////////////////////////////////////////////////////////////////////////
461
462static proc mid(ideal l)
463{
464  int i,j,id;
465  int mid=0;
466  for(i=size(l);i>=1;i--)
467  {
468    for(j=i-1;j>=1;j--)
469    {
470      id=int(l[i]-l[j]);
471      id=max(intvec(id,-id));
472      mid=max(intvec(id,mid));
473    }
474  }
475  return(mid);
476}
477///////////////////////////////////////////////////////////////////////////////
478
479static proc decmide(matrix M,ideal eM0,list bM0)
480{
481  matrix M0=jet(M,0);
482
483  dbprint(printlevel-voice+2,
484    "//computing basis U of generalized eigenspaces of M0...");
485  int t=timer;
486  int i,j;
487  matrix U,M0e;
488  matrix E=freemodule(nrows(M));
489  for(i=ncols(eM0);i>=1;i--)
490  {
491    M0e=E;
492    for(j=max(bM0[i]);j>=1;j--)
493    {
494      M0e=M0e*(M0-eM0[i]*E);
495    }
496    U=syz(M0e)+U;
497  }
498  dbprint(printlevel-voice+2,"//...U computed ["+string(timer-t)+" secs, "
499    +string(memory(1))+" bytes]");
500
501  dbprint(printlevel-voice+2,"//transforming M to U...");
502  t=timer;
503  list daU=detadj(U);
504  daU[2]=(1/number(daU[1]))*daU[2];
505  M=daU[2]*M*U;
506  dbprint(printlevel-voice+2,"//...M transformed ["+string(timer-t)+" secs, "
507    +string(memory(1))+" bytes]");
508
509  dbprint(printlevel-voice+2,
510    "//computing integer differences of eigenvalues of M0...");
511  t=timer;
512  int k;
513  intvec ideM0;
514  ideM0[ncols(eM0)]=0;
515  for(i=ncols(eM0);i>=1;i--)
516  {
517    for(j=ncols(eM0);j>=1;j--)
518    {
519      k=int(eM0[i]-eM0[j]);
520      if(k)
521      {
522        if(k>0)
523        {
524          ideM0[i]=max(intvec(k,ideM0[i]));
525        }
526        else
527        {
528          ideM0[j]=max(intvec(-k,ideM0[j]));
529        }
530      }
531    }
532  }
533  for(i,k=size(bM0),nrows(M);i>=1;i--)
534  {
535    for(j=sum(bM0[i]);j>=1;j--)
536    {
537      ideM0[k]=ideM0[i];
538      k--;
539    }
540  }
541  dbprint(printlevel-voice+2,
542    "//...integer differences of eigenvalues of M0 computed ["+string(timer-t)
543    +" secs, "+string(memory(1))+" bytes]");
544
545  dbprint(printlevel-voice+2,"//transforming M...");
546  t=timer;
547  for(i=nrows(M);i>=1;i--)
548  {
549    if(!ideM0[i])
550    {
551      M[i,i]=M[i,i]+1;
552    }
553    for(j=ncols(M);j>=1;j--)
554    {
555      if(ideM0[i]&&!ideM0[j])
556      {
557        M[i,j]=M[i,j]*var(1);
558      }
559      else
560      {
561        if(!ideM0[i]&&ideM0[j])
562        {
563          M[i,j]=M[i,j]/var(1);
564        }
565      }
566    }
567  }
568  dbprint(printlevel-voice+2,"//...M transformed ["+string(timer-t)+" secs, "
569    +string(memory(1))+" bytes]");
570
571  return(M);
572}
573///////////////////////////////////////////////////////////////////////////////
574
575static proc nonqhmonodromy(poly f,int mu,int opt)
576{
577  pcvinit();
578
579  dbprint(printlevel-voice+2,"//computing kappa, xi and u with "+
580    "u*f^kappa=(matrix(jacob(f))*xi)[1,1]...");
581  list jl=jacoblift(f);
582  def kappa,xi,u=jl[1..3];
583  dbprint(printlevel-voice+2,"//...kappa, xi and u computed");
584  dbprint(printlevel-voice+2,"//kappa="+string(kappa));
585  if(kappa==1)
586  {
587    dbprint(printlevel-voice+2,
588      "//f quasihomogenous with respect to suitable coordinates");
589  }
590  else
591  {
592    dbprint(printlevel-voice+2,
593      "//f not quasihomogenous for any choice of coordinates");
594  }
595  dbprint(printlevel-voice+2,"//xi=");
596  dbprint(printlevel-voice+2,xi);
597  dbprint(printlevel-voice+2,"//u="+string(u));
598
599  int K,N,prevN;
600  list e,P1,P2,Pe,V1,V2,Ve,Vnablae;
601
602  dbprint(printlevel-voice+2,"//increasing K and N...");
603  K,N,P1,P2,Pe,V1,V2,Ve=incK(f,mu,K,1,N,e,P1,P2,Pe,V1,V2,Ve);
604  dbprint(printlevel-voice+2,"//...K and N increased");
605
606  dbprint(printlevel-voice+2,"//computing C{f}-basis e of Brieskorn lattice "
607    +"H''=Omega^(n+1)/df^dOmega^(n-1)...");
608  int t=timer;
609  e=pcvcv2p(quotV(V1+V2,N),0,N);
610  dbprint(printlevel-voice+2,"//...e computed ["+string(timer-t)+" secs, "
611    +string(memory(1))+" bytes]");
612
613  dbprint(printlevel-voice+2,"//e=");
614  dbprint(printlevel-voice+2,e);
615
616  Pe=e;
617  Ve=pcvp2cv(Pe,0,N);
618
619  if(kappa==1)
620  {
621    dbprint(printlevel-voice+2,
622      "//computing 0-jet M of e-matrix of t*nabla...");
623    matrix M=list(MK(f,mu,kappa,xi,u,K,N,prevN,e,V1,V2,Ve,Vnablae))[1];
624    dbprint(printlevel-voice+2,"//...M computed");
625  }
626  else
627  {
628    dbprint(printlevel-voice+2,
629      "//computing transformation matrix U to simple pole...");
630
631    dbprint(printlevel-voice+2,"//computing t*nabla-stable lattice...");
632    matrix M,prevU;
633    matrix U=freemodule(mu)*var(1)^((mu-1)*(kappa-1));
634    int i;
635    dbprint(printlevel-voice+2,"//comparing with previous lattice...");
636    t=timer;
637    for(i=mu-1;i>=1&&size(reduce(U,std(prevU)))>0;i--)
638    {
639      dbprint(printlevel-voice+2,"//...compared with previous lattice ["
640        +string(timer-t)+" secs, "+string(memory(1))+" bytes]");
641
642      dbprint(printlevel-voice+2,"//increasing K and N...");
643      K,N,P1,P2,Pe,V1,V2,Ve=incK(f,mu,K,kappa-1,N,e,P1,P2,Pe,V1,V2,Ve);
644      dbprint(printlevel-voice+2,"//...K and N increased");
645
646      dbprint(printlevel-voice+2,
647        "//computing (K-1)-jet M of e-matrix of t^kappa*nabla...");
648      M,prevN,Vnablae=MK(f,mu,kappa,xi,u,K,N,prevN,e,V1,V2,Ve,Vnablae);
649      dbprint(printlevel-voice+2,"//...M computed");
650
651      prevU=U;
652
653      dbprint(printlevel-voice+2,"//enlarging lattice...");
654      t=timer;
655      U=interred(jet(module(U)+module(var(1)*diff(U,var(1)))+
656        module(mdivp(M*U,var(1)^(kappa-1))),(kappa-1)*(mu-1)));
657      dbprint(printlevel-voice+2,"//...lattice enlarged ["+string(timer-t)
658        +" secs, "+string(memory(1))+" bytes]");
659
660      dbprint(printlevel-voice+2,"//comparing with previous lattice...");
661      t=timer;
662    }
663    dbprint(printlevel-voice+2,"//...compared with previous lattice ["
664      +string(timer-t)+" secs, "+string(memory(1))+" bytes]");
665    dbprint(printlevel-voice+2,"//...t*nabla-stable lattice computed");
666
667    if(ncols(U)>nrows(U))
668    {
669      dbprint(printlevel-voice+2,
670        "//computing C{f}-basis of t*nabla-stable lattice...");
671      t=timer;
672      U=minbase(U);
673      dbprint(printlevel-voice+2,
674        "//...C{f}-basis of t*nabla-stable lattice computed ["+string(timer-t)
675        +" secs, "+string(memory(1))+" bytes]");
676    }
677
678    U=mdivp(U,var(1)^pcvmindeg(U));
679
680    dbprint(printlevel-voice+2,"//...U computed");
681
682    dbprint(printlevel-voice+2,
683      "//computing determinant and adjoint matrix of U...");
684    list daU=detadj(U);
685    poly p=var(1)^min(intvec(pcvmindeg(daU[2]),pcvmindeg(daU[1])));
686    daU[1]=daU[1]/p;
687    daU[2]=mdivp(daU[2],p);
688    dbprint(printlevel-voice+2,
689      "//...determinant and adjoint matrix of U computed");
690
691    if(K<kappa+pcvmindeg(daU[1]))
692    {
693      dbprint(printlevel-voice+2,"//increasing K and N...");
694      K,N,P1,P2,Pe,V1,V2,Ve=
695        incK(f,mu,K,kappa+pcvmindeg(daU[1])-K,N,e,P1,P2,Pe,V1,V2,Ve);
696      dbprint(printlevel-voice+2,"//...K and N increased");
697
698      dbprint(printlevel-voice+2,"//computing M...");
699      M,prevN,Vnablae=MK(f,mu,kappa,xi,u,K,N,prevN,e,V1,V2,Ve,Vnablae);
700      dbprint(printlevel-voice+2,"//...M computed");
701    }
702
703    dbprint(printlevel-voice+2,"//transforming M/t^kappa to simple pole...");
704    t=timer;
705    M=mdivp(daU[2]*(var(1)^kappa*diff(U,var(1))+M*U),
706      leadcoef(daU[1])*var(1)^(kappa+pcvmindeg(daU[1])-1));
707    dbprint(printlevel-voice+2,"//...M/t^kappa transformed to simple pole ["
708      +string(timer-t)+" secs, "+string(memory(1))+" bytes]");
709  }
710
711  if(opt==0)
712  {
713    dbprint(printlevel-voice+2,
714      "//...computing maximal integer difference delta of eigenvalues of M0");
715    t=timer;
716    list jd=jordan(M);
717    def eM0,bM0=jd[1..2];
718    int delta=mid(eM0);
719    dbprint(printlevel-voice+2,"//...delta computed ["+string(timer-t)
720      +" secs, "+string(memory(1))+" bytes]");
721
722    dbprint(printlevel-voice+2,"//delta="+string(delta));
723
724    if(delta>0)
725    {
726      dbprint(printlevel-voice+2,"//increasing K and N...");
727      if(kappa==1)
728      {
729        K,N,P1,P2,Pe,V1,V2,Ve=incK(f,mu,K,1+delta-K,N,e,P1,P2,Pe,V1,V2,Ve);
730      }
731      else
732      {
733        K,N,P1,P2,Pe,V1,V2,Ve=
734          incK(f,mu,K,kappa+pcvmindeg(daU[1])+delta-K,N,e,P1,P2,Pe,V1,V2,Ve);
735      }
736      dbprint(printlevel-voice+2,"//...K and N increased");
737
738      dbprint(printlevel-voice+2,"//computing M...");
739      M,prevN,Vnablae=MK(f,mu,kappa,xi,u,K,N,prevN,e,V1,V2,Ve,Vnablae);
740      dbprint(printlevel-voice+2,"//...M computed");
741
742      if(kappa>1)
743      {
744        dbprint(printlevel-voice+2,
745          "//transforming M/t^kappa to simple pole...");
746        t=timer;
747        M=mdivp(invunit(daU[1]/var(1)^pcvmindeg(daU[1]),delta)*
748          daU[2]*(var(1)^kappa*diff(U,var(1))+M*U),
749          var(1)^(kappa+pcvmindeg(daU[1])-1));
750        dbprint(printlevel-voice+2,
751          "//...M/t^kappa transformed to simple pole ["+string(timer-t)
752          +" secs, "+string(memory(1))+" bytes]");
753      }
754
755      dbprint(printlevel-voice+2,"//decreasing delta...");
756      M=decmide(M,eM0,bM0);
757      delta--;
758      dbprint(printlevel-voice+2,"//delta="+string(delta));
759
760      while(delta>0)
761      {
762        jd=jordan(M);
763        eM0,bM0=jd[1..2];
764        M=decmide(M,eM0,bM0);
765        delta--;
766        dbprint(printlevel-voice+2,"//delta="+string(delta));
767      }
768      dbprint(printlevel-voice+2,"//...delta decreased");
769    }
770  }
771
772  dbprint(printlevel-voice+2,"//computing 0-jet M0 of M...");
773  matrix M0=jet(M,0);
774  dbprint(printlevel-voice+2,"//...M0 computed");
775
776  return(M0);
777}
778///////////////////////////////////////////////////////////////////////////////
779
780static proc qhmonodromy(poly f,intvec w)
781{
782  dbprint(printlevel-voice+2,"//computing basis e of Milnor algebra...");
783  int t=timer;
784  ideal e=kbase(std(jacob(f)));
785  dbprint(printlevel-voice+2,"//...e computed ["+string(timer-t)+" secs, "
786    +string(memory(1))+" bytes]");
787
788  dbprint(printlevel-voice+2,
789    "//computing Milnor number mu and quasihomogeneous degree d...");
790  int mu,d=size(e),(transpose(leadexp(f))*w)[1];
791  dbprint(printlevel-voice+2,"...mu and d computed");
792
793  dbprint(printlevel-voice+2,"//computing te-matrix M of t*nabla...");
794  matrix M[mu][mu];
795  int i;
796  for(i=mu;i>=1;i--)
797  {
798    M[i,i]=number((transpose(leadexp(e[i])+1)*w)[1])/d;
799  }
800  dbprint(printlevel-voice+2,"//...M computed");
801
802  return(M);
803}
804///////////////////////////////////////////////////////////////////////////////
805
806proc monodromy(poly f, list #)
807"USAGE:   monodromy(f[,opt]); f poly, opt int
808ASSUME:  The polynomial f in a series ring (local ordering) defines
809         an isolated hypersurface singularity.
810RETURN:  The procedure returns a residue matrix M of the meromorphic
811         Gauss-Manin connection of the singularity defined by f
812         or an empty matrix if the assumptions are not fulfilled.
813         If opt=0 (default), exp(2*pi*i*M) is a monodromy matrix of f,
814         else, only the characteristic polynomial of exp(2*pi*i*M) coincides
815         with the characteristic polynomial of the monodromy of f.
816THEORY:  The procedure uses an algorithm by Brieskorn (See E. Brieskorn,
817         manuscipta math. 2 (1970), 103-161) to compute a connection matrix of
818         the meromorphic Gauss-Manin connection up to arbitrarily high order,
819         and an algorithm of Gerard and Levelt (See R. Gerard, A.H.M. Levelt,
820         Ann. Inst. Fourier, Grenoble 23,1 (1973), pp. 157-195) to transform
821         it to a simple pole.
822DISPLAY: The procedure displays more comments for higher printlevel.
823EXAMPLE: example monodromy; shows an example.
824"
825{
826  int opt;
827  if(size(#)>0)
828  {
829    if(typeof(opt)=="int")
830    {
831      opt=#[1];
832    }
833    else
834    {
835      print("\\second parameter no int");
836      return();
837    }
838   
839  }
840
841  dbprint(printlevel-voice+2,"//basering="+string(basering));
842
843  int i;
844  for(i=nvars(basering);i>=1;i--)
845  {
846    if(1<var(i))
847    {
848      i=-1;
849    }
850  }
851
852  if(i<0)
853  {
854    print("//no series ring (local ordering)");
855
856    matrix M[1][0];
857    return(M);
858  }
859  else
860  {
861    dbprint(printlevel-voice+2,"//f="+string(f));
862
863    dbprint(printlevel-voice+2,"//computing milnor number mu of f...");
864    int t=timer;
865    int mu=milnor(f);
866    dbprint(printlevel-voice+2,"//...mu computed ["+string(timer-t)+" secs, "
867      +string(memory(1))+" bytes]");
868
869    dbprint(printlevel-voice+2,"//mu="+string(mu));
870
871    if(mu<=0)
872    {
873      if(mu==0)
874      {
875        print("//no singularity");
876      }
877      else
878      {
879        print("//non isolated singularity");
880      }
881
882      matrix M[1][0];
883      return(M);
884    }
885    else
886    {
887      dbprint(printlevel-voice+2,"//computing weight vector w...");
888      intvec w=qhweight(f);
889      dbprint(printlevel-voice+2,"//...w computed");
890
891      dbprint(printlevel-voice+2,"//w="+string(w));
892
893      if(w==0)
894      {
895        dbprint(printlevel-voice+2,
896          "//f not quasihomogeneous with respect to given coordinates");
897        return(nonqhmonodromy(f,mu,opt));
898      }
899      else
900      {
901        dbprint(printlevel-voice+2,
902          "//f quasihomogeneous with respect to given coordinates");
903        return(qhmonodromy(f,w));
904      }
905    }
906  }
907}
908example
909{ "EXAMPLE:"; echo=2;
910  ring R=0,(x,y),ds;
911  poly f=x2y2+x6+y6;
912  matrix M=monodromy(f);
913  print(M);
914}
915///////////////////////////////////////////////////////////////////////////////
916
917proc H''basis(poly f)
918"USAGE:   H''basis(f); f poly
919ASSUME:  The polynomial f in a series ring (local ordering) defines
920         an isolated hypersurface singularity.
921RETURN:  The procedure returns a list of representatives of a C{f}-basis of the
922         Brieskorn lattice H''=Omega^(n+1)/df^dOmega^(n-1).
923THEORY:  H'' is a free C{f}-module of rank milnor(f).
924DISPLAY: The procedure displays more comments for higher printlevel.
925EXAMPLE: example H''basis; shows an example.
926"
927{
928  pcvinit();
929
930  dbprint(printlevel-voice+2,"//basering="+string(basering));
931
932  int i;
933  for(i=nvars(basering);i>=1;i--)
934  {
935    if(1<var(i))
936    {
937      i=-1;
938    }
939  }
940
941  if(i<0)
942  {
943    print("//no series ring (local ordering)");
944
945    return(list());
946  }
947  else
948  {
949    dbprint(printlevel-voice+2,"//f="+string(f));
950
951    dbprint(printlevel-voice+2,"//computing milnor number mu of f...");
952    int t=timer;
953    int mu=milnor(f);
954    dbprint(printlevel-voice+2,"//...mu computed ["+string(timer-t)+" secs, "
955      +string(memory(1))+" bytes]");
956
957    dbprint(printlevel-voice+2,"//mu="+string(mu));
958
959    if(mu<=0)
960    {
961      if(mu==0)
962      {
963        print("//no singularity");
964      }
965      else
966      {
967        print("//non isolated singularity");
968      }
969
970      return(list());
971    }
972    else
973    {
974      dbprint(printlevel-voice+2,"//computing kappa, xi and u with "+
975        "u*f^kappa=(matrix(jacob(f))*xi)[1,1]...");
976      list jl=jacoblift(f);
977      def kappa,xi,u=jl[1..3];
978      dbprint(printlevel-voice+2,"//...kappa, xi and u computed");
979      dbprint(printlevel-voice+2,"//kappa="+string(kappa));
980      if(kappa==1)
981      {
982        dbprint(printlevel-voice+2,
983          "//f quasihomogenous with respect to suitable coordinates");
984      }
985      else
986      {
987        dbprint(printlevel-voice+2,
988          "//f not quasihomogenous for any choice of coordinates");
989      }
990      dbprint(printlevel-voice+2,"//xi=");
991      dbprint(printlevel-voice+2,xi);
992      dbprint(printlevel-voice+2,"//u="+string(u));
993
994      int K,N,prevN;
995      list e,P1,P2,Pe,V1,V2,Ve,Vnablae;
996
997      dbprint(printlevel-voice+2,"//increasing K and N...");
998      K,N,P1,P2,Pe,V1,V2,Ve=incK(f,mu,K,1,N,e,P1,P2,Pe,V1,V2,Ve);
999      dbprint(printlevel-voice+2,"//...K and N increased");
1000
1001      dbprint(printlevel-voice+2,
1002        "//computing C{f}-basis e of Brieskorn lattice "
1003        +"H''=Omega^(n+1)/df^dOmega^(n-1)...");
1004      t=timer;
1005      e=pcvcv2p(quotV(V1+V2,N),0,N);
1006      dbprint(printlevel-voice+2,"//...e computed ["+string(timer-t)+" secs, "
1007        +string(memory(1))+" bytes]");
1008
1009      dbprint(printlevel-voice+2,"//e=");
1010      dbprint(printlevel-voice+2,e);
1011
1012      return(e);
1013    }
1014  }
1015}
1016example
1017{ "EXAMPLE:"; echo=2;
1018  ring R=0,(x,y),ds;
1019  poly f=x2y2+x6+y6;
1020  H''basis(f);
1021}
1022///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.