source: git/Singular/LIB/homolog.lib @ 01cc3b2

fieker-DuValspielwiese
Last change on this file since 01cc3b2 was 0a0355d, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: from 2-0 branch git-svn-id: file:///usr/local/Singular/svn/trunk@7602 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 52.6 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: homolog.lib,v 1.16 2004-11-04 16:51:28 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         or the 2 commands: list L=mres(M,2);  Ext_R(k,L[2]);
469EXAMPLE: example Ext_R; shows an example
470"
471{
472  // In case M is known to be a SB, set attrib(M,"isSB",1); in order to
473  // avoid  unnecessary SB computations
474
475//------------ initialisation -------------------------------------------------
476  module m1,m2,ret,ret0;
477  matrix ker,kb;
478  list L1,L2,L3,L,resl,K;
479  int k,max,ii,t1,t2,di;
480  int s = size(v);
481  intvec v1 = sort(v)[1];
482  max = v1[s];                 // the maximum integer occurring in intvec v
483  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
484  // --------------- Variante mit sres
485  for( ii=1; ii<=size(#); ii++ )
486  {
487    t2=1; // return a list if t2=1
488    if( typeof(#[ii])=="string" )
489    {
490      if ( #[ii]=="sres" ) { t1=1; t2=0; } // use sres instead of mres if t1=1
491    }
492  }
493//----------------- compute resolution of coker(M) ----------------------------
494  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
495  if( t1==1 )
496  {
497    if( attrib(M,"isSB")==0 ) { M=std(M); }
498    resl = sres(M,max+1);
499  }
500  else { resl = mres(M,max+1); }
501  for( ii=1; ii<=s; ii++ )
502  {
503//-----------------  apply Hom(_,R) at k-th place -----------------------------
504    k=v[ii];
505    if( k<0 )                   // Ext^k=0 for negative k
506    {
507      dbprint(p-1,"// Ext^i=0 for i<0!");
508      ret    = gen(1);
509      ret0   = std(ret);
510      L1[ii] = ret;
511      L2[ii] = ret0;
512      L3[ii] = matrix(kbase(ret0));
513      di=dim(ret0);
514      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
515      if (di==0)
516      {
517        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(ret0)));
518      }
519      dbprint(p,"");
520    }
521    else
522    {
523      m2 = transpose(resl[k+1]);
524      if( k==0 ) { m1=0*gen(nrows(m2)); }
525      else { m1 = transpose(resl[k]); }
526//----------------- presentation of ker(m2)/im(m1) ----------------------------
527      ker = syz(m2);
528      ret = modulo(ker,m1);
529      dbprint(p-1,"// Computing Ext^"+string(k)+":",
530         "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of M,",
531         "// then F"+string(k)+"*-->F"+string(k+1)+"* is given by:",m2,
532         "// and F"+string(k-1)+"*-->F"+string(k)+"* is given by:",m1,"");
533      ret0 = std(ret);
534      di=dim(ret0);
535      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
536      if (di==0)
537      {
538        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(ret0)));
539      }
540      dbprint(p,"");
541      if( t2 )
542      {
543         if( vdim(ret0)>=0 )
544         {
545            kb = kbase(ret0);
546            if ( size(ker)!=0 ) { kb = matrix(ker)*kb; }
547            dbprint(p-1,
548              "// columns of matrix are kbase of Ext^"+string(k)+" in F"
549              +string(k)+"*:",kb,"");
550            L3[ii] = kb;
551         }
552         L2[ii] = ret0;
553      }
554      L1[ii] = ret;
555    }
556  }
557  if( t2 )
558  {
559     if( s>1 ) { L = L1,L2,L3; return(L); }
560     else { L = ret,ret0,kb; return(L); }
561  }
562  else
563  {
564     if( s>1 ) { return(L1); }
565     else { return(ret); }
566  }
567}
568example
569{"EXAMPLE:";     echo=2;
570  int p      = printlevel;
571  printlevel = 1;
572  ring r     = 0,(x,y,z),dp;
573  ideal i    = x2y,y2z,z3x;
574  module E   = Ext_R(1,i);    //computes Ext^1(r/i,r)
575  is_zero(E);
576
577  qring R    = std(x2+yz);
578  intvec v   = 0,2;
579  printlevel = 2;             //shows what is going on
580  ideal i    = x,y,z;         //computes Ext^i(r/(x,y,z),r/(x2+yz)), i=0,2
581  list L     = Ext_R(v,i,1);  //over the qring R=r/(x2+yz), std and kbase
582  printlevel = p;
583}
584///////////////////////////////////////////////////////////////////////////////
585
586proc Ext (intvec v, module M, module N, list #)
587"USAGE:   Ext(v,M,N[,any]);  v int resp. intvec, M,N modules
588COMPUTE: A presentation of Ext^k(M',N'); for k=v[1],v[2],... where
589         M'=coker(M) and N'=coker(N). Let
590@example
591       0 <-- M' <-- F0 <-M-- F1 <-- F2 <--... ,
592       0 <-- N' <-- G0 <--N- G1
593@end example
594         be a free resolution of M', resp. a presentation of N'. Consider
595         the commutative diagram
596@example
597           0                  0                  0
598           |^                 |^                 |^
599   --> Hom(Fk-1,N') -Ak-> Hom(Fk,N') -Ak+1-> Hom(Fk+1,N')
600           |^                 |^                 |^
601   --> Hom(Fk-1,G0) -Ak-> Hom(Fk,G0) -Ak+1-> Hom(Fk+1,G0)
602                              |^                 |^
603                              |C                 |B
604                          Hom(Fk,G1) ------> Hom(Fk+1,G1)
605
606      (Ak,Ak+1 induced by M and B,C induced by N).
607@end example
608         Let K=modulo(Ak+1,B), J=module(Ak)+module(C) and Ext=modulo(K,J),
609         then we have exact sequences
610@example
611    R^p --K-> Hom(Fk,G0) --Ak+1-> Hom(Fk+1,G0)/im(B),
612
613    R^q -Ext-> R^p --K-> Hom(Fk,G0)/(im(Ak)+im(C)).
614@end example
615         Hence, Ext presents Ext^k(M',N').
616RETURN:  - module Ext, a presentation of Ext^k(M',N') if v is of type int@*
617         - a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.@*
618         - In case of a third argument of any type return a list l:
619@format
620             l[1] = module Ext/list of Ext^k
621             l[2] = SB of Ext/list of SB of Ext^k
622             l[3] = matrix/list of matrices, each representing a kbase of Ext^k
623                       (if finite dimensional)
624@end format
625DISPLAY: printlevel >=0: dimension, vdim of Ext^k for each k  (default).
626@*       printlevel >=1: matrices Ak, Ak+1 and kbase of Ext^k in Hom(Fk,G0)
627         (if finite dimensional)
628NOTE:    In order to compute Ext^k(M,N) use the command Ext(k,syz(M),syz(N));
629         or: list P=mres(M,2); list Q=mres(N,2); Ext(k,P[2],Q[2]);
630EXAMPLE: example Ext; shows an example
631"
632{
633//---------- initialisation ---------------------------------------------------
634  int k,max,ii,l,row,col,di;
635  module A,B,C,D,M1,M2,N1,ker,imag,extMN,extMN0;
636  matrix kb;
637  list L1,L2,L3,L,resM,K;
638  ideal  test1;
639  intmat Be;
640  int s = size(v);
641  intvec v1 = sort(v)[1];
642  max = v1[s];                    // the maximum integer occurring in intvec v
643  int p = printlevel-voice+3;     // p=printlevel+1 (default: p=1)
644//---------- test: coker(N)=basering, coker(N)=0 ? ----------------------------
645  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
646  N1 = std(N);
647  if( size(N1)==0 )      //coker(N)=basering, in this case proc Ext_R is faster
648  { printlevel=printlevel+1;
649    if( size(#)==0 )
650    { def E = Ext_R(v,M);
651      printlevel=printlevel-1;
652      return(E);
653    }
654    else
655    { def E = Ext_R(v,M,#[1]);
656       printlevel=printlevel-1;
657       return(E);
658     }
659  }
660  if( dim(N1)==-1 )                          //coker(N)=0, all Ext-groups are 0
661  { dbprint(p-1,"2nd module presents 0, hence Ext^k=0, for all k");
662    for( ii=1; ii<=s; ii++ )
663    { k=v[ii];
664      extMN    = gen(1);
665      extMN0   = std(extMN);
666      L1[ii] = extMN;
667      L2[ii] = extMN0;
668      L3[ii] = matrix(kbase(extMN0));
669      di=dim(extMN0);
670      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
671      if (di==0)
672      {
673        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(extMN0)));
674      }
675      dbprint(p,"");
676    }
677  }
678  else
679  {
680    if( size(N1) < size(N) ) { N=N1;}
681    row = nrows(N);
682//---------- resolution of M -------------------------------------------------
683    resM = mres(M,max+1);
684    for( ii=1; ii<=s; ii++ )
685    { k=v[ii];
686      if( k<0  )                                   // Ext^k is 0 for negative k
687      { dbprint(p-1,"// Ext^k=0 for k<0!");
688        extMN    = gen(1);
689        extMN0   = std(extMN);
690        L1[ii] = extMN;
691        L2[ii] = extMN0;
692        L3[ii] = matrix(kbase(extMN0));
693        di=dim(extMN0);
694        dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
695        if (di==0)
696        {
697          dbprint(p,"// vdim of Ext^"+string(k)+":       "
698                  +string(vdim(extMN0)));
699        }
700        dbprint(p,"");
701      }
702      else
703      { M2 = resM[k+1];
704        if( k==0 ) { M1=0; }
705        else { M1 = resM[k]; }
706        col = nrows(M2);
707        D = kohom(N,col);
708//---------- computing homology ----------------------------------------------
709        imag  = kontrahom(M1,row);
710        A     = kontrahom(M2,row);
711        B     = kohom(N,ncols(M2));
712        ker   = modulo(A,B);
713        imag  = imag,D;
714        extMN = modulo(ker,imag);
715        dbprint(p-1,"// Computing Ext^"+string(k)+
716         " (help Ext; gives an explanation):",
717      "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of coker(M),",
718      "// and 0<--coker(N)<--G0<--G1 a presentation of coker(N),",
719      "// then Hom(F"+string(k)+",G0)-->Hom(F"+string(k+1)+
720      ",G0) is given by:",A,
721      "// and Hom(F"+string(k-1)+",G0) + Hom(F"+string(k)+",G1)-->Hom(F"
722                                      +string(k)+",G0) is given by:",imag,"");
723        extMN0 = std(extMN);
724        di=dim(extMN0);
725        dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
726        if (di==0)
727        {
728          dbprint(p,"// vdim of Ext^"+string(k)+":       "
729                  +string(vdim(extMN0)));
730        }
731        dbprint(p,"");
732
733//---------- more information -------------------------------------------------
734        if( size(#)>0 )
735        { if( vdim(extMN0) >= 0 )
736          { kb = kbase(extMN0);
737            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
738            dbprint(p-1,"// columns of matrix are kbase of Ext^"+
739                     string(k)+" in Hom(F"+string(k)+",G0)",kb,"");
740            if( p>0 )
741            { for (l=1;l<=ncols(kb);l=l+1)
742              {
743                "// element",l,"of kbase of Ext^"+string(k)+" in Hom(F"+string(k)+",G0)";
744                "// as matrix: F"+string(k)+"-->G0";
745                print(matrix(ideal(kb[l]),row,col));
746              }
747              "";
748            }
749            L3[ii] = matrix(kb);
750          }
751          L2[ii] = extMN0;
752        }
753        L1[ii] = extMN;
754      }
755    }
756  }
757  if( size(#) )
758  {  if( s>1 ) { L = L1,L2,L3; return(L); }
759     else { L = extMN,extMN0,matrix(kb); return(L); }
760  }
761  else
762  {  if( s>1 ) { return(L1); }
763     else { return(extMN); }
764  }
765}
766example
767{"EXAMPLE:";   echo=2;
768  int p      = printlevel;
769  printlevel = 1;
770  ring r     = 0,(x,y),dp;
771  ideal i    = x2-y3;
772  ideal j    = x2-y5;
773  list E     = Ext(0..2,i,j);    // Ext^k(r/i,r/j) for k=0,1,2 over r
774  qring R    = std(i);
775  ideal j    = fetch(r,j);
776  module M   = [-x,y],[-y2,x];
777  printlevel = 2;
778  module E1  = Ext(1,M,j);       // Ext^1(R^2/M,R/j) over R=r/i
779  list l     = Ext(4,M,M,1);     // Ext^4(R^2/M,R^2/M) over R=r/i
780  printlevel = p;
781}
782///////////////////////////////////////////////////////////////////////////////
783
784proc Hom (module M, module N, list #)
785"USAGE:   Hom(M,N,[any]);  M,N=modules
786COMPUTE: A presentation of Hom(M',N'), M'=coker(M), N'=coker(N) as follows:
787         let
788@example
789   F1 --M-> F0 -->M' --> 0,    G1 --N-> G0 --> N' --> 0
790@end example
791         be presentations of M' and N'. Consider
792@example
793                                  0               0
794                                  |^              |^
795       0 --> Hom(M',N') ----> Hom(F0,N') ----> Hom(F1,N')
796                                  |^              |^
797  (A:  induced by M)          Hom(F0,G0) --A-> Hom(F1,G0)
798                                  |^              |^
799  (B,C:induced by N)              |C              |B
800                              Hom(F0,G1) ----> Hom(F1,G1)
801
802@end example
803         Let D=modulo(A,B) and Hom=modulo(D,C), then we have exact sequences
804@example
805   R^p  --D-> Hom(F0,G0) --A-> Hom(F1,G0)/im(B),
806
807 R^q -Hom-> R^p --D-> Hom(F0,G0)/im(C) --A-> Hom(F1,G0)/im(B).
808@end example
809         Hence Hom presents Hom(M',N')
810RETURN:  module Hom, a presentation of Hom(M',N'), resp., in case of
811         3 arguments, a list l (of size <=3):
812@format
813           - l[1] = Hom
814           - l[2] = SB of Hom
815           - l[3] = kbase of coker(Hom) (if finite dimensional, not 0),
816                    represented by elements in Hom(F0,G0) via mapping D
817@end format
818DISPLAY: printlevel >=0: (affine) dimension of Hom  (default)
819@*       printlevel >=1: D and C and kbase of coker(Hom) in Hom(F0,G0)
820@*       printlevel >=2: elements of kbase of coker(Hom) as matrix :F0-->G0
821NOTE:    DISPLAY is as described only for a direct call of 'Hom'. Calling 'Hom'
822         from another proc has the same effect as decreasing printlevel by 1.
823EXAMPLE: example Hom;  shows examples
824"
825{
826//---------- initialisation ---------------------------------------------------
827  int l,p,di;
828  matrix kb;
829  module A,B,C,D,homMN,homMN0;
830  list L;
831//---------- computation of Hom -----------------------------------------------
832  B = kohom(N,ncols(M));
833  A = kontrahom(M,nrows(N));
834  C = kohom(N,nrows(M));
835  D = modulo(A,B);
836  homMN = modulo(D,C);
837  homMN0= std(homMN);
838  p = printlevel-voice+3;       // p=printlevel+1 (default: p=1)
839  di= dim(homMN0);
840  dbprint(p,"// dimension of Hom:  "+string(di));
841  if (di==0)
842  {
843    dbprint(p,"// vdim of Hom:       "+string(vdim(homMN0)));
844  }
845  dbprint(p,"");
846  dbprint(p-1,
847   "// given  F1 --M-> F0 -->M'--> 0 and  G1 --N-> G0 -->N'--> 0,",
848   "// show D = ker( Hom(F0,G0) --> Hom(F1,G0)/im(Hom(F1,G1)->Hom(F1,G0)) )",D,
849   "// show C = im ( Hom(F0,G1) --> Hom(F0,G0) )",C,"");
850//---------- extra output if size(#)>0 ----------------------------------------
851  if( size(#)>0 )
852  {
853     if( vdim(homMN0)>0 )
854     {
855        kb = kbase(homMN0);
856        kb = matrix(D)*kb;
857        if( p>2 )
858        {
859          for (l=1;l<=ncols(kb);l=l+1)
860          {
861            "// element",l,"of kbase of Hom in Hom(F0,G0) as matrix: F0-->G0:";
862             print(matrix(ideal(kb[l]),nrows(N),nrows(M)));
863          }
864        }
865        else
866        {
867          dbprint(p-1,"// columns of matrix are kbase of Hom in Hom(F0,G0)",
868                      kb);
869        }
870        L=homMN,homMN0,kb;
871        return(L);
872     }
873     L=homMN,homMN0;
874     return(L);
875  }
876  return(homMN);
877}
878example
879{"EXAMPLE:";  echo = 2;
880  int p     = printlevel;
881  printlevel= 1;   //in 'example proc' printlevel has to be increased by 1
882  ring r    = 0,(x,y),dp;
883  ideal i   = x2-y3,xy;
884  qring q   = std(i);
885  ideal i   = fetch(r,i);
886  module M  = [-x,y],[-y2,x],[x3];
887  module H  = Hom(M,i);
888  print(H);
889
890  printlevel= 2;
891  list L    = Hom(M,i,1);"";
892
893  printlevel=1;
894  ring s    = 3,(x,y,z),(c,dp);
895  ideal i   = jacob(ideal(x2+y5+z4));
896  qring rq=std(i);
897  matrix M[2][2]=xy,x3,5y,4z,x2;
898  matrix N[3][2]=x2,x,y3,3xz,x2z,z;
899  print(M);
900  print(N);
901  list l=Hom(M,N,1);
902  printlevel = p;
903}
904///////////////////////////////////////////////////////////////////////////////
905
906proc homology (matrix A,matrix B,module M,module N,list #)
907"USAGE:   homology(A,B,M,N);
908COMPUTE: Let M and N be submodules of R^m and R^n presenting M'=R^m/M, N'=R^n/N
909         (R=basering) and let A,B matrices inducing maps
910@example
911    R^k --A--> R^m --B--> R^n.
912@end example
913         Compute a presentation of the module
914@example
915    ker(B)/im(A) := ker(M'/im(A) --B--> N'/im(BM)+im(BA)).
916@end example
917         If B induces a map M'-->N' (i.e BM=0) and if im(A) is contained in
918         ker(B) (that is, BA=0) then ker(B)/im(A) is the homology of the
919         complex
920@example
921    R^k--A-->M'--B-->N'.
922@end example
923RETURN:  module H, a presentation of ker(B)/im(A).
924NOTE:    homology returns a free module of rank m if ker(B)=im(A).
925EXAMPLE: example homology; shows examples
926"
927{
928  module ker,ima;
929  ker = modulo(B,N);
930  ima = A,M;
931  return(modulo(ker,ima));
932}
933example
934{"EXAMPLE";    echo=2;
935  ring r;
936  ideal id=maxideal(4);
937  qring qr=std(id);
938  module N=maxideal(3)*freemodule(2);
939  module M=maxideal(2)*freemodule(2);
940  module B=[2x,0],[x,y],[z2,y];
941  module A=M;
942  module H=homology(A,B,M,N);
943  H=std(H);
944  // dimension of homology:
945  dim(H);
946  // vector space dimension:
947  vdim(H);
948
949  ring s=0,x,ds;
950  qring qs=std(x4);
951  module A=[x];
952  module B=A;
953  module M=[x3];
954  module N=M;
955  homology(A,B,M,N);
956}
957//////////////////////////////////////////////////////////////////////////////
958
959proc kernel (matrix A,module M,module N)
960"USAGE:   kernel(A,M,N);
961COMPUTE: Let M and N be submodules of R^m and R^n, presenting M'=R^m/M,
962         N'=R^n/N (R=basering), and let A:R^m-->R^n be a matrix inducing a
963         map A':M'-->N'. Then kernel(A,M,N); computes a presentation K of
964         ker(A') as in the commutative diagram:
965@example
966          ker(A') --->  M' --A'--> N'
967             |^         |^         |^
968             |          |          |
969             R^r  ---> R^m --A--> R^n
970             |^         |^         |^
971             |K         |M         |N
972             |          |          |
973             R^s  ---> R^p -----> R^q
974@end example
975RETURN:  module K, a presentation of ker(A':coker(M)->coker(N)).
976EXAMPLE: example kernel; shows examples.
977"
978{
979  module M1 = modulo(A,N);
980  return(modulo(M1,M));
981}
982example
983{"EXAMPLE";    echo=2;
984  ring r;
985  module N=[2x,x],[0,y];
986  module M=maxideal(1)*freemodule(2);
987  matrix A[2][2]=2x,0,x,y,z2,y;
988  module K=kernel(A,M,N);
989  // dimension of kernel:
990  dim(std(K));
991  // vector space dimension of kernel:
992  vdim(std(K));
993  print(K);
994}
995///////////////////////////////////////////////////////////////////////////////
996
997proc kohom (matrix M, int j)
998"USAGE:   kohom(A,k); A=matrix, k=integer
999RETURN:  matrix Hom(R^k,A), i.e. let A be a matrix defining a map F1-->F2
1000         of free R-modules, then the matrix of Hom(R^k,F1)-->Hom(R^k,F2)
1001         is computed (R=basering).
1002EXAMPLE: example kohom; shows an example.
1003"
1004{
1005  if (j==1)
1006  { return(M);}
1007  if (j>1)
1008  { return(tensor(M,diag(1,j))); }
1009  else { return(0);}
1010}
1011example
1012{"EXAMPLE:";   echo=2;
1013  ring r;
1014  matrix n[2][3]=x,y,5,z,77,33;
1015  print(kohom(n,3));
1016}
1017///////////////////////////////////////////////////////////////////////////////
1018
1019proc kontrahom (matrix M, int j)
1020"USAGE:   kontrahom(A,k); A=matrix, k=integer
1021RETURN:  matrix Hom(A,R^k), i.e. let A be a matrix defining a map F1-->F2 of
1022         free  R-modules, then the matrix of Hom(F2,R^k)-->Hom(F1,R^k) is
1023         computed (R=basering).
1024EXAMPLE: example kontrahom; shows an example.
1025"
1026{
1027if (j==1)
1028  { return(transpose(M));}
1029  if (j>1)
1030  { return(transpose(tensor(diag(1,j),M)));}
1031  else { return(0);}
1032}
1033example
1034{"EXAMPLE:";  echo=2;
1035  ring r;
1036  matrix n[2][3]=x,y,5,z,77,33;
1037  print(kontrahom(n,3));
1038}
1039///////////////////////////////////////////////////////////////////////////////
1040///////////////////////////////////////////////////////////////////////////////
1041
1042proc tensorMod(module Phi, module Psi)
1043"USAGE:   tensorMod(M,N);  M,N modules
1044COMPUTE: presentation matrix A of the tensor product T of the modules
1045         M'=coker(M), N'=coker(N): if matrix(M) defines a map M: R^r-->R^s and
1046         matrix(N) defines a map N: R^p-->R^q, then A defines a presentation
1047@example
1048         R^(sp+rq) --A-> R^(sq)  --> T --> 0 .
1049@end example
1050RETURN:  matrix A satisfying coker(A) = tensorprod(coker(M),coker(N)) .
1051EXAMPLE: example tensorMod; shows an example.
1052"
1053{
1054   int s=nrows(Phi);
1055   int q=nrows(Psi);
1056   matrix A=tensor(unitmat(s),Psi);
1057   matrix B=tensor(Phi,unitmat(q));
1058   matrix R=concat(A,B);
1059   return(R);
1060}
1061example
1062{"EXAMPLE:";  echo=2;
1063  ring A=0,(x,y,z),dp;
1064  matrix M[3][3]=1,2,3,4,5,6,7,8,9;
1065  matrix N[2][2]=x,y,0,z;
1066  print(M);
1067  print(N);
1068  print(tensorMod(M,N));
1069}
1070///////////////////////////////////////////////////////////////////////////////
1071proc Tor(intvec v, module M, module N, list #)
1072"USAGE:   Tor(v,M,N[,any]);  v int resp. intvec, M,N modules
1073COMPUTE: a presentation of Tor_k(M',N'), for k=v[1],v[2],... , where
1074         M'=coker(M) and N'=coker(N): let
1075@example
1076       0 <-- M' <-- G0 <-M-- G1
1077       0 <-- N' <-- F0 <--N- F1 <-- F2 <--...
1078@end example
1079         be a presentation of M', resp. a free resolution of N', and consider
1080         the commutative diagram
1081@example
1082          0                    0                    0
1083          |^                   |^                   |^
1084  Tensor(M',Fk+1) -Ak+1-> Tensor(M',Fk) -Ak-> Tensor(M',Fk-1)
1085          |^                   |^                   |^
1086  Tensor(G0,Fk+1) -Ak+1-> Tensor(G0,Fk) -Ak-> Tensor(G0,Fk-1)
1087                               |^                   |^
1088                               |C                   |B
1089                          Tensor(G1,Fk) ----> Tensor(G1,Fk-1)
1090
1091       (Ak,Ak+1 induced by N and B,C induced by M).
1092@end example
1093         Let K=modulo(Ak,B), J=module(C)+module(Ak+1) and Tor=modulo(K,J),
1094         then we have exact sequences
1095@example
1096    R^p  --K-> Tensor(G0,Fk) --Ak-> Tensor(G0,Fk-1)/im(B),
1097
1098    R^q -Tor-> R^p --K-> Tensor(G0,Fk)/(im(C)+im(Ak+1)).
1099@end example
1100         Hence, Tor presents Tor_k(M',N').
1101RETURN:  - if v is of type int: module Tor, a presentation of Tor_k(M',N');@*
1102         - if v is of type intvec: a list of Tor_k(M',N') (k=v[1],v[2],...);@*
1103         - in case of a third argument of any type: list l with
1104@format
1105     l[1] = module Tor/list of Tor_k(M',N'),
1106     l[2] = SB of Tor/list of SB of Tor_k(M',N'),
1107     l[3] = matrix/list of matrices, each representing a kbase of Tor_k(M',N')
1108                (if finite dimensional), or 0.
1109@end format
1110DISPLAY: printlevel >=0: (affine) dimension of Tor_k for each k  (default).
1111@*       printlevel >=1: matrices Ak, Ak+1 and kbase of Tor_k in Tensor(G0,Fk)
1112         (if finite dimensional).
1113NOTE:    In order to compute Tor_k(M,N) use the command Tor(k,syz(M),syz(N));
1114         or: list P=mres(M,2); list Q=mres(N,2); Tor(k,P[2],Q[2]);
1115EXAMPLE: example Tor; shows an example
1116{
1117//---------- initialisation ---------------------------------------------------
1118  int k,max,ii,l,row,col,di;
1119  module A,B,C,D,N1,N2,M1,ker,imag,Im,Im1,Im2,f,torMN,torMN0;
1120  matrix kb;
1121  list L1,L2,L3,L,resN,K;
1122  ideal  test1;
1123  intmat Be;
1124  int s = size(v);
1125  intvec v1 = sort(v)[1];
1126  max = v1[s];                   // maximum integer occurring in intvec v
1127  int p = printlevel-voice+3;    // p=printlevel+1 (default: p=1)
1128
1129//---------- test: coker(M)=basering, coker(M)=0 ? ----------------------------
1130  if( max<0 ) { dbprint(p,"// Tor_i=0 for i<0!"); return([1]); }
1131  M1 = std(M);
1132
1133  if( size(M1)==0 or size(N)==0 )  // coker(M)=basering ==> Tor_i=0 for i>0
1134  {
1135    dbprint(p-1,"// one of the modules M',N' is free, hence Tor_i=0 for i<>0");
1136    for( ii=1; ii<=s; ii++ )
1137    {
1138      k=v[ii];
1139      if (k==0) { torMN=module(tensorMod(M1,N)); }
1140      else { torMN  = gen(1); }
1141      torMN0 = std(torMN);
1142      L1[ii] = torMN;
1143      L2[ii] = torMN0;
1144      L3[ii] = matrix(kbase(torMN0));
1145      di=dim(torMN0);
1146      dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1147      if (di==0)
1148      {
1149        dbprint(p,"// vdim of Tor_"+string(k)+":       "
1150                  +string(vdim(torMN0)));
1151      }
1152      dbprint(p,"");
1153    }
1154
1155    if( size(#) )
1156    {  if( s>1 ) { L = L1,L2,L3; return(L); }
1157       else { L = torMN,torMN0,L3[1]; return(L); }
1158    }
1159    else
1160    {  if( s>1 ) { return(L1); }
1161       else { return(torMN); }
1162    }
1163  }
1164
1165  if( dim(M1)==-1 )              // coker(M)=0, all Tor's are 0
1166  { dbprint(p-1,"2nd module presents 0, hence Tor_k=0, for all k");
1167    for( ii=1; ii<=s; ii++ )
1168    { k=v[ii];
1169      torMN    = gen(1);
1170      torMN0   = std(torMN);
1171      L1[ii] = torMN;
1172      L2[ii] = torMN0;
1173      L3[ii] = matrix(kbase(torMN0));
1174      di=dim(torMN0);
1175      dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1176      if (di==0)
1177      {
1178        dbprint(p,"// vdim of Tor_"+string(k)+":       "
1179                  +string(vdim(torMN0)));
1180      }
1181      dbprint(p,"");
1182    }
1183  }
1184  else
1185  {
1186    if( size(M1) < size(M) ) { M=M1;}
1187    row = nrows(M);
1188//---------- resolution of N -------------------------------------------------
1189    resN = mres(N,max+1);
1190    for( ii=1; ii<=s; ii++ )
1191    { k=v[ii];
1192      if( k<0  )                             // Tor_k is 0 for negative k
1193      { dbprint(p-1,"// Tor_k=0 for k<0!");
1194        torMN  = gen(1);
1195        torMN0 = std(torMN);
1196        L1[ii] = torMN;
1197        L2[ii] = torMN0;
1198        L3[ii] = matrix(kbase(torMN0));
1199        di=dim(torMN0);
1200        dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1201        if (di==0)
1202        {
1203          dbprint(p,"// vdim of Tor_"+string(k)+":       "
1204                    +string(vdim(torMN0)));
1205        }
1206        dbprint(p,"");
1207      }
1208      else
1209      {
1210        N2 = resN[k+1];
1211        if( k==0 ) { torMN=module(tensorMod(M,N)); }
1212        else
1213        {
1214          N1 = resN[k];
1215          col = ncols(N1);
1216
1217//---------- computing homology ----------------------------------------------
1218          imag = tensor(unitmat(nrows(N1)),M);
1219          f = tensor(matrix(N1),unitmat(row));
1220          Im1 = tensor(unitmat(col),M);
1221          Im2 = tensor(matrix(N2),unitmat(row));
1222          ker = modulo(f,imag);
1223          Im  = Im2,Im1;
1224          torMN = modulo(ker,Im);
1225          dbprint(p-1,"// Computing Tor_"+string(k)+
1226                      " (help Tor; gives an explanation):",
1227          "// Let 0 <- coker(M) <- G0 <-M- G1 be the present. of coker(M),",
1228          "// and 0 <- coker(N) <- F0 <-N- F1 <- F2 <- ... a resolution of",
1229          "// coker(N), then Tensor(G0,F"+string(k)+")-->Tensor(G0,F"+
1230          string(k-1)+") is given by:",f,
1231          "// and Tensor(G0,F"+string(k+1)+") + Tensor(G1,F"+string(k)+
1232          ")-->Tensor(G0,F"+string(k)+") is given by:",Im,"");
1233        }
1234
1235        torMN0 = std(torMN);
1236        di=dim(torMN0);
1237        dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1238        if (di==0)
1239        {
1240          dbprint(p,"// vdim of Tor_"+string(k)+":       "
1241                    +string(vdim(torMN0)));
1242        }
1243        dbprint(p,"");
1244
1245//---------- more information -------------------------------------------------
1246        if( size(#)>0 )
1247        { if( vdim(torMN0) >= 0 )
1248          { kb = kbase(torMN0);
1249            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
1250            dbprint(p-1,"// columns of matrix are kbase of Tor_"+
1251                     string(k)+" in Tensor(G0,F"+string(k)+")",kb,"");
1252            L3[ii] = matrix(kb);
1253          }
1254          L2[ii] = torMN0;
1255        }
1256        L1[ii] = torMN;
1257      }
1258    }
1259  }
1260  if( size(#) )
1261  {  if( s>1 ) { L = L1,L2,L3; return(L); }
1262     else { L = torMN,torMN0,matrix(kb); return(L); }
1263  }
1264  else
1265  {  if( s>1 ) { return(L1); }
1266     else { return(torMN); }
1267  }
1268}
1269example
1270{"EXAMPLE:";   echo=2;
1271  int p      = printlevel;
1272  printlevel = 1;
1273  ring r     = 0,(x,y),dp;
1274  ideal i    = x2,y;
1275  ideal j    = x;
1276  list E     = Tor(0..2,i,j);    // Tor_k(r/i,r/j) for k=0,1,2 over r
1277
1278  qring R    = std(i);
1279  ideal j    = fetch(r,j);
1280  module M   = [x,0],[0,x];
1281  printlevel = 2;
1282  module E1  = Tor(1,M,j);       // Tor_1(R^2/M,R/j) over R=r/i
1283
1284  list l     = Tor(3,M,M,1);     // Tor_3(R^2/M,R^2/M) over R=r/i
1285  printlevel = p;
1286}
1287///////////////////////////////////////////////////////////////////////////////
1288proc fitting(module M, int n)
1289"USAGE:   fitting (M,n);  M module, n int
1290RETURN:  ideal, (standard basis of) n-th Fitting ideal of M'=coker(M).
1291EXAMPLE: example fitting; shows an example
1292"
1293{
1294  n=nrows(M)-n;
1295  if(n<=0){return(ideal(1));}
1296  if((n>nrows(M))||(n>ncols(M))){return(ideal(0));}
1297  return(std(minor(M,n)));
1298}
1299example
1300{"EXAMPLE:";   echo=2;
1301  ring R=0,x(0..4),dp;
1302  matrix M[2][4]=x(0),x(1),x(2),x(3),x(1),x(2),x(3),x(4);
1303  print(M);
1304  fitting(M,-1);
1305  fitting(M,0);
1306  fitting(M,1);
1307  fitting(M,2);
1308}
1309///////////////////////////////////////////////////////////////////////////////
1310proc isLocallyFree(matrix S, int r)
1311"USAGE:   isLocallyFree(M,r);  M module, r int
1312RETURN:  1 if M'=coker(M) is locally free of constant rank r;@*
1313         0 if this is not the case.
1314EXAMPLE: example isLocallyFree; shows an example.
1315"
1316{
1317   ideal F=fitting(S,r);
1318   ideal G=fitting(S,r-1);
1319   if((deg(F[1])==0)&&(size(G)==0)){return(1);}
1320   return(0);
1321}
1322example
1323{"EXAMPLE:";   echo=2;
1324  ring R=0,(x,y,z),dp;
1325  matrix M[2][3];     // the presentation matrix
1326  M=x-1,y-1,z,y-1,x-2,x;
1327  ideal I=fitting(M,0); // 0-th Fitting ideal of coker(M)
1328  qring Q=I;
1329  matrix M=fetch(R,M);
1330  isLocallyFree(M,1); // as R/I-module, coker(M) is locally free of rk 1
1331  isLocallyFree(M,0);
1332}
1333///////////////////////////////////////////////////////////////////////////////
1334proc flatteningStrat (module M)
1335"USAGE:   flatteningStrat(M);  M module
1336RETURN:  list of ideals.
1337         The list entries L[1],...,L[r] describe the flattening stratification
1338         of M'=coker(M): setting L[0]=0, L[r+1]=1, the flattening
1339         stratification is given by the open sets Spec(A/V(L[i-1])) \ V(L[i]),
1340         i=1,...,r+1  (A = basering).
1341NOTE:    for more information see the book 'A Singular Introduction to
1342         Commutative Algebra' (by Greuel/Pfister, Springer 2002).
1343EXAMPLE: example flatteningStrat; shows an example
1344"
1345{
1346   list l;
1347   int v,w;
1348   ideal F;
1349   while(1)
1350   {
1351      F=interred(fitting(M,w));
1352      if(F[1]==1){return(l);}
1353      if(size(F)!=0){v++;l[v]=F;}
1354      w++;
1355   }
1356   return(l);
1357}
1358example
1359{"EXAMPLE:";   echo=2;
1360  ring A = 0,x(0..4),dp;
1361  // presentation matrix:
1362  matrix M[2][4] = x(0),x(1),x(2),x(3),x(1),x(2),x(3),x(4);
1363  list L = flatteningStrat(M);
1364  L;
1365}
1366///////////////////////////////////////////////////////////////////////////////
1367proc isFlat(module M)
1368"USAGE:   isFlat(M);  M module
1369RETURN:  1 if M'=coker(M) is flat;@*
1370         0 if this is not the case.
1371EXAMPLE: example isFlat; shows an example.
1372"
1373{
1374   if (size(ideal(M))==0) {return(1);}
1375   int w;
1376   ideal F=fitting(M,0);
1377   while(size(F)==0)
1378   {
1379      w++;
1380      F=fitting(M,w);
1381   }
1382   if (deg(std(F)[1])==0) {return(1);}
1383   return(0);
1384}
1385example
1386{"EXAMPLE:";   echo=2;
1387  ring A = 0,(x,y),dp;
1388  matrix M[3][3] = x-1,y,x,x,x+1,y,x2,xy+x+1,x2+y;
1389  print(M);
1390  isFlat(M);             // coker(M) is not flat over A=Q[x,y]
1391
1392  qring B = std(x2+x-y);   // the ring B=Q[x,y]/<x2+x-y>
1393  matrix M = fetch(A,M);
1394  isFlat(M);             // coker(M) is flat over B
1395
1396  setring A;
1397  qring C = std(x2+x+y);   // the ring C=Q[x,y]/<x2+x+y>
1398  matrix M = fetch(A,M);
1399  isFlat(M);             // coker(M) is not flat over C
1400}
1401///////////////////////////////////////////////////////////////////////////////
1402proc flatLocus(module M)
1403"USAGE:   flatLocus(M);  M module
1404RETURN:  ideal I, s.th. complement of V(I) is flat locus of coker(M).
1405NOTE:    computation is based on Fitting ideals;@*
1406         output is not radical (in general)
1407EXAMPLE: example flatLocus; shows an example
1408"
1409{
1410   if (size(ideal(M))==0) {return(ideal(1));}
1411   int v,w;
1412   ideal F=fitting(M,0);
1413   while(size(F)==0)
1414   {
1415      w++;
1416      F=fitting(M,w);
1417   }
1418   if(typeof(basering)=="qring")
1419   {
1420      for(v=w+1;v<=nrows(M);v++)
1421      {
1422         F=F+intersect(fitting(M,v),quotient(ideal(0),fitting(M,v-1)));
1423      }
1424   }
1425   return(interred(F));
1426}
1427example
1428{"EXAMPLE:";   echo=2;
1429  ring R=0,(x,y,z),dp;
1430  matrix M[2][3]=x,y,z,0,x3,z3;
1431  ideal I=flatLocus(M); // coker(M) is flat outside V(x,yz)
1432  I;                    // computed ideal not radical
1433  ideal J=radical(I);
1434  J;
1435
1436  qring r=std(J);
1437  matrix M=fetch(r,M);
1438  flatLocus(M);         // coker(M) is flat over Spec(Q[x,y,z]/<x,yz>)
1439
1440  isFlat(M);            // flatness test
1441}
1442///////////////////////////////////////////////////////////////////////////////
1443proc isReg(ideal I, module N)
1444"USAGE:   isReg(I,M);  I ideal, M module
1445RETURN:  1 if given (ordered) list of generators for I is coker(M)-sequence;@*
1446         0 if this is not the case.
1447EXAMPLE: example isReg; shows an example.
1448"
1449{
1450  int n=nrows(N);
1451  int i;
1452  while(i<ncols(I))
1453  {
1454     i++;
1455     N=std(N);
1456     if(size(reduce(quotient(N,I[i]),N))!=0){return(0);}
1457     N=N+I[i]*freemodule(n);
1458  }
1459  if (size(reduce(freemodule(n),std(N)))==0){return(0);}
1460  return(1);
1461}
1462example
1463{"EXAMPLE:";   echo=2;
1464  ring R = 0,(x,y,z),dp;
1465  ideal I = x*(y-1),y,z*(y-1);
1466  isReg(I,0);             // given list of generators is Q[x,y,z]-sequence
1467
1468  I = x*(y-1),z*(y-1),y;  // change sorting of generators
1469  isReg(I,0);
1470
1471  ring r = 0,(x,y,z),ds;  // local ring
1472  ideal I=fetch(R,I);
1473  isReg(I,0);             // result independent of sorting of generators
1474}
1475
1476///////////////////////////////////////////////////////////////////////////////
1477// the following static procedures are used by KoszulHomology:
1478//  * binom_int  (binomial coeff. as integer, or -1 if too large)
1479//  * basisNumber
1480//  * basisElement
1481//  * KoszulMap
1482// for details, see 'A Singular Introduction to Commutative Algebra' (by
1483// Greuel/Pfister, Springer 2002), Chapter 7
1484
1485static proc binom_int(int n, int p)
1486{
1487  string s = binomial(n,p);
1488  if (size(s) > 9) { return(-1); } // result too large for integer
1489  execute("int a ="+s);
1490  return(a);
1491}
1492
1493static proc basisNumber(int n,intvec v)
1494{
1495  int p=size(v);
1496  if(p==1){return(v[1]);}
1497  int j=n-1;
1498  int b;
1499  while(j>=n-v[1]+1)
1500  {
1501    b=b+binom_int(j,p-1);
1502    j--;
1503  }
1504  intvec w=v-v[1];
1505  w=w[2..size(w)];
1506  b=b+basisNumber(n-v[1],w);
1507  return(b);
1508}
1509
1510static proc basisElement(int n,int p,int N)
1511{
1512  if(p==1){return(N);}
1513  int s,R;
1514  while(R<N)
1515  {
1516     s++;
1517     R=R+binom_int(n-s,p-1);
1518  }
1519  R=N-R+binom_int(n-s,p-1);
1520  intvec v=basisElement(n-s,p-1,R);
1521  intvec w=s,v+s;
1522  return(w);
1523}
1524
1525proc KoszulMap(ideal x,int p)
1526{
1527  int n=size(x);
1528  int a=binom_int(n,p-1);
1529  int b=binom_int(n,p);
1530
1531  matrix M[a][b];
1532
1533  if(p==1){M=x;return(M);}
1534  int j,k;
1535  intvec v,w;
1536  for(j=1;j<=b;j++)
1537  {
1538    v=basisElement(n,p,j);
1539    w=v[2..p];
1540    M[basisNumber(n,w),j]=x[v[1]];
1541    for(k=2;k<p;k++)
1542    {
1543      w=v[1..k-1],v[k+1..p];
1544      M[basisNumber(n,w),j]=(-1)^(k-1)*x[v[k]];
1545    }
1546    w=v[1..p-1];
1547    M[basisNumber(n,w),j]=(-1)^(p-1)*x[v[p]];
1548  }
1549  return(M);
1550}
1551///////////////////////////////////////////////////////////////////////////////
1552
1553proc KoszulHomology(ideal x, module M, int p)
1554"USAGE:   KoszulHomology(I,M,p);  I ideal, M module, p int
1555COMPUTE: A presentation of the p-th Koszul homology module H_p(f_1,...,f_k;M'),
1556         where M'=coker(M) and f_1,...,f_k are the given (ordered list
1557         of non-zero generators of the) ideal I.
1558         The computed presentation is minimized via prune.
1559         In particular, if H_p(f_1,...,f_k;M')=0 then the return value is 0.
1560RETURN:  module H, s.th. coker(H) = H_p(f_1,...,f_k;M').
1561NOTE:    size of input ideal has to be <= 20.
1562EXAMPLE: example KoszulHomology; shows an example.
1563{
1564  x=simplify(x,2);
1565  int n     = size(x);
1566  if (n==0)
1567  {
1568    ERROR("// KoszulHomology only for non-zero ideals");
1569  }
1570  if (n>20)
1571  {
1572    ERROR("// too many generators in input ideal");
1573  }
1574  if (p>n)
1575  {
1576    module hom=0;
1577    return(hom);
1578  }
1579
1580  int a     = binom_int(n,p-1);  // n over p-1 independent of char(basering)
1581  int b     = binom_int(n,p);
1582
1583  matrix N  = matrix(M);
1584  module ker= freemodule(nrows(N));
1585  if(p!=0)
1586  {
1587    module im= tensor(unitmat(a),N);
1588    module f = tensor(KoszulMap(x,p),unitmat(nrows(N)));
1589    ker      = modulo(f,im);
1590  }
1591  module im1 = tensor(unitmat(b),N);
1592  module im2 = tensor(KoszulMap(x,p+1),unitmat(nrows(N)));
1593  module hom = modulo(ker,im1+im2);
1594  hom        = prune(hom);
1595  return(hom);
1596}
1597example
1598{"EXAMPLE:";   echo=2;
1599  ring R=0,x(1..3),dp;
1600  ideal x=maxideal(1);
1601  module M=0;
1602  KoszulHomology(x,M,0);  // H_0(x,R), x=(x_1,x_2,x_3)
1603
1604  KoszulHomology(x,M,1);  // H_1(x,R), x=(x_1,x_2,x_3)
1605
1606  qring S=std(x(1)*x(2));
1607  module M=0;
1608  ideal x=maxideal(1);
1609  KoszulHomology(x,M,1);
1610
1611  KoszulHomology(x,M,2);
1612}
1613///////////////////////////////////////////////////////////////////////////////
1614proc depth(module M,ideal #)
1615"USAGE:   depth(M,[I]);  M module, I ideal
1616RETURN:  int,
1617         - if called with 1 argument: the depth of M'=coker(M) w.r.t. the
1618         maxideal in the basering (which is then assumed to be local)@*
1619         - if called with 2 arguments: the depth of M'=coker(M) w.r.t. the
1620         ideal I.
1621NOTE:    procedure makes use of KoszulHomology.
1622EXAMPLE: example depth; shows an example.
1623"
1624{
1625  ideal m=maxideal(1);
1626  int n=size(m);
1627  int i;
1628
1629  if (size(#)==0)
1630  {
1631    // depth(M') over local basering
1632    while(i<n)
1633    {
1634      i++;
1635      if(size(KoszulHomology(m,M,i))==0){return(n-i+1);}
1636    }
1637    return(0);
1638  }
1639
1640  ideal I=simplify(#,2);
1641  while(i<size(I))
1642  {
1643    i++;
1644    if(size(KoszulHomology(I,M,i))==0){return(size(I)-i+1);}
1645  }
1646  return(0);
1647}
1648example
1649{"EXAMPLE:";   echo=2;
1650  ring R=0,(x,y,z),dp;
1651  ideal I=x2,xy,yz;
1652  module M=0;
1653  depth(M,I);   // depth(<x2,xy,yz>,Q[x,y,z])
1654  ring r=0,(x,y,z),ds;  // local ring
1655  matrix M[2][2]=x,xy,1+yz,0;
1656  print(M);
1657  depth(M);     // depth(maxideal,coker(M))
1658  ideal I=x;
1659  depth(M,I);   // depth(<x>,coker(M))
1660  I=x+z;
1661  depth(M,I);   // depth(<x+z>,coker(M))
1662}
1663
1664///////////////////////////////////////////////////////////////////////////////
1665proc isCM(module M)
1666"USAGE:   isCM(M);  M module
1667RETURN:  1 if M'=coker(M) is Cohen-Macaulay;@*
1668         0 if this is not the case.
1669ASSUME:  basering is local.
1670EXAMPLE: example isCM; shows an example
1671"
1672{
1673  // test if basering is local:
1674  ideal m=maxideal(1);
1675  int i;
1676  poly f=1;
1677  for (i=1; i<=size(m); i++)
1678  {
1679    f=f+m[i];
1680  }
1681  if (ord(f)>0)
1682  {
1683    print("// basering must be local -- result has no meaning");
1684    return(0);
1685  }
1686
1687  return(depth(M)==dim(std(Ann(M))));
1688}
1689example
1690{"EXAMPLE:";   echo=2;
1691  ring R=0,(x,y,z),ds;  // local ring R = Q[x,y,z]_<x,y,z>
1692  module M=xz,yz,z2;
1693  isCM(M);             // test if R/<xz,yz,z2> is Cohen-Macaulay
1694
1695  M=x2+y2,z7;          // test if R/<x2+y2,z7> is Cohen-Macaulay
1696  isCM(M);
1697}
1698
Note: See TracBrowser for help on using the repository browser.