source: git/Singular/LIB/sheafcoh.lib @ 1b2216

spielwiese
Last change on this file since 1b2216 was 1e1ec4, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Updated LIBs according to master add: new LIBs from master fix: updated LIBs due to minpoly/(de)numerator changes fix: -> $Id$ fix: Fixing wrong rebase of SW on master (LIBs)
  • Property mode set to 100644
File size: 37.9 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Commutative Algebra";
4info="
5LIBRARY:  sheafcoh.lib   Procedures for Computing Sheaf Cohomology
6AUTHORS:  Wolfram Decker, decker@mathematik.uni-kl.de
7          Christoph Lossen,  lossen@mathematik.uni-kl.de
8          Gerhard Pfister,  pfister@mathematik.uni-kl.de
9          Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.uni-kl.de}
10
11PROCEDURES:
12 truncate(phi,d);        truncation of coker(phi) at d
13 truncateFast(phi,d);    truncation of coker(phi) at d (fast+ experimental)
14 CM_regularity(M);       Castelnuovo-Mumford regularity of coker(M)
15 sheafCohBGG(M,l,h);     cohomology of sheaf associated to coker(M)
16 sheafCohBGG2(M,l,h);    cohomology of sheaf associated to coker(M), experimental version
17 sheafCoh(M,l,h);        cohomology of sheaf associated to coker(M)
18 dimH(i,M,d);            compute h^i(F(d)), F sheaf associated to coker(M)
19 dimGradedPart()
20
21 displayCohom(B,l,h,n);  display intmat as Betti diagram (with zero rows)
22
23KEYWORDS: sheaf cohomology
24";
25
26///////////////////////////////////////////////////////////////////////////////
27LIB "matrix.lib";
28LIB "nctools.lib";
29LIB "homolog.lib";
30
31
32///////////////////////////////////////////////////////////////////////////////
33static proc jacobM(matrix M)
34{
35   int n=nvars(basering);
36   matrix B=transpose(diff(M,var(1)));
37   int i;
38   for(i=2;i<=n;i++)
39   {
40     B=concat(B,transpose(diff(M,var(i))));
41   }
42   return(transpose(B));
43}
44
45
46///////////////////////////////////////////////////////////////////////////////
47/**
48  let M = { w_1, ..., w_k }, k = size(M) == ncols(M), n = nvars(currRing).
49  assuming that nrows(M) <= m*n;
50  computes transpose(M) * transpose( var(1) I_m | ... | var(n) I_m ) :== transpose(module{f_1, ... f_k}),
51  where f_i = \sum_{j=1}^{m} (w_i, v_j) gen(j),  (w_i, v_j) is a `scalar` multiplication.
52  that is, if w_i = (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m) then
53
54  (a^1_1, ... a^1_m) | (a^2_1, ..., a^2_m) | ... | (a^n_1, ..., a^n_m)
55*  var_1  ... var_1  |  var_2  ...  var_2  | ... |  var_n  ...  var(n)
56*  gen_1  ... gen_m  |  gen_1  ...  gen_m  | ... |  gen_1  ...  gen_m
57+ =>
58  f_i =
59
60   a^1_1 * var(1) * gen(1) + ... + a^1_m * var(1) * gen(m) +
61   a^2_1 * var(2) * gen(1) + ... + a^2_m * var(2) * gen(m) +
62                             ...
63   a^n_1 * var(n) * gen(1) + ... + a^n_m * var(n) * gen(m);
64
65   NOTE: for every f_i we run only ONCE along w_i saving partial sums into a temporary array of polys of size m
66
67*/
68static proc TensorModuleMult(int m, module M)
69{
70  return( system("tensorModuleMult", m, M) ); // trick!
71
72  int n = nvars(basering);
73  int k = ncols(M);
74
75  int g, cc, vv;
76
77  poly h;
78
79  module Temp; // = {f_1, ..., f_k }
80
81  intvec exp;
82  vector pTempSum, w;
83
84  for( int i = k; i > 0; i-- ) // for every w \in M
85  {
86    pTempSum[m] = 0;
87    w = M[i];
88
89    while(w != 0) // for each term of w...
90    {
91      exp = leadexp(w);
92      g = exp[n+1]; // module component!
93      h = w[g];
94
95      w = w - h * gen(g);
96
97      cc = g % m;
98
99      if( cc == 0)
100      {
101        cc = m;
102      }
103
104      vv = 1 + (g - cc) / m;
105
106      pTempSum = pTempSum + h * var(vv) * gen(cc);
107    }
108
109    Temp[i] = pTempSum;
110  }
111
112  Temp = transpose(Temp);
113
114  return(Temp);
115}
116
117///////////////////////////////////////////////////////////////////////////////
118static proc max(int i,int j)
119{
120  if(i>j){return(i);}
121  return(j);
122}
123
124///////////////////////////////////////////////////////////////////////////////
125static proc min(int i,int j)
126{
127  if(i>j){return(j);}
128  return(i);
129}
130
131
132///////////////////////////////////////////////////////////////////////////////
133proc truncate(module phi, int d)
134"USAGE:   truncate(M,d);  M module, d int
135ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
136         vector as an attribute
137RETURN:  module
138NOTE:    Output is a presentation matrix for the truncation of coker(M)
139         at degree d.
140EXAMPLE: example truncate; shows an example
141KEYWORDS: truncated module
142"
143{
144  if ( typeof(attrib(phi,"isHomog"))=="string" ) {
145    if (size(phi)==0) {
146      // assign weights 0 to generators of R^n (n=nrows(phi))
147      intvec v;
148      v[nrows(phi)]=0;
149      attrib(phi,"isHomog",v);
150    }
151    else {
152      ERROR("No admissible degree vector assigned");
153    }
154  }
155  else {
156   intvec v=attrib(phi,"isHomog");
157  }
158  int i,m,dummy;
159  int s = nrows(phi);
160  module L; // TOO BIG!!!
161  for (i=1; i<=s; i++) {
162    if (d>v[i]) {
163      L = L+maxideal(d-v[i])*gen(i);
164    }
165    else {
166      L = L+gen(i);
167    }
168  }
169  L = modulo(L,phi);
170  L = minbase(prune(L));
171  if (size(L)==0) {return(L);}
172
173  // it only remains to set the degrees for L:
174  // ------------------------------------------
175  m = v[1];
176  for(i=2; i<=size(v); i++) {  if(v[i]<m) { m = v[i]; } }
177  dummy = homog(L);
178  intvec vv = attrib(L,"isHomog");
179  if (d>m) { vv = vv+d; }
180  else     { vv = vv+m; }
181  attrib(L,"isHomog",vv);
182  return(L);
183}
184example
185{"EXAMPLE:";
186   echo = 2;
187   ring R=0,(x,y,z),dp;
188   module M=maxideal(3);
189   homog(M);
190   // compute presentation matrix for truncated module (R/<x,y,z>^3)_(>=2)
191   module M2=truncate(M,2);
192   print(M2);
193   dimGradedPart(M2,1);
194   dimGradedPart(M2,2);
195   // this should coincide with:
196   dimGradedPart(M,2);
197   // shift grading by 1:
198   intvec v=1;
199   attrib(M,"isHomog",v);
200   M2=truncate(M,2);
201   print(M2);
202   dimGradedPart(M2,3);
203}
204
205///////////////////////////////////////////////////////////////////////////////
206
207proc truncateFast(module M, int d)
208"USAGE:   truncateFast(M,d);  M module, d int
209ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
210         vector as an attribute 'isHomog'
211RETURN:  module
212NOTE:    Output is a presentation matrix for the truncation of coker(M)
213         at d.
214         Fast + experimental version. M shoud be a SB!
215DISPLAY: If @code{printlevel}>=1, step-by step timings will be printed.
216         If @code{printlevel}>=2 we add progress debug messages
217         if @code{printlevel}>=3, even all intermediate results...
218EXAMPLE: example truncateFast; shows an example
219KEYWORDS: truncated module
220"
221{
222//  int PL = printlevel + 1;
223  int PL = printlevel - voice + 2;
224
225  dbprint(PL-1, "// truncateFast(M: "+ string(nrows(M)) + " x " + string(ncols(M)) +", " + string(d) + "):");
226  dbprint(PL-2, M);
227
228  intvec save = option(get);
229  if( PL >= 2 )
230  {
231    option(prot);
232    option(mem);
233  }
234
235  int tTruncateBegin=timer;
236
237  if (attrib(M,"isSB")!=1)
238  {
239    ERROR("M must be a standard basis!");
240  }
241
242  dbprint(PL-1, "// M is a SB! ");
243
244  if ( typeof(attrib(M,"isHomog"))=="string" ) {
245    if (size(M)==0) {
246      // assign weights 0 to generators of R^n (n=nrows(M))
247      intvec v;
248      v[nrows(M)]=0;
249      attrib(M,"isHomog",v);
250    }
251    else {
252      ERROR("No admissible degree vector assigned");
253    }
254  }
255  else {
256   intvec v=attrib(M,"isHomog");
257  }
258
259  dbprint(PL-1, "// weighting(M): ["+ string(v) + "]");
260
261  int i,m,dummy;
262  int s = nrows(M);
263
264   int tKBaseBegin = timer;
265    module L = kbase(M, d); // TODO: check whether this is always correct!?!
266
267
268    dbprint(PL-1, "// L = kbase(M,d): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
269    dbprint(PL-2, L);
270    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");
271
272   int tModuloBegin = timer;
273    L = modulo(L,M);
274
275    dbprint(PL-1, "// L = modulo(L,M): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
276    dbprint(PL-2, L);
277    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");
278
279   int tPruneBegin = timer;
280    L = prune(L);
281
282    dbprint(PL-1, "// L = prune(L): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
283    dbprint(PL-2, L);
284    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");
285
286   int tPruneEnd = timer;
287    L = minbase(L);
288   int tMinBaseEnd = timer;
289
290    dbprint(PL-1, "// L = minbase(L): "+string(nrows(L)) + " x " + string(ncols(L)) +"");
291    dbprint(PL-2, L);
292    dbprint(PL-1, "// weighting(L): ["+ string(attrib(L, "isHomog")) + "]");
293
294
295
296
297  if (size(L)!=0)
298  {
299
300    // it only remains to set the degrees for L:
301    // ------------------------------------------
302    m = v[1];
303    for(i=2; i<=size(v); i++) {  if(v[i]<m) { m = v[i]; } }
304    dummy = homog(L);
305    intvec vv = attrib(L,"isHomog");
306    if (d>m) { vv = vv+d; }
307    else     { vv = vv+m; }
308    attrib(L,"isHomog",vv);
309  }
310
311  int tTruncateEnd=timer;
312
313  dbprint(PL-1, "// corrected weighting(L): ["+ string(attrib(L, "isHomog")) + "]");
314
315
316  if(PL > 0)
317  {
318  "
319  -------------- TIMINGS --------------
320  Trunc        Time: ", tTruncateEnd - tTruncateBegin, "
321    :: Before .Time: ", tKBaseBegin - tTruncateBegin, "
322    :: kBase   Time: ", tModuloBegin - tKBaseBegin, "
323    :: Modulo  Time: ", tPruneBegin - tModuloBegin, "
324    :: Prune   Time: ", tPruneEnd - tPruneBegin, "
325    :: Minbase Time: ", tMinBaseEnd - tPruneEnd, "
326    :: After  .Time: ", tTruncateEnd - tMinBaseEnd;
327  }
328
329  option(set, save);
330
331  return(L);
332}
333example
334{"EXAMPLE:";
335   echo = 2;
336   ring R=0,(x,y,z,u,v),dp;
337   module M=maxideal(3);
338   homog(M);
339   // compute presentation matrix for truncated module (R/<x,y,z,u>^3)_(>=2)
340   int t=timer;
341   module M2t=truncate(M,2);
342   t = timer - t;
343   "// Simple truncate: ", t;
344   t=timer;
345   module M2=truncateFast(std(M),2);
346   t = timer - t;
347   "// Fast truncate: ", t;
348   print(M2);
349   "// Check: M2t == M2?: ", size(NF(M2, std(M2t))) + size(NF(M2t, std(M2)));
350
351   dimGradedPart(M2,1);
352   dimGradedPart(M2,2);
353   // this should coincide with:
354   dimGradedPart(M,2);
355   // shift grading by 1:
356   intvec v=1;
357   attrib(M,"isHomog",v);
358   t=timer;
359   M2t=truncate(M,2);
360   t = timer - t;
361   "// Simple truncate: ", t;
362   t=timer;
363   M2=truncateFast(std(M),2);
364   t = timer - t;
365   "// Fast truncate: ", t;
366   print(M2);
367   "// Check: M2t == M2?: ", size(NF(M2, std(M2t))) + size(NF(M2t, std(M2))); //?
368   dimGradedPart(M2,3);
369}
370
371///////////////////////////////////////////////////////////////////////////////
372
373
374
375
376
377proc dimGradedPart(module phi, int d)
378"USAGE:   dimGradedPart(M,d);  M module, d int
379ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
380         vector as an attribute
381RETURN:  int
382NOTE:    Output is the vector space dimension of the graded part of degree d
383         of coker(M).
384EXAMPLE: example dimGradedPart; shows an example
385KEYWORDS: graded module, graded piece
386"
387{
388  if ( typeof(attrib(phi,"isHomog"))=="string" ) {
389    if (size(phi)==0) {
390      // assign weights 0 to generators of R^n (n=nrows(phi))
391      intvec v;
392      v[nrows(phi)]=0;
393    }
394    else { ERROR("No admissible degree vector assigned"); }
395  }
396  else {
397    intvec v=attrib(phi,"isHomog");
398  }
399  int s = nrows(phi);
400  int i,m,dummy;
401  module L,LL;
402  for (i=1; i<=s; i++) {
403    if (d>v[i]) {
404      L = L+maxideal(d-v[i])*gen(i);
405      LL = LL+maxideal(d+1-v[i])*gen(i);
406    }
407    else {
408      L = L+gen(i);
409      if (d==v[i]) {
410        LL = LL+maxideal(1)*gen(i);
411      }
412      else {
413        LL = LL+gen(i);
414      }
415    }
416  }
417  LL=LL,phi;
418  L = modulo(L,LL);
419  L = std(prune(L));
420  if (size(L)==0) {return(0);}
421  return(vdim(L));
422}
423example
424{"EXAMPLE:";
425   echo = 2;
426   ring R=0,(x,y,z),dp;
427   module M=maxideal(3);
428   // assign compatible weight vector (here: 0)
429   homog(M);
430   // compute dimension of graded pieces of R/<x,y,z>^3 :
431   dimGradedPart(M,0);
432   dimGradedPart(M,1);
433   dimGradedPart(M,2);
434   dimGradedPart(M,3);
435   // shift grading:
436   attrib(M,"isHomog",intvec(2));
437   dimGradedPart(M,2);
438}
439
440///////////////////////////////////////////////////////////////////////////////
441
442proc CM_regularity (module M)
443"USAGE:   CM_regularity(M);    M module
444ASSUME:   @code{M} is graded, and it comes assigned with an admissible degree
445         vector as an attribute
446RETURN:  integer, the Castelnuovo-Mumford regularity of coker(M)
447NOTE:    procedure calls mres
448EXAMPLE: example CM_regularity; shows an example
449KEYWORDS: Castelnuovo-Mumford regularity
450"
451{
452  if ( typeof(attrib(M,"isHomog"))=="string" ) {
453    if (size(M)==0) {
454      // assign weights 0 to generators of R^n (n=nrows(M))
455      intvec v;
456      v[nrows(M)]=0;
457      attrib(M,"isHomog",v);
458    }
459    else {
460      ERROR("No admissible degree vector assigned");
461    }
462  }
463
464  if( attrib(CM_regularity,"Algorithm") == "minres_res" )
465  {
466    def L = minres( res(M,0) ); // let's try it out!
467  } else
468  {
469    def L = mres(M,0);
470  }
471
472  intmat BeL = betti(L);
473  int r = nrows(module(matrix(BeL)));  // last non-zero row
474  if (typeof(attrib(BeL,"rowShift"))!="string") {
475    int shift = attrib(BeL,"rowShift");
476  }
477  return(r+shift-1);
478}
479example
480{"EXAMPLE:";
481   echo = 2;
482   ring R=0,(x,y,z,u),dp;
483   resolution T1=mres(maxideal(1),0);
484   module M=T1[3];
485   intvec v=2,2,2,2,2,2;
486   attrib(M,"isHomog",v);
487   CM_regularity(M);
488}
489
490///////////////////////////////////////////////////////////////////////////////
491proc sheafCohBGG(module M,int l,int h)
492"USAGE:   sheafCohBGG(M,l,h);    M module, l,h int
493ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
494         vector as an attribute, @code{h>=l}, and the basering has @code{n+1}
495         variables.
496RETURN:  intmat, cohomology of twists of the coherent sheaf F on P^n
497         associated to coker(M). The range of twists is determined by @code{l},
498         @code{h}.
499DISPLAY: The intmat is displayed in a diagram of the following form: @*
500  @format
501                l            l+1                      h
502  ----------------------------------------------------------
503      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
504           ...............................................
505      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
506      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
507  ----------------------------------------------------------
508    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
509  @end format
510         A @code{'-'} in the diagram refers to a zero entry; a @code{'*'}
511         refers to a negative entry (= dimension not yet determined).
512         refers to a not computed dimension. @*
513NOTE:    This procedure is based on the Bernstein-Gel'fand-Gel'fand
514         correspondence and on Tate resolution ( see [Eisenbud, Floystad,
515         Schreyer: Sheaf cohomology and free resolutions over exterior
516         algebras, Trans AMS 355 (2003)] ).@*
517         @code{sheafCohBGG(M,l,h)} does not compute all values in the above
518         table. To determine all values of @code{h^i(F(d))}, @code{d=l..h},
519         use @code{sheafCohBGG(M,l-n,h+n)}.
520SEE ALSO: sheafCoh, dimH
521EXAMPLE: example sheafCohBGG; shows an example
522"
523{
524  int i,j,k,row,col;
525  if( typeof(attrib(M,"isHomog"))!="intvec" )
526  {
527     if (size(M)==0) { attrib(M,"isHomog",0); }
528     else { ERROR("No admissible degree vector assigned"); }
529  }
530  int n=nvars(basering)-1;
531  int ell=l+n;
532  def R=basering;
533  int reg = CM_regularity(M);
534  int bound=max(reg+1,h-1);
535  module MT=truncate(M,bound);
536  int m=nrows(MT);
537  MT=transpose(jacobM(MT));
538  MT=syz(MT);
539  matrix ML[n+1][1]=maxideal(1);
540  matrix S=transpose(outer(ML,unitmat(m)));
541  matrix SS=transpose(S*MT);
542  //--- to the exterior algebra
543  def AR = Exterior();
544  setring AR;
545  intvec saveopt=option(get);
546  option(redSB);
547  option(redTail);
548  module EM=imap(R,SS);
549  intvec w;
550  //--- here we are with our matrix
551  int bound1=max(1,bound-ell+1);
552  for (i=1; i<=nrows(EM); i++)
553  {
554     w[i]=-bound-1;
555  }
556  attrib(EM,"isHomog",w);
557  resolution RE=mres(EM,bound1);
558  intmat Betti=betti(RE);
559  k=ncols(Betti);
560  row=nrows(Betti);
561  int shift=attrib(Betti,"rowShift")+(k+ell-1);
562  intmat newBetti[n+1][h-l+1];
563  for (j=1; j<=row; j++)
564  {
565    for (i=l; i<=h; i++)
566    {
567      if ((k+1-j-i+ell-shift>0) and (j+i-ell+shift>=1))
568      {
569        newBetti[n+2-shift-j,i-l+1]=Betti[j,k+1-j-i+ell-shift];
570      }
571      else { newBetti[n+2-shift-j,i-l+1]=-1; }
572    }
573  }
574  for (j=2; j<=n+1; j++)
575  {
576    for (i=1; i<j; i++)
577    {
578      newBetti[j,i]=-1;
579    }
580  }
581  int d=k-h+ell-1;
582  for (j=1; j<=n; j++)
583  {
584    for (i=h-l+1; i>=k+j; i--)
585    {
586      newBetti[j,i]=-1;
587    }
588  }
589  displayCohom(newBetti,l,h,n);
590  option(set,saveopt);
591  setring R;
592  return(newBetti);
593}
594example
595{"EXAMPLE:";
596   echo = 2;
597   // cohomology of structure sheaf on P^4:
598   //-------------------------------------------
599   ring r=0,x(1..5),dp;
600   module M=0;
601   def A=sheafCohBGG(M,-9,4);
602   // cohomology of cotangential bundle on P^3:
603   //-------------------------------------------
604   ring R=0,(x,y,z,u),dp;
605   resolution T1=mres(maxideal(1),0);
606   module M=T1[3];
607   intvec v=2,2,2,2,2,2;
608   attrib(M,"isHomog",v);
609   def B=sheafCohBGG(M,-8,4);
610}
611
612
613///////////////////////////////////////////////////////////////////////////////
614static proc showResult( def R, int l, int h )
615{
616  int PL = 1; // printlevel - voice + 2;
617// int PL = printlevel + 1;
618
619  intmat Betti;
620  if(typeof(R)=="resolution")
621  {
622    Betti = betti(R);
623  } else
624  {
625    if(typeof(R)!="intmat")
626    {
627      ERROR("Wrong input!!!");
628    }
629
630    Betti = R;
631  }
632
633  int n=nvars(basering)-1;
634  int ell = l + n;
635
636  int k     = ncols(Betti);
637  int row   = nrows(Betti);
638  int shift = attrib(Betti,"rowShift") + (k + ell - 1);
639
640  int iWTH = h-l+1;
641
642  int d = k - h + ell - 1;
643
644  if( PL > 1 )
645  {
646    "// l: ", l;
647    "// h: ", h;
648    "// n: ", n;
649    "// ell: ", ell;
650    "// k: ", k;
651    "// row: ", row;
652    "// shift: ", shift;
653    "// iWTH: ", iWTH;
654    "// d: ", d;
655  }
656
657  intmat newBetti[ n + 1 ][ iWTH ];
658  int i, j;
659
660  for (j=1; j<=row; j++) {
661    for (i=l; i<=h; i++) {
662      if( (n+2-shift-j)>0 ) {
663
664        if (  (k+1-j-i+ell-shift>0) and (j+i-ell+shift>=1)) {
665          newBetti[n+2-shift-j,i-l+1]=Betti[j,k+1-j-i+ell-shift];
666        }
667        else { newBetti[n+2-shift-j,i-l+1]=-1; }
668
669      }
670    }
671  }
672
673  for (j=2; j<=n+1; j++) {
674    for (i=1; i<min(j, iWTH); i++) {
675      newBetti[j,i]= -1;
676    }
677  }
678
679  for (j=1; j<=n; j++) {
680    for (i=iWTH; i>=k+j; i--) {
681      newBetti[j,i]=0; // -1;
682    }
683  }
684
685  if( PL > 0 )
686  {
687    "Cohomology table:";
688    displayCohom(newBetti, l, h, n);
689  }
690
691  return(newBetti);
692}
693///////////////////////////////////////////////////////////////////////////////
694
695
696
697///////////////////////////////////////////////////////////////////////////////
698proc sheafCohBGG2(module M,int l,int h)
699"USAGE:   sheafCohBGG2(M,l,h);    M module, l,h int
700ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
701         vector as an attribute, @code{h>=l}, and the basering has @code{n+1}
702         variables.
703RETURN:  intmat, cohomology of twists of the coherent sheaf F on P^n
704         associated to coker(M). The range of twists is determined by @code{l},
705         @code{h}.
706DISPLAY: The intmat is displayed in a diagram of the following form: @*
707  @format
708                l            l+1                      h
709  ----------------------------------------------------------
710      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
711           ...............................................
712      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
713      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
714  ----------------------------------------------------------
715    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
716  @end format
717         A @code{'-'} in the diagram refers to a zero entry; a @code{'*'}
718         refers to a negative entry (= dimension not yet determined).
719         refers to a not computed dimension. @*
720         If @code{printlevel}>=1, step-by step timings will be printed.
721         If @code{printlevel}>=2 we add progress debug messages
722         if @code{printlevel}>=3, even all intermediate results...
723NOTE:    This procedure is based on the Bernstein-Gel'fand-Gel'fand
724         correspondence and on Tate resolution ( see [Eisenbud, Floystad,
725         Schreyer: Sheaf cohomology and free resolutions over exterior
726         algebras, Trans AMS 355 (2003)] ).@*
727         @code{sheafCohBGG(M,l,h)} does not compute all values in the above
728         table. To determine all values of @code{h^i(F(d))}, @code{d=l..h},
729         use @code{sheafCohBGG(M,l-n,h+n)}.
730         Experimental version. Should require less memory.
731SEE ALSO: sheafCohBGG
732EXAMPLE: example sheafCohBGG2; shows an example
733"
734{
735  int PL = printlevel - voice + 2;
736//  int PL = printlevel;
737
738  dbprint(PL-1, "// sheafCohBGG2(M: "+ string(nrows(M)) + " x " + string(ncols(M)) +", " + string(l) + ", " + string(h) + "):");
739  dbprint(PL-2, M);
740
741  intvec save = option(get);
742
743  if( PL >= 2 )
744  {
745    option(prot);
746    option(mem);
747  }
748
749  def isCoker = attrib(M, "isCoker");
750  if( typeof(isCoker) == "int" )
751  {
752    if( isCoker > 0 )
753    {
754      dbprint(PL-1, "We are going to assume that M is given by coker matrix (that is, M is not a submodule presentation!)");
755    }
756  }
757
758  int i,j,k,row,col;
759
760  if( typeof(attrib(M,"isHomog"))!="intvec" )
761  {
762     if (size(M)==0) { attrib(M,"isHomog",0); }
763     else { ERROR("No admissible degree vector assigned"); }
764  }
765
766  dbprint(PL-1, "// weighting(M): ["+ string(attrib(M, "isHomog")) + "]");
767
768  option(redSB); option(redTail);
769
770  def R=basering;
771
772  int n = nvars(R) - 1;
773  int ell = l + n;
774
775
776/////////////////////////////////////////////////////////////////////////////
777// computations
778
779  int tBegin=timer;
780  int reg   = CM_regularity(M);
781  int tCMEnd = timer;
782
783  dbprint(PL-1, "// CM_reg(M): "+ string(reg));
784
785  int bound = max(reg + 1, h - 1);
786
787  dbprint(PL-1, "// bound: "+ string(bound));
788
789  ///////////////////////////////////////////////////////////////
790  int tSTDBegin=timer;
791  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!)
792  int tSTDEnd = timer;
793
794  dbprint(PL-1, "// M = std(M: "+string(nrows(M)) + " x " + string(ncols(M))  + ")");
795  dbprint(PL-2, M);
796  dbprint(PL-1, "// weighting(M): ["+ string(attrib(M, "isHomog")) + "]");
797
798
799  printlevel = printlevel + 1;
800  int tTruncateBegin=timer;
801  module MT = truncateFast(M, bound);
802  int tTruncateEnd=timer;
803  printlevel = printlevel - 1;
804
805  dbprint(PL-1, "// MT = truncateFast(M: "+string(nrows(MT)) + " x " + string(ncols(MT)) +", " + string(bound) + ")");
806  dbprint(PL-2, MT);
807  dbprint(PL-1, "// weighting(MT): ["+ string(attrib(MT, "isHomog")) + "]");
808
809  int m=nrows(MT);
810
811  ///////////////////////////////////////////////////////////////
812  int tTransposeJacobBegin=timer;
813  MT = jacob(MT); // ! :(
814  int tTransposeJacobEnd=timer;
815
816  dbprint(PL-1, "// MT = jacob(MT: "+string(nrows(MT)) + " x " + string(ncols(MT))  + ")");
817  dbprint(PL-2, MT);
818  dbprint(PL-1, "// weighting(MT): ["+ string(attrib(MT, "isHomog")) + "]");
819
820  int tSyzBegin=timer;
821  MT = syz(MT);
822  int tSyzEnd=timer;
823
824  dbprint(PL-1, "// MT = syz(MT: "+string(nrows(MT)) + " x " + string(ncols(MT))  + ")");
825  dbprint(PL-2, MT);
826  dbprint(PL-1, "// weighting(MT): ["+ string(attrib(MT, "isHomog")) + "]");
827
828  int tMatrixOppBegin=timer;
829  module SS = TensorModuleMult(m, MT);
830  int tMatrixOppEnd=timer;
831
832  dbprint(PL-1, "// SS = TensorModuleMult("+ string(m)+ ", MT: "+string(nrows(MT)) + " x " + string(ncols(MT))  + ")");
833  dbprint(PL-2, SS);
834  dbprint(PL-1, "// weighting(SS): ["+ string(attrib(SS, "isHomog")) + "]");
835
836  //--- to the exterior algebra
837  def AR = Exterior(); setring AR;
838
839  dbprint(PL-1, "// Test: var(1) * var(1): "+ string(var(1) * var(1)));
840
841  int maxbound = max(1, bound - ell + 1);
842//  int maxbound = max(1, bound - l + 1); // As In M2!!!
843
844  dbprint(PL-1, "// maxbound: "+ string(maxbound));
845
846  //--- here we are with our matrix
847  module EM=imap(R,SS);
848  intvec w;
849  for (i=1; i<=nrows(EM); i++)
850  {
851     w[i]=-bound-1;
852  }
853
854  attrib(EM,"isHomog",w);
855
856  ///////////////////////////////////////////////////////////////
857
858  dbprint(PL-1, "// EM: "+string(nrows(EM)) + " x " + string(ncols(EM))  + ")");
859  dbprint(PL-2, EM);
860  dbprint(PL-1, "// weighting(EM): ["+ string(attrib(EM, "isHomog")) + "]");
861
862  int tResulutionBegin=timer;
863  resolution RE = nres(EM, maxbound); // TODO: Plural computes one too many syzygies...?!
864  int tMinResBegin=timer;
865  RE = minres(RE);
866  int tBettiBegin=timer;
867  intmat Betti = betti(RE); // betti(RE, 1);?
868  int tResulutionEnd=timer;
869
870  int tEnd = tResulutionEnd;
871
872  if( PL > 0 )
873  {
874    //    list L = RE; // TODO: size(L/RE) is wrong!
875    "
876        ----      RESULTS  ----
877        Tate Resolution:
878    ";
879    RE;
880    "Betti numbers for Tate resolution (diagonal cohomology table):";
881    print(Betti, "betti"); // Diagonal form!
882  }
883
884//  printlevel = printlevel + 1;
885  Betti = showResult(Betti, l, h ); // Show usual form of cohomology table
886//  printlevel = printlevel - 1;
887
888  if(PL > 0)
889  {
890  "
891      ----      TIMINGS     -------
892      Trunc      Time: ", tTruncateEnd - tTruncateBegin, "
893      Reg        Time: ", tCMEnd - tBegin, "
894      kStd       Time: ", tSTDEnd - tSTDBegin, "
895      Jacob      Time: ", tTransposeJacobEnd - tTransposeJacobBegin, "
896      Syz        Time: ", tSyzEnd - tSyzBegin, "
897      Mat        Time: ", tMatrixOppEnd - tMatrixOppBegin, "
898      ------------------------------
899      Res        Time: ", tResulutionEnd - tResulutionBegin, "
900      :: NRes    Time: ", tMinResBegin - tResulutionBegin, "
901      :: MinRes .Time: ", tBettiBegin - tMinResBegin, "
902      :: Betti  .Time: ", tResulutionEnd - tBettiBegin, "
903      ---------------------------------------------------------
904      Total Time: ", tEnd - tBegin, "
905      ---------------------------------------------------------
906      ";
907  }
908
909  setring R;
910
911  option(set, save);
912
913  return(Betti);
914}
915example
916{"EXAMPLE:";
917   echo = 2;
918   int pl = printlevel;
919   int l,h, t;
920
921   //-------------------------------------------
922   // cohomology of structure sheaf on P^4:
923   //-------------------------------------------
924   ring r=32001,x(1..5),dp;
925
926   module M= getStructureSheaf(); // OO_P^4
927
928   l = -12; h = 12; // range of twists: l..h
929
930   printlevel = 0;
931   //////////////////////////////////////////////
932   t = timer;
933
934   def A = sheafCoh(M, l, h); // global Ext method:
935
936   "Time: ", timer - t;
937   //////////////////////////////////////////////
938   t = timer;
939
940   A = sheafCohBGG(M, l, h);  // BGG method (without optimization):
941
942   "Time: ", timer - t;
943   //////////////////////////////////////////////
944   t = timer;
945
946   A = sheafCohBGG2(M, l, h); // BGG method (with optimization)
947
948   "Time: ", timer - t;
949   //////////////////////////////////////////////
950   printlevel = pl;
951
952   kill A, r;
953
954   //-------------------------------------------
955   // cohomology of cotangential bundle on P^3:
956   //-------------------------------------------
957   ring R=32001,(x,y,z,u),dp;
958
959   module M = getCotangentialBundle();
960
961   l = -12; h = 11; // range of twists: l..h
962
963
964   //////////////////////////////////////////////
965   printlevel = 0;
966   t = timer;
967
968   def B = sheafCoh(M, l, h); // global Ext method:
969
970   "Time: ", timer - t;
971   //////////////////////////////////////////////
972   t = timer;
973
974   B = sheafCohBGG(M, l, h);  // BGG method (without optimization):
975
976   "Time: ", timer - t;
977   //////////////////////////////////////////////
978   t = timer;
979
980   B = sheafCohBGG2(M, l, h); // BGG method (with optimization)
981
982   "Time: ", timer - t;
983   //////////////////////////////////////////////
984   printlevel = pl;
985}
986
987
988///////////////////////////////////////////////////////////////////////////////
989
990proc dimH(int i,module M,int d)
991"USAGE:   dimH(i,M,d);    M module, i,d int
992ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
993         vector as an attribute, @code{h>=l}, and the basering @code{S} has
994         @code{n+1} variables.
995RETURN:  int,  vector space dimension of @math{H^i(F(d))} for F the coherent
996         sheaf on P^n associated to coker(M).
997NOTE:    The procedure is based on local duality as described in [Eisenbud:
998         Computing cohomology. In Vasconcelos: Computational methods in
999         commutative algebra and algebraic geometry. Springer (1998)].
1000SEE ALSO: sheafCoh, sheafCohBGG
1001EXAMPLE: example dimH; shows an example
1002"
1003{
1004  if( typeof(attrib(M,"isHomog"))=="string" )
1005  {
1006    if (size(M)==0)
1007    {
1008      // assign weights 0 to generators of R^n (n=nrows(M))
1009      intvec v;
1010      v[nrows(M)]=0;
1011      attrib(M,"isHomog",v);
1012    }
1013    else
1014    {
1015      ERROR("No admissible degree vector assigned");
1016    }
1017  }
1018  int Result;
1019  int n=nvars(basering)-1;
1020  if ((i>0) and (i<=n)) {
1021    list L=Ext_R(n-i,M,1)[2];
1022    def N=L[1];
1023    return(dimGradedPart(N,-n-1-d));
1024  }
1025  else
1026  {
1027    if (i==0)
1028    {
1029      list L=Ext_R(intvec(n+1,n+2),M,1)[2];
1030      def N0=L[2];
1031      def N1=L[1];
1032      Result=dimGradedPart(M,d) - dimGradedPart(N0,-n-1-d)
1033                                - dimGradedPart(N1,-n-1-d);
1034      return(Result);
1035    }
1036    else {
1037      return(0);
1038    }
1039  }
1040}
1041example
1042{"EXAMPLE:";
1043   echo = 2;
1044   ring R=0,(x,y,z,u),dp;
1045   resolution T1=mres(maxideal(1),0);
1046   module M=T1[3];
1047   intvec v=2,2,2,2,2,2;
1048   attrib(M,"isHomog",v);
1049   dimH(0,M,2);
1050   dimH(1,M,0);
1051   dimH(2,M,1);
1052   dimH(3,M,-5);
1053}
1054
1055
1056///////////////////////////////////////////////////////////////////////////////
1057
1058proc sheafCoh(module M,int l,int h,list #)
1059"USAGE:   sheafCoh(M,l,h);    M module, l,h int
1060ASSUME:  @code{M} is graded, and it comes assigned with an admissible degree
1061         vector as an attribute, @code{h>=l}. The basering @code{S} has
1062         @code{n+1} variables.
1063RETURN:  intmat, cohomology of twists of the coherent sheaf F on P^n
1064         associated to coker(M). The range of twists is determined by @code{l},
1065         @code{h}.
1066DISPLAY: The intmat is displayed in a diagram of the following form: @*
1067  @format
1068                l            l+1                      h
1069  ----------------------------------------------------------
1070      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
1071           ...............................................
1072      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
1073      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
1074  ----------------------------------------------------------
1075    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
1076  @end format
1077         A @code{'-'} in the diagram refers to a zero entry.
1078NOTE:    The procedure is based on local duality as described in [Eisenbud:
1079         Computing cohomology. In Vasconcelos: Computational methods in
1080         commutative algebra and algebraic geometry. Springer (1998)].@*
1081         By default, the procedure uses @code{mres} to compute the Ext
1082         modules. If called with the additional parameter @code{\"sres\"},
1083         the @code{sres} command is used instead.
1084SEE ALSO: dimH, sheafCohBGG
1085EXAMPLE: example sheafCoh; shows an example
1086"
1087{
1088  int use_sres;
1089  if( typeof(attrib(M,"isHomog"))!="intvec" )
1090  {
1091     if (size(M)==0) { attrib(M,"isHomog",0); }
1092     else { ERROR("No admissible degree vector assigned"); }
1093  }
1094  if (size(#)>0)
1095  {
1096    if (#[1]=="sres") { use_sres=1; }
1097  }
1098  int i,j;
1099  module N,N0,N1;
1100  int n=nvars(basering)-1;
1101  intvec v=0..n+1;
1102  int col=h-l+1;
1103  intmat newBetti[n+1][col];
1104  if (use_sres) { list L=Ext_R(v,M,1,"sres")[2]; }
1105  else          { list L=Ext_R(v,M,1)[2]; }
1106  for (i=l; i<=h; i++)
1107  {
1108    N0=L[n+2];
1109    N1=L[n+1];
1110    newBetti[n+1,i-l+1]=dimGradedPart(M,i) - dimGradedPart(N0,-n-1-i)
1111                             - dimGradedPart(N0,-n-1-i);
1112  }
1113  for (j=1; j<=n; j++)
1114  {
1115     N=L[j];
1116     attrib(N,"isSB",1);
1117     if (dim(N)>=0) {
1118       for (i=l; i<=h; i++)
1119       {
1120         newBetti[j,i-l+1]=dimGradedPart(N,-n-1-i);
1121       }
1122     }
1123  }
1124  displayCohom(newBetti,l,h,n);
1125  return(newBetti);
1126}
1127example
1128{"EXAMPLE:";
1129   echo = 2;
1130   //
1131   // cohomology of structure sheaf on P^4:
1132   //-------------------------------------------
1133   ring r=0,x(1..5),dp;
1134   module M=0;
1135   def A=sheafCoh(0,-7,2);
1136   //
1137   // cohomology of cotangential bundle on P^3:
1138   //-------------------------------------------
1139   ring R=0,(x,y,z,u),dp;
1140   resolution T1=mres(maxideal(1),0);
1141   module M=T1[3];
1142   intvec v=2,2,2,2,2,2;
1143   attrib(M,"isHomog",v);
1144   def B=sheafCoh(M,-6,2);
1145}
1146
1147///////////////////////////////////////////////////////////////////////////////
1148proc displayCohom (intmat data, int l, int h, int n)
1149"USAGE:   displayCohom(data,l,h,n);  data intmat, l,h,n int
1150ASSUME:  @code{h>=l}, @code{data} is the return value of
1151         @code{sheafCoh(M,l,h)} or of @code{sheafCohBGG(M,l,h)}, and the
1152         basering has @code{n+1} variables.
1153RETURN:  none
1154NOTE:    The intmat is displayed in a diagram of the following form: @*
1155  @format
1156                l            l+1                      h
1157  ----------------------------------------------------------
1158      n:     h^n(F(l))    h^n(F(l+1))   ......    h^n(F(h))
1159           ...............................................
1160      1:     h^1(F(l))    h^1(F(l+1))   ......    h^1(F(h))
1161      0:     h^0(F(l))    h^0(F(l+1))   ......    h^0(F(h))
1162  ----------------------------------------------------------
1163    chi:     chi(F(l))    chi(F(l+1))   ......    chi(F(h))
1164  @end format
1165         where @code{F} refers to the associated sheaf of @code{M} on P^n.@*
1166         A @code{'-'} in the diagram refers to a zero entry,  a @code{'*'}
1167         refers to a negative entry (= dimension not yet determined).
1168"
1169{
1170  int i,j,k,dat,maxL;
1171  intvec notSumCol;
1172  notSumCol[h-l+1]=0;
1173  string s;
1174  maxL=4;
1175  for (i=1;i<=nrows(data);i++)
1176  {
1177    for (j=1;j<=ncols(data);j++)
1178    {
1179      if (size(string(data[i,j]))>=maxL-1)
1180      {
1181        maxL=size(string(data[i,j]))+2;
1182      }
1183    }
1184  }
1185  string Row="    ";
1186  string Row1="----";
1187  for (i=l; i<=h; i++) {
1188    for (j=1; j<=maxL-size(string(i)); j++)
1189    {
1190      Row=Row+" ";
1191    }
1192    Row=Row+string(i);
1193    for (j=1; j<=maxL; j++) { Row1 = Row1+"-"; }
1194  }
1195  print(Row);
1196  print(Row1);
1197  for (j=1; j<=n+1; j++)
1198  {
1199    s = string(n+1-j);
1200    Row = "";
1201    for(k=1; k<4-size(s); k++) { Row = Row+" "; }
1202    Row = Row + s+":";
1203    for (i=0; i<=h-l; i++)
1204    {
1205      dat = data[j,i+1];
1206      if (dat>0) { s = string(dat); }
1207      else
1208      {
1209        if (dat==0) { s="-"; }
1210        else        { s="*"; notSumCol[i+1]=1; }
1211      }
1212      for(k=1; k<=maxL-size(s); k++) { Row = Row+" "; }
1213      Row = Row + s;
1214    }
1215    print(Row);
1216  }
1217  print(Row1);
1218  Row="chi:";
1219  for (i=0; i<=h-l; i++)
1220  {
1221    dat = 0;
1222    if (notSumCol[i+1]==0)
1223    {
1224      for (j=0; j<=n; j++) { dat = dat + (-1)^j * data[n+1-j,i+1]; }
1225      s = string(dat);
1226    }
1227    else { s="*"; }
1228    for (k=1; k<=maxL-size(s); k++) { Row = Row+" "; }
1229    Row = Row + s;
1230  }
1231  print(Row);
1232}
1233///////////////////////////////////////////////////////////////////////////////
1234
1235proc getStructureSheaf(list #)
1236{
1237
1238  if( size(#) == 0 )
1239  {
1240    module M = 0;
1241    intvec v = 0;
1242    attrib(M,"isHomog",v);
1243//    homog(M);
1244
1245    attrib(M, "isCoker", 1);
1246
1247//     attrib(M);
1248    return(M);
1249  }
1250
1251  if( typeof(#[1]) == "ideal")
1252  {
1253    ideal I = #[1];
1254
1255    if( size(#) == 2 )
1256    {
1257      if( typeof(#[2]) == "int" )
1258      {
1259        if( #[2] != 0 )
1260        {
1261          qring @@@@QQ = std(I);
1262
1263          module M = getStructureSheaf();
1264
1265          export M;
1266
1267//          keepring @@@@QQ; // This is a bad idea... :(?
1268          return (@@@@QQ);
1269        }
1270      }
1271    }
1272
1273/*
1274    // This seems to be wrong!!!
1275    module M = I * gen(1);
1276    homog(M);
1277
1278    M = modulo(gen(1), module(I * gen(1))); // basering^1 / I
1279
1280    homog(M);
1281
1282    attrib(M, "isCoker", 1);
1283
1284    attrib(M);
1285    return(M);
1286*/
1287  }
1288
1289  ERROR("Wrong argument");
1290
1291}
1292example
1293{"EXAMPLE:";
1294   echo = 2; int pl = printlevel;
1295   printlevel = voice;
1296
1297
1298   ////////////////////////////////////////////////////////////////////////////////
1299   ring r;
1300   module M = getStructureSheaf();
1301   "Basering: ";
1302   basering;
1303   "Module: ", string(M), ", grading is given by weights: ", attrib(M, "isHomog");
1304
1305   def A=sheafCohBGG2(M,-9,9);
1306   print(A);
1307
1308   ////////////////////////////////////////////////////////////////////////////////
1309   setring r;
1310   module M = getStructureSheaf(ideal(var(1)), 0);
1311
1312   "Basering: ";
1313   basering;
1314   "Module: ", string(M), ", grading is given by weights: ", attrib(M, "isHomog");
1315
1316   def A=sheafCohBGG2(M,-9,9);
1317   print(A);
1318
1319   ////////////////////////////////////////////////////////////////////////////////
1320   setring r;
1321   def Q = getStructureSheaf(ideal(var(1)), 1); // returns a new ring!
1322   setring Q; // M was exported in the new ring!
1323
1324   "Basering: ";
1325   basering;
1326   "Module: ", string(M), ", grading is given by weights: ", attrib(M, "isHomog");
1327
1328   def A=sheafCohBGG2(M,-9,9);
1329   print(A);
1330
1331   printlevel = pl;
1332}
1333
1334
1335proc getCotangentialBundle()
1336{
1337  resolution T1=mres(maxideal(1),3);
1338  module M=T1[3];
1339//  attrib(M,"isHomog");
1340//  homog(M);
1341  attrib(M, "isCoker", 1);
1342  // attrib(M);
1343  return (M);
1344}
1345
1346proc getIdealSheafPullback(ideal I, ideal pi)
1347{
1348  def save = basering;
1349  map P = save, pi;
1350  return( P(I) );
1351}
1352
1353// TODO: set attributes!
1354
1355
1356proc getIdealSheaf(ideal I)
1357{
1358  int i = homog(I);
1359  resolution FI = mres(I,2); // Syz + grading...
1360  module M = FI[2];
1361  attrib(M, "isCoker", 1);
1362  //  attrib(M);
1363  return(M);
1364}
1365
1366
1367
1368
1369
1370/*
1371Examples:
1372---------
1373 LIB "sheafcoh.lib";
1374
1375 ring S = 32003, x(0..4), dp;
1376 module MI=maxideal(1);
1377 attrib(MI,"isHomog",intvec(-1));
1378 resolution kos = nres(MI,0);
1379 print(betti(kos),"betti");
1380 LIB "random.lib";
1381 matrix alpha0 = random(32002,10,3);
1382 module pres = module(alpha0)+kos[3];
1383 attrib(pres,"isHomog",intvec(1,1,1,1,1,1,1,1,1,1));
1384 resolution fcokernel = mres(pres,0);
1385 print(betti(fcokernel),"betti");
1386 module dir = transpose(pres);
1387 attrib(dir,"isHomog",intvec(-1,-1,-1,-2,-2,-2,
1388                             -2,-2,-2,-2,-2,-2,-2));
1389 resolution fdir = mres(dir,2);
1390 print(betti(fdir),"betti");
1391 ideal I = groebner(flatten(fdir[2]));
1392 resolution FI = mres(I,0);
1393 print(betti(FI),"betti");
1394 module F=FI[2];
1395 int t=timer;
1396 def A1=sheafCoh(F,-8,8);
1397 timer-t;
1398 t=timer;
1399 def A2=sheafCohBGG(F,-8,8);
1400 timer-t;
1401
1402 LIB "sheafcoh.lib";
1403 LIB "random.lib";
1404 ring S = 32003, x(0..4), dp;
1405 resolution kos = nres(maxideal(1),0);
1406 betti(kos);
1407 matrix kos5 = kos[5];
1408 matrix tphi = transpose(dsum(kos5,kos5));
1409 matrix kos3 = kos[3];
1410 matrix psi = dsum(kos3,kos3);
1411 matrix beta1 = random(32002,20,2);
1412 matrix tbeta1tilde = transpose(psi*beta1);
1413 matrix tbeta0 = lift(tphi,tbeta1tilde);
1414 matrix kos4 = kos[4];
1415 matrix tkos4pluskos4 = transpose(dsum(kos4,kos4));
1416 matrix tgammamin1 = random(32002,20,1);
1417 matrix tgamma0 = tkos4pluskos4*tgammamin1;
1418 matrix talpha0 = concat(tbeta0,tgamma0);
1419 matrix zero[20][1];
1420 matrix tpsi = transpose(psi);
1421 matrix tpresg = concat(tpsi,zero);
1422 matrix pres = module(transpose(talpha0))
1423                    + module(transpose(tpresg));
1424 module dir = transpose(pres);
1425 dir = prune(dir);
1426 homog(dir);
1427 intvec deg_dir = attrib(dir,"isHomog");
1428 attrib(dir,"isHomog",deg_dir-2);        // set degrees
1429 resolution fdir = mres(prune(dir),2);
1430 print(betti(fdir),"betti");
1431 ideal I = groebner(flatten(fdir[2]));
1432 resolution FI = mres(I,0);
1433
1434 module F=FI[2];
1435 def A1=sheafCoh(F,-5,7);
1436 def A2=sheafCohBGG(F,-5,7);
1437
1438*/
Note: See TracBrowser for help on using the repository browser.