source: git/Singular/LIB/homolog.lib @ a890be

fieker-DuValspielwiese
Last change on this file since a890be was b54a393, checked in by Hans Schönemann <hannes@…>, 14 years ago
*pfister: sagbi fixes git-svn-id: file:///usr/local/Singular/svn/trunk@12381 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 60.1 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Commutative Algebra";
4info="
5LIBRARY:  homolog.lib   Procedures for Homological Algebra
6AUTHORS:  Gert-Martin Greuel, greuel@mathematik.uni-kl.de,
7@*        Bernd Martin,  martin@math.tu-cottbus.de
8@*        Christoph Lossen,  lossen@mathematik.uni-kl.de
9
10PROCEDURES:
11 canonMap(id);          the kernel and the cokernel of the canonical map
12 cup(M);                cup: Ext^1(M',M') x Ext^1() --> Ext^2()
13 cupproduct(M,N,P,p,q); cup: Ext^p(M',N') x Ext^q(N',P') --> Ext^p+q(M',P')
14 depth(I,M);            depth(I,M'),    I ideal, M module, M'=coker(M)
15 Ext_R(k,M);            Ext^k(M',R),    M module, R basering, M'=coker(M)
16 Ext(k,M,N);            Ext^k(M',N'),   M,N modules, M'=coker(M), N'=coker(N)
17 fitting(M,n);          n-th Fitting ideal of M'=coker(M), M module, n int
18 flatteningStrat(M);    Flattening stratification of M'=coker(M), M module
19 Hom(M,N);              Hom(M',N'),     M,N modules, M'=coker(M), N'=coker(N)
20 homology(A,B,M,N);     ker(B)/im(A),   homology of complex R^k--A->M'--B->N'
21 isCM(M);               test if coker(M) is Cohen-Macaulay, M module
22 isFlat(M);             test if coker(M) is flat, M module
23 isLocallyFree(M,r);    test if coker(M) is locally free of constant rank r
24 isReg(I,M);            test if I is coker(M)-sequence, I ideal, M module
25 hom_kernel(A,M,N);     ker(M'--A->N')  M,N modules, A matrix
26 kohom(A,k);            Hom(R^k,A),     A matrix over basering R
27 kontrahom(A,k);        Hom(A,R^k),     A matrix over basering R
28 KoszulHomology(I,M,n); n-th Koszul homology H_n(I,coker(M)), I=ideal
29 tensorMod(M,N);        Tensor product of modules M'=coker(M), N'=coker(N)
30 Tor(k,M,N);            Tor_k(M',N'),   M,N modules, M'=coker(M), N'=coker(N)
31";
32
33LIB "general.lib";
34LIB "deform.lib";
35LIB "matrix.lib";
36LIB "poly.lib";
37LIB "primdec.lib";
38///////////////////////////////////////////////////////////////////////////////
39
40proc cup (module M,list #)
41"USAGE:   cup(M,[,any,any]);  M=module
42COMPUTE: cup-product Ext^1(M',M') x Ext^1(M',M') ---> Ext^2(M',M'), where
43         M':=R^m/M, if M in R^m, R basering (i.e. M':=coker(matrix(M))).
44@*       If called with >= 2 arguments: compute symmetrized cup-product
45ASSUME:  all Ext's  are finite dimensional
46RETURN:  - if called with 1 argument: matrix, the columns of the output present
47         the coordinates of b_i&b_j with respect to a kbase of Ext^2, where
48         b_1,b_2,... is a kbase of Ext^1 and & denotes cup product;@*
49         - if called with 2 arguments: matrix, the columns of the output
50         present the coordinates of (1/2)(b_i&b_j + b_j&b_i) with respect to
51         a kbase of Ext^2;
52         - if called with 3 arguments: list,
53@format
54      L[1] = matrix see above (symmetric case, for >=2 arguments)
55      L[2] = matrix of kbase of Ext^1
56      L[3] = matrix of kbase of Ext^2
57@end format
58NOTE:    printlevel >=1;  shows what is going on.
59         printlevel >=2;  shows result in another representation.
60@*       For computing cupproduct of M itself, apply proc to syz(M)!
61EXAMPLE: example cup; shows examples
62"
63{
64//---------- initialization ---------------------------------------------------
65   int    i,j,k,f0,f1,f2,f3,e1,e2;
66   module M1,M2,A,B,C,ker,ima,ext1,ext2,ext10,ext20;
67   matrix cup[1][0];
68   matrix kb1,lift1,kb2,mA,mB,mC;
69   ideal  tes1,tes2,null;
70   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
71//-----------------------------------------------------------------------------
72//take a resolution of M<--F(0)<--- ...  <---F(3)
73//apply Hom(-,M) and compute the Ext's
74//-----------------------------------------------------------------------------
75   list resM = nres(M,3);
76   M1 = resM[2];
77   M2 = resM[3];
78   f0 = nrows(M);
79   f1 = ncols(M);
80   f2 = ncols(M1);
81   f3 = ncols(M2);
82   tes1 = simplify(ideal(M),10);
83   tes2=simplify(ideal(M1),10);
84   if ((tes1[1]*tes2[1]==0) or (tes1[1]==1) or (tes2[1]==1))
85   {
86      dbprint(p,"// Ext == 0 , hence 'cup' is the zero-map");
87      return(@cup);
88   }
89//------ compute Ext^1 --------------------------------------------------------
90   B     = kohom(M,f2);
91   A     = kontrahom(M1,f0);
92   C     = intersect(A,B);
93   C     = reduce(C,std(null));C = simplify(C,10);
94   ker   = lift(A,C)+syz(A);
95   ima   = kohom(M,f1);
96   ima   = ima + kontrahom(M,f0);
97   ext1  = modulo(ker,ima);
98   ext10 = std(ext1);
99   e1    = vdim(ext10);
100   dbprint(p-1,"// vdim (Ext^1) = "+string(e1));
101   if (e1 < 0)
102   {
103     "// Ext^1 not of finite dimension";
104     return(cup);
105   }
106   kb1 = kbase(ext10);
107   kb1 = matrix(ker)*kb1;
108   dbprint(p-1,"// kbase of Ext^1(M,M)",
109              "//  - the columns present the kbase elements in Hom(F(1),F(0))",
110              "//  - F(*) a free resolution of M",kb1);
111
112//------ compute the liftings of Ext^1 ----------------------------------------
113   mC = matrix(A)*kb1;
114   lift1 =lift(B,mC);
115   dbprint(p-1,"// lift kbase of Ext^1:",
116    "//  - the columns present liftings of kbase elements into Hom(F(2),F(1))",
117    "//  - F(*) a free resolution of M ",lift1);
118
119//------ compute Ext^2  -------------------------------------------------------
120   B    = kohom(M,f3);
121   A    = kontrahom(M2,f0);
122   C    = intersect(A,B);
123   C    = reduce(C,std(null));C = simplify(C,10);
124   ker  = lift(A,C)+syz(A);
125   ima  = kohom(M,f2);
126   ima  = ima + kontrahom(M1,f0);
127   ext2 = modulo(ker,ima);
128   ext20= std(ext2);
129   e2   = vdim(ext20);
130   if (e2<0)
131   {
132     "// Ext^2 not of finite dimension";
133     return(cup);
134   }
135   dbprint(p-1,"// vdim (Ext^2) = "+string(e2));
136   kb2 = kbase(ext20);
137   kb2 = matrix(ker)*kb2;
138   dbprint(p-1,"// kbase of Ext^2(M,M)",
139             "//  - the columns present the kbase elements in Hom(F(2),F(0))",
140             "//  - F(*) is a free resolution of M ",kb2);
141//-------  compute: cup-products of base-elements -----------------------------
142   for (i=1;i<=e1;i=i+1)
143   {
144     for (j=1;j<=e1;j=j+1)
145     {
146       mA = matrix(ideal(lift1[j]),f1,f2);
147       mB = matrix(ideal(kb1[i]),f0,f1);
148       mC = mB*mA;
149       if (size(#)==0)
150       {                                          //non symmestric
151         mC = matrix(ideal(mC),f0*f2,1);
152         cup= concat(cup,mC);
153       }
154       else                                       //symmetric version
155       {
156         if (j>=i)
157         {
158           if (j>i)
159           {
160             mA = matrix(ideal(lift1[i]),f1,f2);
161             mB = matrix(ideal(kb1[j]),f0,f1);
162             mC = mC+mB*mA;mC=(1/2)*mC;
163           }
164           mC = matrix(ideal(mC),f0*f2,1);
165           cup= concat(cup,mC);
166         }
167       }
168     }
169   }
170   dbprint(p-1,"// matrix of cup-products (in Ext^2)",cup,
171               "////// end level 2 //////");
172//------- compute: presentation of base-elements -----------------------------
173   cup = lift(ker,cup);
174   cup = lift_kbase(cup,ext20);
175   if( p>2 )
176   {
177     "// the associated matrices of the bilinear mapping 'cup' ";
178     "// corresponding to the kbase elements of Ext^2(M,M) are shown,";
179     "//  i.e. the rows of the final matrix are written as matrix of";
180     "//  a bilinear form on Ext^1 x Ext^1";
181     matrix BL[e1][e1];
182     for (k=1;k<=e2;k=k+1)
183     {
184       "//-----component "+string(k)+":";
185       for (i=1;i<=e1;i=i+1)
186       {
187         for (j=1;j<=e1;j=j+1)
188         {
189           if (size(#)==0) { BL[i,j]=cup[k,j+e1*(i-1)]; }
190           else
191           {
192             if (i<=j)
193             {
194               BL[i,j]=cup[k,j+e1*(i-1)-binomial(i,2)];
195               BL[j,i]=BL[i,j];
196             }
197           }
198         }
199       }
200       print(BL);
201     }
202     "////// end level 3 //////";
203   }
204   if (size(#)>2) { return(cup,kb1,kb2);}
205   else {return(cup);}
206}
207example
208{"EXAMPLE";  echo=2;
209  int p      = printlevel;
210  ring  rr   = 32003,(x,y,z),(dp,C);
211  ideal  I   = x4+y3+z2;
212  qring  o   = std(I);
213  module M   = [x,y,0,z],[y2,-x3,z,0],[z,0,-y,-x3],[0,z,x,-y2];
214  print(cup(M));
215  print(cup(M,1));
216  // 2nd EXAMPLE  (shows what is going on)
217  printlevel = 3;
218  ring   r   = 0,(x,y),(dp,C);
219  ideal  i   = x2-y3;
220  qring  q   = std(i);
221  module M   = [-x,y],[-y2,x];
222  print(cup(M));
223  printlevel = p;
224}
225///////////////////////////////////////////////////////////////////////////////
226
227proc cupproduct (module M,N,P,int p,q,list #)
228"USAGE:   cupproduct(M,N,P,p,q[,any]);  M,N,P modules, p,q integers
229COMPUTE: cup-product Ext^p(M',N') x Ext^q(N',P') ---> Ext^(p+q)(M',P'),
230         where M':=R^m/M, if M in R^m, R basering (i.e. M':=coker(matrix(M)))
231ASSUME:  all Ext's  are of finite dimension
232RETURN:  - if called with 5 arguments: matrix of the associated linear map
233         Ext^p (tensor) Ext^q --> Ext^(p+q),  i.e. the columns of <matrix>
234         present the coordinates of the cup products (b_i & c_j) with respect
235         to a kbase of Ext^p+q  (b_i resp. c_j are the choosen bases of Ext^p,
236         resp. Ext^q).@*
237         - if called with 6 arguments: list L,
238@format
239      L[1] = matrix (see above)
240      L[2] = matrix of kbase of Ext^p(M',N')
241      L[3] = matrix of kbase of Ext^q(N',P')
242      L[4] = matrix of kbase of Ext^p+q(N',P')
243@end format
244NOTE:    printlevel >=1;  shows what is going on.
245         printlevel >=2;  shows the result in another representation.@*
246         For computing the cupproduct of M,N itself, apply proc to syz(M),
247         syz(N)!
248EXAMPLE: example cupproduct; shows examples
249"
250{
251//---------- initialization ---------------------------------------------------
252   int    e1,e2,e3,i,j,k,f0,f1,f2;
253   module M1,M2,N1,N2,P1,P2,A,B,C,ker,ima,extMN,extMN0,extMP,
254          extMP0,extNP,extNP0;
255   matrix cup[1][0];
256   matrix kbMN,kbMP,kbNP,lift1,mA,mB,mC;
257   ideal  test1,test2,null;
258   int pp = printlevel-voice+3;  // pp=printlevel+1 (default: p=1)
259//-----------------------------------------------------------------------------
260//compute resolutions of M and N
261//                     M<--F(0)<--- ...  <---F(p+q+1)
262//                     N<--G(0)<--- ...  <---G(q+1)
263//-----------------------------------------------------------------------------
264   list resM = nres(M,p+q+1);
265   M1 = resM[p];
266   M2 = resM[p+1];
267   list resN = nres(N,q+1);
268   N1 = resN[q];
269   N2 = resN[q+1];
270   P1 = resM[p+q];
271   P2 = resM[p+q+1];
272//-------test: Ext==0?---------------------------------------------------------
273   test1 = simplify(ideal(M1),10);
274   test2 = simplify(ideal(N),10);
275   if (test1[1]==0) { dbprint(pp,"//Ext(M,N)=0");return(cup); }
276   test1 = simplify(ideal(N1),10);
277   test2 = simplify(ideal(P),10);
278   if (test1[1]==0) { dbprint(pp,"//Ext(N,P)=0");return(cup); }
279   test1 = simplify(ideal(P1),10);
280   if (test1[1]==0) { dbprint(pp,"//Ext(M,P)=0");return(cup); }
281 //------ compute kbases of Ext's ---------------------------------------------
282 //------ Ext(M,N)
283   test1 = simplify(ideal(M2),10);
284   if (test1[1]==0) { ker = freemodule(ncols(M1)*nrows(N));}
285   else
286   {
287     A   = kontrahom(M2,nrows(N));
288     B   = kohom(N,ncols(M2));
289     C   = intersect(A,B);
290     C   = reduce(C,std(ideal(0)));C=simplify(C,10);
291     ker = lift(A,C)+syz(A);
292   }
293   ima   = kohom(N,ncols(M1));
294   A     = kontrahom(M1,nrows(N));
295   ima   = ima+A;
296   extMN = modulo(ker,ima);
297   extMN0= std(extMN);
298   e1    = vdim(extMN0);
299   dbprint(pp-1,"// vdim Ext(M,N) = "+string(e1));
300   if (e1 < 0)
301   {
302     "// Ext(M,N) not of finite dimension";
303     return(cup);
304   }
305   kbMN  = kbase(extMN0);
306   kbMN = matrix(ker)*kbMN;
307   dbprint(pp-1,"// kbase of Ext^p(M,N)",
308          "//  - the columns present the kbase elements in Hom(F(p),G(0))",
309          "//  - F(*),G(*) are free resolutions of M and N",kbMN);
310//------- Ext(N,P)
311   test1 = simplify(ideal(N2),10);
312   if (test1[1]==0) {  ker = freemodule(ncols(N1)*nrows(P)); }
313   else
314   {
315     A = kontrahom(N2,nrows(P));
316     B = kohom(P,ncols(N2));
317     C = intersect(A,B);
318     C = reduce(C,std(ideal(0)));C=simplify(C,10);
319     ker = lift(A,C)+syz(A);
320   }
321   ima   = kohom(P,ncols(N1));
322   A     = kontrahom(N1,nrows(P));
323   ima   = ima+A;
324   extNP = modulo(ker,ima);
325   extNP0= std(extNP);
326   e2    = vdim(extNP0);
327   dbprint(pp-1,"// vdim Ext(N,P) = "+string(e2));
328   if (e2 < 0)
329   {
330     "// Ext(N,P) not of finite dimension";
331     return(cup);
332   }
333   kbNP  = kbase(extNP0);
334   kbNP  = matrix(ker)*kbNP;
335   dbprint(pp-1,"// kbase of Ext(N,P):",kbNP,
336          "// kbase of Ext^q(N,P)",
337          "//  - the columns present the kbase elements in Hom(G(q),H(0))",
338          "//  - G(*),H(*) are free resolutions of N and P",kbNP);
339
340//------ Ext(M,P)
341   test1 = simplify(ideal(P2),10);
342   if (test1[1]==0) { ker = freemodule(ncols(P1)*nrows(P)); }
343   else
344   {
345     A = kontrahom(P2,nrows(P));
346     B = kohom(P,ncols(P2));
347     C = intersect(A,B);
348     C = reduce(C,std(ideal(0)));C=simplify(C,10);
349     ker = lift(A,C)+syz(A);
350   }
351   ima   = kohom(P,ncols(P1));
352   A     = kontrahom(P1,nrows(P));
353   ima   = ima+A;
354   extMP = modulo(ker,ima);
355   extMP0= std(extMP);
356   e3    = vdim(extMP0);
357   dbprint(pp-1,"// vdim Ext(M,P) = "+string(e3));
358   if (e3 < 0)
359   {
360     "// Ext(M,P) not of finite dimension";
361     return(cup);
362   }
363   kbMP  = kbase(extMP0);
364   kbMP  = matrix(ker)*kbMP;
365   dbprint(pp-1,"// kbase of Ext^p+q(M,P)",
366          "//  - the columns present the kbase elements in Hom(F(p+q),H(0))",
367          "//  - F(*),H(*) are free resolutions of M and P",kbMP);
368//----- lift kbase of Ext(M,N) ------------------------------------------------
369   lift1 = kbMN;
370   for (i=1;i<=q;i=i+1)
371   {
372     mA = kontrahom(resM[p+i],nrows(resN[i]));
373     mB = kohom(resN[i],ncols(resM[p+i]));
374     lift1 = lift(mB,mA*lift1);
375   }
376   dbprint(pp-1,"// lifting of kbase of Ext^p(M,N)",
377   "//  - the columns present liftings of kbase elements"+
378   " in Hom(F(p+q),G(q))",lift1);
379//-------  compute: cup-products of base-elements -----------------------------
380   f0 = nrows(P);
381   f1 = ncols(N1);
382   f2 = ncols(resM[p+q]);
383   for (i=1;i<=e1;i=i+1)
384   {
385     for (j=1;j<=e2;j=j+1)
386     {
387       mA = matrix(ideal(lift1[j]),f1,f2);
388       mB = matrix(ideal(kbMP[i]),f0,f1);
389       mC = mB*mA;
390       mC = matrix(ideal(mC),f0*f2,1);
391       cup= concat(cup,mC);
392     }
393   }
394   dbprint(pp-1,"// matrix of cup-products (in Ext^p+q)",cup,
395                "////// end level 2 //////");
396//------- compute: presentation of base-elements -----------------------------
397   cup = lift(ker,cup);
398   cup = lift_kbase(cup,extMP0);
399//------- special output ------------------------------------------------------
400   if (pp>2)
401   {
402     "// the associated matrices of the bilinear mapping 'cup' ";
403     "// corresponding to the kbase elements of Ext^p+q(M,P) are shown,";
404     "//  i.e. the rows of the final matrix are written as matrix of";
405     "//  a bilinear form on Ext^p x Ext^q";
406     matrix BL[e1][e2];
407     for (k=1;k<=e3;k=k+1)
408     {
409       "//----component "+string(k)+":";
410       for (i=1;i<=e1;i=i+1)
411       {
412         for (j=1;j<=e2;j=j+1)
413         {
414           BL[i,j]=cup[k,j+e1*(i-1)];
415         }
416        }
417        print(BL);
418      }
419     "////// end level 3 //////";
420   }
421   if (size(#)) { return(cup,kbMN,kbNP,kbMP);}
422   else         { return(cup); }
423}
424example
425{"EXAMPLE";  echo=2;
426  int p      = printlevel;
427  ring  rr   = 32003,(x,y,z),(dp,C);
428  ideal  I   = x4+y3+z2;
429  qring  o   = std(I);
430  module M   = [x,y,0,z],[y2,-x3,z,0],[z,0,-y,-x3],[0,z,x,-y2];
431  print(cupproduct(M,M,M,1,3));
432  printlevel = 3;
433  list l     = (cupproduct(M,M,M,1,3,"any"));
434  show(l[1]);show(l[2]);
435  printlevel = p;
436}
437///////////////////////////////////////////////////////////////////////////////
438
439proc Ext_R (intvec v, module M, list #)
440"USAGE:   Ext_R(v,M[,p]);  v int resp. intvec , M module, p int
441COMPUTE: A presentation of Ext^k(M',R); for k=v[1],v[2],..., M'=coker(M).
442         Let
443@example
444  0 <-- M' <-- F0 <-M-- F1 <-- F2 <-- ...
445@end example
446         be a free resolution of M'. If
447@example
448        0 --> F0* -A1-> F1* -A2-> F2* -A3-> ...
449@end example
450         is the dual sequence, Fi*=Hom(Fi,R), then Ext^k = ker(Ak+1)/im(Ak)
451         is presented as in the following exact sequences:
452@example
453    R^p --syz(Ak+1)-> Fk* ---Ak+1---->  Fk+1* ,
454    R^q ----Ext^k---> R^p --syz(Ak+1)-> Fk*/im(Ak).
455@end example
456         Hence, Ext^k=modulo(syz(Ak+1),Ak) presents Ext^k(M',R).
457RETURN:  - module Ext, a presentation of Ext^k(M',R) if v is of type int@*
458         - a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.@*
459         - In case of a third argument of type int return a list l:
460@format
461     l[1] = module Ext^k resp. list of Ext^k
462     l[2] = SB of Ext^k resp. list of SB of Ext^k
463     l[3] = matrix resp. list of matrices, each representing a kbase of Ext^k
464              (if finite dimensional)
465@end format
466DISPLAY: printlevel >=0: (affine) dimension of Ext^k for each k  (default)
467         printlevel >=1: Ak, Ak+1 and kbase of Ext^k in Fk*
468NOTE:    In order to compute Ext^k(M,R) use the command Ext_R(k,syz(M));@*
469         By default, the procedure uses the @code{mres} command. If called
470         with the additional parameter @code{\"sres\"}, the @code{sres} command
471         is used instead.@*
472         If the attribute @code{\"isHomog\"} has been set for the input module, it
473         is also set for the returned module (accordingly).
474EXAMPLE: example Ext_R; shows an example
475"
476{
477  // In case M is known to be a SB, set attrib(M,"isSB",1); in order to
478  // avoid  unnecessary SB computations
479
480  //------------ check for weight vector (graded case) -----------------------
481  int withWeight;
482  intvec weightM,weightR,ww;
483  if ( typeof(attrib(M,"isHomog"))!="string" )
484  {
485    weightM=attrib(M,"isHomog");
486    withWeight=1;
487  }
488
489  //------------ initialisation ----------------------------------------------
490  module m1t,m1,m2,m2t,ret,ret0,ker;
491  vector leadCol;
492  matrix kb;
493  module G;
494  list L1,L2,L3,L,K;
495  resolution resm2;
496  int j,k,max,ii,t1,t2,di,leadComp,shift;
497  intvec A1,A2,A3;
498  int s = size(v);
499  intvec v1 = sort(v)[1];
500  max = v1[s];                 // the maximum integer occurring in intvec v
501  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
502  // --------------- Variante mit sres
503  for( ii=1; ii<=size(#); ii++ )
504  {
505    if (typeof(#[ii])=="int") { // return a list if t2=1
506      t2=1;
507    }
508    else {
509      if (typeof(#[ii])=="string") {
510        // NOTE: at this writing, sres does not return weights
511        if ( #[ii]=="sres" ) { t1=1; } // use sres instead of mres if t1=1
512      }
513    }
514  }
515
516//----------------- compute resolution of coker(M) ----------------------------
517  if( max<0 )
518  {
519    dbprint(p,"// Ext^i=0 for i<0!");
520    module Result=[1];
521    if (withWeight==1) { attrib(Result,"isHomog",intvec(0)); }
522    if (s==1) {
523      if (t2==0) { return(Result); }
524      else       { return( list(Result,Result,matrix(0)) ); }
525    }
526    list Out, KBOut;
527    for (j=1;j<=s;j++) {
528      Out[j] = Result;
529      KBOut[j] = matrix(0);
530    }
531    if (t2==0) { return(Out); }
532    else       { return( list(Out,Out,KBOut) ); }
533  }
534  if( t1==1 )
535  { // compute resolution via sres command
536    if( attrib(M,"isSB")==0 ) {
537      if (size(M)==0) { attrib(M,"isSB")=1; }
538      else { M=std(M); }
539    }
540    list resl = sres(M,max+1);
541    if (withWeight) {
542      // ****
543      // **** at this writing, sres does not return weights, we have to
544      // **** go through the resolution to compute them
545      // ****
546      attrib(resl,"isHomog",weightM);  // weightM = weights of M
547      G=resl[1];
548      attrib(G,"isHomog",weightM);
549      resl[1]=G;
550      weightR=weightM;
551
552      for (j=2; j<=size(resl); j++) {
553        if (size(G)!=0) {
554          ww=0;
555          for (k=1; k<=ncols(G); k++) {
556            if (size(G[k])==0) { ww[k]=0; }
557            else {
558              leadCol = leadmonom(G[k]);
559              leadComp = nrows(leadCol);
560              ww[k] = deg(leadCol)+weightR[leadComp];
561            }
562          }
563          G=resl[j];
564          attrib(G,"isHomog",ww);
565          resl[j]=G;
566          weightR=ww;
567        }
568      }
569    }
570  }
571  else { list resl = mres(M,max+1);
572    if ((withWeight) and (size(M)==0)) {
573      // ***** At this writing: special treatment for zero module needed
574      G=resl[1];
575      attrib(G,"isHomog",weightM);
576      resl[1]=G;
577    }
578  }
579  for( ii=1; ii<=s; ii++ )
580  {
581//-----------------  apply Hom(_,R) at k-th place -----------------------------
582    k=v[ii];
583    dbprint(p-1,"// Computing Ext^"+string(k)+":");
584    if( k<0 )                   // Ext^k=0 for negative k
585    {
586      dbprint(p-1,"// Ext^i=0 for i<0!");
587      ret    = gen(1);
588      ret0   = std(ret);
589      if (withWeight==1) {
590        attrib(ret,"isHomog",intvec(0));
591        attrib(ret0,"isHomog",intvec(0));
592      }
593      L1[ii] = ret;
594      L2[ii] = ret0;
595      L3[ii] = matrix(kbase(ret0));
596      di=dim(ret0);
597      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
598      if (di==0)
599      {
600        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(ret0)));
601      }
602      dbprint(p,"");
603    }
604    else
605    {
606      m2t = resl[k+1];
607      m2 = transpose(m2t);
608      if ((typeof(attrib(m2t,"isHomog"))!="string" ) && (withWeight))
609      {
610        // ------------- compute weights for dual -----------------------------
611        weightR=attrib(m2t,"isHomog");
612        // --------------------------------------------------------------------
613        // *** set correct weights (at this writing, shift in resolution
614        // *** is not considered when defining the weights for the
615        // *** modules in the resolution):
616        A1=attrib(M,"isHomog");
617        A2=attrib(resl[1],"isHomog");
618        shift=A1[1]-A2[1];
619        for (j=1; j<=size(weightR); j++) { weightR[j]=weightR[j]+shift; }
620        attrib(m2t,"isHomog",weightR);
621        // --------------------------------------------------------------------
622        ww=0;
623        for (j=1; j<=nrows(m2); j++) {
624          if (size(m2t[j])==0) { ww[j]=0; }
625          else {
626            leadCol = leadmonom(m2t[j]);
627            leadComp = nrows(leadCol);
628            ww[j] = deg(leadCol)+weightR[leadComp];
629          }
630        }
631        attrib(m2,"isHomog",-ww);  // dualize --> negative weights
632        // --------------------------------------------------------------------
633        // *** the following should be replaced by the syz command,
634        // *** but syz forgets weights.....
635        resm2 = nres(m2,2);
636        ker = resm2[2];
637        if ((size(ker)>0) and (size(m2)>0)) {
638          // ------------------------------------------------------------------
639          // *** set correct weights (at this writing, shift in resolution
640          // *** is not considered when defining the weights for the
641          // *** modules in the resolution):
642          A1=attrib(resm2,"isHomog");
643          A2=attrib(resm2[1],"isHomog");
644          A3=attrib(ker,"isHomog");
645          shift=A1[1]-A2[1];
646          for (j=1; j<=size(A3); j++) { A3[j]=A3[j]+shift; }
647          // *** set correct weights where weights are undetermined due to
648          // *** zero columns in m2  (read weights from m2t)
649          for (j=1; j<=ncols(m2); j++) {
650            if (size(m2[j])==0) { A3[j]=-weightR[j]; }
651          }
652          attrib(ker,"isHomog",A3);
653          // ------------------------------------------------------------------
654        }
655      }
656      else {
657        ker = syz(m2);
658      }
659
660      if( k==0 ) { matrix MMM1[ncols(m2)][1];
661        m1=MMM1;
662      }
663      else { // k>0
664        m1t = resl[k];
665        m1 = transpose(resl[k]);
666        if ((typeof(attrib(m1t,"isHomog"))!="string" ) && (withWeight)) {
667        // ------------- compute weights for dual -----------------------------
668          weightR=attrib(resl[k],"isHomog");
669          // ------------------------------------------------------------------
670          // *** set correct weights (at this writing, shift in resolution
671          // *** is not considered when defining the weights for the
672          // *** modules in the resolution):
673          A1=attrib(M,"isHomog");
674          A2=attrib(resl[1],"isHomog");
675          shift=A1[1]-A2[1];
676          for (j=1; j<=size(weightR); j++) { weightR[j]=weightR[j]+shift; }
677          attrib(m1t,"isHomog",weightR);
678          // ------------------------------------------------------------------
679          ww=0;
680          for (j=1; j<=nrows(m1); j++) {
681            if (size(m1t[j])==0) { ww[j]=0; }
682            else {
683              leadCol = leadmonom(m1t[j]);
684              leadComp = nrows(leadCol);
685              ww[j] = deg(leadCol)+weightR[leadComp];
686            }
687          }
688          attrib(m1,"isHomog",-ww);  // dualize --> negative weights
689        }
690      }
691      //----------------- presentation of ker(m2)/im(m1) ---------------------
692      if ((k==0) and (size(M)==0)) {
693        ret = M;
694        if (withWeight) { attrib(ret,"isHomog",-weightM); }
695      }
696      else {
697        ret = modulo(ker,m1);
698      }
699      dbprint(p-1,
700         "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of M,",
701         "// then F"+string(k)+"*-->F"+string(k+1)+"* is given by:",m2,
702         "// and F"+string(k-1)+"*-->F"+string(k)+"* is given by:",m1,"");
703      ret0 = std(ret);
704
705      di=dim(ret0);
706      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
707      if (di==0)
708      {
709        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(ret0)));
710      }
711      dbprint(p,"");
712      if( t2 )
713      {
714         if( vdim(ret0)>=0 )
715         {
716            kb = kbase(ret0);
717            if ( size(ker)!=0 ) { kb = matrix(ker)*kb; }
718            dbprint(p-1,
719              "// columns of matrix are kbase of Ext^"+string(k)+" in F"
720              +string(k)+"*:",kb,"");
721            L3[ii] = kb;
722         }
723         L2[ii] = ret0;
724      }
725      L1[ii] = ret;
726    }
727  }
728  if( t2 )
729  {
730     if( s>1 ) { L = L1,L2,L3; return(L); }
731     else { L = ret,ret0,kb; return(L); }
732  }
733  else
734  {
735     if( s>1 ) { return(L1); }
736     else { return(ret); }
737  }
738}
739example
740{"EXAMPLE:";     echo=2;
741  int p      = printlevel;
742  printlevel = 1;
743  ring r     = 0,(x,y,z),dp;
744  ideal i    = x2y,y2z,z3x;
745  module E   = Ext_R(1,i);    //computes Ext^1(r/i,r)
746  is_zero(E);
747
748  qring R    = std(x2+yz);
749  intvec v   = 0,2;
750  printlevel = 2;             //shows what is going on
751  ideal i    = x,y,z;         //computes Ext^i(r/(x,y,z),r/(x2+yz)), i=0,2
752  list L     = Ext_R(v,i,1);  //over the qring R=r/(x2+yz), std and kbase
753  printlevel = p;
754}
755///////////////////////////////////////////////////////////////////////////////
756
757proc Ext (intvec v, module M, module N, list #)
758"USAGE:   Ext(v,M,N[,any]);  v int resp. intvec, M,N modules
759COMPUTE: A presentation of Ext^k(M',N'); for k=v[1],v[2],... where
760         M'=coker(M) and N'=coker(N). Let
761@example
762       0 <-- M' <-- F0 <-M-- F1 <-- F2 <--... ,
763       0 <-- N' <-- G0 <--N- G1
764@end example
765         be a free resolution of M', resp. a presentation of N'. Consider
766         the commutative diagram
767@example
768           0                  0                  0
769           |^                 |^                 |^
770   --> Hom(Fk-1,N') -Ak-> Hom(Fk,N') -Ak+1-> Hom(Fk+1,N')
771           |^                 |^                 |^
772   --> Hom(Fk-1,G0) -Ak-> Hom(Fk,G0) -Ak+1-> Hom(Fk+1,G0)
773                              |^                 |^
774                              |C                 |B
775                          Hom(Fk,G1) ------> Hom(Fk+1,G1)
776
777      (Ak,Ak+1 induced by M and B,C induced by N).
778@end example
779         Let K=modulo(Ak+1,B), J=module(Ak)+module(C) and Ext=modulo(K,J),
780         then we have exact sequences
781@example
782    R^p --K-> Hom(Fk,G0) --Ak+1-> Hom(Fk+1,G0)/im(B),
783
784    R^q -Ext-> R^p --K-> Hom(Fk,G0)/(im(Ak)+im(C)).
785@end example
786         Hence, Ext presents Ext^k(M',N').
787RETURN:  - module Ext, a presentation of Ext^k(M',N') if v is of type int@*
788         - a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.@*
789         - In case of a third argument of any type return a list l:
790@format
791             l[1] = module Ext/list of Ext^k
792             l[2] = SB of Ext/list of SB of Ext^k
793             l[3] = matrix/list of matrices, each representing a kbase of Ext^k
794                       (if finite dimensional)
795@end format
796DISPLAY: printlevel >=0: dimension, vdim of Ext^k for each k  (default).
797@*       printlevel >=1: matrices Ak, Ak+1 and kbase of Ext^k in Hom(Fk,G0)
798         (if finite dimensional)
799NOTE:    In order to compute Ext^k(M,N) use the command Ext(k,syz(M),syz(N));
800         or: list P=mres(M,2); list Q=mres(N,2); Ext(k,P[2],Q[2]);
801EXAMPLE: example Ext; shows an example
802"
803{
804//---------- initialisation ---------------------------------------------------
805  int k,max,ii,l,row,col,di;
806  module A,B,C,D,M1,M2,N1,ker,imag,extMN,extMN0;
807  matrix kb;
808  list L1,L2,L3,L,resM,K;
809  ideal  test1;
810  intmat Be;
811  int s = size(v);
812  intvec v1 = sort(v)[1];
813  max = v1[s];                    // the maximum integer occurring in intvec v
814  int p = printlevel-voice+3;     // p=printlevel+1 (default: p=1)
815//---------- test: coker(N)=basering, coker(N)=0 ? ----------------------------
816  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
817  N1 = std(N);
818  if( size(N1)==0 )      //coker(N)=basering, in this case proc Ext_R is faster
819  { printlevel=printlevel+1;
820    if( size(#)==0 )
821    { def E = Ext_R(v,M);
822      printlevel=printlevel-1;
823      return(E);
824    }
825    else
826    { def E = Ext_R(v,M,#[1]);
827       printlevel=printlevel-1;
828       return(E);
829     }
830  }
831  if( dim(N1)==-1 )                          //coker(N)=0, all Ext-groups are 0
832  { dbprint(p-1,"2nd module presents 0, hence Ext^k=0, for all k");
833    for( ii=1; ii<=s; ii++ )
834    { k=v[ii];
835      extMN    = gen(1);
836      extMN0   = std(extMN);
837      L1[ii] = extMN;
838      L2[ii] = extMN0;
839      L3[ii] = matrix(kbase(extMN0));
840      di=dim(extMN0);
841      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
842      if (di==0)
843      {
844        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(extMN0)));
845      }
846      dbprint(p,"");
847    }
848  }
849  else
850  {
851    if( size(N1) < size(N) ) { N=N1;}
852    row = nrows(N);
853//---------- resolution of M -------------------------------------------------
854    resM = mres(M,max+1);
855    for( ii=1; ii<=s; ii++ )
856    { k=v[ii];
857      if( k<0  )                                   // Ext^k is 0 for negative k
858      { dbprint(p-1,"// Ext^k=0 for k<0!");
859        extMN    = gen(1);
860        extMN0   = std(extMN);
861        L1[ii] = extMN;
862        L2[ii] = extMN0;
863        L3[ii] = matrix(kbase(extMN0));
864        di=dim(extMN0);
865        dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
866        if (di==0)
867        {
868          dbprint(p,"// vdim of Ext^"+string(k)+":       "
869                  +string(vdim(extMN0)));
870        }
871        dbprint(p,"");
872      }
873      else
874      { M2 = resM[k+1];
875        if( k==0 ) { M1=0; }
876        else { M1 = resM[k]; }
877        col = nrows(M2);
878        D = kohom(N,col);
879//---------- computing homology ----------------------------------------------
880        imag  = kontrahom(M1,row);
881        A     = kontrahom(M2,row);
882        B     = kohom(N,ncols(M2));
883        ker   = modulo(A,B);
884        imag  = imag,D;
885        extMN = modulo(ker,imag);
886        dbprint(p-1,"// Computing Ext^"+string(k)+
887         " (help Ext; gives an explanation):",
888      "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of coker(M),",
889      "// and 0<--coker(N)<--G0<--G1 a presentation of coker(N),",
890      "// then Hom(F"+string(k)+",G0)-->Hom(F"+string(k+1)+
891      ",G0) is given by:",A,
892      "// and Hom(F"+string(k-1)+",G0) + Hom(F"+string(k)+",G1)-->Hom(F"
893                                      +string(k)+",G0) is given by:",imag,"");
894        extMN0 = std(extMN);
895        di=dim(extMN0);
896        dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
897        if (di==0)
898        {
899          dbprint(p,"// vdim of Ext^"+string(k)+":       "
900                  +string(vdim(extMN0)));
901        }
902        dbprint(p,"");
903
904//---------- more information -------------------------------------------------
905        if( size(#)>0 )
906        { if( vdim(extMN0) >= 0 )
907          { kb = kbase(extMN0);
908            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
909            dbprint(p-1,"// columns of matrix are kbase of Ext^"+
910                     string(k)+" in Hom(F"+string(k)+",G0)",kb,"");
911            if( p>0 )
912            { for (l=1;l<=ncols(kb);l=l+1)
913              {
914                "// element",l,"of kbase of Ext^"+string(k)+" in Hom(F"+string(k)+",G0)";
915                "// as matrix: F"+string(k)+"-->G0";
916                print(matrix(ideal(kb[l]),row,col));
917              }
918              "";
919            }
920            L3[ii] = matrix(kb);
921          }
922          L2[ii] = extMN0;
923        }
924        L1[ii] = extMN;
925      }
926    }
927  }
928  if( size(#) )
929  {  if( s>1 ) { L = L1,L2,L3; return(L); }
930     else { L = extMN,extMN0,matrix(kb); return(L); }
931  }
932  else
933  {  if( s>1 ) { return(L1); }
934     else { return(extMN); }
935  }
936}
937example
938{"EXAMPLE:";   echo=2;
939  int p      = printlevel;
940  printlevel = 1;
941  ring r     = 0,(x,y),dp;
942  ideal i    = x2-y3;
943  ideal j    = x2-y5;
944  list E     = Ext(0..2,i,j);    // Ext^k(r/i,r/j) for k=0,1,2 over r
945  qring R    = std(i);
946  ideal j    = fetch(r,j);
947  module M   = [-x,y],[-y2,x];
948  printlevel = 2;
949  module E1  = Ext(1,M,j);       // Ext^1(R^2/M,R/j) over R=r/i
950  list l     = Ext(4,M,M,1);     // Ext^4(R^2/M,R^2/M) over R=r/i
951  printlevel = p;
952}
953///////////////////////////////////////////////////////////////////////////////
954
955proc Hom (module M, module N, list #)
956"USAGE:   Hom(M,N,[any]);  M,N=modules
957COMPUTE: A presentation of Hom(M',N'), M'=coker(M), N'=coker(N) as follows:
958         let
959@example
960   F1 --M-> F0 -->M' --> 0,    G1 --N-> G0 --> N' --> 0
961@end example
962         be presentations of M' and N'. Consider
963@example
964                                  0               0
965                                  |^              |^
966       0 --> Hom(M',N') ----> Hom(F0,N') ----> Hom(F1,N')
967                                  |^              |^
968  (A:  induced by M)          Hom(F0,G0) --A-> Hom(F1,G0)
969                                  |^              |^
970  (B,C:induced by N)              |C              |B
971                              Hom(F0,G1) ----> Hom(F1,G1)
972
973@end example
974         Let D=modulo(A,B) and Hom=modulo(D,C), then we have exact sequences
975@example
976   R^p  --D-> Hom(F0,G0) --A-> Hom(F1,G0)/im(B),
977
978 R^q -Hom-> R^p --D-> Hom(F0,G0)/im(C) --A-> Hom(F1,G0)/im(B).
979@end example
980         Hence Hom presents Hom(M',N')
981RETURN:  module Hom, a presentation of Hom(M',N'), resp., in case of
982         3 arguments, a list l (of size <=3):
983@format
984           - l[1] = Hom
985           - l[2] = SB of Hom
986           - l[3] = kbase of coker(Hom) (if finite dimensional, not 0),
987                    represented by elements in Hom(F0,G0) via mapping D
988@end format
989DISPLAY: printlevel >=0: (affine) dimension of Hom  (default)
990@*       printlevel >=1: D and C and kbase of coker(Hom) in Hom(F0,G0)
991@*       printlevel >=2: elements of kbase of coker(Hom) as matrix :F0-->G0
992NOTE:    DISPLAY is as described only for a direct call of 'Hom'. Calling 'Hom'
993         from another proc has the same effect as decreasing printlevel by 1.
994EXAMPLE: example Hom;  shows examples
995"
996{
997//---------- initialisation ---------------------------------------------------
998  int l,p,di;
999  matrix kb;
1000  module A,B,C,D,homMN,homMN0;
1001  list L;
1002//---------- computation of Hom -----------------------------------------------
1003  B = kohom(N,ncols(M));
1004  A = kontrahom(M,nrows(N));
1005  C = kohom(N,nrows(M));
1006  D = modulo(A,B);
1007  homMN = modulo(D,C);
1008  homMN0= std(homMN);
1009  p = printlevel-voice+3;       // p=printlevel+1 (default: p=1)
1010  di= dim(homMN0);
1011  dbprint(p,"// dimension of Hom:  "+string(di));
1012  if (di==0)
1013  {
1014    dbprint(p,"// vdim of Hom:       "+string(vdim(homMN0)));
1015  }
1016  dbprint(p,"");
1017  dbprint(p-1,
1018   "// given  F1 --M-> F0 -->M'--> 0 and  G1 --N-> G0 -->N'--> 0,",
1019   "// show D = ker( Hom(F0,G0) --> Hom(F1,G0)/im(Hom(F1,G1)->Hom(F1,G0)) )",D,
1020   "// show C = im ( Hom(F0,G1) --> Hom(F0,G0) )",C,"");
1021//---------- extra output if size(#)>0 ----------------------------------------
1022  if( size(#)>0 )
1023  {
1024     if( vdim(homMN0)>0 )
1025     {
1026        kb = kbase(homMN0);
1027        kb = matrix(D)*kb;
1028        if( p>2 )
1029        {
1030          for (l=1;l<=ncols(kb);l=l+1)
1031          {
1032            "// element",l,"of kbase of Hom in Hom(F0,G0) as matrix: F0-->G0:";
1033             print(matrix(ideal(kb[l]),nrows(N),nrows(M)));
1034          }
1035        }
1036        else
1037        {
1038          dbprint(p-1,"// columns of matrix are kbase of Hom in Hom(F0,G0)",
1039                      kb);
1040        }
1041        L=homMN,homMN0,kb;
1042        return(L);
1043     }
1044     L=homMN,homMN0;
1045     return(L);
1046  }
1047  return(homMN);
1048}
1049example
1050{"EXAMPLE:";  echo = 2;
1051  int p     = printlevel;
1052  printlevel= 1;   //in 'example proc' printlevel has to be increased by 1
1053  ring r    = 0,(x,y),dp;
1054  ideal i   = x2-y3,xy;
1055  qring q   = std(i);
1056  ideal i   = fetch(r,i);
1057  module M  = [-x,y],[-y2,x],[x3];
1058  module H  = Hom(M,i);
1059  print(H);
1060
1061  printlevel= 2;
1062  list L    = Hom(M,i,1);"";
1063
1064  printlevel=1;
1065  ring s    = 3,(x,y,z),(c,dp);
1066  ideal i   = jacob(ideal(x2+y5+z4));
1067  qring rq=std(i);
1068  matrix M[2][2]=xy,x3,5y,4z,x2;
1069  matrix N[3][2]=x2,x,y3,3xz,x2z,z;
1070  print(M);
1071  print(N);
1072  list l=Hom(M,N,1);
1073  printlevel = p;
1074}
1075///////////////////////////////////////////////////////////////////////////////
1076
1077proc homology (matrix A,matrix B,module M,module N,list #)
1078"USAGE:   homology(A,B,M,N);
1079COMPUTE: Let M and N be submodules of R^m and R^n presenting M'=R^m/M, N'=R^n/N
1080         (R=basering) and let A,B matrices inducing maps
1081@example
1082    R^k --A--> R^m --B--> R^n.
1083@end example
1084         Compute a presentation of the module
1085@example
1086    ker(B)/im(A) := ker(M'/im(A) --B--> N'/im(BM)+im(BA)).
1087@end example
1088         If B induces a map M'-->N' (i.e BM=0) and if im(A) is contained in
1089         ker(B) (that is, BA=0) then ker(B)/im(A) is the homology of the
1090         complex
1091@example
1092    R^k--A-->M'--B-->N'.
1093@end example
1094RETURN:  module H, a presentation of ker(B)/im(A).
1095NOTE:    homology returns a free module of rank m if ker(B)=im(A).
1096EXAMPLE: example homology; shows examples
1097"
1098{
1099  module ker,ima;
1100  ker = modulo(B,N);
1101  ima = A,M;
1102  return(modulo(ker,ima));
1103}
1104example
1105{"EXAMPLE";    echo=2;
1106  ring r;
1107  ideal id=maxideal(4);
1108  qring qr=std(id);
1109  module N=maxideal(3)*freemodule(2);
1110  module M=maxideal(2)*freemodule(2);
1111  module B=[2x,0],[x,y],[z2,y];
1112  module A=M;
1113  module H=homology(A,B,M,N);
1114  H=std(H);
1115  // dimension of homology:
1116  dim(H);
1117  // vector space dimension:
1118  vdim(H);
1119
1120  ring s=0,x,ds;
1121  qring qs=std(x4);
1122  module A=[x];
1123  module B=A;
1124  module M=[x3];
1125  module N=M;
1126  homology(A,B,M,N);
1127}
1128//////////////////////////////////////////////////////////////////////////////
1129
1130proc hom_kernel (matrix A,module M,module N)
1131"USAGE:   hom_kernel(A,M,N);
1132COMPUTE: Let M and N be submodules of R^m and R^n, presenting M'=R^m/M,
1133         N'=R^n/N (R=basering), and let A:R^m-->R^n be a matrix inducing a
1134         map A':M'-->N'. Then ker(A,M,N); computes a presentation K of
1135         ker(A') as in the commutative diagram:
1136@example
1137          ker(A') --->  M' --A'--> N'
1138             |^         |^         |^
1139             |          |          |
1140             R^r  ---> R^m --A--> R^n
1141             |^         |^         |^
1142             |K         |M         |N
1143             |          |          |
1144             R^s  ---> R^p -----> R^q
1145@end example
1146RETURN:  module K, a presentation of ker(A':coker(M)->coker(N)).
1147EXAMPLE: example hom_kernel; shows examples.
1148"
1149{
1150  module M1 = modulo(A,N);
1151  return(modulo(M1,M));
1152}
1153example
1154{"EXAMPLE";    echo=2;
1155  ring r;
1156  module N=[2x,x],[0,y];
1157  module M=maxideal(1)*freemodule(2);
1158  matrix A[2][3]=2x,0,x,y,z2,y;
1159  module K=hom_kernel(A,M,N);
1160  // dimension of kernel:
1161  dim(std(K));
1162  // vector space dimension of kernel:
1163  vdim(std(K));
1164  print(K);
1165}
1166///////////////////////////////////////////////////////////////////////////////
1167
1168proc kohom (matrix M, int j)
1169"USAGE:   kohom(A,k); A=matrix, k=integer
1170RETURN:  matrix Hom(R^k,A), i.e. let A be a matrix defining a map F1-->F2
1171         of free R-modules, then the matrix of Hom(R^k,F1)-->Hom(R^k,F2)
1172         is computed (R=basering).
1173EXAMPLE: example kohom; shows an example.
1174"
1175{
1176  if (j==1)
1177  { return(M);}
1178  if (j>1)
1179  { return(tensor(M,diag(1,j))); }
1180  else { return(0);}
1181}
1182example
1183{"EXAMPLE:";   echo=2;
1184  ring r;
1185  matrix n[2][3]=x,y,5,z,77,33;
1186  print(kohom(n,3));
1187}
1188///////////////////////////////////////////////////////////////////////////////
1189
1190proc kontrahom (matrix M, int j)
1191"USAGE:   kontrahom(A,k); A=matrix, k=integer
1192RETURN:  matrix Hom(A,R^k), i.e. let A be a matrix defining a map F1-->F2 of
1193         free  R-modules, then the matrix of Hom(F2,R^k)-->Hom(F1,R^k) is
1194         computed (R=basering).
1195EXAMPLE: example kontrahom; shows an example.
1196"
1197{
1198if (j==1)
1199  { return(transpose(M));}
1200  if (j>1)
1201  { return(transpose(tensor(diag(1,j),M)));}
1202  else { return(0);}
1203}
1204example
1205{"EXAMPLE:";  echo=2;
1206  ring r;
1207  matrix n[2][3]=x,y,5,z,77,33;
1208  print(kontrahom(n,3));
1209}
1210///////////////////////////////////////////////////////////////////////////////
1211///////////////////////////////////////////////////////////////////////////////
1212
1213proc tensorMod(module Phi, module Psi)
1214"USAGE:   tensorMod(M,N);  M,N modules
1215COMPUTE: presentation matrix A of the tensor product T of the modules
1216         M'=coker(M), N'=coker(N): if matrix(M) defines a map M: R^r-->R^s and
1217         matrix(N) defines a map N: R^p-->R^q, then A defines a presentation
1218@example
1219         R^(sp+rq) --A-> R^(sq)  --> T --> 0 .
1220@end example
1221RETURN:  matrix A satisfying coker(A) = tensorprod(coker(M),coker(N)) .
1222EXAMPLE: example tensorMod; shows an example.
1223"
1224{
1225   int s=nrows(Phi);
1226   int q=nrows(Psi);
1227   matrix A=tensor(unitmat(s),Psi);
1228   matrix B=tensor(Phi,unitmat(q));
1229   matrix R=concat(A,B);
1230   return(R);
1231}
1232example
1233{"EXAMPLE:";  echo=2;
1234  ring A=0,(x,y,z),dp;
1235  matrix M[3][3]=1,2,3,4,5,6,7,8,9;
1236  matrix N[2][2]=x,y,0,z;
1237  print(M);
1238  print(N);
1239  print(tensorMod(M,N));
1240}
1241///////////////////////////////////////////////////////////////////////////////
1242proc Tor(intvec v, module M, module N, list #)
1243"USAGE:   Tor(v,M,N[,any]);  v int resp. intvec, M,N modules
1244COMPUTE: a presentation of Tor_k(M',N'), for k=v[1],v[2],... , where
1245         M'=coker(M) and N'=coker(N): let
1246@example
1247       0 <-- M' <-- G0 <-M-- G1
1248       0 <-- N' <-- F0 <--N- F1 <-- F2 <--...
1249@end example
1250         be a presentation of M', resp. a free resolution of N', and consider
1251         the commutative diagram
1252@example
1253          0                    0                    0
1254          |^                   |^                   |^
1255  Tensor(M',Fk+1) -Ak+1-> Tensor(M',Fk) -Ak-> Tensor(M',Fk-1)
1256          |^                   |^                   |^
1257  Tensor(G0,Fk+1) -Ak+1-> Tensor(G0,Fk) -Ak-> Tensor(G0,Fk-1)
1258                               |^                   |^
1259                               |C                   |B
1260                          Tensor(G1,Fk) ----> Tensor(G1,Fk-1)
1261
1262       (Ak,Ak+1 induced by N and B,C induced by M).
1263@end example
1264         Let K=modulo(Ak,B), J=module(C)+module(Ak+1) and Tor=modulo(K,J),
1265         then we have exact sequences
1266@example
1267    R^p  --K-> Tensor(G0,Fk) --Ak-> Tensor(G0,Fk-1)/im(B),
1268
1269    R^q -Tor-> R^p --K-> Tensor(G0,Fk)/(im(C)+im(Ak+1)).
1270@end example
1271         Hence, Tor presents Tor_k(M',N').
1272RETURN:  - if v is of type int: module Tor, a presentation of Tor_k(M',N');@*
1273         - if v is of type intvec: a list of Tor_k(M',N') (k=v[1],v[2],...);@*
1274         - in case of a third argument of any type: list l with
1275@format
1276     l[1] = module Tor/list of Tor_k(M',N'),
1277     l[2] = SB of Tor/list of SB of Tor_k(M',N'),
1278     l[3] = matrix/list of matrices, each representing a kbase of Tor_k(M',N')
1279                (if finite dimensional), or 0.
1280@end format
1281DISPLAY: printlevel >=0: (affine) dimension of Tor_k for each k  (default).
1282@*       printlevel >=1: matrices Ak, Ak+1 and kbase of Tor_k in Tensor(G0,Fk)
1283         (if finite dimensional).
1284NOTE:    In order to compute Tor_k(M,N) use the command Tor(k,syz(M),syz(N));
1285         or: list P=mres(M,2); list Q=mres(N,2); Tor(k,P[2],Q[2]);
1286EXAMPLE: example Tor; shows an example
1287{
1288//---------- initialisation ---------------------------------------------------
1289  int k,max,ii,l,row,col,di;
1290  module A,B,C,D,N1,N2,M1,ker,imag,Im,Im1,Im2,f,torMN,torMN0;
1291  matrix kb;
1292  list L1,L2,L3,L,resN,K;
1293  ideal  test1;
1294  intmat Be;
1295  int s = size(v);
1296  intvec v1 = sort(v)[1];
1297  max = v1[s];                   // maximum integer occurring in intvec v
1298  int p = printlevel-voice+3;    // p=printlevel+1 (default: p=1)
1299
1300//---------- test: coker(M)=basering, coker(M)=0 ? ----------------------------
1301  if( max<0 ) { dbprint(p,"// Tor_i=0 for i<0!"); return([1]); }
1302  M1 = std(M);
1303
1304  if( size(M1)==0 or size(N)==0 )  // coker(M)=basering ==> Tor_i=0 for i>0
1305  {
1306    dbprint(p-1,"// one of the modules M',N' is free, hence Tor_i=0 for i<>0");
1307    for( ii=1; ii<=s; ii++ )
1308    {
1309      k=v[ii];
1310      if (k==0) { torMN=module(tensorMod(M1,N)); }
1311      else { torMN  = gen(1); }
1312      torMN0 = std(torMN);
1313      L1[ii] = torMN;
1314      L2[ii] = torMN0;
1315      L3[ii] = matrix(kbase(torMN0));
1316      di=dim(torMN0);
1317      dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1318      if (di==0)
1319      {
1320        dbprint(p,"// vdim of Tor_"+string(k)+":       "
1321                  +string(vdim(torMN0)));
1322      }
1323      dbprint(p,"");
1324    }
1325
1326    if( size(#) )
1327    {  if( s>1 ) { L = L1,L2,L3; return(L); }
1328       else { L = torMN,torMN0,L3[1]; return(L); }
1329    }
1330    else
1331    {  if( s>1 ) { return(L1); }
1332       else { return(torMN); }
1333    }
1334  }
1335
1336  if( dim(M1)==-1 )              // coker(M)=0, all Tor's are 0
1337  { dbprint(p-1,"2nd module presents 0, hence Tor_k=0, for all k");
1338    for( ii=1; ii<=s; ii++ )
1339    { k=v[ii];
1340      torMN    = gen(1);
1341      torMN0   = std(torMN);
1342      L1[ii] = torMN;
1343      L2[ii] = torMN0;
1344      L3[ii] = matrix(kbase(torMN0));
1345      di=dim(torMN0);
1346      dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1347      if (di==0)
1348      {
1349        dbprint(p,"// vdim of Tor_"+string(k)+":       "
1350                  +string(vdim(torMN0)));
1351      }
1352      dbprint(p,"");
1353    }
1354  }
1355  else
1356  {
1357    if( size(M1) < size(M) ) { M=M1;}
1358    row = nrows(M);
1359//---------- resolution of N -------------------------------------------------
1360    resN = mres(N,max+1);
1361    for( ii=1; ii<=s; ii++ )
1362    { k=v[ii];
1363      if( k<0  )                             // Tor_k is 0 for negative k
1364      { dbprint(p-1,"// Tor_k=0 for k<0!");
1365        torMN  = gen(1);
1366        torMN0 = std(torMN);
1367        L1[ii] = torMN;
1368        L2[ii] = torMN0;
1369        L3[ii] = matrix(kbase(torMN0));
1370        di=dim(torMN0);
1371        dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1372        if (di==0)
1373        {
1374          dbprint(p,"// vdim of Tor_"+string(k)+":       "
1375                    +string(vdim(torMN0)));
1376        }
1377        dbprint(p,"");
1378      }
1379      else
1380      {
1381        N2 = resN[k+1];
1382        if( k==0 ) { torMN=module(tensorMod(M,N)); }
1383        else
1384        {
1385          N1 = resN[k];
1386          col = ncols(N1);
1387
1388//---------- computing homology ----------------------------------------------
1389          imag = tensor(unitmat(nrows(N1)),M);
1390          f = tensor(matrix(N1),unitmat(row));
1391          Im1 = tensor(unitmat(col),M);
1392          Im2 = tensor(matrix(N2),unitmat(row));
1393          ker = modulo(f,imag);
1394          Im  = Im2,Im1;
1395          torMN = modulo(ker,Im);
1396          dbprint(p-1,"// Computing Tor_"+string(k)+
1397                      " (help Tor; gives an explanation):",
1398          "// Let 0 <- coker(M) <- G0 <-M- G1 be the present. of coker(M),",
1399          "// and 0 <- coker(N) <- F0 <-N- F1 <- F2 <- ... a resolution of",
1400          "// coker(N), then Tensor(G0,F"+string(k)+")-->Tensor(G0,F"+
1401          string(k-1)+") is given by:",f,
1402          "// and Tensor(G0,F"+string(k+1)+") + Tensor(G1,F"+string(k)+
1403          ")-->Tensor(G0,F"+string(k)+") is given by:",Im,"");
1404        }
1405
1406        torMN0 = std(torMN);
1407        di=dim(torMN0);
1408        dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1409        if (di==0)
1410        {
1411          dbprint(p,"// vdim of Tor_"+string(k)+":       "
1412                    +string(vdim(torMN0)));
1413        }
1414        dbprint(p,"");
1415
1416//---------- more information -------------------------------------------------
1417        if( size(#)>0 )
1418        { if( vdim(torMN0) >= 0 )
1419          { kb = kbase(torMN0);
1420            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
1421            dbprint(p-1,"// columns of matrix are kbase of Tor_"+
1422                     string(k)+" in Tensor(G0,F"+string(k)+")",kb,"");
1423            L3[ii] = matrix(kb);
1424          }
1425          L2[ii] = torMN0;
1426        }
1427        L1[ii] = torMN;
1428      }
1429    }
1430  }
1431  if( size(#) )
1432  {  if( s>1 ) { L = L1,L2,L3; return(L); }
1433     else { L = torMN,torMN0,matrix(kb); return(L); }
1434  }
1435  else
1436  {  if( s>1 ) { return(L1); }
1437     else { return(torMN); }
1438  }
1439}
1440example
1441{"EXAMPLE:";   echo=2;
1442  int p      = printlevel;
1443  printlevel = 1;
1444  ring r     = 0,(x,y),dp;
1445  ideal i    = x2,y;
1446  ideal j    = x;
1447  list E     = Tor(0..2,i,j);    // Tor_k(r/i,r/j) for k=0,1,2 over r
1448
1449  qring R    = std(i);
1450  ideal j    = fetch(r,j);
1451  module M   = [x,0],[0,x];
1452  printlevel = 2;
1453  module E1  = Tor(1,M,j);       // Tor_1(R^2/M,R/j) over R=r/i
1454
1455  list l     = Tor(3,M,M,1);     // Tor_3(R^2/M,R^2/M) over R=r/i
1456  printlevel = p;
1457}
1458///////////////////////////////////////////////////////////////////////////////
1459proc fitting(module M, int n)
1460"USAGE:   fitting (M,n);  M module, n int
1461RETURN:  ideal, (standard basis of) n-th Fitting ideal of M'=coker(M).
1462EXAMPLE: example fitting; shows an example
1463"
1464{
1465  n=nrows(M)-n;
1466  if(n<=0){return(ideal(1));}
1467  if((n>nrows(M))||(n>ncols(M))){return(ideal(0));}
1468  return(std(minor(M,n)));
1469}
1470example
1471{"EXAMPLE:";   echo=2;
1472  ring R=0,x(0..4),dp;
1473  matrix M[2][4]=x(0),x(1),x(2),x(3),x(1),x(2),x(3),x(4);
1474  print(M);
1475  fitting(M,-1);
1476  fitting(M,0);
1477  fitting(M,1);
1478  fitting(M,2);
1479}
1480///////////////////////////////////////////////////////////////////////////////
1481proc isLocallyFree(matrix S, int r)
1482"USAGE:   isLocallyFree(M,r);  M module, r int
1483RETURN:  1 if M'=coker(M) is locally free of constant rank r;@*
1484         0 if this is not the case.
1485EXAMPLE: example isLocallyFree; shows an example.
1486"
1487{
1488   ideal F=fitting(S,r);
1489   ideal G=fitting(S,r-1);
1490   if((deg(F[1])==0)&&(size(G)==0)){return(1);}
1491   return(0);
1492}
1493example
1494{"EXAMPLE:";   echo=2;
1495  ring R=0,(x,y,z),dp;
1496  matrix M[2][3];     // the presentation matrix
1497  M=x-1,y-1,z,y-1,x-2,x;
1498  ideal I=fitting(M,0); // 0-th Fitting ideal of coker(M)
1499  qring Q=I;
1500  matrix M=fetch(R,M);
1501  isLocallyFree(M,1); // as R/I-module, coker(M) is locally free of rk 1
1502  isLocallyFree(M,0);
1503}
1504///////////////////////////////////////////////////////////////////////////////
1505proc flatteningStrat (module M)
1506"USAGE:   flatteningStrat(M);  M module
1507RETURN:  list of ideals.
1508         The list entries L[1],...,L[r] describe the flattening stratification
1509         of M'=coker(M): setting L[0]=0, L[r+1]=1, the flattening
1510         stratification is given by the open sets Spec(A/V(L[i-1])) \ V(L[i]),
1511         i=1,...,r+1  (A = basering).
1512NOTE:    for more information see the book 'A Singular Introduction to
1513         Commutative Algebra' (by Greuel/Pfister, Springer 2002).
1514EXAMPLE: example flatteningStrat; shows an example
1515"
1516{
1517   list l;
1518   int v,w;
1519   ideal F;
1520   while(1)
1521   {
1522      F=interred(fitting(M,w));
1523      if(F[1]==1){return(l);}
1524      if(size(F)!=0){v++;l[v]=F;}
1525      w++;
1526   }
1527   return(l);
1528}
1529example
1530{"EXAMPLE:";   echo=2;
1531  ring A = 0,x(0..4),dp;
1532  // presentation matrix:
1533  matrix M[2][4] = x(0),x(1),x(2),x(3),x(1),x(2),x(3),x(4);
1534  list L = flatteningStrat(M);
1535  L;
1536}
1537///////////////////////////////////////////////////////////////////////////////
1538proc isFlat(module M)
1539"USAGE:   isFlat(M);  M module
1540RETURN:  1 if M'=coker(M) is flat;@*
1541         0 if this is not the case.
1542EXAMPLE: example isFlat; shows an example.
1543"
1544{
1545   if (size(ideal(M))==0) {return(1);}
1546   int w;
1547   ideal F=fitting(M,0);
1548   while(size(F)==0)
1549   {
1550      w++;
1551      F=fitting(M,w);
1552   }
1553   if (deg(std(F)[1])==0) {return(1);}
1554   return(0);
1555}
1556example
1557{"EXAMPLE:";   echo=2;
1558  ring A = 0,(x,y),dp;
1559  matrix M[3][3] = x-1,y,x,x,x+1,y,x2,xy+x+1,x2+y;
1560  print(M);
1561  isFlat(M);             // coker(M) is not flat over A=Q[x,y]
1562
1563  qring B = std(x2+x-y);   // the ring B=Q[x,y]/<x2+x-y>
1564  matrix M = fetch(A,M);
1565  isFlat(M);             // coker(M) is flat over B
1566
1567  setring A;
1568  qring C = std(x2+x+y);   // the ring C=Q[x,y]/<x2+x+y>
1569  matrix M = fetch(A,M);
1570  isFlat(M);             // coker(M) is not flat over C
1571}
1572///////////////////////////////////////////////////////////////////////////////
1573proc flatLocus(module M)
1574"USAGE:   flatLocus(M);  M module
1575RETURN:  ideal I, s.th. complement of V(I) is flat locus of coker(M).
1576NOTE:    computation is based on Fitting ideals;@*
1577         output is not radical (in general)
1578EXAMPLE: example flatLocus; shows an example
1579"
1580{
1581   if (size(ideal(M))==0) {return(ideal(1));}
1582   int v,w;
1583   ideal F=fitting(M,0);
1584   while(size(F)==0)
1585   {
1586      w++;
1587      F=fitting(M,w);
1588   }
1589   if(typeof(basering)=="qring")
1590   {
1591      for(v=w+1;v<=nrows(M);v++)
1592      {
1593         F=F+intersect(fitting(M,v),quotient(ideal(0),fitting(M,v-1)));
1594      }
1595   }
1596   return(interred(F));
1597}
1598example
1599{"EXAMPLE:";   echo=2;
1600  ring R=0,(x,y,z),dp;
1601  matrix M[2][3]=x,y,z,0,x3,z3;
1602  ideal I=flatLocus(M); // coker(M) is flat outside V(x,yz)
1603  I;                    // computed ideal not radical
1604  ideal J=radical(I);
1605  J;
1606
1607  qring r=std(J);
1608  matrix M=fetch(r,M);
1609  flatLocus(M);         // coker(M) is flat over Spec(Q[x,y,z]/<x,yz>)
1610
1611  isFlat(M);            // flatness test
1612}
1613///////////////////////////////////////////////////////////////////////////////
1614proc isReg(ideal I, module N)
1615"USAGE:   isReg(I,M);  I ideal, M module
1616RETURN:  1 if given (ordered) list of generators for I is coker(M)-sequence;@*
1617         0 if this is not the case.
1618EXAMPLE: example isReg; shows an example.
1619"
1620{
1621  int n=nrows(N);
1622  int i;
1623  while(i<ncols(I))
1624  {
1625     i++;
1626     N=std(N);
1627     if(size(reduce(quotient(N,I[i]),N))!=0){return(0);}
1628     N=N+I[i]*freemodule(n);
1629  }
1630  if (size(reduce(freemodule(n),std(N)))==0){return(0);}
1631  return(1);
1632}
1633example
1634{"EXAMPLE:";   echo=2;
1635  ring R = 0,(x,y,z),dp;
1636  ideal I = x*(y-1),y,z*(y-1);
1637  isReg(I,0);             // given list of generators is Q[x,y,z]-sequence
1638
1639  I = x*(y-1),z*(y-1),y;  // change sorting of generators
1640  isReg(I,0);
1641
1642  ring r = 0,(x,y,z),ds;  // local ring
1643  ideal I=fetch(R,I);
1644  isReg(I,0);             // result independent of sorting of generators
1645}
1646
1647///////////////////////////////////////////////////////////////////////////////
1648// the following static procedures are used by KoszulHomology:
1649//  * binom_int  (binomial coeff. as integer, or -1 if too large)
1650//  * basisNumber
1651//  * basisElement
1652//  * KoszulMap
1653// for details, see 'A Singular Introduction to Commutative Algebra' (by
1654// Greuel/Pfister, Springer 2002), Chapter 7
1655
1656static proc binom_int(int n, int p)
1657{
1658  bigint s = binomial(n,p);
1659  int a=int(s);
1660  if ((s!=0)&&(a==0)) { return(-1); }
1661  return(a);
1662}
1663
1664static proc basisNumber(int n,intvec v)
1665{
1666  int p=size(v);
1667  if(p==1){return(v[1]);}
1668  int j=n-1;
1669  int b;
1670  while(j>=n-v[1]+1)
1671  {
1672    b=b+binom_int(j,p-1);
1673    j--;
1674  }
1675  intvec w=v-v[1];
1676  w=w[2..size(w)];
1677  b=b+basisNumber(n-v[1],w);
1678  return(b);
1679}
1680
1681static proc basisElement(int n,int p,int N)
1682{
1683  if(p==1){return(N);}
1684  int s,R;
1685  while(R<N)
1686  {
1687     s++;
1688     R=R+binom_int(n-s,p-1);
1689  }
1690  R=N-R+binom_int(n-s,p-1);
1691  intvec v=basisElement(n-s,p-1,R);
1692  intvec w=s,v+s;
1693  return(w);
1694}
1695
1696proc KoszulMap(ideal x,int p)
1697{
1698  int n=size(x);
1699  int a=binom_int(n,p-1);
1700  int b=binom_int(n,p);
1701
1702  matrix M[a][b];
1703
1704  if(p==1){M=x;return(M);}
1705  int j,k;
1706  intvec v,w;
1707  for(j=1;j<=b;j++)
1708  {
1709    v=basisElement(n,p,j);
1710    w=v[2..p];
1711    M[basisNumber(n,w),j]=x[v[1]];
1712    for(k=2;k<p;k++)
1713    {
1714      w=v[1..k-1],v[k+1..p];
1715      M[basisNumber(n,w),j]=(-1)^(k-1)*x[v[k]];
1716    }
1717    w=v[1..p-1];
1718    M[basisNumber(n,w),j]=(-1)^(p-1)*x[v[p]];
1719  }
1720  return(M);
1721}
1722///////////////////////////////////////////////////////////////////////////////
1723
1724proc KoszulHomology(ideal x, module M, int p)
1725"USAGE:   KoszulHomology(I,M,p);  I ideal, M module, p int
1726COMPUTE: A presentation of the p-th Koszul homology module H_p(f_1,...,f_k;M'),
1727         where M'=coker(M) and f_1,...,f_k are the given (ordered list
1728         of non-zero) generators of the ideal I.
1729         The computed presentation is minimized via prune.
1730         In particular, if H_p(f_1,...,f_k;M')=0 then the return value is 0.
1731RETURN:  module H, s.th. coker(H) = H_p(f_1,...,f_k;M').
1732NOTE:    size of input ideal has to be <= 20.
1733EXAMPLE: example KoszulHomology; shows an example.
1734{
1735  x=simplify(x,2);
1736  int n     = size(x);
1737  if (n==0)
1738  {
1739    ERROR("// KoszulHomology only for non-zero ideals");
1740  }
1741  if (n>20)
1742  {
1743    ERROR("// too many generators in input ideal");
1744  }
1745  if (p>n)
1746  {
1747    module hom=0;
1748    return(hom);
1749  }
1750
1751  int a     = binom_int(n,p-1);  // n over p-1 independent of char(basering)
1752  int b     = binom_int(n,p);
1753
1754  matrix N  = matrix(M);
1755  module ker= freemodule(nrows(N));
1756  if(p!=0)
1757  {
1758    module im= tensor(unitmat(a),N);
1759    module f = tensor(KoszulMap(x,p),unitmat(nrows(N)));
1760    ker      = modulo(f,im);
1761  }
1762  module im1 = tensor(unitmat(b),N);
1763  module im2 = tensor(KoszulMap(x,p+1),unitmat(nrows(N)));
1764  module hom = modulo(ker,im1+im2);
1765  hom        = prune(hom);
1766  return(hom);
1767}
1768example
1769{"EXAMPLE:";   echo=2;
1770  ring R=0,x(1..3),dp;
1771  ideal x=maxideal(1);
1772  module M=0;
1773  KoszulHomology(x,M,0);  // H_0(x,R), x=(x_1,x_2,x_3)
1774
1775  KoszulHomology(x,M,1);  // H_1(x,R), x=(x_1,x_2,x_3)
1776
1777  qring S=std(x(1)*x(2));
1778  module M=0;
1779  ideal x=maxideal(1);
1780  KoszulHomology(x,M,1);
1781
1782  KoszulHomology(x,M,2);
1783}
1784///////////////////////////////////////////////////////////////////////////////
1785proc depth(module M,ideal #)
1786"USAGE:   depth(M,[I]);  M module, I ideal
1787RETURN:  int,
1788         - if called with 1 argument: the depth of M'=coker(M) w.r.t. the
1789         maxideal in the basering (which is then assumed to be local)@*
1790         - if called with 2 arguments: the depth of M'=coker(M) w.r.t. the
1791         ideal I.
1792NOTE:    procedure makes use of KoszulHomology.
1793EXAMPLE: example depth; shows an example.
1794"
1795{
1796  ideal m=maxideal(1);
1797  int n=size(m);
1798  int i;
1799
1800  if (size(#)==0)
1801  {
1802    // depth(M') over local basering
1803    while(i<n)
1804    {
1805      i++;
1806      if(size(KoszulHomology(m,M,i))==0){return(n-i+1);}
1807    }
1808    return(0);
1809  }
1810
1811  ideal I=simplify(#,2);
1812  while(i<size(I))
1813  {
1814    i++;
1815    if(size(KoszulHomology(I,M,i))==0){return(size(I)-i+1);}
1816  }
1817  return(0);
1818}
1819example
1820{"EXAMPLE:";   echo=2;
1821  ring R=0,(x,y,z),dp;
1822  ideal I=x2,xy,yz;
1823  module M=0;
1824  depth(M,I);   // depth(<x2,xy,yz>,Q[x,y,z])
1825  ring r=0,(x,y,z),ds;  // local ring
1826  matrix M[2][2]=x,xy,1+yz,0;
1827  print(M);
1828  depth(M);     // depth(maxideal,coker(M))
1829  ideal I=x;
1830  depth(M,I);   // depth(<x>,coker(M))
1831  I=x+z;
1832  depth(M,I);   // depth(<x+z>,coker(M))
1833}
1834
1835///////////////////////////////////////////////////////////////////////////////
1836proc isCM(module M)
1837"USAGE:   isCM(M);  M module
1838RETURN:  1 if M'=coker(M) is Cohen-Macaulay;@*
1839         0 if this is not the case.
1840ASSUME:  basering is local.
1841EXAMPLE: example isCM; shows an example
1842"
1843{
1844  // test if basering is local:
1845  ideal m=maxideal(1);
1846  int i;
1847  poly f=1;
1848  for (i=1; i<=size(m); i++)
1849  {
1850    f=f+m[i];
1851  }
1852  if (ord(f)>0)
1853  {
1854    print("// basering must be local -- result has no meaning");
1855    return(0);
1856  }
1857
1858  return(depth(M)==dim(std(Ann(M))));
1859}
1860example
1861{"EXAMPLE:";   echo=2;
1862  ring R=0,(x,y,z),ds;  // local ring R = Q[x,y,z]_<x,y,z>
1863  module M=xz,yz,z2;
1864  isCM(M);             // test if R/<xz,yz,z2> is Cohen-Macaulay
1865
1866  M=x2+y2,z7;          // test if R/<x2+y2,z7> is Cohen-Macaulay
1867  isCM(M);
1868}
1869
1870proc canonMap(list #)
1871"USAGE:  canonMap(id); id= ideal/module,
1872RETURN:  a list L, the kernel in two different representations and
1873@*       the cokernel of the canonical map
1874@*       M ---> Ext^c_R(Ext^c_R(M,R),R) given by presentations
1875@*       Here M is the R-module (R=basering) given by the presentation
1876@*       defined by id, i.e. M=R/id resp. M=R^n/id
1877@*       c is the codimension of M
1878@*       L[1] is the preimage of the kernel in R resp. R^n
1879@*       L[2] is a presentation of the kernel
1880@*       L[3] is a presentation of the cokernel
1881EXAMPLE: example canonMap; shows an example
1882"
1883{
1884   module M=#[1];
1885   int c=nvars(basering)-dim(std(M));
1886   if(c==0)
1887   {
1888      module K=syz(transpose(M));
1889      module Ke=syz(transpose(K));
1890      module Co=modulo(syz(transpose(syz(K))),transpose(K));
1891    }
1892   else
1893   {
1894      int i;
1895      resolution F=mres(M,c+1);
1896      module K=syz(transpose(F[c+1]));
1897      K=simplify(reduce(K,std(transpose(F[c]))),2);
1898      module A=modulo(K,transpose(F[c]));
1899      resolution G=nres(A,c+1);
1900      for(i=1;i<=c;i++)
1901      {
1902         K=lift(transpose(F[c-i+1]),K*G[i]);
1903      }
1904      module Ke=modulo(transpose(K),transpose(G[c]));
1905      module Co=modulo(syz(transpose(G[c+1])),transpose(K)+transpose(G[c]));
1906   }
1907   return(list(Ke,prune(modulo(Ke,M)),prune(Co)));
1908
1909}
1910example
1911{ "EXAMPLE:"; echo = 2;
1912   ring s=0,(x,y),dp;
1913   ideal i = x,y;
1914   canonMap(i);
1915   ring R = 0,(x,y,z,w),dp;
1916   ideal I1 = x,y;
1917   ideal I2 = z,w;
1918   ideal I = intersect(I1,I2);
1919   canonMap(I);
1920   module M = syz(I);
1921   canonMap(M);
1922   ring S = 0,(x,y,z,t),Wp(3,4,5,1);
1923   ideal I = x-t3,y-t4,z-t5;
1924   ideal J = eliminate(I,t);
1925   ring T = 0,(x,y,z),Wp(3,4,5);
1926   ideal p = imap(S,J);
1927   ideal p2 = p^2;
1928   canonMap(p2);
1929}
1930
Note: See TracBrowser for help on using the repository browser.