source: git/Singular/LIB/sheafcoh.lib @ 91c7251

///////////////////////////////////////////////////////////////////////////
version="version sheafcoh.lib 4.0.0.0 Jun_2013 "; // $Id$
category="Commutative Algebra";
info="
LIBRARY:  sheafcoh.lib   Procedures for Computing Sheaf Cohomology
AUTHORS:  Wolfram Decker, decker@mathematik.uni-kl.de
          Christoph Lossen,  lossen@mathematik.uni-kl.de
          Gerhard Pfister,  pfister@mathematik.uni-kl.de
          Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.uni-kl.de}

PROCEDURES:
 truncate(phi,d);        truncation of coker(phi) at d
 truncateFast(phi,d);    truncation of coker(phi) at d (fast+ experimental)
 CM_regularity(M);       Castelnuovo-Mumford regularity of coker(M)
 sheafCohBGG(M,l,h);     cohomology of sheaf associated to coker(M)
 sheafCohBGG2(M,l,h);    cohomology of sheaf associated to coker(M), experimental version
 sheafCoh(M,l,h);        cohomology of sheaf associated to coker(M)
 dimH(i,M,d);            compute h^i(F(d)), F sheaf associated to coker(M)
 dimGradedPart()

 displayCohom(B,l,h,n);  display intmat as Betti diagram (with zero rows)

KEYWORDS: sheaf cohomology
";

///////////////////////////////////////////////////////////////////////////////
LIB "matrix.lib";
LIB "nctools.lib";
LIB "homolog.lib";


///////////////////////////////////////////////////////////////////////////////
static proc jacobM(matrix M)
{
   int n=nvars(basering);
   matrix B=transpose(diff(M,var(1)));
   int i;
   for(i=2;i<=n;i++)
   {
     B=concat(B,transpose(diff(M,var(i))));
   }
   return(transpose(B));
}


///////////////////////////////////////////////////////////////////////////////
/**
  let M = { w_1, ..., w_k }, k = size(M) == ncols(M), n = nvars(currRing).
  assuming that nrows(M) <= m*n;
  computes transpose(M) * transpose( var(1) I_m | ... | var(n) I_m ) :== transpose(module{f_1, ... f_k}),
  where f_i = \sum_{j=1}^{m} (w_i, v_j) gen(j),  (w_i, v_j) is a `scalar` multiplication.
  that is, if w_i = (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m) then

  (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m)
*  var_1  ... var_1  |  var_2  ...  var_2  | ... |  var_n  ...  var(n)
*  gen_1  ... gen_m  |  gen_1  ...  gen_m  | ... |  gen_1  ...  gen_m
+ =>
  f_i =

   a^1_1 * var(1) * gen(1) + ... + a^1_m * var(1) * gen(m) +
   a^2_1 * var(2) * gen(1) + ... + a^2_m * var(2) * gen(m) +
                             ...
   a^n_1 * var(n) * gen(1) + ... + a^n_m * var(n) * gen(m);

   NOTE: for every f_i we run only ONCE along w_i saving partial sums into a temporary array of polys of size m

*/
static proc TensorModuleMult(int m, module M)
{
  return( system("tensorModuleMult", m, M) ); // trick!

  int n = nvars(basering);
  int k = ncols(M);

  int g, cc, vv;

  poly h;

  module Temp; // = {f_1, ..., f_k }

  intvec exp;
  vector pTempSum, w;

  for( int i = k; i > 0; i-- ) // for every w \in M
  {
    pTempSum[m] = 0;
    w = M[i];

    while(w != 0) // for each term of w...
    {
      exp = leadexp(w);
      g = exp[n+1]; // module component!
      h = w[g];

      w = w - h * gen(g);

      cc = g % m;

      if( cc == 0)
      {
        cc = m;
      }

      vv = 1 + (g - cc) / m;

      pTempSum = pTempSum + h * var(vv) * gen(cc);
    }

    Temp[i] = pTempSum;
  }

  Temp = transpose(Temp);

  return(Temp);
}

///////////////////////////////////////////////////////////////////////////////
static proc max(int i,int j)
{
  if(i>j){return(i);}
  return(j);
}

///////////////////////////////////////////////////////////////////////////////
static proc min(int i,int j)
{
  if(i>j){return(j);}
  return(i);
}


///////////////////////////////////////////////////////////////////////////////
proc truncate(module phi, int d)
"USAGE:   truncate(M,d);  M module, d int
ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute
RETURN:  module
NOTE:    Output is a presentation matrix for the truncation of coker(M)
         at degree d.
EXAMPLE: example truncate; shows an example
KEYWORDS: truncated module
"
{
  if ( typeof(attrib(phi,"isHomog"))=="string" ) {
    if (size(phi)==0) {
      // assign weights 0 to generators of R^n (n=nrows(phi))
      intvec v;
      v[nrows(phi)]=0;
      attrib(phi,"isHomog",v);
    }
    else {
      ERROR("No admissible degree vector assigned");
    }
  }
  else {
   intvec v=attrib(phi,"isHomog");
  }
  int i,m,dummy;
  int s = nrows(phi);
  module L; // TOO BIG!!!
  for (i=1; i<=s; i++) {
    if (d>v[i]) {
      L = L+maxideal(d-v[i])*gen(i);
    }
    else {
      L = L+gen(i);
    }
  }
  L = modulo(L,phi);
  L = minbase(prune(L));
  if (size(L)==0) {return(L);}

  // it only remains to set the degrees for L:
  // ------------------------------------------
  m = v[1];
  for(i=2; i<=size(v); i++) {  if(v[i]<m) { m = v[i]; } }
  dummy = homog(L);
  intvec vv = attrib(L,"isHomog");
  if (d>m) { vv = vv+d; }
  else     { vv = vv+m; }
  attrib(L,"isHomog",vv);
  return(L);
}
example
{"EXAMPLE:";
   echo = 2;
   ring R=0,(x,y,z),dp;
   module M=maxideal(3);
   homog(M);
   // compute presentation matrix for truncated module (R/<x,y,z>^3)_(>=2)
   module M2=truncate(M,2);
   print(M2);
   dimGradedPart(M2,1);
   dimGradedPart(M2,2);
   // this should coincide with:
   dimGradedPart(M,2);
   // shift grading by 1:
   intvec v=1;
   attrib(M,"isHomog",v);
   M2=truncate(M,2);
   print(M2);
   dimGradedPart(M2,3);
}

