source: git/Singular/LIB/spcurve.lib @ 66d68c

spielwiese
Last change on this file since 66d68c was 66d68c, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13499 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 30.8 KB
Line 
1/////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Singularities";
4info="
5LIBRARY: spcurve.lib    Deformations and Invariants of CM-codim 2 Singularities
6AUTHOR:  Anne Fruehbis-Krueger, anne@mathematik.uni-kl.de
7
8PROCEDURES:
9 isCMcod2(i);            presentation matrix of the ideal i, if i is CM
10 CMtype(i);              Cohen-Macaulay type of the ideal i
11 matrixT1(M,n);          1st order deformation T1 in matrix description
12 semiCMcod2(M,T1);       semiuniversal deformation of maximal minors of M
13 discr(sem,n);           discriminant of semiuniversal deformation
14 qhmatrix(M);            weights if M is quasihomogeneous
15 relweight(N,W,a);       relative matrix weight of N w.r.t. weights (W,a)
16 posweight(M,T1,i);      deformation of coker(M) of non-negative weight
17 KSpencerKernel(M);      kernel of the Kodaira-Spencer map
18";
19
20LIB "elim.lib";
21LIB "homolog.lib";
22LIB "inout.lib";
23LIB "poly.lib";
24/////////////////////////////////////////////////////////////////////////////
25
26proc isCMcod2(ideal kurve)
27"USAGE:   isCMcod2(i);   i an ideal
28RETURN:  presentation matrix of i, if i is Cohen-Macaulay of codimension 2 @*
29         a zero matrix otherwise
30EXAMPLE: example isCMcod2; shows an example"
31{
32  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
33//---------------------------------------------------------------------------
34// Compute a minimal free resolution of the ideal and check if the
35// resolution has the expected structure
36//---------------------------------------------------------------------------
37  list kurveres=mres(kurve,0);
38  matrix M=kurveres[2];
39  if ((size(kurveres)>3) &&
40         ((size(kurveres[3])>1) ||
41          ((size(kurveres[3])<=1) && (kurveres[3][1,1]!=0))))
42  {
43    dbprint(p,"//not Cohen-Macaulay, codim 2");
44    matrix ret=0;
45    return(ret);
46  }
47  return(M);
48}
49example
50{ "EXAMPLE:"; echo=2;
51  ring r=32003,(x,y,z),ds;
52  ideal i=xz,yz,x^3-y^4;
53  print(isCMcod2(i));
54}
55/////////////////////////////////////////////////////////////////////////////
56
57proc CMtype(ideal kurve)
58"USAGE:   CMtype(i);  i an ideal, CM of codimension 2
59RETURN:  Cohen-Macaulay type of i (integer)
60         (-1, if i is not Cohen-Macaulay of codimension 2)
61EXAMPLE: example CMtype; shows an example"
62{
63  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
64  int gt = -1;
65//---------------------------------------------------------------------------
66// Compute a minimal free resolution of the ideal and check if the
67// resolution has the expected structure
68//---------------------------------------------------------------------------
69  list kurveres;
70  kurveres=mres(kurve,0);
71  if ((size(kurveres)>3) &&
72         ((size(kurveres[3])>1) ||
73          ((size(kurveres[3])<=1) && (kurveres[3][1,1]!=0))))
74  {
75    dbprint(p,"//not Cohen-Macaulay, codim 2");
76    return(gt);
77  }
78//---------------------------------------------------------------------------
79// Return the Cohen-Macaulay type of i
80//---------------------------------------------------------------------------
81  matrix M = matrix(kurveres[2]);
82  gt = ncols(M);
83  return(gt);
84}
85example
86{ "EXAMPLE:"; echo=2;
87  ring r=32003,(x,y,z),ds;
88  ideal i=xy,xz,yz;
89  CMtype(i);
90}
91/////////////////////////////////////////////////////////////////////////////
92
93proc matrixT1(matrix M ,int n)
94"USAGE:   matrixT1(M,n);  M matrix, n integer
95ASSUME:  M is a presentation matrix of an ideal i, CM of codimension 2;
96         consider i as a family of ideals in a ring in the first n
97         variables where the remaining variables are considered as
98         parameters
99RETURN:  list consisting of the k x (k+1) matrix M and a module K_M such that
100         T1=Mat(k,k+1;R)/K_M is the space of first order deformations of i
101EXAMPLE: example matrixT1; shows an example"
102{
103  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
104//--------------------------------------------------------------------------
105// Initialization and sanity checks
106//--------------------------------------------------------------------------
107  int nr=nrows(M);
108  int nc=ncols(M);
109  if ( nr < nc )
110  {
111    M=transpose(M);
112    int temp=nc;
113    nc=nr;
114    nr=temp;
115    int tra=1;
116  }
117  if ( nr != (nc+1) )
118  {
119    ERROR("not a k x (k+1) matrix");
120  }
121//---------------------------------------------------------------------------
122// Construct the denominator - step by step
123// step 1: initialization
124//---------------------------------------------------------------------------
125  int gt=nc;
126  int i,j;
127  ideal m = M;
128  ideal dx;
129  ideal rv;
130  ideal lv;
131  matrix R[gt][gt]=0;
132  matrix L[gt+1][gt+1]=0;
133  matrix T1[n+gt*gt+(gt+1)*(gt+1)][gt*(gt+1)] = 0;
134//---------------------------------------------------------------------------
135// step 2: the derivatives of the matrix are generators of the denominator
136//---------------------------------------------------------------------------
137  for( i=1; i<= n; i++ )
138  {
139    dx=diff(m,var(i));
140    T1[i,1..gt*(gt+1)] = dx;
141  }
142//---------------------------------------------------------------------------
143// step 3: M*R is a generator as well
144//---------------------------------------------------------------------------
145  for( i=1; i <= gt; i++ )
146  {
147    for ( j=1 ; j <= gt ; j++ )
148    {
149      R[i,j]=1;
150      rv = M * R;
151      T1[n+(i-1)*gt+j,1..gt*(gt+1)] = rv;
152      R[i,j]=0;
153    }
154  }
155//---------------------------------------------------------------------------
156// step 4: so is L*M
157//---------------------------------------------------------------------------
158  for( i=1; i <= (gt+1); i++)
159  {
160    for( j=1 ; j <= (gt+1);j++ )
161    {
162      L[i,j]=1;
163      lv = L * M;
164      T1[n+gt*gt+(i-1)*(gt+1)+j,1..gt*(gt+1)] = lv;
165      L[i,j]=0;
166    }
167  }
168//---------------------------------------------------------------------------
169// Compute the vectorspace basis of T1
170//---------------------------------------------------------------------------
171  module t1 = module(transpose(T1));
172  list result=M,t1;
173  return(result);
174}
175example
176{ "EXAMPLE:"; echo = 2;
177  ring r=32003,(x(1),x(2),x(3)),ds;
178  ideal curve=x(1)*x(2),x(1)*x(3),x(2)*x(3);
179  matrix M=isCMcod2(curve);
180  matrixT1(M,3);
181}
182/////////////////////////////////////////////////////////////////////////////
183
184proc semiCMcod2(matrix M, module t1,list #)
185"USAGE:   semiCMcod2(M,t1[,s]); M matrix, t1 module, s any
186ASSUME:  M is a presentation matrix of an ideal i, CM of codimension 2,
187         and t1 is a presentation of the space of first order deformations
188         of i ((M,t1) as returned by the procedure matrixT1)
189RETURN:  new ring in which the ideal semi describing the semiuniversal
190         deformation of i;
191         if the optional third argument is given, the perturbation matrix
192         of the semiuniversal deformation is returned instead of the ideal.
193NOTE:    The current basering should not contain any variables named
194         A(j) where j is some integer!
195EXAMPLE: example semiCMcod2; shows an example"
196{
197  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
198//---------------------------------------------------------------------------
199// Initialization
200//---------------------------------------------------------------------------
201  module t1erz=kbase(std(t1));
202  int tau=vdim(t1);
203  int gt=ncols(M);
204  int i;
205  def r=basering;
206  if(size(M)!=gt*(gt+1))
207  {
208    gt=gt-1;
209  }
210  for(i=1; i<=size(t1erz); i++)
211  {
212    if(rvar(A(i)))
213    {
214      int jj=-1;
215      break;
216    }
217  }
218  if (defined(jj)>1)
219  {
220    if (jj==-1)
221    {
222      ERROR("Your ring contains a variable T(i)!");
223    }
224  }
225//---------------------------------------------------------------------------
226// Definition of the new ring and the image of M and t1 in the new ring
227//---------------------------------------------------------------------------
228  ring rtemp=0,(A(1..tau)),dp;
229  def rneu=r+rtemp;
230  setring rneu;
231  matrix M=imap(r,M);
232  ideal m=M;
233  module t1erz=imap(r,t1erz);
234//---------------------------------------------------------------------------
235// Construction of the presentation matrix of the versal deformation
236//---------------------------------------------------------------------------
237  matrix N=matrix(m);
238  matrix Mtemp[gt*(gt+1)][1];
239  for( i=1; i<=tau; i++)
240  {
241    Mtemp=t1erz[i];
242    N=N+A(i)*transpose(Mtemp);
243  }
244  ideal n=N;
245  matrix O[gt+1][gt]=n;
246//---------------------------------------------------------------------------
247// Construction of the return value
248//---------------------------------------------------------------------------
249  if(size(#)>0)
250  {
251    matrix semi=O;
252  }
253  else
254  {
255    ideal semi=minor(O,gt);
256  }
257  export semi;
258  return(rneu);
259}
260example
261{ "EXAMPLE:"; echo=2;
262  ring r=32003,(x(1),x(2),x(3)),ds;
263  ideal curve=x(1)*x(2),x(1)*x(3),x(2)*x(3);
264  matrix M=isCMcod2(curve);
265  list l=matrixT1(M,3);
266  def rneu=semiCMcod2(l[1],std(l[2]));
267  setring rneu;
268  semi;
269}
270/////////////////////////////////////////////////////////////////////////////
271
272proc discr(ideal kurve, int n)
273"USAGE:   discr(sem,n);  sem ideal, n integer
274ASSUME:  sem is the versal deformation of an ideal of codimension 2. @*
275         The first n variables of the ring are treated as variables
276         all the others as parameters.
277RETURN:  ideal describing the discriminant
278NOTE:    This is not a powerful algorithm!
279EXAMPLE: example discr; shows an example"
280{
281  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
282//---------------------------------------------------------------------------
283// some sanity checks and initialization
284//---------------------------------------------------------------------------
285  int i;
286  ideal sem=std(kurve);
287  ideal semdiff;
288  ideal J2;
289  int ncol=ncols(matrix(sem));
290  matrix Jacob[n][ncol];
291//---------------------------------------------------------------------------
292// compute the Jacobian matrix
293//---------------------------------------------------------------------------
294  for (i=1; i<=n; i++)
295  {
296    semdiff=diff(sem,var(i));
297    Jacob[i,1..ncol]=semdiff;
298  }
299//---------------------------------------------------------------------------
300// eliminate the first n variables in the ideal generated by
301// the versal deformation and the 2x2 minors of the Jacobian
302//---------------------------------------------------------------------------
303  semdiff=minor(Jacob,2);
304  J2=sem,semdiff;
305  J2=std(J2);
306  poly eli=1;
307  for(i=1; i<=n; i++)
308  {
309    eli=eli*var(i);
310  }
311  ideal dis=eliminate(J2,eli);
312  return(dis);
313}
314example
315{ "EXAMPLE:"; echo=2;
316  ring r=32003,(x(1),x(2),x(3)),ds;
317  ideal curve=x(1)*x(2),x(1)*x(3),x(2)*x(3);
318  matrix M=isCMcod2(curve);
319  list l=matrixT1(M,3);
320  def rneu=semiCMcod2(l[1],std(l[2]));
321  setring rneu;
322  discr(semi,3);
323}
324/////////////////////////////////////////////////////////////////////////////
325
326proc qhmatrix(matrix M)
327"USAGE:   qhmatrix(M);   M a k x (k+1) matrix
328RETURN:  list, consisting of an integer vector containing the weights of
329         the variables of the basering and an integer matrix giving the
330         weights of the entries of M, if M is quasihomogeneous;
331         zero integer vector and zero integer matrix, if M is not
332         quasihomogeneous, i.e. does not allow row and column weights
333EXAMPLE: example qhmatrix; shows an example"
334{
335  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
336//---------------------------------------------------------------------------
337// Initialization and sanity checks
338//---------------------------------------------------------------------------
339  def r=basering;
340  int i,j,temp;
341  int tra=0;
342  int nr=nrows(M);
343  int nc=ncols(M);
344  if ( nr > nc )
345  {
346    M=transpose(M);
347    temp=nc;
348    nc=nr;
349    nr=temp;
350    tra=1;
351  }
352  if ( nc != (nr+1) )
353  {
354    ERROR("not a k x (k+1) matrix");
355  }
356  ideal m=minor(M,nr);
357//---------------------------------------------------------------------------
358// get the weight using the fact that the matrix is quasihomogeneous, if
359// its maximal minors are, and check, whether M is really quasihomogeneous
360//---------------------------------------------------------------------------
361  intvec a=weight(m);
362  string tempstr="ring rneu=" + charstr(r) + ",(" + varstr(r) + "),Ws(" + string(a) + ");";
363  execute(tempstr);
364  def M=imap(r,M);
365  int difset=0;
366  list l;
367  int dif;
368  int donttest=0;
369  int comprow=0;
370  intmat W[nr][nc];
371//---------------------------------------------------------------------------
372// find a row not containing a 0
373//---------------------------------------------------------------------------
374  for(i=1; i<=nr; i++)
375  {
376    if(comprow==0)
377    {
378      comprow=i;
379      for(j=1; j<=nc; j++)
380      {
381        if(M[i,j]==0)
382        {
383          comprow=0;
384          break;
385        }
386      }
387    }
388  }
389//---------------------------------------------------------------------------
390// get the weights of the comprow'th row or use emergency exit
391//---------------------------------------------------------------------------
392  if(comprow==0)
393  {
394    intvec v=0;
395    intmat V=0
396    list ret=v,V;
397    return(ret);
398  }
399  else
400  {
401    for(j=1; j<=nc; j++)
402    {
403      l[j]=deg(lead(M[comprow,j]));
404    }
405  }
406//---------------------------------------------------------------------------
407// do the checks
408//---------------------------------------------------------------------------
409  for(i=1; i<=nr; i++)
410  {
411    if ( i==comprow )
412    {
413// this row should not be tested against itself
414      donttest=1;
415    }
416    else
417    {
418// initialize the difference of the rows, but ignore 0-entries
419      if (M[i,1]!=0)
420      {
421        dif=deg(lead(M[i,1]))-l[1];
422        difset=1;
423      }
424      else
425      {
426        list memo;
427        memo[1]=1;
428      }
429    }
430// check column by column
431    for(j=1; j<=nc; j++)
432    {
433      if(M[i,j]==0)
434      {
435        if(defined(memo)!=0)
436        {
437          memo[size(memo)+1]=j;
438        }
439        else
440        {
441          list memo;
442          memo[1]=j;
443        }
444      }
445      temp=deg(lead(M[i,j]));
446      if((difset!=1) && (donttest!=1) && (M[i,j]!=0))
447      {
448// initialize the difference of the rows, if necessary - still ignore 0s
449        dif=deg(lead(M[i,j]))-l[j];
450        difset=1;
451      }
452// is M[i,j] quasihomogeneous - else emergency exit
453      if(M[i,j]!=jet(M[i,j],temp,a)-jet(M[i,j],temp-1,a))
454      {
455        intvec v=0;
456        intmat V=0;
457        list ret=v,V;
458        return(ret);
459      }
460      if(donttest!=1)
461      {
462// check row and column weights - else emergency exit
463        if(((temp-l[j])!=dif) && (M[i,j]!=0) && (difset==1))
464        {
465          intvec v=0;
466          intmat V=0;
467          list ret=v,V;
468          return(ret);
469        }
470      }
471// set the weight matrix entry
472      W[i,j]=temp;
473    }
474// clean up the 0's we left out
475    if((difset==1) && (defined(memo)!=0))
476    {
477      for(j=1; j<=size(memo); j++)
478      {
479        W[i,memo[j]]=dif+l[memo[j]];
480      }
481      kill memo;
482    }
483    donttest=0;
484  }
485//---------------------------------------------------------------------------
486// transpose, if M was transposed during initialization, and return the list
487//---------------------------------------------------------------------------
488  if ( tra==1 )
489  {
490    W=transpose(W);
491  }
492  setring r;
493  list ret=a,W;
494  return(ret);
495}
496example
497{ "EXAMPLE:"; echo=2;
498  ring r=0,(x,y,z),ds;
499  matrix M[3][2]=z,0,y,x,x^3,y;
500  qhmatrix(M);
501  pmat(M);
502}
503/////////////////////////////////////////////////////////////////////////////
504
505proc relweight(matrix N, intmat W, intvec a)
506"USAGE:   relweight(N,W,a); N matrix, W intmat, a intvec
507ASSUME:  N is a non-zero matrix
508         W is an integer matrix of the same size as N
509         a is an integer vector giving the weights of the variables
510RETURN:  integer, max(a-weighted order(N_ij) - W_ij | all entries ij) @*
511         string \"ERROR\" if sizes do not match
512EXAMPLE: example relweight; shows an example
513"
514{
515  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
516//---------------------------------------------------------------------------
517// Initialization and sanity checks
518//---------------------------------------------------------------------------
519  if ((size(N)!=size(W)) || (ncols(N)!=ncols(W)))
520  {
521    ERROR("matrix size does not match");
522  }
523  if (size(a)!=nvars(basering))
524  {
525    ERROR("length of weight vector != number of variables");
526  }
527  int i,j,temp;
528  def r=basering;
529//---------------------------------------------------------------------------
530// Comparision entry by entry
531//---------------------------------------------------------------------------
532  for(i=1; i<=nrows(N); i++)
533  {
534    for(j=1; j<=ncols(N); j++)
535    {
536      if (N[i,j]!=0)
537      {
538        temp=mindeg1(N[i,j],a)-W[i,j];
539        if (defined(ret))
540        {
541          if(temp > ret)
542          {
543            ret=temp;
544          }
545        }
546        else
547        {
548          int ret=temp;
549        }
550      }
551    }
552  }
553  return(ret);
554}
555example
556{ "EXAMPLE:"; echo=2;
557  ring r=32003,(x,y,z),ds;
558  matrix N[2][3]=z,0,y,x,x^3,y;
559  intmat W[2][3]=1,1,1,1,1,1;
560  intvec a=1,1,1;
561  relweight(N,W,a);
562}
563/////////////////////////////////////////////////////////////////////////////
564
565proc posweight(matrix M, module t1, int choose, list #)
566"USAGE:   posweight(M,t1,n[,s]); M matrix, t1 module, n int, s string @*
567         n=0 : all deformations of non-negative weight @*
568         n=1 : only non-constant deformations of non-negative weight @*
569         n=2 : all deformations of positive weight @*
570ASSUME:  M is a presentation matrix of a Cohen-Macaulay codimension 2
571         ideal and t1 is its T1 space in matrix notation
572RETURN:  new ring containing a list posw, consisting of a presentation
573         matrix describing the deformation given by the generators of T1
574         of non-negative/positive weight and the weight vector for the new
575         variables
576NOTE:   The current basering should not contain any variables named
577         T(i) where i is some integer!
578EXAMPLE: example posweight; shows an example"
579{
580//---------------------------------------------------------------------------
581// Initialization and sanity checks
582//---------------------------------------------------------------------------
583  if (size(#)>0)
584  {
585    if (typeof(#[1])=="string")
586    {
587      string newname=#[1];
588    }
589  }
590  if (attrib(t1,"isSB"))
591  {
592    module t1erz=kbase(t1);
593    int tau=vdim(t1);
594  }
595  else
596  { module t1erz=kbase(std(t1));
597    int tau=vdim(std(t1));
598  }
599  for(int i=1; i<=size(t1erz); i++)
600  {
601    if(rvar(T(i)))
602    {
603      int jj=-1;
604      break;
605    }
606  }
607  kill i;
608  if (defined(jj))
609  {
610    if (jj==-1)
611    {
612      ERROR("Your ring contains a variable T(i)!");
613    }
614  }
615  int pw=0;
616  int i;
617  def r=basering;
618  list l=qhmatrix(M);
619  int gt=ncols(M);
620  if(size(M)!=gt*(gt+1))
621  {
622    gt=gt-1;
623  }
624  matrix erzmat[gt+1][gt];
625  list erz;
626  if ((size(l[1])==1) && (l[1][1]==0) && (size(l[2])==1) && (l[2][1,1]==0))
627  {
628    ERROR("Internal Error: Problem determining the weights.");
629  }
630//---------------------------------------------------------------------------
631// Find the generators of T1 of non-negative weight
632//---------------------------------------------------------------------------
633  int relw;
634  list rlw;
635  for(i=1; i<=tau; i++)
636  {
637    erzmat=t1erz[i];
638    kill relw;
639    def relw=relweight(erzmat,l[2],l[1]);
640    if(typeof(relw)=="int")
641    {
642      if (((choose==0) && (relw>=0))
643           || ((choose==1) && (relw>=0) && (CMtype(minor(M+erzmat,gt))==gt))
644           || ((choose==2) && (relw > 0)))
645      {
646        pw++;
647        rlw[pw]=relw;
648        erz[pw]=erzmat;
649      }
650    }
651    else
652    {
653      ERROR("Internal Error: Problem determining relative weight.");
654    }
655  }
656//---------------------------------------------------------------------------
657// Definition of the new ring and the image of M and erz in the new ring
658//---------------------------------------------------------------------------
659  if(size(rlw)==0)
660  {
661    ERROR("Internal Error: Problem determining relative weight.");
662  }
663  intvec iv=rlw[1..size(rlw)];
664  ring rtemp=0,(T(1..pw)),dp;
665  def rneu=r+rtemp;
666  setring rneu;
667  matrix M=imap(r,M);
668  ideal m=M;
669// we cannot imap erz, if its size=0
670  if(pw==0)
671  {
672    list erz1;
673  }
674  else
675  {
676    list erz1=imap(r,erz);
677  }
678//---------------------------------------------------------------------------
679// Construction of the presentation matrix of the deformation
680//---------------------------------------------------------------------------
681  matrix N=matrix(m);
682  ideal mtemp;
683  matrix Mtemp[gt*(gt+1)][1];
684  for( i=1; i<=pw; i++)
685  {
686    mtemp=erz1[i];
687    Mtemp=mtemp;
688    N=N+T(i)*transpose(Mtemp);
689  }
690  ideal n=N;
691  matrix O[gt+1][gt]=n;
692//---------------------------------------------------------------------------
693// Keep the matrix and return the ring in which it lives
694//---------------------------------------------------------------------------
695  list posw=O,iv;
696  export posw;
697  return(rneu);
698}
699example
700{ "EXAMPLE:"; echo=2;
701  ring r=32003,(x(1),x(2),x(3)),ds;
702  ideal curve=(x(3)-x(1)^2)*x(3),(x(3)-x(1)^2)*x(2),x(2)^2-x(1)^7*x(3);
703  matrix M=isCMcod2(curve);
704  list l=matrixT1(M,3);
705  def rneu=posweight(l[1],std(l[2]),0);
706  setring rneu;
707  pmat(posw[1]);
708  posw[2];
709}
710/////////////////////////////////////////////////////////////////////////////
711
712proc KSpencerKernel(matrix M,list #)
713"USAGE:     KSpencerKernel(M[,s][,v]);  M matrix, s string, v intvec @*
714           optional parameters (please specify in this order, if both are
715           present):
716           *  s = first of the names of the new rings
717              e.g. \"R\" leads to ring names R and R1
718           *  v of size n(n+1) leads to the following module ordering @*
719              gen(v[1]) > gen(v[2]) > ... > gen(v[n(n+1)]) where the matrix
720              entry ij corresponds to gen((i-1)*n+j)
721ASSUME:    M is a quasihomogeneous n x (n+1) matrix where the n minors define
722           an isolated space curve singularity
723RETURN:    new ring containing the coefficient matrix KS representing
724           the kernel of the Kodaira-Spencer map of the family of
725           non-negative deformations having the given singularity as
726           special fibre
727NOTE:      * the initial basering should not contain variables with name
728             e(i) or T(i), since those variable names will internally be
729             used by the script
730           * setting an intvec with 5 entries and name watchProgress
731             shows the progress of the computations: @*
732             watchProgress[1]>0  => option(prot) in groebner commands @*
733             watchProgress[2]>0  => trace output for highcorner @*
734             watchProgress[3]>0  => output of deformed matrix @*
735             watchProgress[4]>0  => result of elimination step @*
736             watchProgress[4]>1  => trace output of multiplications with xyz
737                                    and subsequent reductions @*
738             watchProgress[5]>0  => matrix representing the kernel using print
739EXAMPLE:   example KSpencerKernel; shows an example"
740{
741  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
742//---------------------------------------------------------------------------
743// Initialization and sanity checks
744//---------------------------------------------------------------------------
745  intvec optvec=option(get);
746  if (size(#)>0)
747  {
748    if (typeof(#[1])=="string")
749    {
750      string newname=#[1];
751    }
752    if (typeof(#[1])=="intvec")
753    {
754      intvec desiredorder=#[1];
755    }
756    if (size(#)>1)
757    {
758      if (typeof(#[2])=="intvec")
759      {
760        intvec desiredorder=#[2];
761      }
762    }
763  }
764  if (defined(watchProgress))
765  {
766    if ((typeof(watchProgress)!="intvec") || (size(watchProgress)<5))
767    {
768      "watchProgress should be an intvec with at least 5 entries";
769      "ignoring watchProgress";
770      def kksave=watchProgress;
771      kill watchProgress;
772    }
773  }
774  option(redTail);
775  if (nvars(basering) != 3 )
776  {
777    ERROR("It should be a curve in 3 space");
778  }
779//---------------------------------------------------------------------------
780// change to a basering with the correct weihted order
781//---------------------------------------------------------------------------
782  def rt=basering;
783  list wl=qhmatrix(M);
784  if ((size(wl)!=2) || ((wl[1]==0) && (wl[2]==0)))
785  {
786    ERROR("The matrix was not n x (n+1) or not quasihomogenous");
787  }
788  string ringre=" ring r=" + charstr(rt) + ",(x,y,z), Ws(" + string(wl[1]) + ");";
789  execute(ringre);
790  matrix M=imap(rt,M);
791  int ne=size(M);
792  if (defined(desiredorder)>1)
793  {
794    intvec iv;
795    for(int i=1;i<=size(desiredorder);i++)
796    {
797      iv[desiredorder[i]]=i;
798    }
799  }
800  else
801  {
802    intvec iv=1..ne;
803  }
804  list l=matrixT1(M,3);
805  if (dim(std(l[2])) != 0)
806  {
807    ERROR("The matrix does not define an isolated space curve singularity");
808  }
809  module t1qh=l[2];
810//--------------------------------------------------------------------------
811// Passing to a new ring with extra variables e(i) corresponding to
812// the module generators gen(i) for weighted standard basis computation
813// accepting weights for the gen(i)
814//--------------------------------------------------------------------------
815  int jj=0;
816  for(int i=1; i<=ne; i++)
817  {
818    if(rvar(e(i)))
819    {
820      jj=-1;
821    }
822  }
823  if (jj==-1)
824  {
825    ERROR("Your ring contains a variable e(i)!");
826  }
827  if(defined(desiredorder)>1)
828  {
829    ringre="ring re=" + charstr(r) +",(e(1.." + string(ne) + "),"+
830                        varstr(basering) + "),Ws(";
831    intvec tempiv=intvec(wl[2]);
832    for(i=1;i<=ne;i++)
833    {
834      ringre=ringre + string((-1)*tempiv[desiredorder[i]]) + ",";
835    }
836    ringre= ringre  + string(wl[1]) + ");";
837  }
838  else
839  {
840    ringre="ring re=" + charstr(r) +",(e(1.." + string(ne) + "),"+ varstr(basering)
841                        + "),Ws(" + string((-1)*intvec(wl[2])) + ","
842                        + string(wl[1]) + ");";
843  }
844  execute(ringre);
845  module temp=imap(r,t1qh);
846  ideal t1qh=mod2id(temp,iv);
847  if (defined(watchProgress))
848  {
849    if (watchProgress[1]!=0)
850    {
851      option(prot);
852      "Protocol output of the groebner computation (quasihomogenous case)";
853    }
854  }
855  ideal t1qhs=std(t1qh);
856  if (defined(watchProgress))
857  {
858    if (watchProgress[1]!=0)
859    {
860      "groebner computation finished";
861      option(noprot);
862    }
863  }
864  ideal t1qhsl=lead(t1qhs);
865  module mo=id2mod(t1qhsl,iv);
866//--------------------------------------------------------------------------
867// Return to the initial ring to compute the kbase and noether there
868// (in the new ring t1qh is of course not of dimension 0 but of dimension 3
869// so we have to go back)
870//--------------------------------------------------------------------------
871  setring r;
872  module mo=imap(re,mo);
873  attrib(mo,"isSB",1);                // mo is monomial ==> SB
874  attrib(mo,"isHomog",intvec(wl[2])); // highcorner has to respect the weights
875  vector noe=highcorner(mo);
876  if (defined(watchProgress))
877  {
878    if (watchProgress[2]!=0)
879    {
880      "weights corresponding to the entries of the matrix:";
881      wl;
882      "leading term of the groebner basis (quasihomogeneous case)";
883      mo;
884      "noether";
885      noe;
886    }
887  }
888//--------------------------------------------------------------------------
889// Define the family of curves with the same quasihomogeneous initial
890// matrix M, compute T1 and pass again to the ring with the variables e(i)
891//--------------------------------------------------------------------------
892  def rneu=posweight(M,mo,2);
893  setring rneu;
894  list li=posw;
895  if (size(li)<=1)
896  {
897    ERROR("Internal Error: Problem determining perturbations of weight > 0.")
898  }
899  if (defined(watchProgress))
900  {
901    if(watchProgress[3]!=0)
902    {
903      "perturbed matrix and weights of the perturbations:";
904      li;
905    }
906  }
907  list li2=matrixT1(li[1],3);
908  module Mpert=transpose(matrix(ideal(li2[1])));
909  module t1pert=li2[2];
910  int nv=nvars(rneu)-nvars(r);
911  ring rtemp=0,(T(1..nv)),wp(li[2]);
912  def reneu=re+rtemp;
913  setring reneu;
914  module noe=matrix(imap(r,noe));
915  ideal noet=mod2id(noe,iv);
916  module temp=imap(rneu,t1pert);
917  ideal t1pert=mod2id(temp,iv);
918//--------------------------------------------------------------------------
919// Compute the standard basis and select those generators with leading term
920// divisible by some T(i)
921//--------------------------------------------------------------------------
922  noether=noet[size(noet)];
923  if (defined(watchProgress))
924  {
925    if (watchProgress[1]!=0)
926    {
927      "protocol output of the groebner command (perturbed case)";
928      option(prot);
929    }
930  }
931  ideal t1perts=std(t1pert);
932  noether=noet[size(noet)];
933  t1perts=interred(t1perts);
934  if (defined(Debug))
935  {
936    if (watchProgress[1]!=0)
937    {
938      "groebner computation finished (perturbed case)";
939      option(noprot);
940    }
941  }
942  ideal templ=lead(t1perts);
943  for(int j=1;j<=nv;j++)
944  {
945    templ=subst(templ,T(j),0);
946  }
947  ideal mx;
948  ideal mt;
949  for(j=1;j<=size(t1perts);j++)
950  {
951    if(templ[j]!=0)
952    {
953      mx=mx,t1perts[j];
954    }
955    else
956    {
957      mt=mt,t1perts[j];
958    }
959  }
960//--------------------------------------------------------------------------
961// multiply by the initial ring variables to shift the generators with
962// leading term divisible by some T(i) and reduce afterwards
963//--------------------------------------------------------------------------
964                       // This is obviously no SB, but we have to reduce by
965  attrib(mx,"isSB",1); // it and setting isSB suppresses error messages
966  noether=noet[size(noet)];
967  ideal ker_gen=reduce(mt,mx);
968  ideal ovar=var(ne+1),var(ne+2),var(ne+3);
969  j=1;
970  noether=noet[size(noet)];
971  if (defined(watchProgress))
972  {
973    if (watchProgress[4]!=0)
974    {
975      "generators of the kernel as a C[T]{x} module:";
976      mt;
977      "noether:";
978      noether;
979    }
980  }
981  int zeros;
982  templ=ker_gen;
983  while(zeros==0)
984  {
985    zeros=1;
986    templ=templ*ovar;
987    templ=reduce(templ,mx);
988    if(defined(watchProgress))
989    {
990      if(watchProgress[4]>1)
991      {
992        templ;
993      }
994    }
995    if (size(templ)!= 0)
996    {
997      zeros=0;
998      ker_gen=ker_gen,templ;
999    }
1000  }
1001//-------------------------------------------------------------------------
1002// kill zero entries, keep only one of identical entries
1003//-------------------------------------------------------------------------
1004  ovar=var(1);
1005  for(i=2;i<=ne;i++)
1006  {
1007    ovar=ovar,var(i);
1008  }
1009  ker_gen=ker_gen,ovar^2;
1010  noether=noet[size(noet)];
1011  ker_gen=simplify(ker_gen,10);
1012//-------------------------------------------------------------------------
1013// interreduce ker_gen as a k[T]-module
1014//-------------------------------------------------------------------------
1015  intvec mgen=1..(ne+3);
1016  ideal Mpert=mod2id(imap(rneu,Mpert),iv);
1017  templ=0;
1018  for(i=1;i<=nv;i++)
1019  {
1020    templ[i]=diff(Mpert[size(Mpert)],T(i));
1021  }
1022  templ=templ,ovar^2;
1023  list retl=subrInterred(templ,ker_gen,mgen);
1024// Build up the matrix representing L
1025  module retlm=transpose(retl[2]);
1026  for(i=1;i<=size(retl[1]);i++)
1027  {
1028    if(reduce(retl[1][1,i],std(ovar^2))==0)
1029    {
1030      retlm[i]=0;
1031    }
1032  }
1033  retlm=simplify(transpose(simplify(transpose(retlm),10)),10);
1034  if(defined(watchProgress))
1035  {
1036    if(watchProgress[5]>0)
1037    {
1038      print(retlm);
1039    }
1040  }
1041  ker_gen=retl[3];
1042// we define ret=i(L),(delta_j(t_k))_jk
1043  list ret=id2mod(ker_gen,iv),matrix(retlm);
1044// cleanups - define what we previously killed
1045  if(defined(kksave)>1)
1046  {
1047    def watchProgress=kksave;
1048    export watch Progress;
1049  }
1050  option(set,optvec);
1051  def KS=ret[2];
1052  export KS;
1053  return(reneu);
1054}
1055example
1056{ "EXAMPLE:"; echo=2;
1057  ring r=0,(x,y,z),ds;
1058  matrix M[3][2]=z-x^7,0,y^2,z,x^9,y;
1059  def rneu=KSpencerKernel(M,"ar");
1060  setring rneu;
1061  basering;
1062  print(KS);
1063}
1064///////////////////////////////////////////////////////////////////////////
1065
Note: See TracBrowser for help on using the repository browser.