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

fieker-DuValspielwiese
Last change on this file since edaa83 was edaa83, checked in by Hans Schoenemann <hannes@…>, 9 years ago
unified ring handling: typeof(qring) is now ring
  • Property mode set to 100644
File size: 60.4 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="version homolog.lib 4.0.0.0 Jun_2013 "; // $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,module N,module P,int p,int 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 chosen 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      M=std(M);
537    }
538    list resl = sres(M,max+1);
539    if (withWeight) {
540      // ****
541      // **** at this writing, sres does not return weights, we have to
542      // **** go through the resolution to compute them
543      // ****
544      attrib(resl,"isHomog",weightM);  // weightM = weights of M
545      G=resl[1];
546      attrib(G,"isHomog",weightM);
547      resl[1]=G;
548      weightR=weightM;
549
550      for (j=2; j<=size(resl); j++) {
551        if (size(G)!=0) {
552          ww=0;
553          for (k=1; k<=ncols(G); k++) {
554            if (size(G[k])==0) { ww[k]=0; }
555            else {
556              leadCol = leadmonom(G[k]);
557              leadComp = nrows(leadCol);
558              ww[k] = deg(leadCol)+weightR[leadComp];
559            }
560          }
561          G=resl[j];
562          attrib(G,"isHomog",ww);
563          resl[j]=G;
564          weightR=ww;
565        }
566      }
567    }
568  }
569  else { list resl = mres(M,max+1);
570    if ((withWeight) and (size(M)==0)) {
571      // ***** At this writing: special treatment for zero module needed
572      G=resl[1];
573      attrib(G,"isHomog",weightM);
574      resl[1]=G;
575    }
576  }
577  for( ii=1; ii<=s; ii++ )
578  {
579//-----------------  apply Hom(_,R) at k-th place -----------------------------
580    k=v[ii];
581    dbprint(p-1,"// Computing Ext^"+string(k)+":");
582    if( k<0 )                   // Ext^k=0 for negative k
583    {
584      dbprint(p-1,"// Ext^i=0 for i<0!");
585      ret    = gen(1);
586      ret0   = std(ret);
587      if (withWeight==1) {
588        attrib(ret,"isHomog",intvec(0));
589        attrib(ret0,"isHomog",intvec(0));
590      }
591      L1[ii] = ret;
592      L2[ii] = ret0;
593      L3[ii] = matrix(kbase(ret0));
594      di=dim(ret0);
595      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
596      if (di==0)
597      {
598        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(ret0)));
599      }
600      dbprint(p,"");
601    }
602    else
603    {
604      m2t = resl[k+1];
605      m2 = transpose(m2t);
606      if ((typeof(attrib(m2t,"isHomog"))!="string" ) && (withWeight))
607      {
608        // ------------- compute weights for dual -----------------------------
609        weightR=attrib(m2t,"isHomog");
610        // --------------------------------------------------------------------
611        // *** set correct weights (at this writing, shift in resolution
612        // *** is not considered when defining the weights for the
613        // *** modules in the resolution):
614        A1=attrib(M,"isHomog");
615        A2=attrib(resl[1],"isHomog");
616        shift=A1[1]-A2[1];
617        for (j=1; j<=size(weightR); j++) { weightR[j]=weightR[j]+shift; }
618        attrib(m2t,"isHomog",weightR);
619        // --------------------------------------------------------------------
620        ww=0;
621        for (j=1; j<=nrows(m2); j++) {
622          if (size(m2t[j])==0) { ww[j]=0; }
623          else {
624            leadCol = leadmonom(m2t[j]);
625            leadComp = nrows(leadCol);
626            ww[j] = deg(leadCol)+weightR[leadComp];
627          }
628        }
629        attrib(m2,"isHomog",-ww);  // dualize --> negative weights
630        // --------------------------------------------------------------------
631        // *** the following should be replaced by the syz command,
632        // *** but syz forgets weights.....
633        resm2 = nres(m2,2);
634        ker = resm2[2];
635        if ((size(ker)>0) and (size(m2)>0)) {
636          // ------------------------------------------------------------------
637          // *** set correct weights (at this writing, shift in resolution
638          // *** is not considered when defining the weights for the
639          // *** modules in the resolution):
640          A1=attrib(resm2,"isHomog");
641          A2=attrib(resm2[1],"isHomog");
642          A3=attrib(ker,"isHomog");
643          shift=A1[1]-A2[1];
644          for (j=1; j<=size(A3); j++) { A3[j]=A3[j]+shift; }
645          // *** set correct weights where weights are undetermined due to
646          // *** zero columns in m2  (read weights from m2t)
647          for (j=1; j<=ncols(m2); j++) {
648            if (size(m2[j])==0) { A3[j]=-weightR[j]; }
649          }
650          attrib(ker,"isHomog",A3);
651          // ------------------------------------------------------------------
652        }
653      }
654      else {
655        ker = syz(m2);
656      }
657
658      if( k==0 ) { matrix MMM1[ncols(m2)][1];
659        m1=MMM1;
660      }
661      else { // k>0
662        m1t = resl[k];
663        m1 = transpose(resl[k]);
664        if ((typeof(attrib(m1t,"isHomog"))!="string" ) && (withWeight)) {
665        // ------------- compute weights for dual -----------------------------
666          weightR=attrib(resl[k],"isHomog");
667          // ------------------------------------------------------------------
668          // *** set correct weights (at this writing, shift in resolution
669          // *** is not considered when defining the weights for the
670          // *** modules in the resolution):
671          A1=attrib(M,"isHomog");
672          A2=attrib(resl[1],"isHomog");
673          shift=A1[1]-A2[1];
674          for (j=1; j<=size(weightR); j++) { weightR[j]=weightR[j]+shift; }
675          attrib(m1t,"isHomog",weightR);
676          // ------------------------------------------------------------------
677          ww=0;
678          for (j=1; j<=nrows(m1); j++) {
679            if (size(m1t[j])==0) { ww[j]=0; }
680            else {
681              leadCol = leadmonom(m1t[j]);
682              leadComp = nrows(leadCol);
683              ww[j] = deg(leadCol)+weightR[leadComp];
684            }
685          }
686          attrib(m1,"isHomog",-ww);  // dualize --> negative weights
687        }
688      }
689      //----------------- presentation of ker(m2)/im(m1) ---------------------
690      if ((k==0) and (size(M)==0)) {
691        ret = M;
692        if (withWeight) { attrib(ret,"isHomog",-weightM); }
693      }
694      else {
695        ret = modulo(ker,m1);
696      }
697      dbprint(p-1,
698         "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of M,",
699         "// then F"+string(k)+"*-->F"+string(k+1)+"* is given by:",m2,
700         "// and F"+string(k-1)+"*-->F"+string(k)+"* is given by:",m1,"");
701      ret0 = std(ret);
702
703      di=dim(ret0);
704      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
705      if (di==0)
706      {
707        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(ret0)));
708      }
709      dbprint(p,"");
710      if( t2 )
711      {
712         if( vdim(ret0)>=0 )
713         {
714            kb = kbase(ret0);
715            if ( size(ker)!=0 ) { kb = matrix(ker)*kb; }
716            dbprint(p-1,
717              "// columns of matrix are kbase of Ext^"+string(k)+" in F"
718              +string(k)+"*:",kb,"");
719            L3[ii] = kb;
720         }
721         L2[ii] = ret0;
722      }
723      L1[ii] = ret;
724    }
725  }
726  if( t2 )
727  {
728     if( s>1 ) { L = L1,L2,L3; return(L); }
729     else { L = ret,ret0,kb; return(L); }
730  }
731  else
732  {
733     if( s>1 ) { return(L1); }
734     else { return(ret); }
735  }
736}
737example
738{"EXAMPLE:";     echo=2;
739  int p      = printlevel;
740  printlevel = 1;
741  ring r     = 0,(x,y,z),dp;
742  ideal i    = x2y,y2z,z3x;
743  module E   = Ext_R(1,i);    //computes Ext^1(r/i,r)
744  is_zero(E);
745
746  qring R    = std(x2+yz);
747  intvec v   = 0,2;
748  printlevel = 2;             //shows what is going on
749  ideal i    = x,y,z;         //computes Ext^i(r/(x,y,z),r/(x2+yz)), i=0,2
750  list L     = Ext_R(v,i,1);  //over the qring R=r/(x2+yz), std and kbase
751  printlevel = p;
752}
753///////////////////////////////////////////////////////////////////////////////
754
755proc Ext (intvec v, module M, module N, list #)
756"USAGE:   Ext(v,M,N[,any]);  v int resp. intvec, M,N modules
757COMPUTE: A presentation of Ext^k(M',N'); for k=v[1],v[2],... where
758         M'=coker(M) and N'=coker(N). Let
759@example
760       0 <-- M' <-- F0 <-M-- F1 <-- F2 <--... ,
761       0 <-- N' <-- G0 <--N- G1
762@end example
763         be a free resolution of M', resp. a presentation of N'. Consider
764         the commutative diagram
765@example
766           0                  0                  0
767           |^                 |^                 |^
768   --> Hom(Fk-1,N') -Ak-> Hom(Fk,N') -Ak+1-> Hom(Fk+1,N')
769           |^                 |^                 |^
770   --> Hom(Fk-1,G0) -Ak-> Hom(Fk,G0) -Ak+1-> Hom(Fk+1,G0)
771                              |^                 |^
772                              |C                 |B
773                          Hom(Fk,G1) ------> Hom(Fk+1,G1)
774
775      (Ak,Ak+1 induced by M and B,C induced by N).
776@end example
777         Let K=modulo(Ak+1,B), J=module(Ak)+module(C) and Ext=modulo(K,J),
778         then we have exact sequences
779@example
780    R^p --K-> Hom(Fk,G0) --Ak+1-> Hom(Fk+1,G0)/im(B),
781
782    R^q -Ext-> R^p --K-> Hom(Fk,G0)/(im(Ak)+im(C)).
783@end example
784         Hence, Ext presents Ext^k(M',N').
785RETURN:  - module Ext, a presentation of Ext^k(M',N') if v is of type int@*
786         - a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.@*
787         - In case of a third argument of any type return a list l:
788@format
789             l[1] = module Ext/list of Ext^k
790             l[2] = SB of Ext/list of SB of Ext^k
791             l[3] = matrix/list of matrices, each representing a kbase of Ext^k
792                       (if finite dimensional)
793@end format
794DISPLAY: printlevel >=0: dimension, vdim of Ext^k for each k  (default).
795@*       printlevel >=1: matrices Ak, Ak+1 and kbase of Ext^k in Hom(Fk,G0)
796         (if finite dimensional)
797NOTE:    In order to compute Ext^k(M,N) use the command Ext(k,syz(M),syz(N));
798         or: list P=mres(M,2); list Q=mres(N,2); Ext(k,P[2],Q[2]);
799EXAMPLE: example Ext; shows an example
800"
801{
802//---------- initialisation ---------------------------------------------------
803  int k,max,ii,l,row,col,di;
804  module A,B,C,D,M1,M2,N1,ker,imag,extMN,extMN0;
805  matrix kb;
806  list L1,L2,L3,L,resM,K;
807  ideal  test1;
808  intmat Be;
809  int s = size(v);
810  max = Max(v);                   // the maximum integer occurring in intvec v
811  int p = printlevel-voice+3;     // p=printlevel+1 (default: p=1)
812//---------- test: coker(N)=basering, coker(N)=0 ? ----------------------------
813  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
814  N1 = std(N);
815  if( size(N1)==0 )      //coker(N)=basering, in this case proc Ext_R is faster
816  { printlevel=printlevel+1;
817    if( size(#)==0 )
818    { def E = Ext_R(v,M);
819      printlevel=printlevel-1;
820      return(E);
821    }
822    else
823    { def E = Ext_R(v,M,#[1]);
824       printlevel=printlevel-1;
825       return(E);
826     }
827  }
828  if( dim(N1)==-1 )                          //coker(N)=0, all Ext-groups are 0
829  { dbprint(p-1,"2nd module presents 0, hence Ext^k=0, for all k");
830    for( ii=1; ii<=s; ii++ )
831    { k=v[ii];
832      extMN    = gen(1);
833      extMN0   = std(extMN);
834      L1[ii] = extMN;
835      L2[ii] = extMN0;
836      L3[ii] = matrix(kbase(extMN0));
837      di=dim(extMN0);
838      dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
839      if (di==0)
840      {
841        dbprint(p,"// vdim of Ext^"+string(k)+":       "+string(vdim(extMN0)));
842      }
843      dbprint(p,"");
844    }
845  }
846  else
847  {
848    if( size(N1) < size(N) ) { N=N1;}
849    row = nrows(N);
850//---------- resolution of M -------------------------------------------------
851    resM = mres(M,max+1);
852    for( ii=1; ii<=s; ii++ )
853    { k=v[ii];
854      if( k<0  )                                   // Ext^k is 0 for negative k
855      { dbprint(p-1,"// Ext^k=0 for k<0!");
856        extMN    = gen(1);
857        extMN0   = std(extMN);
858        L1[ii] = extMN;
859        L2[ii] = extMN0;
860        L3[ii] = matrix(kbase(extMN0));
861        di=dim(extMN0);
862        dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
863        if (di==0)
864        {
865          dbprint(p,"// vdim of Ext^"+string(k)+":       "
866                  +string(vdim(extMN0)));
867        }
868        dbprint(p,"");
869      }
870      else
871      { M2 = resM[k+1];
872        if( k==0 ) { M1=0; }
873        else { M1 = resM[k]; }
874        col = nrows(M2);
875        D = kohom(N,col);
876//---------- computing homology ----------------------------------------------
877        imag  = kontrahom(M1,row);
878        A     = kontrahom(M2,row);
879        B     = kohom(N,ncols(M2));
880        ker   = modulo(A,B);
881        imag  = imag,D;
882        extMN = modulo(ker,imag);
883        dbprint(p-1,"// Computing Ext^"+string(k)+
884         " (help Ext; gives an explanation):",
885      "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of coker(M),",
886      "// and 0<--coker(N)<--G0<--G1 a presentation of coker(N),",
887      "// then Hom(F"+string(k)+",G0)-->Hom(F"+string(k+1)+
888      ",G0) is given by:",A,
889      "// and Hom(F"+string(k-1)+",G0) + Hom(F"+string(k)+",G1)-->Hom(F"
890                                      +string(k)+",G0) is given by:",imag,"");
891        extMN0 = std(extMN);
892        di=dim(extMN0);
893        dbprint(p,"// dimension of Ext^"+string(k)+":  "+string(di));
894        if (di==0)
895        {
896          dbprint(p,"// vdim of Ext^"+string(k)+":       "
897                  +string(vdim(extMN0)));
898        }
899        dbprint(p,"");
900
901//---------- more information -------------------------------------------------
902        if( size(#)>0 )
903        { if( vdim(extMN0) >= 0 )
904          { kb = kbase(extMN0);
905            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
906            dbprint(p-1,"// columns of matrix are kbase of Ext^"+
907                     string(k)+" in Hom(F"+string(k)+",G0)",kb,"");
908            if( p>0 )
909            { for (l=1;l<=ncols(kb);l=l+1)
910              {
911                "// element",l,"of kbase of Ext^"+string(k)+" in Hom(F"+string(k)+",G0)";
912                "// as matrix: F"+string(k)+"-->G0";
913                print(matrix(ideal(kb[l]),row,col));
914              }
915              "";
916            }
917            L3[ii] = matrix(kb);
918          }
919          L2[ii] = extMN0;
920        }
921        L1[ii] = extMN;
922      }
923    }
924  }
925  if( size(#) )
926  {  if( s>1 ) { L = L1,L2,L3; return(L); }
927     else { L = extMN,extMN0,matrix(kb); return(L); }
928  }
929  else
930  {  if( s>1 ) { return(L1); }
931     else { return(extMN); }
932  }
933}
934example
935{"EXAMPLE:";   echo=2;
936  int p      = printlevel;
937  printlevel = 1;
938  ring r     = 0,(x,y),dp;
939  ideal i    = x2-y3;
940  ideal j    = x2-y5;
941  list E     = Ext(0..2,i,j);    // Ext^k(r/i,r/j) for k=0,1,2 over r
942  qring R    = std(i);
943  ideal j    = fetch(r,j);
944  module M   = [-x,y],[-y2,x];
945  printlevel = 2;
946  module E1  = Ext(1,M,j);       // Ext^1(R^2/M,R/j) over R=r/i
947  list l     = Ext(4,M,M,1);     // Ext^4(R^2/M,R^2/M) over R=r/i
948  printlevel = p;
949}
950///////////////////////////////////////////////////////////////////////////////
951
952proc Hom (module M, module N, list #)
953"USAGE:   Hom(M,N,[any]);  M,N=modules
954COMPUTE: A presentation of Hom(M',N'), M'=coker(M), N'=coker(N) as follows:
955         let
956@example
957   F1 --M-> F0 -->M' --> 0,    G1 --N-> G0 --> N' --> 0
958@end example
959         be presentations of M' and N'. Consider
960@example
961                                  0               0
962                                  |^              |^
963       0 --> Hom(M',N') ----> Hom(F0,N') ----> Hom(F1,N')
964                                  |^              |^
965  (A:  induced by M)          Hom(F0,G0) --A-> Hom(F1,G0)
966                                  |^              |^
967  (B,C:induced by N)              |C              |B
968                              Hom(F0,G1) ----> Hom(F1,G1)
969
970@end example
971         Let D=modulo(A,B) and Hom=modulo(D,C), then we have exact sequences
972@example
973   R^p  --D-> Hom(F0,G0) --A-> Hom(F1,G0)/im(B),
974
975 R^q -Hom-> R^p --D-> Hom(F0,G0)/im(C) --A-> Hom(F1,G0)/im(B).
976@end example
977         Hence Hom presents Hom(M',N')
978RETURN:  module Hom, a presentation of Hom(M',N'), resp., in case of
979         3 arguments, a list l (of size <=3):
980@format
981           - l[1] = Hom
982           - l[2] = SB of Hom
983           - l[3] = kbase of coker(Hom) (if finite dimensional, not 0),
984                    represented by elements in Hom(F0,G0) via mapping D
985@end format
986DISPLAY: printlevel >=0: (affine) dimension of Hom  (default)
987@*       printlevel >=1: D and C and kbase of coker(Hom) in Hom(F0,G0)
988@*       printlevel >=2: elements of kbase of coker(Hom) as matrix :F0-->G0
989NOTE:    DISPLAY is as described only for a direct call of 'Hom'. Calling 'Hom'
990         from another proc has the same effect as decreasing printlevel by 1.
991EXAMPLE: example Hom;  shows examples
992"
993{
994//---------- initialisation ---------------------------------------------------
995  int l,p,di;
996  matrix kb;
997  module A,B,C,D,homMN,homMN0;
998  list L;
999//---------- computation of Hom -----------------------------------------------
1000  B = kohom(N,ncols(M));
1001  A = kontrahom(M,nrows(N));
1002  C = kohom(N,nrows(M));
1003  D = modulo(A,B);
1004  homMN = modulo(D,C);
1005  homMN0= std(homMN);
1006  p = printlevel-voice+3;       // p=printlevel+1 (default: p=1)
1007  di= dim(homMN0);
1008  dbprint(p,"// dimension of Hom:  "+string(di));
1009  if (di==0)
1010  {
1011    dbprint(p,"// vdim of Hom:       "+string(vdim(homMN0)));
1012  }
1013  dbprint(p,"");
1014  dbprint(p-1,
1015   "// given  F1 --M-> F0 -->M'--> 0 and  G1 --N-> G0 -->N'--> 0,",
1016   "// show D = ker( Hom(F0,G0) --> Hom(F1,G0)/im(Hom(F1,G1)->Hom(F1,G0)) )",D,
1017   "// show C = im ( Hom(F0,G1) --> Hom(F0,G0) )",C,"");
1018//---------- extra output if size(#)>0 ----------------------------------------
1019  if( size(#)>0 )
1020  {
1021     if( vdim(homMN0)>0 )
1022     {
1023        kb = kbase(homMN0);
1024        kb = matrix(D)*kb;
1025        if( p>2 )
1026        {
1027          for (l=1;l<=ncols(kb);l=l+1)
1028          {
1029            "// element",l,"of kbase of Hom in Hom(F0,G0) as matrix: F0-->G0:";
1030             print(matrix(ideal(kb[l]),nrows(N),nrows(M)));
1031          }
1032        }
1033        else
1034        {
1035          dbprint(p-1,"// columns of matrix are kbase of Hom in Hom(F0,G0)",
1036                      kb);
1037        }
1038        L=homMN,homMN0,kb;
1039        return(L);
1040     }
1041     L=homMN,homMN0;
1042     return(L);
1043  }
1044  return(homMN);
1045}
1046example
1047{"EXAMPLE:";  echo = 2;
1048  int p     = printlevel;
1049  printlevel= 1;   //in 'example proc' printlevel has to be increased by 1
1050  ring r    = 0,(x,y),dp;
1051  ideal i   = x2-y3,xy;
1052  qring q   = std(i);
1053  ideal i   = fetch(r,i);
1054  module M  = [-x,y],[-y2,x],[x3];
1055  module H  = Hom(M,i);
1056  print(H);
1057
1058  printlevel= 2;
1059  list L    = Hom(M,i,1);"";
1060
1061  printlevel=1;
1062  ring s    = 3,(x,y,z),(c,dp);
1063  ideal i   = jacob(ideal(x2+y5+z4));
1064  qring rq=std(i);
1065  matrix M[2][2]=xy,x3,5y,4z,x2;
1066  matrix N[3][2]=x2,x,y3,3xz,x2z,z;
1067  print(M);
1068  print(N);
1069  list l=Hom(M,N,1);
1070  printlevel = p;
1071}
1072///////////////////////////////////////////////////////////////////////////////
1073
1074proc homology (matrix A,matrix B,module M,module N,list #)
1075"USAGE:   homology(A,B,M,N);
1076COMPUTE: Let M and N be submodules of R^m and R^n presenting M'=R^m/M, N'=R^n/N
1077         (R=basering) and let A,B matrices inducing maps
1078@example
1079    R^k --A--> R^m --B--> R^n.
1080@end example
1081         Compute a presentation of the module
1082@example
1083    ker(B)/im(A) := ker(M'/im(A) --B--> N'/im(BM)+im(BA)).
1084@end example
1085         If B induces a map M'-->N' (i.e BM=0) and if im(A) is contained in
1086         ker(B) (that is, BA=0) then ker(B)/im(A) is the homology of the
1087         complex
1088@example
1089    R^k--A-->M'--B-->N'.
1090@end example
1091RETURN:  module H, a presentation of ker(B)/im(A).
1092NOTE:    homology returns a free module of rank m if ker(B)=im(A).
1093EXAMPLE: example homology; shows examples
1094"
1095{
1096  module ker,ima;
1097  ker = modulo(B,N);
1098  ima = A,M;
1099  return(modulo(ker,ima));
1100}
1101example
1102{"EXAMPLE";    echo=2;
1103  ring r;
1104  ideal id=maxideal(4);
1105  qring qr=std(id);
1106  module N=maxideal(3)*freemodule(2);
1107  module M=maxideal(2)*freemodule(2);
1108  module B=[2x,0],[x,y],[z2,y];
1109  module A=M;
1110  module H=homology(A,B,M,N);
1111  H=std(H);
1112  // dimension of homology:
1113  dim(H);
1114  // vector space dimension:
1115  vdim(H);
1116
1117  ring s=0,x,ds;
1118  qring qs=std(x4);
1119  module A=[x];
1120  module B=A;
1121  module M=[x3];
1122  module N=M;
1123  homology(A,B,M,N);
1124}
1125//////////////////////////////////////////////////////////////////////////////
1126
1127proc hom_kernel (matrix A,module M,module N)
1128"USAGE:   hom_kernel(A,M,N);
1129COMPUTE: Let M and N be submodules of R^m and R^n, presenting M'=R^m/M,
1130         N'=R^n/N (R=basering), and let A:R^m-->R^n be a matrix inducing a
1131         map A':M'-->N'. Then ker(A,M,N); computes a presentation K of
1132         ker(A') as in the commutative diagram:
1133@example
1134          ker(A') --->  M' --A'--> N'
1135             |^         |^         |^
1136             |          |          |
1137             R^r  ---> R^m --A--> R^n
1138             |^         |^         |^
1139             |K         |M         |N
1140             |          |          |
1141             R^s  ---> R^p -----> R^q
1142@end example
1143RETURN:  module K, a presentation of ker(A':coker(M)->coker(N)).
1144EXAMPLE: example hom_kernel; shows examples.
1145"
1146{
1147  module M1 = modulo(A,N);
1148  return(modulo(M1,M));
1149}
1150example
1151{"EXAMPLE";    echo=2;
1152  ring r;
1153  module N=[2x,x],[0,y];
1154  module M=maxideal(1)*freemodule(2);
1155  matrix A[2][3]=2x,0,x,y,z2,y;
1156  module K=hom_kernel(A,M,N);
1157  // dimension of kernel:
1158  dim(std(K));
1159  // vector space dimension of kernel:
1160  vdim(std(K));
1161  print(K);
1162}
1163///////////////////////////////////////////////////////////////////////////////
1164
1165proc kohom (matrix M, int j)
1166"USAGE:   kohom(A,k); A=matrix, k=integer
1167RETURN:  matrix Hom(R^k,A), i.e. let A be a matrix defining a map F1-->F2
1168         of free R-modules, then the matrix of Hom(R^k,F1)-->Hom(R^k,F2)
1169         is computed (R=basering).
1170EXAMPLE: example kohom; shows an example.
1171"
1172{
1173  if (j==1)
1174  { return(M);}
1175  if (j>1)
1176  { return(tensor(M,diag(1,j))); }
1177  else { return(0);}
1178}
1179example
1180{"EXAMPLE:";   echo=2;
1181  ring r;
1182  matrix n[2][3]=x,y,5,z,77,33;
1183  print(kohom(n,3));
1184}
1185///////////////////////////////////////////////////////////////////////////////
1186
1187proc kontrahom (matrix M, int j)
1188"USAGE:   kontrahom(A,k); A=matrix, k=integer
1189RETURN:  matrix Hom(A,R^k), i.e. let A be a matrix defining a map F1-->F2 of
1190         free  R-modules, then the matrix of Hom(F2,R^k)-->Hom(F1,R^k) is
1191         computed (R=basering).
1192EXAMPLE: example kontrahom; shows an example.
1193"
1194{
1195if (j==1)
1196  { return(transpose(M));}
1197  if (j>1)
1198  { return(transpose(tensor(diag(1,j),M)));}
1199  else { return(0);}
1200}
1201example
1202{"EXAMPLE:";  echo=2;
1203  ring r;
1204  matrix n[2][3]=x,y,5,z,77,33;
1205  print(kontrahom(n,3));
1206}
1207///////////////////////////////////////////////////////////////////////////////
1208///////////////////////////////////////////////////////////////////////////////
1209
1210proc tensorMod(module Phi, module Psi)
1211"USAGE:   tensorMod(M,N);  M,N modules
1212COMPUTE: presentation matrix A of the tensor product T of the modules
1213         M'=coker(M), N'=coker(N): if matrix(M) defines a map M: R^r-->R^s and
1214         matrix(N) defines a map N: R^p-->R^q, then A defines a presentation
1215@example
1216         R^(sp+rq) --A-> R^(sq)  --> T --> 0 .
1217@end example
1218RETURN:  matrix A satisfying coker(A) = tensorprod(coker(M),coker(N)) .
1219EXAMPLE: example tensorMod; shows an example.
1220"
1221{
1222   int s=nrows(Phi);
1223   int q=nrows(Psi);
1224   matrix A=tensor(unitmat(s),Psi);
1225   matrix B=tensor(Phi,unitmat(q));
1226   matrix R=concat(A,B);
1227   return(R);
1228}
1229example
1230{"EXAMPLE:";  echo=2;
1231  ring A=0,(x,y,z),dp;
1232  matrix M[3][3]=1,2,3,4,5,6,7,8,9;
1233  matrix N[2][2]=x,y,0,z;
1234  print(M);
1235  print(N);
1236  print(tensorMod(M,N));
1237}
1238///////////////////////////////////////////////////////////////////////////////
1239proc Tor(intvec v, module M, module N, list #)
1240"USAGE:   Tor(v,M,N[,any]);  v int resp. intvec, M,N modules
1241COMPUTE: a presentation of Tor_k(M',N'), for k=v[1],v[2],... , where
1242         M'=coker(M) and N'=coker(N): let
1243@example
1244       0 <-- M' <-- G0 <-M-- G1
1245       0 <-- N' <-- F0 <--N- F1 <-- F2 <--...
1246@end example
1247         be a presentation of M', resp. a free resolution of N', and consider
1248         the commutative diagram
1249@example
1250          0                    0                    0
1251          |^                   |^                   |^
1252  Tensor(M',Fk+1) -Ak+1-> Tensor(M',Fk) -Ak-> Tensor(M',Fk-1)
1253          |^                   |^                   |^
1254  Tensor(G0,Fk+1) -Ak+1-> Tensor(G0,Fk) -Ak-> Tensor(G0,Fk-1)
1255                               |^                   |^
1256                               |C                   |B
1257                          Tensor(G1,Fk) ----> Tensor(G1,Fk-1)
1258
1259       (Ak,Ak+1 induced by N and B,C induced by M).
1260@end example
1261         Let K=modulo(Ak,B), J=module(C)+module(Ak+1) and Tor=modulo(K,J),
1262         then we have exact sequences
1263@example
1264    R^p  --K-> Tensor(G0,Fk) --Ak-> Tensor(G0,Fk-1)/im(B),
1265
1266    R^q -Tor-> R^p --K-> Tensor(G0,Fk)/(im(C)+im(Ak+1)).
1267@end example
1268         Hence, Tor presents Tor_k(M',N').
1269RETURN:  - if v is of type int: module Tor, a presentation of Tor_k(M',N');@*
1270         - if v is of type intvec: a list of Tor_k(M',N') (k=v[1],v[2],...);@*
1271         - in case of a third argument of any type: list l with
1272@format
1273     l[1] = module Tor/list of Tor_k(M',N'),
1274     l[2] = SB of Tor/list of SB of Tor_k(M',N'),
1275     l[3] = matrix/list of matrices, each representing a kbase of Tor_k(M',N')
1276                (if finite dimensional), or 0.
1277@end format
1278DISPLAY: printlevel >=0: (affine) dimension of Tor_k for each k  (default).
1279@*       printlevel >=1: matrices Ak, Ak+1 and kbase of Tor_k in Tensor(G0,Fk)
1280         (if finite dimensional).
1281NOTE:    In order to compute Tor_k(M,N) use the command Tor(k,syz(M),syz(N));
1282         or: list P=mres(M,2); list Q=mres(N,2); Tor(k,P[2],Q[2]);
1283EXAMPLE: example Tor; shows an example
1284{
1285//---------- initialisation ---------------------------------------------------
1286  int k,max,ii,l,row,col,di;
1287  module A,B,C,D,N1,N2,M1,ker,imag,Im,Im1,Im2,f,torMN,torMN0;
1288  matrix kb;
1289  list L1,L2,L3,L,resN,K;
1290  ideal  test1;
1291  intmat Be;
1292  int s = size(v);
1293  max = Max(v);                  // maximum integer occurring in intvec v
1294  int p = printlevel-voice+3;    // p=printlevel+1 (default: p=1)
1295
1296//---------- test: coker(M)=basering, coker(M)=0 ? ----------------------------
1297  if( max<0 ) { dbprint(p,"// Tor_i=0 for i<0!"); return([1]); }
1298  M1 = std(M);
1299
1300  if( size(M1)==0 or size(N)==0 )  // coker(M)=basering ==> Tor_i=0 for i>0
1301  {
1302    dbprint(p-1,"// one of the modules M',N' is free, hence Tor_i=0 for i<>0");
1303    for( ii=1; ii<=s; ii++ )
1304    {
1305      k=v[ii];
1306      if (k==0) { torMN=module(tensorMod(M1,N)); }
1307      else { torMN  = gen(1); }
1308      torMN0 = std(torMN);
1309      L1[ii] = torMN;
1310      L2[ii] = torMN0;
1311      L3[ii] = matrix(kbase(torMN0));
1312      di=dim(torMN0);
1313      dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1314      if (di==0)
1315      {
1316        dbprint(p,"// vdim of Tor_"+string(k)+":       "
1317                  +string(vdim(torMN0)));
1318      }
1319      dbprint(p,"");
1320    }
1321
1322    if( size(#) )
1323    {  if( s>1 ) { L = L1,L2,L3; return(L); }
1324       else { L = torMN,torMN0,L3[1]; return(L); }
1325    }
1326    else
1327    {  if( s>1 ) { return(L1); }
1328       else { return(torMN); }
1329    }
1330  }
1331
1332  if( dim(M1)==-1 )              // coker(M)=0, all Tor's are 0
1333  { dbprint(p-1,"2nd module presents 0, hence Tor_k=0, for all k");
1334    for( ii=1; ii<=s; ii++ )
1335    { k=v[ii];
1336      torMN    = gen(1);
1337      torMN0   = std(torMN);
1338      L1[ii] = torMN;
1339      L2[ii] = torMN0;
1340      L3[ii] = matrix(kbase(torMN0));
1341      di=dim(torMN0);
1342      dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1343      if (di==0)
1344      {
1345        dbprint(p,"// vdim of Tor_"+string(k)+":       "
1346                  +string(vdim(torMN0)));
1347      }
1348      dbprint(p,"");
1349    }
1350  }
1351  else
1352  {
1353    if( size(M1) < size(M) ) { M=M1;}
1354    row = nrows(M);
1355//---------- resolution of N -------------------------------------------------
1356    resN = mres(N,max+1);
1357    for( ii=1; ii<=s; ii++ )
1358    { k=v[ii];
1359      if( k<0  )                             // Tor_k is 0 for negative k
1360      { dbprint(p-1,"// Tor_k=0 for k<0!");
1361        torMN  = gen(1);
1362        torMN0 = std(torMN);
1363        L1[ii] = torMN;
1364        L2[ii] = torMN0;
1365        L3[ii] = matrix(kbase(torMN0));
1366        di=dim(torMN0);
1367        dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1368        if (di==0)
1369        {
1370          dbprint(p,"// vdim of Tor_"+string(k)+":       "
1371                    +string(vdim(torMN0)));
1372        }
1373        dbprint(p,"");
1374      }
1375      else
1376      {
1377        N2 = resN[k+1];
1378        if( k==0 ) { torMN=module(tensorMod(M,N)); }
1379        else
1380        {
1381          N1 = resN[k];
1382          col = ncols(N1);
1383
1384//---------- computing homology ----------------------------------------------
1385          imag = tensor(unitmat(nrows(N1)),M);
1386          f = tensor(matrix(N1),unitmat(row));
1387          Im1 = tensor(unitmat(col),M);
1388          Im2 = tensor(matrix(N2),unitmat(row));
1389          ker = modulo(f,imag);
1390          Im  = Im2,Im1;
1391          torMN = modulo(ker,Im);
1392          dbprint(p-1,"// Computing Tor_"+string(k)+
1393                      " (help Tor; gives an explanation):",
1394          "// Let 0 <- coker(M) <- G0 <-M- G1 be the present. of coker(M),",
1395          "// and 0 <- coker(N) <- F0 <-N- F1 <- F2 <- ... a resolution of",
1396          "// coker(N), then Tensor(G0,F"+string(k)+")-->Tensor(G0,F"+
1397          string(k-1)+") is given by:",f,
1398          "// and Tensor(G0,F"+string(k+1)+") + Tensor(G1,F"+string(k)+
1399          ")-->Tensor(G0,F"+string(k)+") is given by:",Im,"");
1400        }
1401
1402        torMN0 = std(torMN);
1403        di=dim(torMN0);
1404        dbprint(p,"// dimension of Tor_"+string(k)+":  "+string(di));
1405        if (di==0)
1406        {
1407          dbprint(p,"// vdim of Tor_"+string(k)+":       "
1408                    +string(vdim(torMN0)));
1409        }
1410        dbprint(p,"");
1411
1412//---------- more information -------------------------------------------------
1413        if( size(#)>0 )
1414        { if( vdim(torMN0) >= 0 )
1415          { kb = kbase(torMN0);
1416            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
1417            dbprint(p-1,"// columns of matrix are kbase of Tor_"+
1418                     string(k)+" in Tensor(G0,F"+string(k)+")",kb,"");
1419            L3[ii] = matrix(kb);
1420          }
1421          L2[ii] = torMN0;
1422        }
1423        L1[ii] = torMN;
1424      }
1425    }
1426  }
1427  if( size(#) )
1428  {  if( s>1 ) { L = L1,L2,L3; return(L); }
1429     else { L = torMN,torMN0,matrix(kb); return(L); }
1430  }
1431  else
1432  {  if( s>1 ) { return(L1); }
1433     else { return(torMN); }
1434  }
1435}
1436example
1437{"EXAMPLE:";   echo=2;
1438  int p      = printlevel;
1439  printlevel = 1;
1440  ring r     = 0,(x,y),dp;
1441  ideal i    = x2,y;
1442  ideal j    = x;
1443  list E     = Tor(0..2,i,j);    // Tor_k(r/i,r/j) for k=0,1,2 over r
1444
1445  qring R    = std(i);
1446  ideal j    = fetch(r,j);
1447  module M   = [x,0],[0,x];
1448  printlevel = 2;
1449  module E1  = Tor(1,M,j);       // Tor_1(R^2/M,R/j) over R=r/i
1450
1451  list l     = Tor(3,M,M,1);     // Tor_3(R^2/M,R^2/M) over R=r/i
1452  printlevel = p;
1453}
1454///////////////////////////////////////////////////////////////////////////////
1455proc fitting(module M, int n)
1456"USAGE:   fitting (M,n);  M module, n int
1457RETURN:  ideal, (standard basis of) n-th Fitting ideal of M'=coker(M).
1458EXAMPLE: example fitting; shows an example
1459"
1460{
1461  n=nrows(M)-n;
1462  if(n<=0){return(ideal(1));}
1463  if((n>nrows(M))||(n>ncols(M))){return(ideal(0));}
1464  return(std(minor(M,n)));
1465}
1466example
1467{"EXAMPLE:";   echo=2;
1468  ring R=0,x(0..4),dp;
1469  matrix M[2][4]=x(0),x(1),x(2),x(3),x(1),x(2),x(3),x(4);
1470  print(M);
1471  fitting(M,-1);
1472  fitting(M,0);
1473  fitting(M,1);
1474  fitting(M,2);
1475}
1476///////////////////////////////////////////////////////////////////////////////
1477proc isLocallyFree(matrix S, int r)
1478"USAGE:   isLocallyFree(M,r);  M module, r int
1479RETURN:  1 if M'=coker(M) is locally free of constant rank r;@*
1480         0 if this is not the case.
1481EXAMPLE: example isLocallyFree; shows an example.
1482"
1483{
1484   ideal F=fitting(S,r);
1485   ideal G=fitting(S,r-1);
1486   if((deg(F[1])==0)&&(size(G)==0)){return(1);}
1487   return(0);
1488}
1489example
1490{"EXAMPLE:";   echo=2;
1491  ring R=0,(x,y,z),dp;
1492  matrix M[2][3];     // the presentation matrix
1493  M=x-1,y-1,z,y-1,x-2,x;
1494  ideal I=fitting(M,0); // 0-th Fitting ideal of coker(M)
1495  qring Q=I;
1496  matrix M=fetch(R,M);
1497  isLocallyFree(M,1); // as R/I-module, coker(M) is locally free of rk 1
1498  isLocallyFree(M,0);
1499}
1500///////////////////////////////////////////////////////////////////////////////
1501proc flatteningStrat (module M)
1502"USAGE:   flatteningStrat(M);  M module
1503RETURN:  list of ideals.
1504         The list entries L[1],...,L[r] describe the flattening stratification
1505         of M'=coker(M): setting L[0]=0, L[r+1]=1, the flattening
1506         stratification is given by the open sets Spec(A/V(L[i-1])) \ V(L[i]),
1507         i=1,...,r+1  (A = basering).
1508NOTE:    for more information see the book 'A Singular Introduction to
1509         Commutative Algebra' (by Greuel/Pfister, Springer 2002).
1510EXAMPLE: example flatteningStrat; shows an example
1511"
1512{
1513   list l;
1514   int v,w;
1515   ideal F;
1516   while(1)
1517   {
1518      F=interred(fitting(M,w));
1519      if(F[1]==1){return(l);}
1520      if(size(F)!=0){v++;l[v]=F;}
1521      w++;
1522   }
1523   return(l);
1524}
1525example
1526{"EXAMPLE:";   echo=2;
1527  ring A = 0,x(0..4),dp;
1528  // presentation matrix:
1529  matrix M[2][4] = x(0),x(1),x(2),x(3),x(1),x(2),x(3),x(4);
1530  list L = flatteningStrat(M);
1531  L;
1532}
1533///////////////////////////////////////////////////////////////////////////////
1534proc isFlat(module M)
1535"USAGE:   isFlat(M);  M module
1536RETURN:  1 if M'=coker(M) is flat;@*
1537         0 if this is not the case.
1538EXAMPLE: example isFlat; shows an example.
1539"
1540{
1541   if (size(ideal(M))==0) {return(1);}
1542   int w;
1543   ideal F=fitting(M,0);
1544   while(size(F)==0)
1545   {
1546      w++;
1547      F=fitting(M,w);
1548   }
1549   if (deg(std(F)[1])==0) {return(1);}
1550   return(0);
1551}
1552example
1553{"EXAMPLE:";   echo=2;
1554  ring A = 0,(x,y),dp;
1555  matrix M[3][3] = x-1,y,x,x,x+1,y,x2,xy+x+1,x2+y;
1556  print(M);
1557  isFlat(M);             // coker(M) is not flat over A=Q[x,y]
1558
1559  qring B = std(x2+x-y);   // the ring B=Q[x,y]/<x2+x-y>
1560  matrix M = fetch(A,M);
1561  isFlat(M);             // coker(M) is flat over B
1562
1563  setring A;
1564  qring C = std(x2+x+y);   // the ring C=Q[x,y]/<x2+x+y>
1565  matrix M = fetch(A,M);
1566  isFlat(M);             // coker(M) is not flat over C
1567}
1568///////////////////////////////////////////////////////////////////////////////
1569proc flatLocus(module M)
1570"USAGE:   flatLocus(M);  M module
1571RETURN:  ideal I, s.th. complement of V(I) is flat locus of coker(M).
1572NOTE:    computation is based on Fitting ideals;@*
1573         output is not radical (in general)
1574EXAMPLE: example flatLocus; shows an example
1575"
1576{
1577   if (size(ideal(M))==0) {return(ideal(1));}
1578   int v,w;
1579   ideal F=fitting(M,0);
1580   while(size(F)==0)
1581   {
1582      w++;
1583      F=fitting(M,w);
1584   }
1585   if(size(ideal(basering))>0)
1586   {
1587      for(v=w+1;v<=nrows(M);v++)
1588      {
1589         F=F+intersect(fitting(M,v),quotient(ideal(0),fitting(M,v-1)));
1590      }
1591   }
1592   return(interred(F));
1593}
1594example
1595{"EXAMPLE:";   echo=2;
1596  ring R=0,(x,y,z),dp;
1597  matrix M[2][3]=x,y,z,0,x3,z3;
1598  ideal I=flatLocus(M); // coker(M) is flat outside V(x,yz)
1599  I;                    // computed ideal not radical
1600  ideal J=radical(I);
1601  J;
1602
1603  qring r=std(J);
1604  matrix M=fetch(r,M);
1605  flatLocus(M);         // coker(M) is flat over Spec(Q[x,y,z]/<x,yz>)
1606
1607  isFlat(M);            // flatness test
1608}
1609///////////////////////////////////////////////////////////////////////////////
1610proc isReg(ideal I, module N)
1611"USAGE:   isReg(I,M);  I ideal, M module
1612RETURN:  1 if given (ordered) list of generators for I is coker(M)-sequence;@*
1613         0 if this is not the case.
1614EXAMPLE: example isReg; shows an example.
1615"
1616{
1617  int n=nrows(N);
1618  int i;
1619  while(i<ncols(I))
1620  {
1621     i++;
1622     N=std(N);
1623     if(size(reduce(quotient(N,I[i]),N))!=0){return(0);}
1624     N=N+I[i]*freemodule(n);
1625  }
1626  if (size(reduce(freemodule(n),std(N)))==0){return(0);}
1627  return(1);
1628}
1629example
1630{"EXAMPLE:";   echo=2;
1631  ring R = 0,(x,y,z),dp;
1632  ideal I = x*(y-1),y,z*(y-1);
1633  isReg(I,0);             // given list of generators is Q[x,y,z]-sequence
1634
1635  I = x*(y-1),z*(y-1),y;  // change sorting of generators
1636  isReg(I,0);
1637
1638  ring r = 0,(x,y,z),ds;  // local ring
1639  ideal I=fetch(R,I);
1640  isReg(I,0);             // result independent of sorting of generators
1641}
1642
1643///////////////////////////////////////////////////////////////////////////////
1644// the following static procedures are used by KoszulHomology:
1645//  * binom_int  (binomial coeff. as integer, or -1 if too large)
1646//  * basisNumber
1647//  * basisElement
1648//  * KoszulMap
1649// for details, see 'A Singular Introduction to Commutative Algebra' (by
1650// Greuel/Pfister, Springer 2002), Chapter 7
1651
1652static proc binom_int(int n, int p)
1653{
1654  bigint s = binomial(n,p);
1655  int a=int(s);
1656  if ((s!=0)&&(a==0)) { return(-1); }
1657  return(a);
1658}
1659
1660static proc basisNumber(int n,intvec v)
1661{
1662  int p=size(v);
1663  if(p==1){return(v[1]);}
1664  int j=n-1;
1665  int b;
1666  while(j>=n-v[1]+1)
1667  {
1668    b=b+binom_int(j,p-1);
1669    j--;
1670  }
1671  intvec w=v-v[1];
1672  w=w[2..size(w)];
1673  b=b+basisNumber(n-v[1],w);
1674  return(b);
1675}
1676
1677static proc basisElement(int n,int p,int N)
1678{
1679  if(p==1){return(N);}
1680  int s,R;
1681  while(R<N)
1682  {
1683     s++;
1684     R=R+binom_int(n-s,p-1);
1685  }
1686  R=N-R+binom_int(n-s,p-1);
1687  intvec v=basisElement(n-s,p-1,R);
1688  intvec w=s,v+s;
1689  return(w);
1690}
1691
1692proc KoszulMap(ideal x,int p)
1693{
1694  int n=size(x);
1695  int a=binom_int(n,p-1);
1696  int b=binom_int(n,p);
1697
1698  matrix M[a][b];
1699
1700  if(p==1){M=x;return(M);}
1701  int j,k;
1702  intvec v,w;
1703  for(j=1;j<=b;j++)
1704  {
1705    v=basisElement(n,p,j);
1706    w=v[2..p];
1707    M[basisNumber(n,w),j]=x[v[1]];
1708    for(k=2;k<p;k++)
1709    {
1710      w=v[1..k-1],v[k+1..p];
1711      M[basisNumber(n,w),j]=(-1)^(k-1)*x[v[k]];
1712    }
1713    w=v[1..p-1];
1714    M[basisNumber(n,w),j]=(-1)^(p-1)*x[v[p]];
1715  }
1716  return(M);
1717}
1718///////////////////////////////////////////////////////////////////////////////
1719
1720proc KoszulHomology(ideal x, module M, int p)
1721"USAGE:   KoszulHomology(I,M,p);  I ideal, M module, p int
1722COMPUTE: A presentation of the p-th Koszul homology module H_p(f_1,...,f_k;M'),
1723         where M'=coker(M) and f_1,...,f_k are the given (ordered list
1724         of non-zero) generators of the ideal I.
1725         The computed presentation is minimized via prune.
1726         In particular, if H_p(f_1,...,f_k;M')=0 then the return value is 0.
1727RETURN:  module H, s.th. coker(H) = H_p(f_1,...,f_k;M').
1728NOTE:    size of input ideal has to be <= 20.
1729EXAMPLE: example KoszulHomology; shows an example.
1730{
1731  x=simplify(x,2);
1732  int n     = size(x);
1733  if (n==0)
1734  {
1735    ERROR("// KoszulHomology only for non-zero ideals");
1736  }
1737  if (n>20)
1738  {
1739    ERROR("// too many generators in input ideal");
1740  }
1741  if (p>n)
1742  {
1743    module hom=0;
1744    return(hom);
1745  }
1746
1747  int a     = binom_int(n,p-1);  // n over p-1 independent of char(basering)
1748  int b     = binom_int(n,p);
1749
1750  matrix N  = matrix(M);
1751  module ker= freemodule(nrows(N));
1752  if(p!=0)
1753  {
1754    module im= tensor(unitmat(a),N);
1755    module f = tensor(KoszulMap(x,p),unitmat(nrows(N)));
1756    ker      = modulo(f,im);
1757  }
1758  module im1 = tensor(unitmat(b),N);
1759  module im2 = tensor(KoszulMap(x,p+1),unitmat(nrows(N)));
1760  module hom = modulo(ker,im1+im2);
1761  hom        = prune(hom);
1762  return(hom);
1763}
1764example
1765{"EXAMPLE:";   echo=2;
1766  ring R=0,x(1..3),dp;
1767  ideal x=maxideal(1);
1768  module M=0;
1769  KoszulHomology(x,M,0);  // H_0(x,R), x=(x_1,x_2,x_3)
1770
1771  KoszulHomology(x,M,1);  // H_1(x,R), x=(x_1,x_2,x_3)
1772
1773  qring S=std(x(1)*x(2));
1774  module M=0;
1775  ideal x=maxideal(1);
1776  KoszulHomology(x,M,1);
1777
1778  KoszulHomology(x,M,2);
1779}
1780///////////////////////////////////////////////////////////////////////////////
1781proc depth(module M,ideal #)
1782"USAGE:   depth(M,[I]);  M module, I ideal
1783RETURN:  int,
1784         - if called with 1 argument: the depth of M'=coker(M) w.r.t. the
1785         maxideal in the basering (which is then assumed to be local)@*
1786         - if called with 2 arguments: the depth of M'=coker(M) w.r.t. the
1787         ideal I.
1788NOTE:    procedure makes use of KoszulHomology.
1789EXAMPLE: example depth; shows an example.
1790"
1791{
1792  ideal m=maxideal(1);
1793  int n=size(m);
1794  int i;
1795
1796  if (size(#)==0)
1797  {
1798    // depth(M') over local basering
1799    while(i<n)
1800    {
1801      i++;
1802      if(size(KoszulHomology(m,M,i))==0){return(n-i+1);}
1803    }
1804    return(0);
1805  }
1806
1807  ideal I=simplify(#,2);
1808  while(i<size(I))
1809  {
1810    i++;
1811    if(size(KoszulHomology(I,M,i))==0){return(size(I)-i+1);}
1812  }
1813  return(0);
1814}
1815example
1816{"EXAMPLE:";   echo=2;
1817  ring R=0,(x,y,z),dp;
1818  ideal I=x2,xy,yz;
1819  module M=0;
1820  depth(M,I);   // depth(<x2,xy,yz>,Q[x,y,z])
1821  ring r=0,(x,y,z),ds;  // local ring
1822  matrix M[2][2]=x,xy,1+yz,0;
1823  print(M);
1824  depth(M);     // depth(maxideal,coker(M))
1825  ideal I=x;
1826  depth(M,I);   // depth(<x>,coker(M))
1827  I=x+z;
1828  depth(M,I);   // depth(<x+z>,coker(M))
1829}
1830
1831///////////////////////////////////////////////////////////////////////////////
1832proc isCM(module M)
1833"USAGE:   isCM(M);  M module
1834RETURN:  1 if M'=coker(M) is Cohen-Macaulay;@*
1835         0 if this is not the case.
1836ASSUME:  basering is local.
1837EXAMPLE: example isCM; shows an example
1838"
1839{
1840  // test if basering is local:
1841  ideal m=maxideal(1);
1842  int i;
1843  poly f=1;
1844  for (i=1; i<=size(m); i++)
1845  {
1846    f=f+m[i];
1847  }
1848  if (ord(f)>0)
1849  {
1850    print("// basering must be local -- result has no meaning");
1851    return(0);
1852  }
1853
1854  return(depth(M)==dim(std(Ann(M))));
1855}
1856example
1857{"EXAMPLE:";   echo=2;
1858  ring R=0,(x,y,z),ds;  // local ring R = Q[x,y,z]_<x,y,z>
1859  module M=xz,yz,z2;
1860  isCM(M);             // test if R/<xz,yz,z2> is Cohen-Macaulay
1861
1862  M=x2+y2,z7;          // test if R/<x2+y2,z7> is Cohen-Macaulay
1863  isCM(M);
1864}
1865
1866proc canonMap(list #)
1867"USAGE:  canonMap(id); id= ideal/module,
1868RETURN:  a list L, the kernel in two different representations and
1869@*       the cokernel of the canonical map
1870@*       M ---> Ext^c_R(Ext^c_R(M,R),R) given by presentations
1871@*       Here M is the R-module (R=basering) given by the presentation
1872@*       defined by id, i.e. M=R/id resp. M=R^n/id
1873@*       c is the codimension of M
1874@*       L[1] is the preimage of the kernel in R resp. R^n
1875@*       L[2] is a presentation of the kernel
1876@*       L[3] is a presentation of the cokernel
1877EXAMPLE: example canonMap; shows an example
1878"
1879{
1880   module M=#[1];
1881   int c=nvars(basering)-dim(std(M));
1882   if(c==0)
1883   {
1884      module K=syz(transpose(M));
1885      module Ke=syz(transpose(K));
1886      module Co=modulo(syz(transpose(syz(K))),transpose(K));
1887    }
1888   else
1889   {
1890      int i;
1891      resolution F=mres(M,c+1);
1892      module K=syz(transpose(F[c+1]));
1893      K=simplify(reduce(K,std(transpose(F[c]))),2);
1894      module A=modulo(K,transpose(F[c]));
1895      resolution G=nres(A,c+1);
1896      for(i=1;i<=c;i++)
1897      {
1898         K=lift(transpose(F[c-i+1]),K*G[i]);
1899      }
1900      module Ke=modulo(transpose(K),transpose(G[c]));
1901      module Co=modulo(syz(transpose(G[c+1])),transpose(K)+transpose(G[c]));
1902   }
1903   return(list(Ke,prune(modulo(Ke,M)),prune(Co)));
1904
1905}
1906example
1907{ "EXAMPLE:"; echo = 2;
1908   ring s=0,(x,y),dp;
1909   ideal i = x,y;
1910   canonMap(i);
1911   ring R = 0,(x,y,z,w),dp;
1912   ideal I1 = x,y;
1913   ideal I2 = z,w;
1914   ideal I = intersect(I1,I2);
1915   canonMap(I);
1916   module M = syz(I);
1917   canonMap(M);
1918   ring S = 0,(x,y,z,t),Wp(3,4,5,1);
1919   ideal I = x-t3,y-t4,z-t5;
1920   ideal J = eliminate(I,t);
1921   ring T = 0,(x,y,z),Wp(3,4,5);
1922   ideal p = imap(S,J);
1923   ideal p2 = p^2;
1924   canonMap(p2);
1925}
1926
1927// taken from qhmoduli.lib
1928static proc Max(def data)
1929"USAGE:   Max(data); intvec/list of integers
1930PURPOSE: find the maximal integer contained in 'data'
1931RETURN:  list
1932ASSUME:  'data' contians only integers and is not empty
1933"
1934{
1935  int i;
1936  int max = data[1];
1937
1938  for(i = size(data); i>1;i--)
1939  {
1940    if(data[i] > max) { max = data[i]; }
1941  }
1942  return(max);
1943}
1944example
1945{"EXAMPLE:";  echo = 2;
1946  Max(list(1,2,3));
1947}
Note: See TracBrowser for help on using the repository browser.