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

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