source: git/Singular/LIB/sheafcoh.lib @ a33befb

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