///////////////////////////////////////////////////////////////////////////////
version="$Id: sheafcoh.lib,v 1.15 2007-10-30 17:17:44 Singular Exp $";
category="Commutative Algebra";
info="
LIBRARY:  sheafcoh.lib   Procedures for Computing Sheaf Cohomology
AUTHORS:  Wolfram Decker, decker@math.uni-sb.de,
@*        Christoph Lossen,  lossen@mathematik.uni-kl.de
@*        Gerhard Pfister,  pfister@mathematik.uni-kl.de

PROCEDURES:
 truncate(phi,d);        truncation of coker(phi) at d
 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)

AUXILIARY PROCEDURES:
 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));
}

// returns transposed Jacobian of a module M
static proc tJacobian(module M)
{
  M = transpose(M);

  int N = nvars(basering);
  int k = ncols(M);
  int r = nrows(M);

  module Result;
  Result[N*k] = 0;

  int i, j;
  int l = 1;

  for( j = 1; j <= N; j++ ) // for all variables
  {
    for( i = 1; i <= k; i++ ) // for every v \in transpose(M)
    {
      Result[l] = diff(M[i], var(j));
      l++;
    }
  }

  return(Result);
}


/**
  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)
{
  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);
}

///////////////////////////////////////////////////////////////////////////////
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 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;
  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 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");
    }
  }
  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;
  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);
  setring R;
  return(newBetti);
  option(noredSB);
  option(noredTail);
}
example
{"EXAMPLE:";
   echo = 2;
   // cohomology of structure sheaf on P^4:
   //-------------------------------------------
   ring r=0,x(1..5),dp;
   module M=0;
   def A=sheafCohBGG(0,-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);
}

///////////////////////////////////////////////////////////////////////////////
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. @*
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 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"); }
  }
  intvec ivOptionsSave = option(get);
  option(redSB); option(redTail);

  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 = tJacobian(MT); // transpose(jacobM(MT));
  MT=syz(MT);

  module SS = TensorModuleMult(m, MT);

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

  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 = minres(nres(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);

  setring R;
  option(set, ivOptionsSave);

  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=sheafCohBGG2(0,-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=sheafCohBGG2(M,-8,4);
}


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

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);
}
///////////////////////////////////////////////////////////////////////////////


/*
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);

*/