source: git/Singular/LIB/homolog.lib @ 6fe9a5

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