///////////////////////////////////////////////////////////////////////////////

proc truncateFast(module M, int d)
"USAGE:   truncateFast(M,d);  M module, d int
ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute 'isHomog'
RETURN:  module
NOTE:    Output is a presentation matrix for the truncation of coker(M)
         at d.
         Fast + experimental version. M shoud be a SB!
DISPLAY: If @code{printlevel}>=1, step-by step timings will be printed.
         If @code{printlevel}>=2 we add progress debug messages
         if @code{printlevel}>=3, even all intermediate results...
EXAMPLE: example truncateFast; shows an example
KEYWORDS: truncated module
"
{
//  int PL = printlevel + 1;
  int PL = printlevel - voice + 2;

  dbprint(PL-1, "// truncateFast(M: "+ string(nrows(M)) + " x " + string(ncols(M)) +", " + string(d) + "):");
  dbprint(PL-2, M);

  intvec save = option(get);
  if( PL >= 2 )
  {
    option(prot);
    option(mem);
  }

  int tTruncateBegin=timer;

  if (attrib(M,"isSB")!=1)
  {
    ERROR("M must be a standard basis!");
  }

  dbprint(PL-1, "// M is a SB! ");

  if ( typeof(attrib(M,"isHomog"))=="string" ) {
    if (size(M)==0) {
      // assign weights 0 to generators of R^n (n=nrows(M))
      intvec v;
      v[nrows(M)]=0;
      attrib(M,"isHomog",v);
    }
    else {
      ERROR("No admissible degree vector assigned");
    }
  }
  else {
   intvec v=attrib(M,"isHomog");
  }

  dbprint(PL-1, "// weighting(M): ["+ string(v) + "]");

  int i,m,dummy;
  int s = nrows(M);

   int tKBaseBegin = timer;
    module L = kbase(M, d); // TODO: check whether this is always correct!?!


    dbprint(PL-1, "// L = kbase(M,d): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
    dbprint(PL-2, L);
    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");

   int tModuloBegin = timer;
    L = modulo(L,M);

    dbprint(PL-1, "// L = modulo(L,M): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
    dbprint(PL-2, L);
    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");

   int tPruneBegin = timer;
    L = prune(L);

    dbprint(PL-1, "// L = prune(L): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
    dbprint(PL-2, L);
    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");

   int tPruneEnd = timer;
    L = minbase(L);
   int tMinBaseEnd = timer;

    dbprint(PL-1, "// L = minbase(L): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
    dbprint(PL-2, L);
    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");




  if (size(L)!=0)
  {

    // it only remains to set the degrees for L:
    // ------------------------------------------
    m = v[1];
    for(i=2; i<=size(v); i++) {  if(v[i]<m) { m = v[i]; } }
    dummy = homog(L);
    intvec vv = attrib(L,"isHomog");
    if (d>m) { vv = vv+d; }
    else     { vv = vv+m; }
    attrib(L,"isHomog",vv);
  }

  int tTruncateEnd=timer;

  dbprint(PL-1, "// corrected weighting(L): ["+ string(attrib(L, "isHomog")) + "]");


  if(PL > 0)
  {
  "
  -------------- TIMINGS --------------
  Trunc        Time: ", tTruncateEnd - tTruncateBegin, "
    :: Before .Time: ", tKBaseBegin - tTruncateBegin, "
    :: kBase   Time: ", tModuloBegin - tKBaseBegin, "
    :: Modulo  Time: ", tPruneBegin - tModuloBegin, "
    :: Prune   Time: ", tPruneEnd - tPruneBegin, "
    :: Minbase Time: ", tMinBaseEnd - tPruneEnd, "
    :: After  .Time: ", tTruncateEnd - tMinBaseEnd;
  }

  option(set, save);

  return(L);
}
example
{"EXAMPLE:";
   echo = 2;
   ring R=0,(x,y,z,u,v),dp;
   module M=maxideal(3);
   homog(M);
   // compute presentation matrix for truncated module (R/<x,y,z,u>^3)_(>=2)
   int t=timer;
   module M2t=truncate(M,2);
   t = timer - t;
   "// Simple truncate: ", t;
   t=timer;
   module M2=truncateFast(std(M),2);
   t = timer - t;
   "// Fast truncate: ", t;
   print(M2);
   "// Check: M2t == M2?: ", size(NF(M2, std(M2t))) + size(NF(M2t, std(M2)));

   dimGradedPart(M2,1);
   dimGradedPart(M2,2);
   // this should coincide with:
   dimGradedPart(M,2);
   // shift grading by 1:
   intvec v=1;
   attrib(M,"isHomog",v);
   t=timer;
   M2t=truncate(M,2);
   t = timer - t;
   "// Simple truncate: ", t;
   t=timer;
   M2=truncateFast(std(M),2);
   t = timer - t;
   "// Fast truncate: ", t;
   print(M2);
   "// Check: M2t == M2?: ", size(NF(M2, std(M2t))) + size(NF(M2t, std(M2))); //?
   dimGradedPart(M2,3);
}

///////////////////////////////////////////////////////////////////////////////





proc dimGradedPart(module phi, int d)
"USAGE:   dimGradedPart(M,d);  M module, d int
ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute
RETURN:  int
NOTE:    Output is the vector space dimension of the graded part of degree d
         of coker(M).
EXAMPLE: example dimGradedPart; shows an example
KEYWORDS: graded module, graded piece
"
{
  if ( typeof(attrib(phi,"isHomog"))=="string" ) {
    if (size(phi)==0) {
      // assign weights 0 to generators of R^n (n=nrows(phi))
      intvec v;
      v[nrows(phi)]=0;
    }
    else { ERROR("No admissible degree vector assigned"); }
  }
  else {
    intvec v=attrib(phi,"isHomog");
  }
  int s = nrows(phi);
  int i,m,dummy;
  module L,LL;
  for (i=1; i<=s; i++) {
    if (d>v[i]) {
      L = L+maxideal(d-v[i])*gen(i);
      LL = LL+maxideal(d+1-v[i])*gen(i);
    }
    else {
      L = L+gen(i);
      if (d==v[i]) {
        LL = LL+maxideal(1)*gen(i);
      }
      else {
        LL = LL+gen(i);
      }
    }
  }
  LL=LL,phi;
  L = modulo(L,LL);
  L = std(prune(L));
  if (size(L)==0) {return(0);}
  return(vdim(L));
}
example
{"EXAMPLE:";
   echo = 2;
   ring R=0,(x,y,z),dp;
   module M=maxideal(3);
   // assign compatible weight vector (here: 0)
   homog(M);
   // compute dimension of graded pieces of R/<x,y,z>^3 :
   dimGradedPart(M,0);
   dimGradedPart(M,1);
   dimGradedPart(M,2);
   dimGradedPart(M,3);
   // shift grading:
   attrib(M,"isHomog",intvec(2));
   dimGradedPart(M,2);
}

///////////////////////////////////////////////////////////////////////////////

proc CM_regularity (module M)
"USAGE:   CM_regularity(M);    M module
ASSUME:   @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute
RETURN:  integer, the Castelnuovo-Mumford regularity of coker(M)
NOTE:    procedure calls mres
EXAMPLE: example CM_regularity; shows an example
KEYWORDS: Castelnuovo-Mumford regularity
"
{
  if ( typeof(attrib(M,"isHomog"))=="string" ) {
    if (size(M)==0) {
      // assign weights 0 to generators of R^n (n=nrows(M))
      intvec v;
      v[nrows(M)]=0;
      attrib(M,"isHomog",v);
    }
    else {
      ERROR("No admissible degree vector assigned");
    }
  }

  if( attrib(CM_regularity,"Algorithm") == "minres_res" )
  {
    def L = minres( res(M,0) ); // let's try it out!
  } else
  {
    def L = mres(M,0);
  }

  intmat BeL = betti(L);
  int r = nrows(module(matrix(BeL)));  // last non-zero row
  if (typeof(attrib(BeL,"rowShift"))!="string") {
    int shift = attrib(BeL,"rowShift");
  }
  return(r+shift-1);
}
example
{"EXAMPLE:";
   echo = 2;
   ring R=0,(x,y,z,u),dp;
   resolution T1=mres(maxideal(1),0);
   module M=T1[3];
   intvec v=2,2,2,2,2,2;
   attrib(M,"isHomog",v);
   CM_regularity(M);
}

///////////////////////////////////////////////////////////////////////////////
proc sheafCohBGG(module M,int l,int h)
"USAGE:   sheafCohBGG(M,l,h);    M module, l,h int
ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute, @code{h>=l}, and the basering has @code{n+1}
         variables.
RETURN:  intmat, cohomology of twists of the coherent sheaf F on P^n
         associated to coker(M). The range of twists is determined by @code{l},
         @code{h}.
DISPLAY: The intmat is displayed in a diagram of the following form: @*
  @format
                l            l+1                      h
  ----------------------------------------------------------
      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
           ...............................................
      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
  ----------------------------------------------------------
    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
  @end format
         A @code{'-'} in the diagram refers to a zero entry; a @code{'*'}
         refers to a negative entry (= dimension not yet determined).
         refers to a not computed dimension. @*
NOTE:    This procedure is based on the Bernstein-Gel'fand-Gel'fand
         correspondence and on Tate resolution ( see [Eisenbud, Floystad,
         Schreyer: Sheaf cohomology and free resolutions over exterior
         algebras, Trans AMS 355 (2003)] ).@*
         @code{sheafCohBGG(M,l,h)} does not compute all values in the above
         table. To determine all values of @code{h^i(F(d))}, @code{d=l..h},
         use @code{sheafCohBGG(M,l-n,h+n)}.
SEE ALSO: sheafCoh, dimH
EXAMPLE: example sheafCohBGG; shows an example
"
{
  int i,j,k,row,col;
  if( typeof(attrib(M,"isHomog"))!="intvec" )
  {
     if (size(M)==0) { attrib(M,"isHomog",0); }
     else { ERROR("No admissible degree vector assigned"); }
  }
  int n=nvars(basering)-1;
  int ell=l+n;
  def R=basering;
  int reg = CM_regularity(M);
  int bound=max(reg+1,h-1);
  module MT=truncate(M,bound);
  int m=nrows(MT);
  MT=transpose(jacobM(MT));
  MT=syz(MT);
  matrix ML[n+1][1]=maxideal(1);
  matrix S=transpose(outer(ML,unitmat(m)));
  matrix SS=transpose(S*MT);
  //--- to the exterior algebra
  def AR = Exterior();
  setring AR;
  intvec saveopt=option(get);
  option(redSB);
  option(redTail);
  module EM=imap(R,SS);
  intvec w;
  //--- here we are with our matrix
  int bound1=max(1,bound-ell+1);
  for (i=1; i<=nrows(EM); i++)
  {
     w[i]=-bound-1;
  }
  attrib(EM,"isHomog",w);
  resolution RE=mres(EM,bound1);
  intmat Betti=betti(RE);
  k=ncols(Betti);
  row=nrows(Betti);
  int shift=attrib(Betti,"rowShift")+(k+ell-1);
  intmat newBetti[n+1][h-l+1];
  for (j=1; j<=row; j++)
  {
    for (i=l; i<=h; i++)
    {
      if ((k+1-j-i+ell-shift>0) and (j+i-ell+shift>=1))
      {
        newBetti[n+2-shift-j,i-l+1]=Betti[j,k+1-j-i+ell-shift];
      }
      else { newBetti[n+2-shift-j,i-l+1]=-1; }
    }
  }
  for (j=2; j<=n+1; j++)
  {
    for (i=1; i<j; i++)
    {
      newBetti[j,i]=-1;
    }
  }
  int d=k-h+ell-1;
  for (j=1; j<=n; j++)
  {
    for (i=h-l+1; i>=k+j; i--)
    {
      newBetti[j,i]=-1;
    }
  }
  displayCohom(newBetti,l,h,n);
  option(set,saveopt);
  setring R;
  return(newBetti);
}
example
{"EXAMPLE:";
   echo = 2;
   // cohomology of structure sheaf on P^4:
   //-------------------------------------------
   ring r=0,x(1..5),dp;
   module M=0;
   def A=sheafCohBGG(M,-9,4);
   // cohomology of cotangential bundle on P^3:
   //-------------------------------------------
   ring R=0,(x,y,z,u),dp;
   resolution T1=mres(maxideal(1),0);
   module M=T1[3];
   intvec v=2,2,2,2,2,2;
   attrib(M,"isHomog",v);
   def B=sheafCohBGG(M,-8,4);
}


///////////////////////////////////////////////////////////////////////////////
static proc showResult( def R, int l, int h )
{
  int PL = 1; // printlevel - voice + 2;
// int PL = printlevel + 1;

  intmat Betti;
  if(typeof(R)=="resolution")
  {
    Betti = betti(R);
  } else
  {
    if(typeof(R)!="intmat")
    {
      ERROR("Wrong input!!!");
    }

    Betti = R;
  }

  int n=nvars(basering)-1;
  int ell = l + n;

  int k     = ncols(Betti);
  int row   = nrows(Betti);
  int shift = attrib(Betti,"rowShift") + (k + ell - 1);

  int iWTH = h-l+1;

  int d = k - h + ell - 1;

  if( PL > 1 )
  {
    "// l: ", l;
    "// h: ", h;
    "// n: ", n;
    "// ell: ", ell;
    "// k: ", k;
    "// row: ", row;
    "// shift: ", shift;
    "// iWTH: ", iWTH;
    "// d: ", d;
  }

  intmat newBetti[ n + 1 ][ iWTH ];
  int i, j;

  for (j=1; j<=row; j++) {
    for (i=l; i<=h; i++) {
      if( (n+2-shift-j)>0 ) {

        if (  (k+1-j-i+ell-shift>0) and (j+i-ell+shift>=1)) {
          newBetti[n+2-shift-j,i-l+1]=Betti[j,k+1-j-i+ell-shift];
        }
        else { newBetti[n+2-shift-j,i-l+1]=-1; }

      }
    }
  }

  for (j=2; j<=n+1; j++) {
    for (i=1; i<min(j, iWTH); i++) {
      newBetti[j,i]= -1;
    }
  }

  for (j=1; j<=n; j++) {
    for (i=iWTH; i>=k+j; i--) {
      newBetti[j,i]=0; // -1;
    }
  }

  if( PL > 0 )
  {
    "Cohomology table:";
    displayCohom(newBetti, l, h, n);
  }

  return(newBetti);
}
///////////////////////////////////////////////////////////////////////////////



///////////////////////////////////////////////////////////////////////////////
proc sheafCohBGG2(module M,int l,int h)
"USAGE:   sheafCohBGG2(M,l,h);    M module, l,h int
ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute, @code{h>=l}, and the basering has @code{n+1}
         variables.
RETURN:  intmat, cohomology of twists of the coherent sheaf F on P^n
         associated to coker(M). The range of twists is determined by @code{l},
         @code{h}.
DISPLAY: The intmat is displayed in a diagram of the following form: @*
  @format
                l            l+1                      h
  ----------------------------------------------------------
      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
           ...............................................
      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
  ----------------------------------------------------------
    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
  @end format
         A @code{'-'} in the diagram refers to a zero entry; a @code{'*'}
         refers to a negative entry (= dimension not yet determined).
         refers to a not computed dimension. @*
         If @code{printlevel}>=1, step-by step timings will be printed.
         If @code{printlevel}>=2 we add progress debug messages
         if @code{printlevel}>=3, even all intermediate results...
NOTE:    This procedure is based on the Bernstein-Gel'fand-Gel'fand
         correspondence and on Tate resolution ( see [Eisenbud, Floystad,
         Schreyer: Sheaf cohomology and free resolutions over exterior
         algebras, Trans AMS 355 (2003)] ).@*
         @code{sheafCohBGG(M,l,h)} does not compute all values in the above
         table. To determine all values of @code{h^i(F(d))}, @code{d=l..h},
         use @code{sheafCohBGG(M,l-n,h+n)}.
         Experimental version. Should require less memory.
SEE ALSO: sheafCohBGG
EXAMPLE: example sheafCohBGG2; shows an example
"
{
  int PL = printlevel - voice + 2;
//  int PL = printlevel;

  dbprint(PL-1, "// sheafCohBGG2(M: "+ string(nrows(M)) + " x " + string(ncols(M)) +", " + string(l) + ", " + string(h) + "):");
  dbprint(PL-2, M);

  intvec save = option(get);

  if( PL >= 2 )
  {
    option(prot);
    option(mem);
  }

  def isCoker = attrib(M, "isCoker");
  if( typeof(isCoker) == "int" )
  {
    if( isCoker > 0 )
    {
      dbprint(PL-1, "We are going to assume that M is given by coker matrix (that is, M is not a submodule presentation!)");
    }
  }

  int i,j,k,row,col;

  if( typeof(attrib(M,"isHomog"))!="intvec" )
  {
     if (size(M)==0) { attrib(M,"isHomog",0); }
     else { ERROR("No admissible degree vector assigned"); }
  }

  dbprint(PL-1, "// weighting(M): ["+ string(attrib(M, "isHomog")) + "]");

  option(redSB); option(redTail);

  def R=basering;

  int n = nvars(R) - 1;
  int ell = l + n;


/////////////////////////////////////////////////////////////////////////////
// computations

  int tBegin=timer;
  int reg   = CM_regularity(M);
  int tCMEnd = timer;

  dbprint(PL-1, "// CM_reg(M): "+ string(reg));

  int bound = max(reg + 1, h - 1);

  dbprint(PL-1, "// bound: "+ string(bound));

  ///////////////////////////////////////////////////////////////
  int tSTDBegin=timer;
  M = std(M); // for kbase! // NOTE: this should be after CM_regularity, since otherwise CM_regularity computes JUST TOOOOOOO LONG sometimes (see Reg_Hard examples!)
  int tSTDEnd = timer;

  dbprint(PL-1, "// M = std(M: "+string(nrows(M)) + " x " + string(ncols(M))  + ")");
  dbprint(PL-2, M);
  dbprint(PL-1, "// weighting(M): ["+ string(attrib(M, "isHomog")) + "]");


  printlevel = printlevel + 1;
  int tTruncateBegin=timer;
  module MT = truncateFast(M, bound);
  int tTruncateEnd=timer;
  printlevel = printlevel - 1;

  dbprint(PL-1, "// MT = truncateFast(M: "+string(nrows(MT)) + " x " + string(ncols(MT)) +", " + string(bound) + ")");
  dbprint(PL-2, MT);
  dbprint(PL-1, "// weighting(MT): ["+ string(attrib(MT, "isHomog")) + "]");

  int m=nrows(MT);

  ///////////////////////////////////////////////////////////////
  int tTransposeJacobBegin=timer;
  MT = jacob(MT); // ! :(
  int tTransposeJacobEnd=timer;

  dbprint(PL-1, "// MT = jacob(MT: "+string(nrows(MT)) + " x " + string(ncols(MT))  + ")");
  dbprint(PL-2, MT);
  dbprint(PL-1, "// weighting(MT): ["+ string(attrib(MT, "isHomog")) + "]");

  int tSyzBegin=timer;
  MT = syz(MT);
  int tSyzEnd=timer;

  dbprint(PL-1, "// MT = syz(MT: "+string(nrows(MT)) + " x " + string(ncols(MT))  + ")");
  dbprint(PL-2, MT);
  dbprint(PL-1, "// weighting(MT): ["+ string(attrib(MT, "isHomog")) + "]");

  int tMatrixOppBegin=timer;
  module SS = TensorModuleMult(m, MT);
  int tMatrixOppEnd=timer;

  dbprint(PL-1, "// SS = TensorModuleMult("+ string(m)+ ", MT: "+string(nrows(MT)) + " x " + string(ncols(MT))  + ")");
  dbprint(PL-2, SS);
  dbprint(PL-1, "// weighting(SS): ["+ string(attrib(SS, "isHomog")) + "]");

  //--- to the exterior algebra
  def AR = Exterior(); setring AR;

  dbprint(PL-1, "// Test: var(1) * var(1): "+ string(var(1) * var(1)));

  int maxbound = max(1, bound - ell + 1);
//  int maxbound = max(1, bound - l + 1); // As In M2!!!

  dbprint(PL-1, "// maxbound: "+ string(maxbound));

  //--- here we are with our matrix
  module EM=imap(R,SS);
  intvec w;
  for (i=1; i<=nrows(EM); i++)
  {
     w[i]=-bound-1;
  }

  attrib(EM,"isHomog",w);

  ///////////////////////////////////////////////////////////////

  dbprint(PL-1, "// EM: "+string(nrows(EM)) + " x " + string(ncols(EM))  + ")");
  dbprint(PL-2, EM);
  dbprint(PL-1, "// weighting(EM): ["+ string(attrib(EM, "isHomog")) + "]");

  int tResulutionBegin=timer;
  resolution RE = nres(EM, maxbound); // TODO: Plural computes one too many syzygies...?!
  int tMinResBegin=timer;
  RE = minres(RE);
  int tBettiBegin=timer;
  intmat Betti = betti(RE); // betti(RE, 1);?
  int tResulutionEnd=timer;

  int tEnd = tResulutionEnd;

  if( PL > 0 )
  {
    //    list L = RE; // TODO: size(L/RE) is wrong!
    "
        ----      RESULTS  ----
        Tate Resolution:
    ";
    RE;
    "Betti numbers for Tate resolution (diagonal cohomology table):";
    print(Betti, "betti"); // Diagonal form!
  }

//  printlevel = printlevel + 1;
  Betti = showResult(Betti, l, h ); // Show usual form of cohomology table
//  printlevel = printlevel - 1;

  if(PL > 0)
  {
  "
      ----      TIMINGS     -------
      Trunc      Time: ", tTruncateEnd - tTruncateBegin, "
      Reg        Time: ", tCMEnd - tBegin, "
      kStd       Time: ", tSTDEnd - tSTDBegin, "
      Jacob      Time: ", tTransposeJacobEnd - tTransposeJacobBegin, "
      Syz        Time: ", tSyzEnd - tSyzBegin, "
      Mat        Time: ", tMatrixOppEnd - tMatrixOppBegin, "
      ------------------------------
      Res        Time: ", tResulutionEnd - tResulutionBegin, "
      :: NRes    Time: ", tMinResBegin - tResulutionBegin, "
      :: MinRes .Time: ", tBettiBegin - tMinResBegin, "
      :: Betti  .Time: ", tResulutionEnd - tBettiBegin, "
      ---------------------------------------------------------
      Total Time: ", tEnd - tBegin, "
      ---------------------------------------------------------
      ";
  }

  setring R;

  option(set, save);

  return(Betti);
}
example
{"EXAMPLE:";
   echo = 2;
   int pl = printlevel;
   int l,h, t;

   //-------------------------------------------
   // cohomology of structure sheaf on P^4:
   //-------------------------------------------
   ring r=32001,x(1..5),dp;

   module M= getStructureSheaf(); // OO_P^4

   l = -12; h = 12; // range of twists: l..h

   printlevel = 0;
   //////////////////////////////////////////////
   t = timer;

   def A = sheafCoh(M, l, h); // global Ext method:

   "Time: ", timer - t;
   //////////////////////////////////////////////
   t = timer;

   A = sheafCohBGG(M, l, h);  // BGG method (without optimization):

   "Time: ", timer - t;
   //////////////////////////////////////////////
   t = timer;

   A = sheafCohBGG2(M, l, h); // BGG method (with optimization)

   "Time: ", timer - t;
   //////////////////////////////////////////////
   printlevel = pl;

   kill A, r;

   //-------------------------------------------
   // cohomology of cotangential bundle on P^3:
   //-------------------------------------------
   ring R=32001,(x,y,z,u),dp;

   module M = getCotangentialBundle();

   l = -12; h = 11; // range of twists: l..h


   //////////////////////////////////////////////
   printlevel = 0;
   t = timer;

   def B = sheafCoh(M, l, h); // global Ext method:

   "Time: ", timer - t;
   //////////////////////////////////////////////
   t = timer;

   B = sheafCohBGG(M, l, h);  // BGG method (without optimization):

   "Time: ", timer - t;
   //////////////////////////////////////////////
   t = timer;

   B = sheafCohBGG2(M, l, h); // BGG method (with optimization)

   "Time: ", timer - t;
   //////////////////////////////////////////////
   printlevel = pl;
}


///////////////////////////////////////////////////////////////////////////////

proc dimH(int i,module M,int d)
"USAGE:   dimH(i,M,d);    M module, i,d int
ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute, @code{h>=l}, and the basering @code{S} has
         @code{n+1} variables.
RETURN:  int,  vector space dimension of @math{H^i(F(d))} for F the coherent
         sheaf on P^n associated to coker(M).
NOTE:    The procedure is based on local duality as described in [Eisenbud:
         Computing cohomology. In Vasconcelos: Computational methods in
         commutative algebra and algebraic geometry. Springer (1998)].
SEE ALSO: sheafCoh, sheafCohBGG
EXAMPLE: example dimH; shows an example
"
{
  if( typeof(attrib(M,"isHomog"))=="string" )
  {
    if (size(M)==0)
    {
      // assign weights 0 to generators of R^n (n=nrows(M))
      intvec v;
      v[nrows(M)]=0;
      attrib(M,"isHomog",v);
    }
    else
    {
      ERROR("No admissible degree vector assigned");
    }
  }
  int Result;
  int n=nvars(basering)-1;
  if ((i>0) and (i<=n)) {
    list L=Ext_R(n-i,M,1)[2];
    def N=L[1];
    return(dimGradedPart(N,-n-1-d));
  }
  else
  {
    if (i==0)
    {
      list L=Ext_R(intvec(n+1,n+2),M,1)[2];
      def N0=L[2];
      def N1=L[1];
      Result=dimGradedPart(M,d) - dimGradedPart(N0,-n-1-d)
                                - dimGradedPart(N1,-n-1-d);
      return(Result);
    }
    else {
      return(0);
    }
  }
}
example
{"EXAMPLE:";
   echo = 2;
   ring R=0,(x,y,z,u),dp;
   resolution T1=mres(maxideal(1),0);
   module M=T1[3];
   intvec v=2,2,2,2,2,2;
   attrib(M,"isHomog",v);
   dimH(0,M,2);
   dimH(1,M,0);
   dimH(2,M,1);
   dimH(3,M,-5);
}


///////////////////////////////////////////////////////////////////////////////

proc sheafCoh(module M,int l,int h,list #)
"USAGE:   sheafCoh(M,l,h);    M module, l,h int
ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
         vector as an attribute, @code{h>=l}. The basering @code{S} has
         @code{n+1} variables.
RETURN:  intmat, cohomology of twists of the coherent sheaf F on P^n
         associated to coker(M). The range of twists is determined by @code{l},
         @code{h}.
DISPLAY: The intmat is displayed in a diagram of the following form: @*
  @format
                l            l+1                      h
  ----------------------------------------------------------
      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
           ...............................................
      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
  ----------------------------------------------------------
    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
  @end format
         A @code{'-'} in the diagram refers to a zero entry.
NOTE:    The procedure is based on local duality as described in [Eisenbud:
         Computing cohomology. In Vasconcelos: Computational methods in
         commutative algebra and algebraic geometry. Springer (1998)].@*
         By default, the procedure uses @code{mres} to compute the Ext
         modules. If called with the additional parameter @code{\"sres\"},
         the @code{sres} command is used instead.
SEE ALSO: dimH, sheafCohBGG
EXAMPLE: example sheafCoh; shows an example
"
{
  int use_sres;
  if( typeof(attrib(M,"isHomog"))!="intvec" )
  {
     if (size(M)==0) { attrib(M,"isHomog",0); }
     else { ERROR("No admissible degree vector assigned"); }
  }
  if (size(#)>0)
  {
    if (#[1]=="sres") { use_sres=1; }
  }
  int i,j;
  module N,N0,N1;
  int n=nvars(basering)-1;
  intvec v=0..n+1;
  int col=h-l+1;
  intmat newBetti[n+1][col];
  if (use_sres) { list L=Ext_R(v,M,1,"sres")[2]; }
  else          { list L=Ext_R(v,M,1)[2]; }
  for (i=l; i<=h; i++)
  {
    N0=L[n+2];
    N1=L[n+1];
    newBetti[n+1,i-l+1]=dimGradedPart(M,i) - dimGradedPart(N0,-n-1-i)
                             - dimGradedPart(N0,-n-1-i);
  }
  for (j=1; j<=n; j++)
  {
     N=L[j];
     attrib(N,"isSB",1);
     if (dim(N)>=0) {
       for (i=l; i<=h; i++)
       {
         newBetti[j,i-l+1]=dimGradedPart(N,-n-1-i);
       }
     }
  }
  displayCohom(newBetti,l,h,n);
  return(newBetti);
}
example
{"EXAMPLE:";
   echo = 2;
   //
   // cohomology of structure sheaf on P^4:
   //-------------------------------------------
   ring r=0,x(1..5),dp;
   module M=0;
   def A=sheafCoh(0,-7,2);
   //
   // cohomology of cotangential bundle on P^3:
   //-------------------------------------------
   ring R=0,(x,y,z,u),dp;
   resolution T1=mres(maxideal(1),0);
   module M=T1[3];
   intvec v=2,2,2,2,2,2;
   attrib(M,"isHomog",v);
   def B=sheafCoh(M,-6,2);
}

///////////////////////////////////////////////////////////////////////////////
proc displayCohom (intmat data, int l, int h, int n)
"USAGE:   displayCohom(data,l,h,n);  data intmat, l,h,n int
ASSUME:  @code{h>=l}, @code{data} is the return value of
         @code{sheafCoh(M,l,h)} or of @code{sheafCohBGG(M,l,h)}, and the
         basering has @code{n+1} variables.
RETURN:  none
NOTE:    The intmat is displayed in a diagram of the following form: @*
  @format
                l            l+1                      h
  ----------------------------------------------------------
      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
           ...............................................
      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
  ----------------------------------------------------------
    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
  @end format
         where @code{F} refers to the associated sheaf of @code{M} on P^n.@*
         A @code{'-'} in the diagram refers to a zero entry,  a @code{'*'}
         refers to a negative entry (= dimension not yet determined).
"
{
  int i,j,k,dat,maxL;
  intvec notSumCol;
  notSumCol[h-l+1]=0;
  string s;
  maxL=4;
  for (i=1;i<=nrows(data);i++)
  {
    for (j=1;j<=ncols(data);j++)
    {
      if (size(string(data[i,j]))>=maxL-1)
      {
        maxL=size(string(data[i,j]))+2;
      }
    }
  }
  string Row="    ";
  string Row1="----";
  for (i=l; i<=h; i++) {
    for (j=1; j<=maxL-size(string(i)); j++)
    {
      Row=Row+" ";
    }
    Row=Row+string(i);
    for (j=1; j<=maxL; j++) { Row1 = Row1+"-"; }
  }
  print(Row);
  print(Row1);
  for (j=1; j<=n+1; j++)
  {
    s = string(n+1-j);
    Row = "";
    for(k=1; k<4-size(s); k++) { Row = Row+" "; }
    Row = Row + s+":";
    for (i=0; i<=h-l; i++)
    {
      dat = data[j,i+1];
      if (dat>0) { s = string(dat); }
      else
      {
        if (dat==0) { s="-"; }
        else        { s="*"; notSumCol[i+1]=1; }
      }
      for(k=1; k<=maxL-size(s); k++) { Row = Row+" "; }
      Row = Row + s;
    }
    print(Row);
  }
  print(Row1);
  Row="chi:";
  for (i=0; i<=h-l; i++)
  {
    dat = 0;
    if (notSumCol[i+1]==0)
    {
      for (j=0; j<=n; j++) { dat = dat + (-1)^j * data[n+1-j,i+1]; }
      s = string(dat);
    }
    else { s="*"; }
    for (k=1; k<=maxL-size(s); k++) { Row = Row+" "; }
    Row = Row + s;
  }
  print(Row);
}
///////////////////////////////////////////////////////////////////////////////

proc getStructureSheaf(list #)
{

  if( size(#) == 0 )
  {
    module M = 0;
    intvec v = 0;
    attrib(M,"isHomog",v);
//    homog(M);

    attrib(M, "isCoker", 1);

//     attrib(M);
    return(M);
  }

  if( typeof(#[1]) == "ideal")
  {
    ideal I = #[1];

    if( size(#) == 2 )
    {
      if( typeof(#[2]) == "int" )
      {
        if( #[2] != 0 )
        {
          qring @@@@QQ = std(I);

          module M = getStructureSheaf();

          export M;

//          keepring @@@@QQ; // This is a bad idea... :(?
          return (@@@@QQ);
        }
      }
    }

/*
    // This seems to be wrong!!!
    module M = I * gen(1);
    homog(M);

    M = modulo(gen(1), module(I * gen(1))); // basering^1 / I

    homog(M);

    attrib(M, "isCoker", 1);

    attrib(M);
    return(M);
*/
  }

  ERROR("Wrong argument");

}
example
{"EXAMPLE:";
   echo = 2; int pl = printlevel;
   printlevel = voice;


   ////////////////////////////////////////////////////////////////////////////////
   ring r;
   module M = getStructureSheaf();
   "Basering: ";
   basering;
   "Module: ", string(M), ", grading is given by weights: ", attrib(M, "isHomog");

   def A=sheafCohBGG2(M,-9,9);
   print(A);

   ////////////////////////////////////////////////////////////////////////////////
   setring r;
   module M = getStructureSheaf(ideal(var(1)), 0);

   "Basering: ";
   basering;
   "Module: ", string(M), ", grading is given by weights: ", attrib(M, "isHomog");

   def A=sheafCohBGG2(M,-9,9);
   print(A);

   ////////////////////////////////////////////////////////////////////////////////
   setring r;
   def Q = getStructureSheaf(ideal(var(1)), 1); // returns a new ring!
   setring Q; // M was exported in the new ring!

   "Basering: ";
   basering;
   "Module: ", string(M), ", grading is given by weights: ", attrib(M, "isHomog");

   def A=sheafCohBGG2(M,-9,9);
   print(A);

   printlevel = pl;
}


proc getCotangentialBundle()
{
  resolution T1=mres(maxideal(1),3);
  module M=T1[3];
//  attrib(M,"isHomog");
//  homog(M);
  attrib(M, "isCoker", 1);
  // attrib(M);
  return (M);
}

proc getIdealSheafPullback(ideal I, ideal pi)
{
  def save = basering;
  map P = save, pi;
  return( P(I) );
}

// TODO: set attributes!


proc getIdealSheaf(ideal I)
{
  int i = homog(I);
  resolution FI = mres(I,2); // Syz + grading...
  module M = FI[2];
  attrib(M, "isCoker", 1);
  //  attrib(M);
  return(M);
}





/*
Examples:
---------
 LIB "sheafcoh.lib";

 ring S = 32003, x(0..4), dp;
 module MI=maxideal(1);
 attrib(MI,"isHomog",intvec(-1));
 resolution kos = nres(MI,0);
 print(betti(kos),"betti");
 LIB "random.lib";
 matrix alpha0 = random(32002,10,3);
 module pres = module(alpha0)+kos[3];
 attrib(pres,"isHomog",intvec(1,1,1,1,1,1,1,1,1,1));
 resolution fcokernel = mres(pres,0);
 print(betti(fcokernel),"betti");
 module dir = transpose(pres);
 attrib(dir,"isHomog",intvec(-1,-1,-1,-2,-2,-2,
                             -2,-2,-2,-2,-2,-2,-2));
 resolution fdir = mres(dir,2);
 print(betti(fdir),"betti");
 ideal I = groebner(flatten(fdir[2]));
 resolution FI = mres(I,0);
 print(betti(FI),"betti");
 module F=FI[2];
 int t=timer;
 def A1=sheafCoh(F,-8,8);
 timer-t;
 t=timer;
 def A2=sheafCohBGG(F,-8,8);
 timer-t;

 LIB "sheafcoh.lib";
 LIB "random.lib";
 ring S = 32003, x(0..4), dp;
 resolution kos = nres(maxideal(1),0);
 betti(kos);
 matrix kos5 = kos[5];
 matrix tphi = transpose(dsum(kos5,kos5));
 matrix kos3 = kos[3];
 matrix psi = dsum(kos3,kos3);
 matrix beta1 = random(32002,20,2);
 matrix tbeta1tilde = transpose(psi*beta1);
 matrix tbeta0 = lift(tphi,tbeta1tilde);
 matrix kos4 = kos[4];
 matrix tkos4pluskos4 = transpose(dsum(kos4,kos4));
 matrix tgammamin1 = random(32002,20,1);
 matrix tgamma0 = tkos4pluskos4*tgammamin1;
 matrix talpha0 = concat(tbeta0,tgamma0);
 matrix zero[20][1];
 matrix tpsi = transpose(psi);
 matrix tpresg = concat(tpsi,zero);
 matrix pres = module(transpose(talpha0))
                    + module(transpose(tpresg));
 module dir = transpose(pres);
 dir = prune(dir);
 homog(dir);
 intvec deg_dir = attrib(dir,"isHomog");
 attrib(dir,"isHomog",deg_dir-2);        // set degrees
 resolution fdir = mres(prune(dir),2);
 print(betti(fdir),"betti");
 ideal I = groebner(flatten(fdir[2]));
 resolution FI = mres(I,0);

 module F=FI[2];
 def A1=sheafCoh(F,-5,7);
 def A2=sheafCohBGG(F,-5,7);

*/