source: git/Singular/LIB/homolog.lib @ 4d68980

// $Id: homolog.lib,v 1.10 1999-12-13 15:33:47 obachman Exp $
//(BM/GMG)
///////////////////////////////////////////////////////////////////////////////

version="$Id: homolog.lib,v 1.10 1999-12-13 15:33:47 obachman Exp $";
info="
LIBRARY:  homolog.lib   PROCEDURES FOR HOMOLOGICAL ALGEBRA

AUTHORS:  Gert-Martin Greuel, email: greuel@mathematik.uni-kl.de, 
          Bernd Martin, email: martin@math.tu-cottbus.de

PROCEDURES:
 cup(M);                cup: Ext^1(M',M') x Ext^1() --> Ext^2()
 cupproduct(M,N,P,p,q); cup: Ext^p(M',N') x Ext^q(N',P') --> Ext^p+q(M',P')
 Ext_R(k,M);            Ext^k(M',R),    M module, R basering, M'=coker(M)
 Ext(k,M,N);            Ext^k(M',N'),   M,N modules, M'=coker(M), N'=coker(N)
 Hom(M,N);              Hom(M',N'),     M,N modules, M'=coker(M), N'=coker(N)
 homology(A,B,M,N)      ker(B)/im(A),   homology of complex R^k--A->M'--B->N'
 kernel(A,M,N);         ker(M'--A->N')  M,N modules, A matrix
 kohom(A,k);            Hom(R^k,A),     A matrix over basering R
 kontrahom(A,k);        Hom(A,R^k),     A matrix over basering R
";

LIB "general.lib";
LIB "deform.lib";
LIB "matrix.lib";
LIB "poly.lib";
///////////////////////////////////////////////////////////////////////////////

proc cup (module M,list #)
"USAGE:   cup(M,[,any,any]);  M=module
COMPUTE: cup-product Ext^1(M',M') x Ext^1(M',M') ---> Ext^2(M',M'),
         where M':=R^m/M, if M in R^m, R basering (i.e. M':=coker(matrix(M)))
         in case of a second argument: symmetrized cup-product
ASSUME:  all Ext's  are finite dimensional
RETURN:  matrix of the associated linear map,
         i.e. the columns of <matrix> present the coordinates of b_i & b_j
         (resp. (1/2)(b_i & b_j + b_j & b_i) in the symmetric version)
         with respect to a kbase of Ext^2,
         (where b_1,b_2,... is a kbase of Ext^1 and & denotes cup product)
         in case of  a third argument return a list:
               L[1] = matrix see above (and symmetric case)
               L[2] = matrix of kbase of Ext^1
               L[3] = matrix of kbase of Ext^2
NOTE:    printlevel >=1;  shows what is going on
         printlevel >=2;  shows result in another representation
         For computing cupproduct of M itself, apply proc to syz(M)!
EXAMPLE: example cup; shows examples
"
{
//---------- initialization ---------------------------------------------------
   int    i,j,k,f0,f1,f2,f3,e1,e2;
   module M1,M2,A,B,C,ker,ima,ext1,ext2,ext10,ext20;
   matrix cup[1][0];
   matrix kb1,lift1,kb2,mA,mB,mC;
   ideal  tes1,tes2,null;
   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//-----------------------------------------------------------------------------
//take a resolution of M<--F(0)<--- ...  <---F(3)
//apply Hom(-,M) and compute the Ext's
//-----------------------------------------------------------------------------
   list resM = nres(M,3);
   M1 = resM[2];
   M2 = resM[3];
   f0 = nrows(M);
   f1 = ncols(M);
   f2 = ncols(M1);
   f3 = ncols(M2);
   tes1 = simplify(ideal(M),10);
   tes2=simplify(ideal(M1),10);
   if ((tes1[1]*tes2[1]==0) or (tes1[1]==1) or (tes2[1]==1))
   {
      dbprint(p,"// Ext == 0 , hence 'cup' is the zero-map");
      return(@cup);
   }
//------ compute Ext^1 --------------------------------------------------------
   B     = kohom(M,f2);
   A     = kontrahom(M1,f0);
   C     = intersect(A,B);
   C     = reduce(C,std(null));C = simplify(C,10);
   ker   = lift(A,C)+syz(A);
   ima   = kohom(M,f1);
   ima   = ima + kontrahom(M,f0);
   ext1  = modulo(ker,ima);
   ext10 = std(ext1);
   e1    = vdim(ext10);
   dbprint(p-1,"// vdim (Ext^1) = "+string(e1));
   if (e1 < 0)
   {
     "// Ext^1 not of finite dimension";
     return(cup);
   }
   kb1 = kbase(ext10);
   kb1 = matrix(ker)*kb1;
   dbprint(p-1,"// kbase of Ext^1(M,M)",
               "//  - the columns present the kbase elements in Hom(F(1),F(0))",
               "//  - F(*) a free resolution of M",kb1);

//------ compute the liftings of Ext^1 ----------------------------------------
   mC = matrix(A)*kb1;
   lift1 =lift(B,mC);
   dbprint(p-1,"// lift kbase of Ext^1:",
    "//  - the columns present liftings of kbase elements into Hom(F(2),F(1))",
    "//  - F(*) a free resolution of M ",lift1);

//------ compute Ext^2  -------------------------------------------------------
   B    = kohom(M,f3);
   A    = kontrahom(M2,f0);
   C    = intersect(A,B);
   C    = reduce(C,std(null));C = simplify(C,10);
   ker  = lift(A,C)+syz(A);
   ima  = kohom(M,f2);
   ima  = ima + kontrahom(M1,f0);
   ext2 = modulo(ker,ima);
   ext20= std(ext2);
   e2   = vdim(ext20);
   if (e2<0)
   {
     "// Ext^2 not of finite dimension";
     return(cup);
   }
   dbprint(p-1,"// vdim (Ext^2) = "+string(e2));
   kb2 = kbase(ext20);
   kb2 = matrix(ker)*kb2;
   dbprint(p-1,"// kbase of Ext^2(M,M)",
               "//  - the columns present the kbase elements in Hom(F(2),F(0))",
               "//  - F(*) is a a free resolution of M ",kb2);
//-------  compute: cup-products of base-elements -----------------------------
   for (i=1;i<=e1;i=i+1)
   {
     for (j=1;j<=e1;j=j+1)
     {
       mA = matrix(ideal(lift1[j]),f1,f2);
       mB = matrix(ideal(kb1[i]),f0,f1);
       mC = mB*mA;
       if (size(#)==0)
       {                                          //non symmestric
         mC = matrix(ideal(mC),f0*f2,1);
         cup= concat(cup,mC);
       }
       else                                       //symmetric version
       {
         if (j>=i)
         {
           if (j>i)
           {
             mA = matrix(ideal(lift1[i]),f1,f2);
             mB = matrix(ideal(kb1[j]),f0,f1);
             mC = mC+mB*mA;mC=(1/2)*mC;
           }
           mC = matrix(ideal(mC),f0*f2,1);
           cup= concat(cup,mC);
         }
       }
     }
   }
   dbprint(p-1,"// matrix of cup-products (in Ext^2)",cup,"////// end level 2 //////");
//------- comptute: presentation of base-elements -----------------------------
   cup = lift(ker,cup);
   cup = lift_kbase(cup,ext20);
   if( p>2 )
   {
     "// the associated matrices of the bilinear mapping 'cup' ";
     "// corresponding to the kbase elements of Ext^2(M,M) are shown,";
     "//  i.e. the rows of the final matrix are written as matrix of";
     "//  a bilinear form on Ext^1 x Ext^1";
     matrix BL[e1][e1];
     for (k=1;k<=e2;k=k+1)
     {
       "//_____"+string(k)+". component:";
       for (i=1;i<=e1;i=i+1)
       {
         for (j=1;j<=e1;j=j+1)
         {
           if (size(#)==0) { BL[i,j]=cup[k,j+e1*(i-1)]; }
           else
           {
             if (i<=j)
             {
               BL[i,j]=cup[k,j+e1*(i-1)-binomial(i,2)];
               BL[j,i]=BL[i,j];
             }
           }
         }
       }
       print(BL);
     }
     "////// end level 3 //////";
   }
   if (size(#)>2) { return(cup,kb1,kb2);}
   else {return(cup);}
}
example
{"EXAMPLE";  echo=2;
  int p      = printlevel;
  ring  rr   = 32003,(x,y,z),(dp,C);
  ideal  I   = x4+y3+z2;
  qring  o   = std(I);
  module M   = [x,y,0,z],[y2,-x3,z,0],[z,0,-y,-x3],[0,z,x,-y2];
  print(cup(M));
  print(cup(M,1));
  // 2nd EXAMPLE  (shows what is going on)
  printlevel = 3;
  ring   r   = 0,(x,y),(dp,C);
  ideal  i   = x2-y3;
  qring  q   = std(i);
  module M   = [-x,y],[-y2,x];
  print(cup(M));
  printlevel = p;
}
///////////////////////////////////////////////////////////////////////////////

proc cupproduct (module M,N,P,int p,q,list #)
"USAGE:   cupproduct(M,N,P,p,q[,any]);  M,N,P modules, p,q integers
COMPUTE: cup-product Ext^p(M',N') x Ext^q(N',P') ---> Ext^p+q(M',P')
         where M':=R^m/M, if M in R^m, R basering (i.e. M':=coker(matrix(M)))
ASSUME:  all Ext's  are of finite dimension
RETURN:  matrix of the associated linear map Ext^p(tensor)Ext^q-->Ext^p+q
         i.e. the columnes of <matrix> present the coordinates of
         the cup products (b_i & c_j) with respect to a kbase of Ext^p+q
         (b_i resp. c_j are choosen bases of Ext^p resp. Ext^q)
         in case of a 6th argument:
            return a list
            L[1] = matrix (see above)
            L[2] = matrix of kbase of Ext^p(M',N')
            L[3] = matrix of kbase of Ext^q(N',P')
            L[4] = matrix of kbase of Ext^p+q(N',P')
NOTE:    printlevel >=1;  shows what is going on
         printlevel >=2;  shows result in another representation
         For computing cupproduct of M,N itself, apply proc to syz(M),syz(N)!
EXAMPLE: example cupproduct; shows examples
"
{
//---------- initialization ---------------------------------------------------
   int    e1,e2,e3,i,j,k,f0,f1,f2;
   module M1,M2,N1,N2,P1,P2,A,B,C,ker,ima,extMN,extMN0,extMP,
          extMP0,extNP,extNP0;
   matrix cup[1][0];
   matrix kbMN,kbMP,kbNP,lift1,mA,mB,mC;
   ideal  test1,test2,null;
   int pp = printlevel-voice+3;  // pp=printlevel+1 (default: p=1)
//-----------------------------------------------------------------------------
//compute resolutions of M and N
//                     M<--F(0)<--- ...  <---F(p+q+1)
//                     N<--G(0)<--- ...  <---G(q+1)
//-----------------------------------------------------------------------------
   list resM = nres(M,p+q+1);
   M1 = resM[p];
   M2 = resM[p+1];
   list resN = nres(N,q+1);
   N1 = resN[q];
   N2 = resN[q+1];
   P1 = resM[p+q];
   P2 = resM[p+q+1];
//-------test: Ext==0?---------------------------------------------------------
   test1 = simplify(ideal(M1),10);
   test2 = simplify(ideal(N),10);
   if (test1[1]==0) { dbprint(pp,"//Ext(M,N)=0");return(cup); }
   test1 = simplify(ideal(N1),10);
   test2 = simplify(ideal(P),10);
   if (test1[1]==0) { dbprint(pp,"//Ext(N,P)=0");return(cup); }
   test1 = simplify(ideal(P1),10);
   if (test1[1]==0) { dbprint(pp,"//Ext(M,P)=0");return(cup); }
 //------ compute kbases of Ext's ---------------------------------------------
 //------ Ext(M,N)
   test1 = simplify(ideal(M2),10);
   if (test1[1]==0) { ker = freemodule(ncols(M1)*nrows(N));}
   else
   {
     A   = kontrahom(M2,nrows(N));
     B   = kohom(N,ncols(M2));
     C   = intersect(A,B);
     C   = reduce(C,std(ideal(0)));C=simplify(C,10);
     ker = lift(A,C)+syz(A);
   }
   ima   = kohom(N,ncols(M1));
   A     = kontrahom(M1,nrows(N));
   ima   = ima+A;
   extMN = modulo(ker,ima);
   extMN0= std(extMN);
   e1    = vdim(extMN0);
   dbprint(pp-1,"// vdim Ext(M,N) = "+string(e1));
   if (e1 < 0)
   {
     "// Ext(M,N) not of finite dimension";
     return(cup);
   }
   kbMN  = kbase(extMN0);
   kbMN = matrix(ker)*kbMN;
   dbprint(pp-1,"// kbase of Ext^p(M,N)",
          "//  - the columns present the kbase elements in Hom(F(p),G(0))",
          "//  - F(*),G(*) are free resolutions of M and N",kbMN);
//------- Ext(N,P)
   test1 = simplify(ideal(N2),10);
   if (test1[1]==0) {  ker = freemodule(ncols(N1)*nrows(P)); }
   else
   {
     A = kontrahom(N2,nrows(P));
     B = kohom(P,ncols(N2));
     C = intersect(A,B);
     C = reduce(C,std(ideal(0)));C=simplify(C,10);
     ker = lift(A,C)+syz(A);
   }
   ima   = kohom(P,ncols(N1));
   A     = kontrahom(N1,nrows(P));
   ima   = ima+A;
   extNP = modulo(ker,ima);
   extNP0= std(extNP);
   e2    = vdim(extNP0);
   dbprint(pp-1,"// vdim Ext(N,P) = "+string(e2));
   if (e2 < 0)
   {
     "// Ext(N,P) not of finite dimension";
     return(cup);
   }
   kbNP  = kbase(extNP0);
   kbNP  = matrix(ker)*kbNP;
   dbprint(pp-1,"// kbase of Ext(N,P):",kbNP,
          "// kbase of Ext^q(N,P)",
          "//  - the columns present the kbase elements in Hom(G(q),H(0))",
          "//  - G(*),H(*) are free resolutions of N and P",kbNP);

//------ Ext(M,P)
   test1 = simplify(ideal(P2),10);
   if (test1[1]==0) { ker = freemodule(ncols(P1)*nrows(P)); }
   else
   {
     A = kontrahom(P2,nrows(P));
     B = kohom(P,ncols(P2));
     C = intersect(A,B);
     C = reduce(C,std(ideal(0)));C=simplify(C,10);
     ker = lift(A,C)+syz(A);
   }
   ima   = kohom(P,ncols(P1));
   A     = kontrahom(P1,nrows(P));
   ima   = ima+A;
   extMP = modulo(ker,ima);
   extMP0= std(extMP);
   e3    = vdim(extMP0);
   dbprint(pp-1,"// vdim Ext(M,P) = "+string(e3));
   if (e3 < 0)
   {
     "// Ext(M,P) not of finite dimension";
     return(cup);
   }
   kbMP  = kbase(extMP0);
   kbMP  = matrix(ker)*kbMP;
   dbprint(pp-1,"// kbase of Ext^p+q(M,P)",
          "//  - the columns present the kbase elements in Hom(F(p+q),H(0))",
          "//  - F(*),H(*) are free resolutions of M and P",kbMP);
//----- lift kbase of Ext(M,N) ------------------------------------------------
   lift1 = kbMN;
   for (i=1;i<=q;i=i+1)
   {
     mA = kontrahom(resM[p+i],nrows(resN[i]));
     mB = kohom(resN[i],ncols(resM[p+i]));
     lift1 = lift(mB,mA*lift1);
   }
   dbprint(pp-1,"// lifting of kbase of Ext^p(M,N)",
   "//  - the columns present the lifting of kbase elements in Hom(F(p+q),G(q))",lift1);
//-------  compute: cup-products of base-elements -----------------------------
   f0 = nrows(P);
   f1 = ncols(N1);
   f2 = ncols(resM[p+q]);
   for (i=1;i<=e1;i=i+1)
   {
     for (j=1;j<=e2;j=j+1)
     {
       mA = matrix(ideal(lift1[j]),f1,f2);
       mB = matrix(ideal(kbMP[i]),f0,f1);
       mC = mB*mA;
       mC = matrix(ideal(mC),f0*f2,1);
       cup= concat(cup,mC);
     }
   }
   dbprint(pp-1,"// matrix of cup-products (in Ext^p+q)",cup,"////// end level 2 //////");
//------- comptute: presentation of base-elements -----------------------------
   cup = lift(ker,cup);
   cup = lift_kbase(cup,extMP0);
//------- special output ------------------------------------------------------
   if (pp>2)
   {
     "// the associated matrices of the bilinear mapping 'cup' ";
     "// corresponding to the kbase elements of Ext^p+q(M,P) are shown,";
     "//  i.e. the rows of the final matrix are written as matrix of";
     "//  a bilinear form on Ext^p x Ext^q";
     matrix BL[e1][e2];
     for (k=1;k<=e3;k=k+1)
     {
       "//----"+string(k)+". component:";
       for (i=1;i<=e1;i=i+1)
       {
         for (j=1;j<=e2;j=j+1)
         {
           BL[i,j]=cup[k,j+e1*(i-1)];
         }
        }
        print(BL);
      }
     "////// end level 3 //////";
   }
   if (size(#)) { return(cup,kbMN,kbNP,kbMP);}
   else         { return(cup); }
}
example
{"EXAMPLE";  echo=2;
  int p      = printlevel;
  ring  rr   = 32003,(x,y,z),(dp,C);
  ideal  I   = x4+y3+z2;
  qring  o   = std(I);
  module M   = [x,y,0,z],[y2,-x3,z,0],[z,0,-y,-x3],[0,z,x,-y2];
  print(cupproduct(M,M,M,1,3));
  printlevel = 3;
  list l     = (cupproduct(M,M,M,1,3,"any"));
  show(l[1]);show(l[2]);
  printlevel = p;
}
///////////////////////////////////////////////////////////////////////////////

proc Ext_R (intvec v, module M, list #)
"USAGE:   Ext_R(v,M[,p]);  v=int/intvec , M=module, p=int
COMPUTE: A presentation of Ext^k(M',R); for k=v[1],v[2],..., M'=coker(M).
         Let ...--> F2 --> F1 --M-> F0-->M'-->0 be a free resolution of M'. If
         0 <-- F0* <-A1-- F1* <-A2-- F2* <-A3--... denotes the dual sequence,
         Fi*=Hom(Fi,R), then Ext^k = ker(Ak)/im(Ak+1) is presented as in the
         following exact sequences:
                 Fk-1* <-Ak-- Fk* <-syz(Ak)-- R^p
                 Fk*/im(Ak+1) <-syz(Ak)-- R^p <-Ext^k-- R^q
         Hence Ext^k=modulo(syz(Ak),Ak+1) presents Ext^k(M',R);
RETURN:  Ext^k, of type module, a presentation of Ext^k(M',R) if v is of type
         int, resp. a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.
         In case of a third argument of type int return a list:
           [1] = module Ext^k/list of Ext^k
           [2] = SB of Ext^k/list of SB of Ext^k
           [3] = matrix/list of matrices, each representing kbase of Ext^k
                 (if finite dimensional)
DISPLAY: printlevel >=0: degree of Ext^k for each k  (default)
         printlevel >=1: Ak, Ak+1 and kbase of Ext^k in Fk*
NOTE:    In order to compute Ext^k(M,R) use the command Ext_R(k,syz(M));
         or the 2 commands: list L=mres(M,2);  Ext_R(k,L[2]);
EXAMPLE: example Ext_R; shows examples
"
{

// In case M is known to be a SB, set attrib(M,"isSB",1); in order to
// avoid  unnecessary SB computations

//------------ initialisation -------------------------------------------------
  module m1,m2,ret,ret0;
  matrix ker,kb;
  list L1,L2,L3,L,resl,K;
  int k,max,ii,t1,t2;
  int s = size(v);
  intvec v1 = sort(v)[1];
  max = v1[s];                 // the maximum integer occuring in intvec v
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
  // --------------- Variante mit sres
  for( ii=1; ii<=size(#); ii++ )
  {
    t2=1; // return a list if t2=1
    if( typeof(#[ii])=="string" )
    {
      if ( #[ii]=="sres" ) { t1=1; t2=0; } // use sres instead of mres if t1=1
    }
  }
//----------------- compute resolution of coker(M) ----------------------------
  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
  if( t1==1 )
  {
    if( attrib(M,"isSB")==0 ) { M=std(M); }
    resl = sres(M,max+1);
  }
  else { resl = mres(M,max+1); }
  for( ii=1; ii<=s; ii++ )
  {
//-----------------  apply Hom(_,R) at k-th place -----------------------------
    k=v[ii];
    if( k<0 )                   // Ext^k=0 for negative k
    {
      dbprint(p-1,"// Ext^i=0 for i<0!");
      ret    = gen(1);
      ret0   = std(ret);
      L1[ii] = ret;
      L2[ii] = ret0;
      L3[ii] = matrix(kbase(ret0));
      dbprint(p,"// degree of Ext^"+string(k)+":");
      if( p>=0 ) { degree(ret0);"";}
    }
    else
    {
      m2 = transpose(resl[k+1]);
      if( k==0 ) { m1=0*gen(nrows(m2)); }
      else { m1 = transpose(resl[k]); }
//----------------- presentation of ker(m2)/im(m1) ----------------------------
      ker = syz(m2);
      ret = modulo(ker,m1);
      dbprint(p-1,"// Computing Ext^"+string(k)+":",
         "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of M,",
         "// then F"+string(k)+"*-->F"+string(k+1)+"* is given by:",m2,
         "// and F"+string(k-1)+"*-->F"+string(k)+"* is given by:",m1,"");
      ret0 = std(ret);
      dbprint(p,"// degree of Ext^"+string(k)+":");
      if( p>0 ) { degree(ret0);"";}
      if( t2 )
      {
         if( vdim(ret0)>=0 )
         {
            kb = kbase(ret0);
            if ( size(ker)!=0 ) { kb = matrix(ker)*kb; }
            dbprint(p-1,
            "// columns of matrix are kbase of Ext^"+string(k)+" in F"+string(k)+"*:",kb,"");
            L3[ii] = kb;
         }
         L2[ii] = ret0;
      }
      L1[ii] = ret;
    }
  }
  if( t2 )
  {
     if( s>1 ) { L = L1,L2,L3; return(L); }
     else { L = ret,ret0,kb; return(L); }
  }
  else
  {
     if( s>1 ) { return(L1); }
     else { return(ret); }
  }
}
example
{"EXAMPLE:";     echo=2;
  int p      = printlevel;
  printlevel = 1;
  ring r     = 0,(x,y,z),dp;
  ideal i    = x2y,y2z,z3x;
  module E   = Ext_R(1,i);       //computes Ext^1(r/i,r)
  is_zero(E);
  module m   = [x],[0,y];
  list L     = Ext_R(2..3,m);    //computes Ext^i(r^2/m,r), i=2,3
  show(L);"";
  qring R    = std(x2+yz);
  intvec v   = 0,2,4;
  printlevel = 2;                //shows what is going on
  ideal i    = x,y,z;            //computes Ext^i(r/(x,y,z),r/(x2+yz)), i=0,2,4
  list L     = Ext_R(v,i,1);     //over the qring R=r/(x2+yz), std and kbase
  printlevel = p;
}
///////////////////////////////////////////////////////////////////////////////

proc Ext (intvec v, module M, module N, list #)
"USAGE:   Ext(v,M,N[,any]);  v=int/intvec, M,N=modules
COMPUTE: A presentation of Ext^k(M',N'); for k=v[1],v[2],... where
         M'=coker(M) and N'=coker(N). Let
         0<--M'<-- F0 <-M-- F1 <-- F2 <--...  resp. 0<--N'<-- G0 <--N- G1 be
         a free resolution of M' resp. a presentations of N'. Consider
@format
                      0                  0                  0
                      |^                 |^                 |^
               --> Hom(Fk-1,N') -Ak-> Hom(Fk,N') -Ak+1-> Hom(Fk+1,N')
                      |^                 |^                 |^
               --> Hom(Fk-1,G0) -Ak-> Hom(Fk,G0) -Ak+1-> Hom(Fk+1,G0)
                                         |^                 |^
                                         |C                 |B
                                      Hom(Fk,G1) -----> Hom(Fk+1,G1)
@end format
         (Ak,Ak+1 induced by M and B,C induced by N).
         Let K=modulo(Ak+1,B), J=module(Ak)+module(C) and Ext=modulo(K,J),
         then we have exact sequences
                  R^p  --K-> Hom(Fk,G0) --Ak+1-> Hom(Fk+1,G0)/im(B)
         R^q -Ext-> R^p --K->Hom(Fk,G0)/im(Ak)+im(C) --Ak+1->Hom(Fk+1,G0)/im(B)
         Hence Ext presents Ext^k(M',N')
RETURN:  
         Ext, of type module, a presentation of Ext^k(M',N') if v is of type
         int, resp. a list of Ext^k (k=v[1],v[2],...) if v is of type intvec.
         In case of a third argument of any type return a list:
           [1] = module Ext/list of Ext^k
           [2] = SB of Ext/list of SB of Ext^k
           [3] = matrix/list of matrices, each representing a kbase of Ext^k
                 (if finite dimensional)
DISPLAY: printlevel >=0: degree of Ext^k for each k  (default)
         printlevel >=1: matrices Ak, Ak+1 and kbase of Ext^k in Hom(Fk,G0)
                            (if finite dimensional)
NOTE:    In order to compute Ext^k(M,N) use the command Ext(k,syz(M),syz(N));
         or: list P=mres(M,2); list Q=mres(N,2); Ext(k,P[2],Q[2]);
EXAMPLE: example Ext;   shows examples
"
{
//---------- initialisation ---------------------------------------------------
  int k,max,ii,l,row,col;
  module A,B,C,D,M1,M2,N1,ker,imag,extMN,extMN0;
  matrix kb;
  list L1,L2,L3,L,resM,K;
  ideal  test1;
  intmat Be;
  int s = size(v);
  intvec v1 = sort(v)[1];
  max = v1[s];                      // the maximum integer occuring in intvec v
  int p = printlevel-voice+3;       // p=printlevel+1 (default: p=1)
//---------- test: coker(N)=basering, coker(N)=0 ? ----------------------------
  if( max<0 ) { dbprint(p,"// Ext^i=0 for i<0!"); return([1]); }
  N1 = std(N);
  if( size(N1)==0 )      //coker(N)=basering, in this case proc Ext_R is faster
  { printlevel=printlevel+1;
    if( size(#)==0 )
    { def E = Ext_R(v,M);
      printlevel=printlevel-1;
      return(E);
    }
    else
    { def E = Ext_R(v,M,#[1]);
       printlevel=printlevel-1;
       return(E);
     }
  }
  if( dim(N1)==-1 )                          //coker(N)=0, all Ext-groups are 0
  { dbprint(p-1,"2nd module presents 0, hence Ext^k=0, for all k");
    for( ii=1; ii<=s; ii++ )
    { k=v[ii];
      extMN    = gen(1);
      extMN0   = std(extMN);
      L1[ii] = extMN;
      L2[ii] = extMN0;
      L3[ii] = matrix(kbase(extMN0));
      if( p>0 ) { "// degree of Ext^"+string(k)+":"; degree(extMN0);""; }
    }
  }
  else
  {
    if( size(N1) < size(N) ) { N=N1;}
    row = nrows(N);
//---------- resolution of M -------------------------------------------------
    resM = mres(M,max+1);
    for( ii=1; ii<=s; ii++ )
    { k=v[ii];
      if( k<0  )                                    // Ext^k is 0 for negative k
      { dbprint(p-1,"// Ext^k=0 for k<0!");
        extMN    = gen(1);
        extMN0   = std(extMN);
        L1[ii] = extMN;
        L2[ii] = extMN0;
        L3[ii] = matrix(kbase(extMN0));
        if( p>0 ) { "// degree of Ext^"+string(k)+":"; degree(extMN0);""; }
      }
      else
      { M2 = resM[k+1];
        if( k==0 ) { M1=0*gen(nrows(M2)); }
        else { M1 = resM[k]; }
        col = ncols(M1);
        D = kohom(N,col);
//---------- computing homology ----------------------------------------------
        imag  = kontrahom(M1,row);
        A     = kontrahom(M2,row);
        B     = kohom(N,ncols(M2));
        ker   = modulo(A,B);
        imag  = imag,D;
        extMN = modulo(ker,imag);
        dbprint(p-1,"// Computing Ext^"+string(k)+" (help Ext; gives an explanation):",
         "// Let 0<--coker(M)<--F0<--F1<--F2<--... be a resolution of coker(M),",
         "// and 0<--coker(N)<--G0<--G1 a presentation of coker(N),",
         "// then Hom(F"+string(k)+",G0)-->Hom(F"+string(k+1)+",G0) is given by:",A,
         "// and Hom(F"+string(k-1)+",G0) + Hom(F"+string(k)+",G1)-->Hom(F"+string(k)+",G0) is given by:",imag,"");
        extMN0 = std(extMN);
        if( p>0 ) { "// degree of Ext^"+string(k)+":"; degree(extMN0);""; }
//---------- more information -------------------------------------------------
        if( size(#)>0 )
        { if( vdim(extMN0) >= 0 )
          { kb = kbase(extMN0);
            if ( size(ker)!=0) { kb = matrix(ker)*kb; }
            dbpri(p-1,"// columns of matrix are kbase of Ext^"+
                     string(k)+" in Hom(F"+string(k)+",G0)",kb,"");
            if( p>0 )
            { for (l=1;l<=ncols(kb);l=l+1)
              {
                "// element",l,"of kbase of Ext^"+string(k)+" in Hom(F"+string(k)+",G0)";
                "// as matrix: F"+string(k)+"-->G0";
                print(matrix(ideal(kb[l]),row,col));
              }
              "";
            }
            L3[ii] = matrix(kb);
          }
          L2[ii] = extMN0;
        }
        L1[ii] = extMN;
      }
    }
  }
  if( size(#) )
  {  if( s>1 ) { L = L1,L2,L3; return(L); }
     else { L = extMN,extMN0,matrix(kb); return(L); }
  }
  else
  {  if( s>1 ) { return(L1); }
     else { return(extMN); }
  }
}
example
{"EXAMPLE:";   echo=2;
  int p      = printlevel;
  printlevel = 1;
  ring r     = 0,(x,y),dp;
  ideal i    = x2-y3;
  ideal j    = x2-y5;
  list E     = Ext(0..2,i,j);    // Ext^k(r/i,r/j) for k=0,1,2 over r
  qring R    = std(i);
  ideal j    = fetch(r,j);
  module M   = [-x,y],[-y2,x];
  printlevel = 2;
  module E1  = Ext(1,M,j);       // Ext^1(R^2/M,R/j) over R=r/i
  list l     = Ext(4,M,M,1);     // Ext^4(R^2/M,R^2/M) over R=r/i
  printlevel = p;
}
////////////////////////////////////////////////////////////////////////////////

proc Hom (module M, module N, list #)
"USAGE:   Hom(M,N,[any]);  M,N=modules
COMPUTE: A presentation of Hom(M',N'), M'=coker(M), N'=coker(N) as follows:
         Let ...-->F1 --M-> F0-->M'-->0 and ...-->G1 --N-> G0-->N'-->0  be
         presentations of M' and N'. Consider
@format
                                         0               0
                                         |^              |^
              0 --> Hom(M',N') ----> Hom(F0,N') ----> Hom(F1,N')
                                         |^              |^
         (A:  induced by M)          Hom(F0,G0) --A-> Hom(F1,G0)
                                         |^              |^
         (B,C:induced by N)              |C              |B
                                     Hom(F0,G1) ----> Hom(F1,G1)

         Let D=modulo(A,B) and Hom=modulo(D,C), then we have exact sequences
              R^p  --D-> Hom(F0,G0) --A-> Hom(F1,G0)/im(B)
              R^q -Hom-> R^p --D-> Hom(F0,G0)/im(C) --A-> Hom(F1,G0)/im(B).
         Hence Hom presents Hom(M',N')
@end format
RETURN:  Hom, of type module, presentation of Hom(M',N') or,
         in case of 3 arguments, a list:
           [1] = Hom
           [2] = SB of Hom
           [3] = kbase of coker(Hom) (if finite dimensional), represented by
                 elements in Hom(F0,G0) via mapping D
DISPLAY: printlevel >=0: degree of Hom  (default)
         printlevel >=1: D and C and kbase of coker(Hom) in Hom(F0,G0)
         printlevel >=2: elements of kbase of coker(Hom) as matrix :F0-->G0
NOTE:    DISPLAY is as described only for a direct call of 'Hom'. Calling 'Hom'
         from another proc has the same effect as decreasing printlevel by 1.
EXAMPLE: example Hom;  shows examples
"
{
//---------- initialisation ---------------------------------------------------
  int l,p;
  matrix kb;
  module A,B,C,D,homMN,homMN0;
  list L;
//---------- computation of Hom -----------------------------------------------
  B = kohom(N,ncols(M));
  A = kontrahom(M,nrows(N));
  C = kohom(N,nrows(M));
  D = modulo(A,B);
  homMN = modulo(D,C);
  homMN0= std(homMN);
  p = printlevel-voice+3;       // p=printlevel+1 (default: p=1)
  if( p>=0 ) { "// degree of Hom:"; degree(homMN0); ""; }
  dbprint(p-1,"// given ...--> F1 --M-> F0 -->M'--> 0 and ...--> G1 --N-> G0 -->N'--> 0,",
         "// show D=ker(Hom(F0,G0) --> Hom(F1,G0)/im(Hom(F1,G1) --> Hom(F1,G0)))",D,
         "// show C=im(Hom(F0,G1) --> Hom(F0,G0))",C,"");
//---------- extra output if size(#)>0 ----------------------------------------
  if( size(#)>0 )
  {  if( vdim(homMN0)>0 )
     {  kb = kbase(homMN0);
        kb = matrix(D)*kb;
        if( p>2 )
        {  for (l=1;l<=ncols(kb);l=l+1)
           {
             "// element",l,"of kbase of Hom in Hom(F0,G0) as matrix: F0-->G0:";
             print(matrix(ideal(kb[l]),nrows(N),nrows(M)));
           }
        }
        else
       { dbprint(p-1,"// columns of matrix are kbase of Hom in Hom(F0,G0)",kb); }
        L=homMN,homMN0,kb; return(L);
     }
     L=homMN,homMN0; return(L);
  }
  return(homMN);
}
example
{"EXAMPLE:";  echo = 2;
  int p     = printlevel;
  printlevel= 1;   //in 'example proc' printlevel has to be increased by 1
  ring r    = 0,(x,y),dp;
  ideal i   = x2-y3,xy;
  qring q   = std(i);
  ideal i   = fetch(r,i);
  module M  = [-x,y],[-y2,x],[x3];
  module H  = Hom(M,i);
  print(H);
  printlevel= 2;
  list L    = Hom(M,i,1);"";
  ring s    = 3,(x,y,z),(c,dp);
  ideal i   = jacob(ideal(x2+y5+z4));
  qring rq=std(i);
  matrix M[2][2]=xy,x3,5y,4z,x2;
  matrix N[3][2]=x2,x,y3,3xz,x2z,z;
  print(M);
  print(N);
  list l=Hom(M,N,1);
  printlevel = p;
}
////////////////////////////////////////////////////////////////////////////////

proc homology (matrix A,matrix B,module M,module N,list #)
"USAGE:   homology(A,B,M,N);
COMPUTE: Let M and N be submodules of R^m and R^n presenting M'=R^m/M, N'=R^n/N
         (R=basering) and let A,B matrices inducing maps R^k--A-->R^m--B-->R^n.
         Compute a presentation R^q --H-> R^m of the module
              ker(B)/im(A) := ker(M'/im(A) --B--> N'/im(BM)+im(BA)).
         If B induces a map M'--B-->N' (i.e BM=0) and if R^k--A-->M'--B-->N' is
         a complex (i.e. BA=0) then ker(B)/im(A) is the homology of this complex
RETURN:  module H, a presentation of ker(B)/im(A)
NOTE:    homology returns a free module of rank m if ker(B)=im(A)
EXAMPLE: example homology; shows examples
"
{
  module ker,ima;
  ker = modulo(B,N);
  ima = A,M;
  return(modulo(ker,ima));
}
example
{"EXAMPLE";    echo=2;
  ring r;
  ideal id=maxideal(4);
  qring qr=std(id);
  module N=maxideal(3)*freemodule(2);
  module M=maxideal(2)*freemodule(2);
  module B=[2x,0],[x,y],[z2,y];
  module A=M;
  degree(std(homology(A,B,M,N)));"";
  ring s=0,x,ds;
  qring qs=std(x4);
  module A=[x];module B=A;
  module M=[x3];module N=M;
  homology(A,B,M,N);
}
//////////////////////////////////////////////////////////////////////////////

proc kernel (matrix A,module M,module N)
"USAGE:   kernel(A,M,N);
COMPUTE: Let M and N be submodules of R^m and R^n presenting M'=R^m/M, N'=R^n/N
         (R=basering) and let A:R^m-->R^n a matrix inducing a map A':M'-->N'.
         Compute a presentation K of ker(A') as in the commutative diagram:
@format
                       ker(A') --->  M' --A'--> N'
                           |^        |^         |^
                           |         |          |
                          R^r  ---> R^m --A--> R^n
                           |^        |^         |^
                           |K        |M         |N
                           |         |          |
                          R^s  ---> R^p -----> R^q
@end format
RETURN:  module K, a presentation of ker(A')
EXAMPLE: example kernel; shows examples
"
{
  module M1 = modulo(A,N);
  return(modulo(M1,M));
}
example
{"EXAMPLE";    echo=2;
  ring r;
  module N=[2x,x],[0,y];
  module M=maxideal(1)*freemodule(2);
  matrix A[2][2]=2x,0,x,y,z2,y;
  module K=kernel(A,M,N);
  degree(std(K));
  print(K);
}
////////////////////////////////////////////////////////////////////////////////

proc kohom (matrix M, int j)
"USAGE:   kohom(A,k); A=matrix, k=integer
RETURN:  matrix Hom(R^k,A) i.e. let A be a matrix defining a map: F1 --> F2 of
         free R-modules, the matrix of Hom(R^k,F1)-->Hom(R^k,F2) is computed
EXAMPLE: example kohom; shows an example
"
{
  if (j==1)
  { return(M);}
  if (j>1)
  { return(outer(M,diag(1,j))); }
  else { return(0);}
}
example
{"EXAMPLE:";   echo=2;
  ring r;
  matrix n[2][3]=x,y,5,z,77,33;
  print(kohom(n,3));
}
/////////////////////////////////////////////////////////////////////////////////

proc kontrahom (matrix M, int j)
"USAGE:   kontrahom(A,k); A=matrix, k=integer
RETURN:  matrix Hom(A,R^k), i.e. let A be a matrix defining a map: F1 --> F2 of
         free  R-modules, the matrix of Hom(F2,R^k)-->Hom(F1,R^k) is computed
EXAMPLE: example kontrahom; shows an example
"
{
if (j==1)
  { return(transpose(M));}
  if (j>1)
  { return(transpose(outer(diag(1,j),M)));}
  else { return(0);}
}
example
{"EXAMPLE:";  echo=2;
  ring r;
  matrix n[2][3]=x,y,5,z,77,33;
  print(kontrahom(n,3));
}
///////////////////////////////////////////////////////////////////////////////