source: git/Singular/LIB/derham.lib @ 76d26c

spielwiese
Last change on this file since 76d26c was 76d26c, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: rename minIntRoot to minIntRootD to avoid name clash
  • Property mode set to 100644
File size: 180.1 KB
RevLine 
[0e8a5a]1///////////////////////////////////////////////////////////////////////////////
[08fa62]2version="version derham.lib 4.0.0.0 Jun_2013 ";
[0e8a5a]3category="Noncommutative";
4info="
5LIBRARY:  derham.lib      Computation of deRham cohomology
6
7AUTHORS:  Cornelia Rottner, rottner@mathematik.uni-kl.de
8
9OVERVIEW:
10A library for computing the de Rham cohomology of complements of complex affine
11varieties.
12
13
14REFERENCES:
15[OT] Oaku, T.; Takayama, N.: Algorithms of D-modules - restriction, tensor product,
16     localzation, and local cohomology groups}, J. Pure Appl. Algebra 156, 267-308
17     (2001)
18[R]  Rottner, C.: Computing de Rham Cohomology,diploma thesis (2012)
19[W1] Walther, U.: Algorithmic computation of local cohomology modules and the local
20     cohomological dimension of algebraic varieties}, J. Pure Appl. Algebra 139,
21     303-321 (1999)
22[W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of
23     Complex Affine Varieties}, J. Symbolic Computation 29, 796-839 (2000)
24[W3] Walther, U.: Computing the cup product structure for complements of complex
25     affine varieties, J. Pure Appl. Algebra 164, 247-273 (2001)
26
27
28PROCEDURES:
29
30deRhamCohomology(list[,opt]); computes the de Rham cohomology
31MVComplex(list);              computes the Mayer-Vietoris complex
32";
33
34LIB "nctools.lib";
35LIB "matrix.lib";
36LIB "qhmoduli.lib";
37LIB "general.lib";
[08fa62]38//LIB "dmod.lib";
[0e8a5a]39LIB "bfun.lib";
40LIB "dmodapp.lib";
41LIB "poly.lib";
42LIB "schreyer.lib";
43LIB "dmodloc.lib";
44
45
46////////////////////////////////////////////////////////////////////////////////////
47
48proc deRhamCohomology(list L,list #)
49"USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices
50        optional list consisting of one up to three strings @*
51        The optional strings may be one of the strings@*
52        -'noCE': compute quasi-isomorphic complexes without using Cartan-Eilenberg
53         resolutionsq@*
54        -'Vdres': compute quasi-isomorphic complexes using Cartan-Eilenberg
55         resolutions; the CE resolutions are computed via V__d-homogenization
56         and without using Schreyer's method @*
57        -'Sres': compute quasi-isomorphic complexes using Cartan-Eilenberg
58         resolutions in the homogenized Weyl algebra via Schreyer's method@*
59        one of the strings@*
60        -'iterativeloc': compute localizations by factorizing the polynomials and
61         sucessive localization of the factors @*
62        -'no iterativeloc': compute localizations by directly localizing the
63         product@*
64        and one of the strings
65        -'onlybounds': computes bounds for the minimal and maximal interger roots
66         of the global b-function
67        -'exactroots' computes the minimal and maximal integer root of the global
68         b-function
69        The default is 'noCE', 'iterativeloc' and 'onlybounds'.
70ASSUME: -The basering must be a polynomial ring over the field of rational numbers@*
71RETURN: list, where the ith entry is the (i-1)st de Rham cohomology group of the
72        complement of the complex affine variety given by the polynomials in L
73EXAMPLE:example deRhamCohomology; shows an example
74"
75{
76  intvec saveoptions=option(get);
77  intvec i1,i2;
78  option(none);
79  int recursiveloc=1;
80  int i,j,nr,nc;
81  def R=basering;
82  poly islcm, forlcm;
83  int n=nvars(R);
84  int le=size(L)+n;
85  string Syzstring="noCE";
86  int onlybounds=1;
87  int diffforms;
88  for (i=1; i<=size(#); i++)
89    {
90      if (#[i]=="Sres")
91        {
92          Syzstring="Sres";
93        }
94      if (#[i]=="Vdres")
95        {
96          Syzstring="Vdres";
97        }
98      if (#[i]=="noiterativeloc")
99        {
100          recursiveloc=0;
101        }
102      if (#[i]=="exactroots")
103        {
104          onlybounds=0;
105        }
106      if (#[i]=="diffforms")
107        {
108          diffforms=1;
109        }
110    }
111  for (i=1; i<=size(L); i++)
112    {
113      if (L[i]==0)
114        {
115          L=delete(L,i);
116          i=i-1;
117        }
118    }
119  if (size(L)==0)
120    {
121      return (list(0));//////////////////////////////////////////////////////////////////stimmt das jetzt?!??????????????????????????????????
122    }
123  for (i=1; i<= size(L); i++)
124    {
125      if (leadcoef(L[i])-L[i]==0)
126        {
127          return(list(1));    ///////////////////////////////////////////////////////////////stimmt das jetzt?!????????????????????????????????????
128        }
129    }
130  if (size(L)==0)
131    {
132      /*the complement of the variety given by the input is the whole space*/
133      return(list(1));
134    }
135  for (i=1; i<=size(L); i++)
136    {
137      if (typeof(L[i])!="poly")
138        {
139          print("The input list must consist of polynomials");
140          return();
141        }
142    }
143  if (size(L)==1 and Syzstring=="noCE")
144    {
145      Syzstring="Sres";
146    }
147  /* 1st step: compute the Mayer-Vietoris Complex and its Fourier transform*/
148  def W=MVComplex(L,recursiveloc);//new ring that contains the MV complex
149  setring W;
150  list fortoVdstrict=MV;
151  if (diffforms==0)
152    {
153      ideal IFourier=var(n+1);
154      for (i=2;i<=n;i++)
155        {
156          IFourier=IFourier,var(n+i);
157        }
158      for (i=1; i<=n;i++)
159        {
160          IFourier=IFourier,-var(i);
161        }
162      map cFourier=W,IFourier;
163      matrix sup;
164      for (i=1; i<=size(MV); i++)
165        {
166          sup=fortoVdstrict[i];
167          /*takes the Fourier transform of the MV complex*/
168          fortoVdstrict[i]=cFourier(sup);
169        }
170    }
171  /* 2nd step: Compute a V_d-strict free complex that is quasi-isomorphic to the
172     complex fortoVdstrict
173     The 1st entry of the list rem will be the quasi-isomorphic complex, the 2nd
174     entry contains the cohomology modules and is needed for the computation of the
175     global b-function*/
176  if (Syzstring=="noCE")
177    {
178      list rem=quasiisomorphicVdComplex(fortoVdstrict,diffforms);
179      list quasiiso=rem[3];
180    }
181  else
182    {
183      list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring,diffforms);
184      if (diffforms==1)
185        {
186          list quasiiso=list(matrix(1,1,1));
187        }
188    }
189  list newcomplex=rem[1];
190////////////////////////////////////////////////////////////////////////////////////
191  /* 3rd step: Compute the  bounds for the minimal and maximal integer root of the
192     global b-function of newcomplex(i.e. compute the lcm of the b-functions of its
193     cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal
194     integer root.
195
196     If we compute only the bounds, we omit additional Groebner basis computations.
197     However this leads to a higher-dimensional truncated complex.
198
199     Note that the  cohomology modules are already contained in rem[2].
200     minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/
201  if (diffforms==1)
202    {
203      list minmaxk=exactGlobalBFunIntegration(rem[2]);
204    }
205  else
206    {
207      if (onlybounds==1)
208        {
209          list minmaxk=globalBFun(rem[2],Syzstring);
210        }
211      else
212        {
213          list minmaxk=exactGlobalBFun(rem[2],Syzstring);
214        }
215    }
216  if (size(minmaxk)==0)
217    {
218      return (0);
219    }
220  ///////////////////////////////////////////////////////////////////////////Bis hierhin angepasst
221  /*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where
222    C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i-2],
223    m^i=newcomplex[3*i-1], d^i=newcomplex[3*i]), using Thm 5.7 in [W1]:
224    The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set
225    (0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and
226    if it is placed in component k it holds that
227    minmaxk[1]-m^i[k]<=deg(P_(i_j))<=minmaxk[2]-m^i[k]*/
228  int k,l;
229  list truncatedcomplex,shorten,upto;
230  for (i=1; i<=size(newcomplex) div 3; i++)
231    {
232      shorten[3*i-1]=list();
233      for (j=1; j<=size(newcomplex[3*i-1]); j++)
234        {
235          /*shorten[3*i-1][j][k]=minmaxk[k]-m^i[j]+1 (for k=1,2) if this value is
236            positive otherwise we will set it to be list();
237.-            we added +1, because we will use a list, where we put in position l
238            polys of degree l+1*/
239          shorten[3*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1);
240          if (diffforms==1)
241            {
242              shorten[3*i-1][j][1]=1;
243            }
244          shorten[3*i-1][j][2]=minmaxk[2]-newcomplex[3*i-1][j]+1;
245          upto[size(upto)+1]=shorten[3*i-1][j][2];
246          if (shorten[3*i-1][j][2]<=0)
247            {
248              shorten[3*i-1][j]=list();
249            }
250          else
251            {
252              if (shorten[3*i-1][j][1]<=0)
253                {
254                  shorten[3*i-1][j][1]=1;
255                }
256            }
257        }
258    }
259  int iupto=Max(upto);//maximal degree +1 of the polynomials we have to consider
260  if (iupto<=0)
261    {
262      return(list(0));
263    }
264  list allpolys;
265  /*allpolys[i] will consist list of all monomials in D(1),...,D(n) of degree i-1*/
266  allpolys[1]=list(1);
267  list minvar;
268  list keepv;
269  minvar[1]=list(1);
270  for (i=1; i<=iupto-1; i++)
271    {
272      allpolys[i+1]=list();
273      minvar[i+1]=list();
274      for (k=1; k<=size(allpolys[i]); k++)
275        {
276          for (j=minvar[i][k]; j<=nvars(W) div 2; j++)
277            {
278              if (diffforms==0)
279                {
280                  allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j);
281                }
282              else
283                {
284                  allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*x(j);
285                }
286              minvar[i+1][size(minvar[i+1])+1]=j;
287            }
288        }
289    }
290  list keepformatrix,sizetruncom,fortrun,fst;
291  int count,stc;
292  intvec v,forin;
293  matrix subm;
294  list keepcount;
295  list passendespoly;
296  /*now we compute the truncation*/
297  for (i=1; i<=size(newcomplex) div 3; i++)
298    {
299      /*truncatedcomplex[2*i-1] will contain all the generators for the truncation
300        of D_n/(x(1),..,x(n))\otimes C[i]*/
301      truncatedcomplex[2*i-1]=list();
302      sizetruncom[2*i-1]=list();
303      sizetruncom[2*i]=list();
304      passendespoly[i]=list();
305      /*truncatedcomplex[2*i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i])
306        ->trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/
307      truncatedcomplex[2*i]=newcomplex[3*i];
308      v=0;count=0;
309      sizetruncom[2*i][1]=0;
310      for (j=1; j<=newcomplex[3*i-2]; j++)
311        {
312          if (size(shorten[3*i-1][j])!=0)
313            {
314              fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]);
315              truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1];
316              for (k=1; k<=size(fortrun[1]); k++)
317                {
318                  for (l=1; l<=size(fortrun[1][k]); l++)
319                    {
320                      passendespoly[i][size(passendespoly[i])+1]=list(fortrun[1][k][l][1],j);
321                    }
322                }
323              count=count+fortrun[2];
324              fst=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1);
325              sizetruncom[2*i-1][size(sizetruncom[2*i-1])+1]=fst;
326              sizetruncom[2*i][size(sizetruncom[2*i])+1]=count;
327              if (v!=0)
328                {
329                  v[size(v)+1]=j;
330                }
331              else
332                {
333                  v[1]=j;
334                }
335            }
336        }
337      if (v!=0)
338        {
339          keepv[i]=v;
340          subm=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i]));
341          truncatedcomplex[2*i]=subm;
342          if (i!=1)
343            {
344              i1=1..nrows(truncatedcomplex[2*(i-1)]);
345              subm=submat(truncatedcomplex[2*(i-1)],i1,v);
346              truncatedcomplex[2*(i-1)]=subm;
347            }
348        }
349      else
350        {
351          keepv[i]=list();
352          truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i]));
353          if (i!=1)
354            {
355              nr=nrows(truncatedcomplex[2*(i-1)]);
356              truncatedcomplex[2*(i-1)]=matrix(0,nr,1);
357            }
358        }
359    }
360  list keeptruncatedcomplex=truncatedcomplex;
361  matrix M;
362  int st,pi,pj;
363  poly ptc;
364  int b,d,ideg,kplus,lplus;
365  int z;
366  poly form,lform,nform;
367  /*computation of the maps*/
368  if (diffforms==1)
369    {
370      def ConvWeyl=makeConverseWeyl(nvars(basering) div 2);
371      setring ConvWeyl;
372      poly form,lform,nform;
373      poly ptc;
374      list truncatedcomplex;
375      matrix M;
376      ideal I=x(1);
377      for (i=2; i<=nvars(basering) div 2; i++)
378        {
379          I=I,var(nvars(basering) div 2 + i);
380        }
381      for (i=1; i<=nvars(basering) div 2; i++)
382        {
383          I=I,var(i);
384        }
385      map transtc=W,I;
386      truncatedcomplex=transtc(truncatedcomplex);
387    }
388  for (i=1; i<size(truncatedcomplex) div 2; i++)
389    {
390      nr=max(1,sizetruncom[2*i][size(sizetruncom[2*i])]);
391      nc=max(1,sizetruncom[2*i+2][size(sizetruncom[2*i+2])]);
392      M=matrix(0,nr,nc);
393      for (k=1; k<=size(truncatedcomplex[2*i-1]);k++)
394        {
395          for (l=1; l<=size(truncatedcomplex[2*(i+1)-1]); l++)
396            {
397              if (size(sizetruncom[2*i])!=1)
398                {
399                  for (j=1; j<=size(truncatedcomplex[2*i-1][k]); j++)
400                    {
401                      for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++)
402                        {
403                          form=truncatedcomplex[2*i-1][k][j][b][1];
404                          form=form*truncatedcomplex[2*i][k,l];
405
406
407                          for (z=1; z<=nvars(basering) div 2; z++)//neu
408                            {//
409                              form=subst(form,var(z),0);//
410                            }//
411
412                          while (form!=0)
413                            {
414                              lform=lead(form);
415                              v=leadexp(lform);
416                              v=v[1..n];
417                              // if (v==(0:n))
418                              //{
419                                  ideg=deg(lform)-sizetruncom[2*(i+1)-1][l][1];
420                                  if (ideg>=0)
421                                    {
422                                      nr=ideg+1;
423                                      st=size(truncatedcomplex[2*(i+1)-1][l][nr]);
424                                      for (d=1; d<=st;d++)
425                                        {
426                                          nc=2*(i+1)-1;
427                                          ptc=truncatedcomplex[nc][l][ideg+1][d][1];
428                                          if (leadmonom(lform)==ptc)
429                                            {
430                                              nr=2*i-1;
431                                              pi=truncatedcomplex[nr][k][j][b][2];
432                                              pi=pi+sizetruncom[2*i][k];
433                                              nc=2*(i+1)-1;
434                                              nr=ideg+1;
435                                              pj=truncatedcomplex[nc][l][nr][d][2];
436                                              pj=pj+sizetruncom[2*(i+1)][l];
437                                              M[pi,pj]=leadcoef(lform);
438                                              break;
439                                            }
440                                        }
441                                    }
442                                  //        }
443
444                              form=form-lform;
445                            }
446                        }
447                    }
448                }
449            }
450        }
451      truncatedcomplex[2*i]=M;
452      truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])];
453    }
454  truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])];
455  if (truncatedcomplex[2*i-1]!=0)
456    {
457      truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1);
458    }
459  if (diffforms==1)
460    {
461      setring W;
462    truncatedcomplex=imap(ConvWeyl,truncatedcomplex);
463  }
464  setring R;
465 list truncatedcomplex=imap(W,truncatedcomplex);
466/*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex,
467  i.e. D^i=C^truncatedcomplex[2*i-1] and d^i=truncatedcomplex[2*i]*/
468 if (diffforms==0)
469   {
470     list derhamhom=findCohomology(truncatedcomplex,le);
471     option(set,saveoptions);
472     return (derhamhom);
473   }
474 list outall=findCohomologyDiffForms(truncatedcomplex,le);
475 setring W;
476 list dimanddiff=imap(R,outall);
477 list alldiffforms=dimanddiff[2];
478 while(size(alldiffforms)<size(passendespoly))
479   {
480     passendespoly=delete(passendespoly,1);
481   }
482 list newdiffforms;
483 matrix Diff;
484 for (i=1; i<=size(alldiffforms); i++)
485   {
486     newdiffforms[i]=list();
487     for (j=1; j<=size(alldiffforms[i]); j++)
488       {
489         Diff=matrix(0,1,newcomplex[3*(i+size(newcomplex) div 3 - size(alldiffforms))-2]);
490         for (k=1; k<=ncols(alldiffforms[i][j]); k++)
491           {
492             if (alldiffforms[i][j][1,k]!=0)
493               {
494                 Diff[1,passendespoly[i][k][2]]=Diff[1,passendespoly[i][k][2]]+alldiffforms[i][j][1,k]*passendespoly[i][k][1];
495               }
496           }
497         newdiffforms[i][j]=Diff;
498       }
499   }
500 list omegacomplex=makeOmega(nvars(W) div 2);
501 list newcomplexmod;
502 for (i=1; i<=size(newcomplex) div 3; i++)
503   {
504     newcomplexmod[2*i-1]=newcomplex[3*i-2];
505     newcomplexmod[2*i]=newcomplex[3*i];
506   }
507 while (size(dimanddiff[1])<size(newcomplexmod) div 2)
508   {
509     newcomplexmod=delete(newcomplexmod,1);
510     newcomplexmod=delete(newcomplexmod,1);
511   }
512 while (size(dimanddiff[1])<size(quasiiso))
513   {
514     quasiiso=delete(quasiiso,1);
515   }
516 while (size(dimanddiff[1])>size(generators))
517   {
518     generators=insert(generators,list());
519   }
520 while (size(dimanddiff[1])>size(quasiiso))
521   {
522     quasiiso=insert(quasiiso,list());
523   }
524 int keepsign;
525 list derhamdiff;
526 list doublecom=makeDoubleComplex(newcomplexmod,omegacomplex,quasiiso,generators);
527 matrix diffform;
528 int stopping;
529 int p;
530 matrix convert;
531 list interim;
532 list correspondingposition;
533 list allforms=list();
534 for (i=1; i<=size(newdiffforms); i++)
535   {
536     derhamdiff[i]=list();
537     allforms[i]=list();
538     for (j=1; j<=size(newdiffforms[i]); j++)
539       {
540         allforms[i][j]=list();
541         keepsign=1;
542         derhamdiff[i][j]=0;
543         diffform=newdiffforms[i][j];//Zeilenform
544         correspondingposition=doublecom[i][1];//needed fpr transformation process
545         interim=transferDiffforms(diffform,correspondingposition);
546         if (size(interim)!=0)
547           {
548             allforms[i][j][size(allforms[i][j])+1]=interim;
549           }
550         stopping=0;
551         p=1;
552         for (k=i; k<=size(newdiffforms); k++)
553           {
554             keepsign=(-1)*keepsign;
555             if (stopping==0)
556               {
557                 if (size(doublecom[k][p][2])==0)
558                   {
559                     stopping=1;
560                   }
561                 else
562                   {
563                     if (size(doublecom[k+1][p][3])!=0)
564                       {
565                         diffform=diffform*doublecom[k][p][2];//Spaltenform
566                         if (diffform!=matrix(0,nrows(diffform),ncols(diffform)))
567                           {
568                              diffform=findPreimage(doublecom[k+1][p][3],transpose(diffform));//Zeilenform
569                             correspondingposition=doublecom[k+1][p+1];//needed for transformation process
570                             interim=transferDiffforms(keepsign*diffform,correspondingposition);
571                             if (size(interim)!=0)
572                               {
573                                 allforms[i][j][size(allforms[i][j])+1]=interim;
574                               }
575                             p=p+1;
576                           }
577                         else
578                           {
579                             stopping=1;
580                           }
581                       }
582                     else
583                       {
584                         stopping=1;
585                       }
586                   }
587               }
588           }
589       }
590   }
591 setring R;
592 list allforms=fetch(W,allforms);
593 option(set,saveoptions);
594 return (allforms);
595}
596
597example
598{ "EXAMPLE:";
599  ring r = 0,(x,y,z),dp;
600  list L=(xy,xz);
601  deRhamCohomology(L);
602}
603
604////////////////////////////////////////////////////////////////////////////////////
605//COMPUTATION OF THE MAYER-VIETORIS COMPLEX
606////////////////////////////////////////////////////////////////////////////////////
607
608proc MVComplex(list L,list #)
609"USAGE:MVComplex(L); L a list of  polynomials
610ASSUME: -Basering is a polynomial ring with n vwariables and rational coefficients
611        -L is a list of non-constant polynomials
612RETURN: ring W: the nth Weyl algebra @*
613        W contains a list MV, which represents the Mayer-Vietrois complex (C^i,d^i) of the
[7fe9f8b]614        polynomials contained in L as follows:@*
615        the C^i are given  by D_n^ncols(C[2*i-1])/im(C[2*i-1]) and the differentials
616        d^i are given by C[2*i]
617EXAMPLE:example MVComplex; shows an example
618"
619{
620  /* We follow algorithm 3.2.5 in [R],if #!=0 we use also  Remark 3.2.6 in [R] for
[0e8a5a]621     an additional iterative localization*/
[7fe9f8b]622  def R=basering;
623  int i;
624  int iterative=1;
625  if (size(#)!=0)
[0e8a5a]626    {
627      iterative=#[1];
628    }
629  for (i=1; i<=size(L); i++)
630    {
631      if (L[i]==0)
632        {
633          print("localization with respect to 0 not possible");
634          return();
635        }
636      if (leadcoef(L[i])-L[i]==0)
637        {
638          print("polynomials must be non-constant");
639          return();
640        }
641    }
642  if (iterative==1)
643    {
644      /*compute the localizations by factorizing the polynomials and iterative
645        localization of the factors */
646      for (i=1; i<=size(L); i++)
647        {
648          L[i]=factorize(L[i],1);
649        }
650    }
651  int r=size(L);
652  int n=nvars(basering);
653  int le=size(L)+n;
654  /*construct the ring Ws*/
655  def W=makeWeyl(n);
656  setring W;
657  list man=ringlist(W);
658  if (n==1)
659    {
660      man[2][1]="x(1)";
661      man[2][2]="D(1)";
662      def Wi=ring(man);
663      setring Wi;
664      kill W;
665      def W=Wi;
666      setring W;
667      list man=ringlist(W);
668    }
669  man[2][size(man[2])+1]="s";;
670  man[3][3]=man[3][2];
671  man[3][2]=list("dp",intvec(1));
672  matrix N=UpOneMatrix(size(man[2]));
673  man[5]=N;
674  matrix M[1][1];
675  man[6]=transpose(concat(transpose(concat(man[6],M)),M));
676  def Ws=ring(man);
677  setring Ws;
678  int j,k,l,c;
679  list L=fetch(R,L);
680  list Cech;
681  ideal J=var(1+n);
682  for (i=2; i<=n; i++)
683    {
684      J=J,var(i+n);
685    }
686  Cech[1]=list(J);
687  list Theta, remminroots;
688  Theta[1]=list(list(list(),1,1));
689  list rem,findminintroot,diffmaps;
690  int minroot,st,sk;
691  intvec k1;
692  poly fred,forfetch;
693  matrix subm;
694  int rmr;
695  if (iterative==0)
696    {/*computation of the modules of the MV complex*/
697      for (i=1; i<=r; i++)
698        {
699          findminintroot=list();
700          Cech[i+1]=list();
701          Theta[i+1]=list();
702          k1=1;
703          for (j=1; j<=i; j++)
704            {
705              k1[size(k1)+1]=size(Theta[j+1]);
706              for (k=1; k<=k1[j]; k++)
707                {
708                  Theta[j+1][size(Theta[j+1])+1]=list(Theta[j][k][1]+list(i));
709                  Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i];
710                  /*We compute the s-parametric annihilator J(s)  and the b-function
711                    of the polynomial L[i] and Cech[i][k] to localize the module
712                    D_n/(D(1),...,D(n))[L[i]^(-1)]\otimes D_n^c/im(Cech[i][k]),
713                    where c=ncols(Cech[i][k]) and the im(Cech[i][k]) is generated by
714                    the rows of the matrix.
715                    If we plug the minimal integer root r(or a smaller integer
716                    value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic to
717                    the above localization*/
718                  rem=SannfsIBM(L[i],Cech[j][k]);
719                  Cech[j+1][size(Cech[j+1])+1]=rem[1];
720                  findminintroot[size(findminintroot)+1]=rem[2];
721                }
722            }
723          /* we compute the minimal root of all b-functions of L[i] computed above,
724             because we want to plug in the same root r in all s-parametric
725             annihilators we computed for L[i]  ->this will ensure  we can compute
726             the maps of the MV complex*/
[76d26c]727          minroot=minIntRootD(findminintroot);
[0e8a5a]728          for (j=1; j<=i; j++)
729            {
730              for (k=1; k<=k1[j]; k++)
731                {
732                  sk=size(Cech[j+1])+1-k;
733                  Cech[j+1][size(Cech[j+1])+1-k]=subst(Cech[j+1][sk],s,minroot);
734                }
735            }
736          remminroots[i]=minroot;
737        }
738      Cech=delete(Cech,1);
739      Theta=delete(Theta,1);
740      list zw;
741      poly reme;
742      /*computation of the maps of the MV complex*/
743      for (i=1; i<r; i++)
744        {
745          diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
746          for (j=1; j<=size(Cech[i]); j++)
747            {
748              for (k=1; k<=size(Cech[i+1]); k++)
749                {
750                  zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
751                  if (zw[2]!=0)
752                    {
753                      rmr=-remminroots[zw[1]];
754                      reme=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
755                      zw[2]=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
756                      diffmaps[i][j,k]=zw[2];
757                    }
758                }
759            }
760        }
761      diffmaps[r]=matrix(0,1,1);
762    }
763  list generators;
764  if (iterative==1)
765    {
766      for (i=1; i<=r;i++)
767        {
768          generators[i]=list();////////////////////////////////////////////////////////////////////
769          Cech[i+1]=list();
770          Theta[i+1]=list();
771          k1=1;
772          for (c=1; c<=size(L[i]); c++)
773            {
774              findminintroot=list();
775              for (j=1; j<=i; j++)
776                {
777                  if (c==1)
778                    {
779                      k1[size(k1)+1]=size(Theta[j+1]);
780                    }
781                  for (k=1; k<=k1[j]; k++)
782                    {
783                      /*We compute the s-parametric annihilator J(s)  und the b-
784                        function of the polynomial L[i][c] and Cech[i][k] to
785                        localize the module D_n/(D(1),...,D(n))[L[i][c]^(-1)]\otimes
786                        D_n^c/im(Cech[i][k]), where c=ncols(Cech[i][k]).
787                        If we plug the minimal integer root r(or a smaller integer
788                        value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic
789                        to the above localization*/
790                      if (c==1)
791                        {
792                          rmr=size(Theta[j+1])+1;
793                          Theta[j+1][rmr]=list(Theta[j][k][1]+list(i));
794                          Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i][c];
795                          rem=SannfsIBM(L[i][c],Cech[j][k]);
796                          Cech[j+1][size(Cech[j+1])+1]=rem[1];
797                          findminintroot[size(findminintroot)+1]=rem[2];
798                        }
799                      else
800                        {
801                          st=size(Theta[j+1])-k1[j]+k;
802                          Theta[j+1][st][2]=Theta[j+1][st][2]*L[i][c];
803                          rem=SannfsIBM(L[i][c],Cech[j+1][size(Cech[j+1])-k1[j]+k]);
804                          Cech[j+1][size(Cech[j+1])-k1[j]+k]=rem[1];
805                          findminintroot[size(findminintroot)+1]=rem[2];
806                        }
807                    }
808                }
809                /* we compute the minimal root of all b-functions of L[i][c]
810                   computed above,because we want to plug in the same root r in all
811                   s-parametric annihilators we computed for L[i]  ->this will
812                   ensure  we can compute the maps of the MV complex*/
[76d26c]813              minroot=minIntRootD(findminintroot);
[0e8a5a]814              for (j=1; j<=i; j++)
815                {
816                  for (k=1; k<=k1[j]; k++)
817                    {
818                      st=size(Cech[j+1])+1-k;
819                      Cech[j+1][st]=subst(Cech[j+1][st],s,minroot);
820                    }
821                }
822              if (c==1)
823                {
824                  remminroots[i]=list();
825                }
826              remminroots[i][c]=minroot;
827            }
828        }
829      Cech=delete(Cech,1);
830      Theta=delete(Theta,1);
831      list zw;
832      poly reme;
833      /*maps of the MV Complex*/
834      for (i=1; i<r; i++)
835        {
836          diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
837          for (j=1; j<=size(Cech[i]); j++)
838            {
839              for (k=1; k<=size(Cech[i+1]); k++)
840                {
841                  zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
842                  if (zw[2]!=0)
843                    {
844                      reme=1;
845                      for (c=1; c<=size(L[zw[1]]);c++)
846                        {
847                          reme=reme*L[zw[1]][c]^(-remminroots[zw[1]][c]);
848                        }
849                      diffmaps[i][j,k]=zw[2]*reme;
850                    }
851                }
852            }
853        }
854      diffmaps[r]=matrix(0,1,1);
855      for (i=1; i<=r; i++)
856        {
857          for (j=1; j<=size(Theta[i]); j++)
858            {
859              generators[i][j]=1;
860              for (c=1; c<=size(Theta[i][j][1]); c++)
861                {
862                  for (k=1; k<=size(L[Theta[i][j][1][c]]); k++)
863                    {
864                      generators[i][j]=generators[i][j]*L[Theta[i][j][1][c]][k]^((-1)*remminroots[Theta[i][j][1][c]][k]);
865                    }
866                }
867            }
868        }
869    }
870  setring W;
871  /*map the modules and maps to the Weyl algebra*/
872  list diffmaps=imap(Ws,diffmaps);
873  list Cechmodules=imap(Ws,Cech);
874  if (iterative==1)
875    {
876      list Theta=imap(Ws,Theta);
877      list generators=imap(Ws,generators);
878    }
879  list Cech;
880  matrix sup;
881  for (i=1; i<=r; i++)
882    {
883      sup=transpose(matrix(Cechmodules[i][1]));
884      Cech[2*i-1]=sup;
885      for (j=2; j<=size(Cechmodules[i]); j++)
886        {
887          sup=transpose(matrix(Cechmodules[i][j]));
888          Cech[2*i-1]=dsum(Cech[2*i-1],sup);
889        }
890      sup=matrix(diffmaps[i]);
891      Cech[2*i]=sup;
892    }
893  list MV=Cech;
894  if (iterative==1)
895    {
896      export Theta;
897      export generators;
898    }
899  export MV;
900
901  return (W);
902}
903
904example
905{ "EXAMPLE:";
906  ring r = 0,(x,y,z),dp;
907  list L=xy,xz;
908  def C=MVComplex(L);
909  setring C;
910  MV;
911}
912
913////////////////////////////////////////////////////////////////////////////////////
914
915static proc SannfsIBM(poly F,ideal myJ)
916"USAGE: SannfsIBM(f,J), F poly, J ideal
917ASSUME: basering is D_n[s], where D_n is the Weyl algebra and s and extra
918        commutative variable@*
919        f is a polynomial in the variables x(1),...,x(n) with rational coefficients
920        @*
921        J is holonomic and f-saturated
922RETURN  AlList of the form (K,g), where K is an ideal and g a univariant polynomial
923        in  the variable s. K is the s-parametric annihilator of F and J and g is
924        the b-function of F and J.
925"
926{
927  /*modified version of the procedure SannfsBM from the library dmod.lib: SannfsBM
928    computes the s-parametric annihilator for J=(x_1,...,x_n)*/
929  /* We use Algorithm 3.1.12 in[R] to compute the s-parametric
930     annihilator. Then we use the s-parametric annihilator to compute the b-function
931     via Algorithm 4.7 in [W1].*/
932  /* We assume that the basering the the nth Weyl algebra D_n. We create the ring
933     D_n[s,t], where t*s=s*t-t*/
934  def save = basering;
935  int N = nvars(basering)-1;
936  int Nnew = N+2;
937  int i,j;
938  string s;
939  list RL = ringlist(basering);
940  list L, Lord;
941  list tmp;
942  intvec iv;
943  L[1] = RL[1];
944  L[4] = RL[4];
945  list Name  = RL[2];
946  Name=delete(Name,size(Name));
947  list RName;
948  RName[1] = "t";
949  RName[2] = "s";
950  list DName;
951 for(i=1;i<=N div 2;i++)
952  {
953    DName[i] = var(N div 2+i);
954    Name=delete(Name,N div 2+1);
955  }
956  tmp[1] = "t";
957  tmp[2] = "s";
958  list NName = tmp +Name+DName;
959  L[2]   = NName;
960  kill NName;
961  tmp[1]  = "lp";
962  iv      = 1,1;
963  tmp[2]  = iv;
964  Lord[1] = tmp;
965  tmp[1]  = "dp";
966  s       = "iv=";
967  for(i=1;i<=Nnew;i++)
[7fe9f8b]968  {
[0e8a5a]969    s = s+"1,";
[7fe9f8b]970  }
[0e8a5a]971  s[size(s)]= ";";
972  execute(s);
973  kill s;
974  tmp[2]    = iv;
975  Lord[2]   = tmp;
976  tmp[1]    = "C";
977  iv        = 0;
978  tmp[2]    = iv;
979  Lord[3]   = tmp;
980  tmp       = 0;
981  L[3]      = Lord;
982  def @R@ = ring(L);
983  setring @R@;
984  matrix @D[Nnew][Nnew];
985  @D[1,2]=t;
986  for(i=1; i<=N div 2; i++)
[7fe9f8b]987  {
[0e8a5a]988    @D[2+i, N div 2+2+i]=1;
989  }
990  def @R = nc_algebra(1,@D);
991  setring @R;
992  kill @R@;
993  /*we start with the computation of the s-parametric annihilator*/
994  poly  F = imap(save,F);
995  ideal myJ=imap(save,myJ);
996  for (i=1; i<=N div 2; i++)
[7fe9f8b]997    {
[0e8a5a]998      myJ=subst(myJ,D(i),D(i)+diff(F,x(i))*t);
[7fe9f8b]999    }
[0e8a5a]1000  ideal I = t*F+s;
1001  I=I,myJ;//the s-parametric annihilator in D_n[s,t]
1002  /*we compute the intersection of I and D_n[s]*/
1003  ideal J = slimgb(I);
1004  ideal K = nselect(J,1);
1005  K = slimgb(K);//the s-parametric annihilator
1006  /*we use K to compute the b-function*/
1007  ideal B=K,F;
1008  B=slimgb(B);
1009  vector p=pIntersect(s,B);
1010  poly f=vec2poly(p,2);
1011  setring save;
1012  poly f=imap(@R,f);
1013  ideal K=imap(@R,K);
1014  return (list(K,f));
1015}
1016
1017////////////////////////////////////////////////////////////////////////////////////
1018//COMPUTATION OF A QUASI-ISOMORPHIC V_D-STRICT FREE COMPLEX
1019////////////////////////////////////////////////////////////////////////////////////
1020
1021static proc quasiisomorphicVdComplex(list L,list #)
1022"USAGE: quasiisomorphicVdComplex(L[,df]); L a list of the form (M_1,f_1,...,M_s,f_s),
1023        where the M_i and f_i are matrices
1024ASSUME: Basering is the Weyl algebra D_n @*
1025        (M_1,f_1,...,M_s,f_s) represents a complex 0->D_n^(r_1)/im(M_1)->
1026        D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i)
1027        is generated by the rows of M_i. In particular it hold:@*
1028        - The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @*
1029        -the image of M_1*f_i is contained in the image of M_(i+1) @*
1030        d is an integer between 1 and n. If no value for d is given, it is assumed
1031        to be n @*
1032        df is an optional int, if df equals 1 a \tilde(V_d)-strict complex
1033        will be computed (instead of a V_d-strict one) (for a definition see [W3])
1034RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @*
1035        L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,i_s,g_s,m_s) such that:@*
1036        -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs
1037         of size i_j@*
1038        -D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0  is a V_d-strict complex
1039         with differentials m_i that is quasi-isomorphic to the complex given by L@*
1040        L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and
1041        the n_i are shift vectors such that:@*
1042        -coker(H_i) is the ith cohomology group of the complex given by L_1@*
1043        -the n_i are the shift vectors of the coker(H_i)
1044THEORY: We follow Proposition 3.2 and Corollary 3.3 in [W3]
1045"
1046{
1047  int tilde;
1048  if (size(#)!=0)
[7fe9f8b]1049    {
[0e8a5a]1050      tilde=#[1];
[7fe9f8b]1051    }
[0e8a5a]1052  def B=basering;
[08fa62]1053  int n=nvars(B) div 2 + 1;//+1 muesste stimmen! bitte kontrollieren!
[0e8a5a]1054  int d=nvars(B) div 2;
1055  int r=size(L) div 2;
1056  int lonc=n+r;
1057  int Kiold=0;
1058  matrix kerold;
1059  // matrix kernew=out[r][2][2];
1060  matrix kernew=diag(1,ncols(L[size(L)-1]));
1061  module mL;
1062  int i;
1063  int k;
1064  matrix testm;
1065  int Kinew=nrows(kernew);
1066  int Jiold=0;
1067  int Jinew=0;
1068  matrix Niold;
1069  matrix Ninew;
1070  list newcomplex;
1071  int Aiold=Kinew;
1072  matrix savediv;
1073  newcomplex[3*lonc-2]=Kinew;
1074  newcomplex[3*lonc-1]=intvec(0:Kinew);
1075  newcomplex[3*lonc]=matrix(0,Kinew,1);
1076  list quasiiso;
1077  quasiiso[lonc]=diag(1,Kinew);
1078  matrix invimage;
1079  matrix keralpha;
1080  intvec v;
1081  int j;
1082  matrix sc;
1083  matrix fnc;
1084  int indk;
1085  int indj;
1086  int Aiold;
1087  list saveres;
1088  matrix Liplus;
1089  for (i=r-1; i>=0; i--)
1090    {
1091      indk=0;
1092      indj=0;
1093      Kiold=Kinew;
1094      kerold=kernew;
1095      if (i!=0)
1096        {
1097          // kernew=divdr(L[2*i],L[2*i+1],1);
1098          kernew=divdr(L[2*i],L[2*i+1]);
1099          mL=slimgb(transpose(L[2*i-1]));
1100          for (k=1; k<=nrows(kernew); k++)
1101            {
1102              testm=reduce(transpose(submat(kernew,k,intvec(1..ncols(kernew)))),mL);
1103              if (testm==matrix(0,nrows(testm),ncols(testm)))
1104                {
1105                  kernew=transpose(deletecol(transpose(kernew),k));
1106                  k=k-1;
1107                }
1108            }
1109          Kinew=nrows(kernew);
1110          if (kernew==matrix(0,nrows(kernew),ncols(kernew)))
1111            {
1112              Kinew=0;
1113              indk=1;
1114            }
1115        }
1116      else
1117        {
1118          Kinew=0;
1119          indk=1;
1120        }
1121      Jiold=Jinew;
1122      Niold=Ninew;
1123      keralpha=transpose(syz(transpose(newcomplex[3*(i+n)+3])));
1124      if (i!=0)
1125        {
1126          invimage=divdr(quasiiso[n+i+1],transpose(concat(transpose(L[2*i]),transpose(L[2*i+1]))));
1127          Ninew=vdStrictIntersect(keralpha,invimage,newcomplex[3*(n+i+1)-1],tilde);//////////////
1128        }
1129      else
1130        {
1131          invimage=divdr(quasiiso[n+i+1],L[2*i+1]);
1132          saveres=vdStrictIntersectPlus(keralpha,invimage,newcomplex[3*(n+i+1)-1],tilde);////////////////////////
1133
1134          ///////////////////BIS HIERHIN VERALLGEMEINERT////////////////////////////////////////////////////////////////////
1135
1136
1137          Ninew=saveres[1];
1138        }
1139      Jinew=nrows(Ninew);
1140      if (Ninew==matrix(0,nrows(Ninew),ncols(Ninew)))
1141        {
1142          Jinew=0;
1143          indk=1;
1144        }
1145      newcomplex[3*(n+i)-2]=Kinew+Jinew;
1146      v=0;
1147      if (indk==0)
1148        {
1149          v=(0:Kinew);
1150          if (indj==0)
1151            {
1152              fnc=transpose(concat(transpose(matrix(0,Kinew,Kiold+Jiold)),transpose(Ninew)));
1153            }
1154          else
1155            {
1156              fnc=matrix(0,Kinew,Kiold+Jiold);
1157            }
1158        }
1159      else
1160        {
1161          if (indj==0)
1162            {
1163              fnc=Ninew;
1164            }
1165          else
1166            {
1167              fnc=matrix(0,1,Kiold+Jiold);
1168              newcomplex[3*(n+i)-2]=1;
1169            }
1170        }
1171      Aiold=Jinew+Kinew;
1172      if (Aiold==0)
1173        {
1174          Aiold=1;
1175        }
1176      newcomplex[3*(n+i)]=fnc;
1177      for (j=1; j<=Jinew; j++)
1178        {
1179          if (tilde==0)
1180            {
1181              v[Kinew+j]=VdDeg(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]);
1182            }
1183          else
1184            {
1185              v[Kinew+j]=VdDegTilde(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]);
1186            }
1187        }
1188      newcomplex[3*(n+i)-1]=v;
1189      if (i==0)
1190        {
1191          quasiiso[n+i]=matrix(0,Jinew,1);
1192        }
1193      else
1194        {
1195          if (indj==0)
1196            {
1197              sc=submat(fnc,intvec(Kinew+1..nrows(fnc)),intvec(1..ncols(fnc)))*quasiiso[n+i+1];
1198              Liplus=transpose(concat(transpose(L[2*i]),transpose(L[2*i+1])));
1199              sc=matrixLift(Liplus,sc);//stimmt das jetzt
1200              sc=submat(sc,intvec(1..nrows(sc)),intvec(1..nrows(L[2*i])));
1201              if (indk==0)
1202                {
1203                  //pi=kernew
1204                  quasiiso[n+i]=transpose(concat(transpose(kernew),transpose(sc)));
1205                }
1206              else
1207                {
1208                  quasiiso[n+i]=sc;
1209                }
1210            }
1211          else
1212            {
1213              if (indk==0)
1214                {
1215                  quasiiso[n+i]=kernew;
1216                }
1217              else
1218                {
1219                  quasiiso[n+i]=matrix(0,1,ncols(kernew));
1220                }
1221            }
1222        }
1223    }
1224  for (i=1; i<=n-1; i++)
[7fe9f8b]1225    {
[0e8a5a]1226      quasiiso[n-i]=list();
1227      if (size(saveres[2][i])!=0)
1228        {
1229          newcomplex[3*(n-i)]=saveres[2][i];
1230          newcomplex[3*(n-i)-2]=nrows(saveres[2][i]);
1231          v=0;
1232          for (j=1; j<=newcomplex[3*(n-i)-2]; j++)
1233            {
1234              if (tilde==0)
1235                {
1236                  v[j]=VdDeg(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(n-i)+2]);
1237                }
1238              else
1239                {
1240                  v[j]=VdDegTilde(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(n-i)+2]);
1241                }
1242            }
1243          newcomplex[3*(n-i)-1]=v;
1244        }
1245      else
1246        {
1247          newcomplex[3*(n-i)]=matrix(0,1,1);
1248          if (newcomplex[3*(n-i)+1]!=0)
1249            {
1250              newcomplex[3*(n-i)]=matrix(0,1,newcomplex[3*(n-i)+1]);
1251            }
1252          newcomplex[3*(n-i)-2]=int(0);
1253          newcomplex[3*(n-i)-1]=intvec(0);
1254        }
[7fe9f8b]1255    }
[0e8a5a]1256  list result;
1257  result[1]=newcomplex;
1258  result[2]=list();
1259  list forsep;
1260  for (i=1; i<=size(L) div 2+1; i++)
1261    {
1262      forsep[2*i]=newcomplex[3*(n+i-1)];
1263      forsep[2*i-1]=matrix(0,1,nrows(forsep[2*i]));
1264    }
1265  forsep=shortExactPieces(forsep);
1266  list listofHis;
1267  matrix forVd;
1268  for (i=1; i<=size(L) div 2; i++)
1269    {
1270      v=0;
1271      listofHis[i]=list(forsep[i+1][1][5]);
1272      forVd=forsep[i+1][2][2];
1273      for (j=1; j<=nrows(forVd); j++)
1274        {
1275          if (tilde==0)
1276            {
1277              v[j]=VdDeg(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)-1]);
1278            }
1279          else
1280            {
1281              v[j]=VdDegTilde(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)-1]);
1282            }
1283        }
1284      listofHis[i][2]=v;
1285    }
1286  result[2]=listofHis;
1287  result[3]=quasiiso;
1288  return(result);
1289}
1290
1291////////////////////////////////////////////////////////////////////////////////////
1292
1293static proc vdStrictIntersect(matrix M, matrix N, intvec v, int tilde)
1294{
1295  def B=basering;
1296  option(returnSB);//                    alternative:erst intersect und dann SB-Berechung mit slimgb
1297  if (tilde==0)
1298    {
1299      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v);
1300    }
1301  else
1302    {
1303      def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v);
1304    }
1305  setring HomWeyl;
1306  matrix M=fetch(B,M);
1307  matrix N=fetch(B,N);
1308  M=nHomogenize(M);
1309  N=nHomogenize(N);
1310  matrix vdintersection=transpose(intersect(transpose(M),transpose(N)));
1311  vdintersection=subst(vdintersection,h,1);
1312  setring B;
1313  matrix vdintersection=fetch(HomWeyl,vdintersection);
1314  option(noreturnSB);
1315  return(vdintersection);
1316}
1317
1318////////////////////////////////////////////////////////////////////////////////////
1319
1320static proc vdStrictIntersectPlus(matrix M, matrix N, intvec v, int tilde)
1321{
1322  def B=basering;
1323  int n=nvars(B) div 2;
1324  matrix vdint=transpose(intersect(transpose(M),transpose(N)));
1325  if (tilde==0)
1326    {
1327      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v);
1328    }
1329  else
1330    {
1331      def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v);
1332    }
1333  setring HomWeyl;
1334  matrix vdint=fetch(B,vdint);
1335  matrix N=fetch(B,N);
1336  vdint=nHomogenize(vdint);
1337  intvec i1;
1338  intvec i2;
1339  int i;
1340  int nr;
1341  int nc;
1342  def ringofSyz=Sres(transpose(vdint),n);////////////////////////////////////////////////////////////////
1343  setring ringofSyz;
1344  matrix vdint=transpose(matrix(RES[2]));
1345  vdint=subst(vdint,h,1);
1346  int logens=ncols(vdint)+1;
1347  int omitemptylist;
1348  matrix zerom;
1349  list rofA;
[08fa62]1350  for (i=3; i<=n+3; i++)////////////////////////////////////////////////////////////////////////////n und si muessen noch definiert werden
[0e8a5a]1351    {
1352      if (size(RES)>=i)
1353        {
1354          zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
1355          if (RES[i]!=zerom)
1356            {
1357              rofA[i-2]=(matrix(RES[i]));
1358              if (i==3)
1359                {
1360                  if (nrows(rofA[i-2])-logens+1!=nrows(vdint))
1361                    {
1362                      //build the resolution
1363                      nr=nrows(vdint)+logens-1;
1364                      nc=ncols(rofA[i-2]);
1365                      rofA[i-2]=matrix(rofA[i-2],nr,nc);
1366                    }
1367
1368                }
1369              if (i!=3)
1370                {
1371                  if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3]))
1372                    {
1373                      nr=nrows(rofA[i-3])+logens-1;
1374                      nc=ncols(rofA[i-2]);
1375                      rofA[i-2]=matrix(rofA[i-2],nr,nc);
1376                    }
1377                }
1378              i1=intvec(logens..nrows(rofA[i-2]));
1379              i2=intvec(1..ncols(rofA[i-2]));
1380              rofA[i-2]=transpose(submat(rofA[i-2],i1,i2));
1381              logens=logens+ncols(rofA[i-2]);
1382              rofA[i-2]=subst(rofA[i-2],h,1);
1383            }
1384          else
1385            {
1386              rofA[i-2]=list();
1387            }
1388        }
1389      else
1390        {
1391          rofA[i-2]=list();
1392        }
1393    }
1394  if(size(rofA[1])==0)
[7fe9f8b]1395    {
[0e8a5a]1396      omitemptylist=1;
[7fe9f8b]1397    }
[0e8a5a]1398  setring B;
1399  vdint=fetch(ringofSyz,vdint);
1400  if (omitemptylist!=1)
1401    {
1402      list rofA=fetch(ringofSyz,rofA);
1403    }
1404  kill HomWeyl;
1405  kill ringofSyz;
1406  return(list(vdint,rofA));
[7fe9f8b]1407}
1408
1409////////////////////////////////////////////////////////////////////////////////////
1410
1411static proc toVdStrictFreeComplex(list L,string Syzstring,list #)
[0e8a5a]1412"USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form
1413        (M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a
[7fe9f8b]1414        string, d an optional integer
1415ASSUME: Basering is the Weyl algebra D_n @*
1416        (M_1,f_1,...,M_s,f_s) represents a complex 0->D_n^(r_1)/im(M_1)->
[0e8a5a]1417        D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i)
[7fe9f8b]1418        is generated by the rows of M_i. In particular it hold:@*
1419        - The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @*
1420        -the image of M_1*f_i is contained in the image of M_(i+1) @*
[0e8a5a]1421        d is an optional integer which indices in the case size(L)=2, whether a
1422        V_d-strict or \tilde(V_d)-strict will be computed@*
[7fe9f8b]1423        Syzstring is either: @*
[0e8a5a]1424        -'Sres' (computes the resolutions and Groebner bases in the homogenized
1425         Weyl algebra using Schreyer's method)@*
[7fe9f8b]1426        or @*
[0e8a5a]1427        -'Vdres' (computes the resolutions via V_d-homogenization and without
[7fe9f8b]1428         Schreyer's method)@*
1429RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @*
1430        L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,i_s,g_s,m_s) such that:@*
[0e8a5a]1431        -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs
[7fe9f8b]1432         of size i_j@*
1433        -D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0  is a V_d-strict complex
1434         with differentials m_i that is quasi-isomorphic to the complex given by L@*
[0e8a5a]1435        L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and
[7fe9f8b]1436        the n_i are shift vectors such that:@*
[0e8a5a]1437        -coker(H_i) is the ith cohomology group of the complex given by L_1@*
[7fe9f8b]1438        -the n_i are the shift vectors of the coker(H_i)
1439THEORY: We follow Algorithm 3.8 in [W2]
1440"
1441{
[0e8a5a]1442  def B=basering;
[7fe9f8b]1443  int n=nvars(B) div 2+2;
1444  int d=nvars(B) div 2;
1445  intvec v;
1446  list out, outall;
1447  int i,j,k,indi,nc,nr;
1448  matrix mem;
1449  intvec i1,i2;
[0e8a5a]1450  int tilde;
[7fe9f8b]1451  if (size(#)!=0)
1452    {
[0e8a5a]1453      for (i=1; i<=size(#); i++)
[7fe9f8b]1454        {
[0e8a5a]1455          if (typeof(#[i])=="int")
1456            {
1457              tilde=#[i];
1458            }
[7fe9f8b]1459        }
1460    }
1461  /* If size(L)=2, our complex consists for only one non-trivial module.
[0e8a5a]1462     Therefore, we just have to compute a V_d-strict resolution of this module.*/
1463  if (size(L)==2)
1464    {
1465      v=(0:ncols(L[1]));
1466      out[3*n-1]=v;
1467      out[3*n-2]=ncols(L[1]);
1468      out[3*n]=L[2];
1469      if (Syzstring=="Vdres")
[7fe9f8b]1470        {
[0e8a5a]1471          /*if Syzstring="Vdres", we compute a V_d-strict Groebner basis of L[1]
1472            using F-homogenization (Prop. 3.9 in [OT]); then we compute the syzygies
1473            and make them V_d-strict using Prop  3.9[OT] and so on*/
1474          out[3*n-3]=VdStrictGB(L[1],d,v);
1475          for (i=n-1; i>=1; i--)
1476            {
1477              out[3*i-2]=nrows(out[3*i]);
1478              v=0;
1479              for (j=1; j<=out[3*i-2]; j++)
1480                {
1481                  mem=submat(out[3*i],j,intvec(1..ncols(out[3*i])));
1482                  v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector
1483                }
1484              out[3*i-1]=v;
1485              if (i!=1)
1486                {
1487                  /*next step in the resolution*/
1488                  out[3*i-3]=transpose(syz(transpose(out[3*i])));
1489                  if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3])))
1490                    {
1491                      /*makes the resolution V_d-strict*/
1492                      out[3*i-3]=VdStrictGB(out[3*i-3],d,out[3*i-1]);
1493                    }
1494                  else
1495                    {
1496                      /*resolution is already computed*/
1497                      out[3*i-3]=matrix(0,1,ncols(out[3*i-3]));
1498                      out[3*i-4]=intvec(0);
1499                      out[3*i-5]=int(0);
1500                      for (j=i-2; j>=1; j--)
1501                        {
1502                          out[3*j]=matrix(0,1,1);
1503                          out[3*j-1]=intvec(0);
1504                          out[3*j-2]=int(0);
1505                        }
1506                      break;
1507                    }
1508                }
1509            }
[7fe9f8b]1510        }
[0e8a5a]1511      else
[7fe9f8b]1512        {
[0e8a5a]1513          /*in the case Syzstring!="Vdres" we compute the resolution in the
1514            homogenized Weyl algebra using Thm 9.10 in[OT]*/
1515          if (tilde==0)
1516            {
1517              def HomWeyl=makeHomogenizedWeyl(d);
1518            }
1519          else
1520            {
1521              def HomWeyl=makeHomogenizedWeylTilde(d);
1522            }
1523          setring HomWeyl;
1524          list L=fetch(B,L);
1525          L[1]=nHomogenize(L[1]);
1526          list out=fetch(B,out);
1527          out[3*n-3]=L[1];
1528          /*computes a ring with a list RES; RES is a V_d-strict resolution of
1529            L[1]*/
1530          def ringofSyz=Sres(transpose(L[1]),d);
1531          setring ringofSyz;
1532          int logens=2;
1533          matrix mem;
1534          list out=fetch(HomWeyl,out);
1535          out[3*n-3]=transpose(matrix(RES[2]));
1536          out[3*n-3]=subst(out[3*n-3],h,1);
1537          for (i=n-1; i>=1; i--)
1538            {
1539              out[3*i-2]=nrows(out[3*i]);
1540              v=0;
1541              for (j=1; j<=out[3*i-2]; j++)
1542                {
1543                  mem=submat(out[3*i],j,intvec(1..ncols(out[3*i])));
1544                  if (tilde==0)
1545                    {
1546                      v[j]=VdDeg(mem,d, out[3*i+2]);
1547                    }
1548                  else
1549                    {
1550                      v[j]=VdDegTilde(mem,d, out[3*i+2]);
1551                    }
1552                }
1553              out[3*i-1]=v;//shift vector such that the resolution RES is V_d-strict
1554              if (i!=1)
1555                {
1556                  indi=0;
1557                  if (size(RES)>=n-i+2)
1558                    {
1559                      nr=nrows(matrix(RES[n-i+2]));
1560                      mem=matrix(0,nr,ncols(matrix(RES[n-i+2])));
1561                      if (matrix(RES[n-i+2])!=mem)
1562                        {
1563                          indi=1;
1564                          out[3*i-3]=(matrix(RES[n-i+2]));
1565                          if (nrows(out[3*i-3])-logens+1!=nrows(out[3*i]))
1566                            {
1567                              mem=out[3*i-3];
1568                              out[3*i-3]=matrix(mem,nrows(mem)+logens-1,ncols(mem));
1569                            }
1570                          mem=out[3*i-3];
1571                          i1=intvec(logens..nrows(mem));
1572                          mem=submat(mem,i1,intvec(1..ncols(mem)));
1573                          out[3*i-3]=transpose(mem);
1574                          out[3*i-3]=subst(out[3*i-3],h,1);
1575                          logens=logens+ncols(out[3*i-3]);
1576                        }
1577                    }
1578                  if(indi==0)
1579                    {
1580                      out[3*i-3]=matrix(0,1,nrows(out[3*i]));
1581                      out[3*i-4]=intvec(0);
1582                      out[3*i-5]=int(0);
1583                      for (j=i-2; j>=1; j--)
1584                        {
1585                          out[3*j]=matrix(0,1,1);
1586                          out[3*j-1]=intvec(0);
1587                          out[3*j-2]=int(0);
1588                        }
1589                      break;
1590                    }
1591                }
1592            }
1593          setring B;
1594          out=fetch(ringofSyz,out);//contains the V_d-strict resolution
1595          kill ringofSyz;
1596        }
1597      outall[1]=out;
1598      outall[2]=list(list(out[3*n-3],out[3*n-1]));
1599      return(outall);
1600    }
1601  /*case size(L)>2: We compute a quasi-isomorphic free complex following Alg 3.8 in
1602    [W2]*/
1603  /* We denote the complex given by L as (C^i,d^i).
1604     We start by computing in the proc shortExaxtPieces representations for the
1605     short exact sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1), where the B^i, Z^i
1606     and H^i are coboundaries, cocycles and cohomology groups, respectively.*/
1607  out=shortExactPieces(L);
[7fe9f8b]1608  list rem;
[0e8a5a]1609  /* shortExactpiecesToVdStrict makes the sequences B^i->Z^i->H^i and
1610     Z^i->C^i->B^(i+1) V_d-strict*/
[7fe9f8b]1611  rem=shortExactPiecesToVdStrict(out,d,Syzstring);
[0e8a5a]1612  /*VdStrictDoubleComplexes computes V_d-strict resolutions over the seqeunces from
1613    proc shortExactPiecesToVdstrict*/
[7fe9f8b]1614  out=VdStrictDoubleComplexes(rem[1],d,Syzstring);
1615  for (i=1;i<=size(out); i++)
[0e8a5a]1616    {
1617      rem[2][i][1]=out[i][1][5][1];
1618      rem[2][i][2]=out[i][1][8][1];
1619    }
1620  /* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences
1621     Z^i->C^i->B^(i+1)) together to obtain a Cartan-Eilenberg resolution of
1622     (C^i,d^i)*/
[7fe9f8b]1623  out=assemblingDoubleComplexes(out);
1624  /*the proc totalComplex takes the total complex of the double complex from the
[0e8a5a]1625    proc assemblingDoubleComplexes*/
[7fe9f8b]1626  out=totalComplex(out);
1627  outall[1]=out;
[0e8a5a]1628  outall[2]=rem[2];//contains the cohomology groups and their shift vectors
[7fe9f8b]1629  return (outall);
1630}
1631
1632////////////////////////////////////////////////////////////////////////////////////
1633
1634
1635static proc sublist(list L,int m,int n)
1636{
1637  list out;
1638  int i; int j;
1639  int count;
1640  for (i=m; i<=n; i++)
1641    {
[0e8a5a]1642      out[size(out)+1]=list();
1643      for (j=1; j<=size(L[i]); j++)
1644        {
1645          count=count+1;
1646          out[size(out)][j]=list(L[i][j],count);
1647        }
[7fe9f8b]1648    }
1649  list o=list(out,count);
1650  return(o);
1651}
1652
1653////////////////////////////////////////////////////////////////////////////////////
1654
[0e8a5a]1655static proc LMSubset(list L,list M, list #)
[7fe9f8b]1656{
1657  int i;
1658  int j=1;
[0e8a5a]1659  if (size(#)==0)
[7fe9f8b]1660    {
[0e8a5a]1661      list position=(M[size(M)],(-1)^(size(L)));
1662    }
1663  else
1664    {
1665      list position=(M[size(M)],1);
[7fe9f8b]1666    }
[0e8a5a]1667  for (i=1; i<=size(L); i++)
1668    {
1669      if (L[i]!=M[j])
1670        {
1671          if (L[i]!=M[i+1] or j!=i)
1672            {
1673              return (L[i],0);
1674            }
1675          else
1676            {
1677              if (size(#)==0)
1678                {
1679                  position=(M[i],(-1)^(i-1));
1680                }
1681              else
1682                {
1683                  position=(M[i],(-1)^(size(L)+1-i));
1684                }
1685              j=j+1;
1686            }
1687        }
1688      j=j+1;
[7fe9f8b]1689
[0e8a5a]1690    }
[7fe9f8b]1691  return (position);
1692}
1693
1694////////////////////////////////////////////////////////////////////////////////////
1695
1696static proc shortExactPieces(list L)
1697{
1698  /*we follow Section 3.3 in [W2]*/
[0e8a5a]1699  /* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex  C=(C^i,d^i)
1700     as in the procedure toVdstrictcomplex*/
[7fe9f8b]1701  matrix  Bnew= divdr(L[2],L[3]);
1702  matrix Bold=Bnew;
1703  matrix Z=divdr(Bnew,L[1]);
1704  list bzh,zcb;
1705  bzh=list(list(),list(),Z,unitmat(ncols(Z)),Z);
1706  zcb=(Z, Bnew, L[1], unitmat(ncols(L[1])), Bnew);
1707  list sep;
1708  /* the list sep will be of size s such that
[0e8a5a]1709     -sep[i]=(sep[i][1],sep[i][2]) is a list of two lists
1710     -sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)->coker(Z^i)
1711      ->coker(H^i) represents the short exact seqeuence B^i(C)->Z^i(C)->H^i(C)
1712     -sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)->coker(C^i)->
1713      coker(B^(i+1)) represents the short exact seqeuence Z^i(C)->C^i->B^(i+1)(C)*/
[7fe9f8b]1714  sep[1]=list(bzh,zcb);
1715  int i;
1716  list out;
1717  for (i=3; i<=size(L)-2; i=i+2)
[0e8a5a]1718    {
1719      /*the proc bzhzcb computes representations for the short exact seqeunces */
1720      out=bzhzcb(Bold, L[i-1] , L[i], L[i+1], L[i+2]);
1721      sep[size(sep)+1]=out[1];
1722      Bold=out[2];
1723    }
[7fe9f8b]1724  bzh=(divdr(L[size(L)-2], L[size(L)-1]),L[size(L)-2], L[size(L)-1]);
1725  bzh[4]=unitmat(ncols(L[size(L)-1]));
1726  bzh[5]=transpose(concat(transpose(L[size(L)-2]),transpose(L[size(L)-1])));
1727  zcb=(L[size(L)-1], unitmat(ncols(L[size(L)-1])), L[size(L)-1],list(),list());
1728  sep[size(sep)+1]=list(bzh,zcb);
1729  return(sep);
1730}
1731
1732////////////////////////////////////////////////////////////////////////////////////
1733
1734static proc bzhzcb (matrix Bold,matrix f0,matrix C1,matrix f1,matrix C2)
1735{
1736  matrix Bnew=divdr(f1,C2);
1737  matrix Z= divdr(Bnew,C1);
1738  matrix lift1= matrixLift(Bnew,f0);
1739  matrix H=transpose(concat(transpose(lift1),transpose(Z)));
1740  list bzh=(Bold, lift1, Z, unitmat(ncols(Z)),H);
1741  list zcb=(Z, Bnew, C1, unitmat(ncols(C1)),Bnew);
1742  list out=(list(bzh, zcb), Bnew);
1743  return(out);
1744}
1745
1746////////////////////////////////////////////////////////////////////////////////////
1747
1748static proc shortExactPiecesToVdStrict(list C,int d,list #)
1749{/* We transform the short exact pieces from procedure shortExactPieces to V_d-
[0e8a5a]1750    strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in
1751    [W2].*/
[7fe9f8b]1752  /* If we compute our Groebner bases in the homogenized Weyl algebra, we already
[0e8a5a]1753     compute some resolutions it omit additional Groebner basis computations later
1754     on.*/
[7fe9f8b]1755  int s =size(C);int i; int j;
1756  string Syzstring="Sres";
[0e8a5a]1757  intvec v=0:ncols(C[s][1][5]);
[7fe9f8b]1758  if (size(#)!=0)
1759    {
[0e8a5a]1760      for (i=1; i<=size(#); i++)
1761        {
1762          if (typeof(#[i])=="string")
1763            {
1764              Syzstring=#[i];
1765            }
1766          if (typeof(#[i])=="intvec")
1767            {
1768               v=#[i];
1769            }
1770        }
[7fe9f8b]1771    }
1772  list out;
1773  list forout;
1774  if (Syzstring=="Vdres")
[0e8a5a]1775    {
1776      out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s));
1777    }
[7fe9f8b]1778  else
[0e8a5a]1779    {
1780      forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s);
1781      list resolutionofA=forout[9];
1782      list resolutionofC=forout[10];
1783      forout=delete(forout,10);
1784      forout=delete(forout,9);
1785      out[s]=list(forout);
1786      for (i=1; i<=size(resolutionofC); i++)
1787        {
1788          out[s][1][5][i+1]=resolutionofC[i];//save the resolutions
1789          out[s][1][1][i+1]=resolutionofA[i];
1790        }
1791    }
[7fe9f8b]1792  out[s][2]=list(list(out[s][1][3][1]));
1793  out[s][2][2]=list(unitmat(ncols(out[s][1][3][1])));
1794  out[s][2][3]=list(out[s][1][3][1]);
1795  out[s][2][4]=list(list());
1796  out[s][2][5]=list(list());
1797  out[s][2][6]=list(out[s][1][7][1]);
1798  out[s][2][7]=list(out[s][2][6][1]);
1799  out[s][2][8]=list(list());
1800  list resolutionofD;
1801  list resolutionofF;
1802  for (i=s-1; i>=2; i--)
[0e8a5a]1803    {
1804      C[i][2][5]=out[i+1][1][1][1];
1805      forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s);
1806      if (Syzstring=="Sres")
1807        {
1808          resolutionofD=forout[3];//save the resolutions
1809          resolutionofF=forout[4];
1810          forout=delete(forout,4);
1811          forout=delete(forout,3);
1812        }
1813      out[i]=forout;
1814      if(Syzstring=="Sres")
1815        {
1816          for (j=2; j<=size(out[i+1][1][1]); j++)
1817            {
1818              out[i][2][5][j]=out[i+1][1][1][j];
1819            }
1820          for (j=1; j<=size(resolutionofD);j++)
1821            {
1822              out[i][1][1][j+1]=resolutionofD[j];
1823              out[i][1][5][j+1]=resolutionofF[j];
1824            }
1825        }
[7fe9f8b]1826    }
1827  out[1]=list(list());//initalize our list
1828  C[1][2][5]=out[2][1][1][1];
1829  /*Compute the last V_d-strict seqeunce*/
1830  if (Syzstring=="Vdres")
1831    {
[0e8a5a]1832      out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
[7fe9f8b]1833    }
[0e8a5a]1834  else
[7fe9f8b]1835    {
[0e8a5a]1836      forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
1837      out[1][2]=delete(forout,9);
1838      list resolutionofA2=forout[9];
1839      for (i=1; i<=size(out[2][1][1]); i++)
1840        {
1841          /*put the modules for the resolutions in the right spot*/
1842          out[1][2][5][i]=out[2][1][1][i];
1843        }
1844      for (i=1; i<=size(resolutionofA2); i++)
1845        {
1846          out[1][2][1][i+1]=resolutionofA2[i];
1847        }
[7fe9f8b]1848    }
1849  out[1][1][3]=list(out[1][2][1][1]);
1850  out[1][1][5]=list(out[1][2][1][1]);
1851  out[1][1][4]=list(unitmat(ncols(out[1][1][3][1])));
1852  out[1][1][7]=list(out[1][2][6][1]);
1853  out[1][1][8]=list(out[1][2][6][1]);
1854  out[1][1][1]=list(list());
1855  out[1][1][2]=list(list());
1856  out[1][1][6]=list(list());
1857  if (Syzstring=="Sres")
1858    {
[0e8a5a]1859      for (i=1; i<=size(out[1][2][1]); i++)
1860        {
1861          out[1][1][3][i]=out[1][2][1][i];
1862          out[1][1][5][i]=out[1][2][1][i];
1863        }
[7fe9f8b]1864    }
1865  list Hi;
[0e8a5a]1866  for (i=1; i<=size(out); i++)
1867    {
1868      Hi[i]=list(out[i][1][5][1],out[i][1][8][1]);
1869    }
[7fe9f8b]1870  list outall;
1871  outall[1]=out;
1872  outall[2]=Hi;
1873  return(outall);
1874}
1875
1876////////////////////////////////////////////////////////////////////////////////////
1877
[0e8a5a]1878static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #)
[7fe9f8b]1879{
1880  /*this is the Algorithm 3.11 in [W2]*/
1881  int omitemptylist;
1882  int lengthofres=si+n-1;
1883  int i,j,logens;
1884  def B=basering;
1885  matrix bi=slimgb(transpose(C[5]));
[0e8a5a]1886  /* Computation of a V_d-strict Groebner basis of C[5]:
1887     -if Syzstring=="Vdres" this is done using the method of weighted homogenization
1888     (Prop. 3.9 [OT])
1889     -else we use the homogenized Weyl algebra for Groebner basis computations
1890     (Prop 9.9 [OT]),
1891     in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit
1892     extra Groebner basis computations later on*/
[7fe9f8b]1893  int nr,nc;
1894  intvec i1,i2;
1895  if (Syzstring=="Vdres")
1896    {
[0e8a5a]1897      if(size(#)==0)
1898        {
1899          matrix J_C=VdStrictGB(C[5],n,list(v));
1900        }
1901      else
1902        {
1903          matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
1904        }
[7fe9f8b]1905    }
1906  else
1907    {
[0e8a5a]1908      if (size(#)==0)
[7fe9f8b]1909        {
[0e8a5a]1910          matrix MC=C[5];
1911          def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v);
1912          setring HomWeyl;
1913          matrix J_C=fetch(B,MC);
1914          J_C=nHomogenize(J_C);
1915          /*computation of V_d-strict resolution of C[5]->needed for proc
1916            VdstrictDoubleComplexes*/
1917          def ringofSyz=Sres(transpose(J_C),lengthofres);
1918          setring ringofSyz;
1919          matrix J_C=transpose(matrix(RES[2]));
1920          J_C=subst(J_C,h,1);
1921          logens=ncols(J_C)+1;
1922          matrix zerom;
1923          list rofC;//will contain resolution of C
1924          for (i=3; i<=n+si+1; i++)
[7fe9f8b]1925            {
[0e8a5a]1926              if (size(RES)>=i)
1927                {
1928                  zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
1929                  if (RES[i]!=zerom)
1930                    {
1931                      rofC[i-2]=(matrix(RES[i]));
[7fe9f8b]1932
[0e8a5a]1933                      if (i==3)
1934                        {
1935                          if (nrows(rofC[i-2])-logens+1!=nrows(J_C))
1936                            {
1937                              //build the resolution
1938                              nr=nrows(J_C)+logens-1;
1939                              nc=ncols(rofC[i-2]);
1940                              rofC[i-2]=matrix(rofC[i-2],nr,nc);
1941                            }
1942
1943                        }
1944                      if (i!=3)
1945                        {
1946                          if (nrows(rofC[i-2])-logens+1!=nrows(rofC[i-3]))
1947                            {
1948                              nr=nrows(rofC[i-3])+logens-1;
1949                              nc=ncols(rofC[i-2]);
1950                              rofC[i-2]=matrix(rofC[i-2],nr,nc);
1951                            }
1952                        }
1953                      i1=intvec(logens..nrows(rofC[i-2]));
1954                      i2=intvec(1..ncols(rofC[i-2]));
1955                      rofC[i-2]=transpose(submat(rofC[i-2],i1,i2));
1956                      logens=logens+ncols(rofC[i-2]);
1957                      rofC[i-2]=subst(rofC[i-2],h,1);
1958                    }
1959                  else
1960                    {
1961                      rofC[i-2]=list();
1962                    }
1963                }
1964              else
1965                {
1966                  rofC[i-2]=list();
1967                }
[7fe9f8b]1968            }
[0e8a5a]1969          if(size(rofC[1])==0)
[7fe9f8b]1970            {
[0e8a5a]1971              omitemptylist=1;
1972            }
1973          setring B;
1974          matrix  J_C=fetch(ringofSyz,J_C);
1975          if (omitemptylist!=1)
1976            {
1977              list rofC=fetch(ringofSyz,rofC);
1978            }
1979          omitemptylist=0;
1980          kill HomWeyl;
1981          kill ringofSyz;
1982        }
1983      else
1984        {
1985          matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
1986        }
[7fe9f8b]1987    }
1988  /* we compute a V_d-strict Groebner basis for C[3]*/
1989  matrix J_A=C[1];
1990  matrix f_CB=C[4];
1991  matrix f_ACB=transpose(concat(transpose(C[2]),transpose(f_CB)));
1992  matrix J_AC=divdr(f_ACB,C[3]);
1993  matrix P=matrixLift(J_AC * prodr(ncols(C[1]),ncols(C[5])) ,J_C);
1994  list storePi;
1995  matrix Pi[1][ncols(J_AC)];
1996  for (i=1; i<=nrows(J_C); i++)
1997    {
[0e8a5a]1998      for (j=1; j<=nrows(J_AC);j++)
1999        {
2000          Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC)));
2001        }
2002      storePi[i]=Pi;
2003      Pi=0;
[7fe9f8b]2004    }
2005  /*we compute the shift vector for C[1]*/
2006  intvec m_a;
2007  list findMin;
[0e8a5a]2008  int comMin;
[7fe9f8b]2009  for (i=1; i<=ncols(J_A); i++)
2010    {
[0e8a5a]2011      for (j=1; j<=size(storePi);j++)
2012        {
2013          if (storePi[j][1,i]!=0)
2014            {
2015              comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v);
2016              comMin=comMin-VdDeg(storePi[j][1,i],n,intvec(0));
2017              findMin[size(findMin)+1]=comMin;
2018            }
2019        }
2020      if (size(findMin)!=0)
2021        {
2022          m_a[i]=Min(findMin);
2023          findMin=list();
2024        }
2025      else
2026        {
2027          m_a[i]=0;
2028        }
[7fe9f8b]2029    }
2030  matrix zero[ncols(J_A)][ncols(J_C)];
2031  matrix g_AB=concat(unitmat(ncols(J_A)),zero);
2032  matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C)))));
2033  intvec m_b=m_a,v;
[0e8a5a]2034  /* computation of a V_d-strict Groebner basis of C[1] (and resolution if
2035     Syzstring=='Vdres') */
[7fe9f8b]2036  if (Syzstring=="Vdres")
[0e8a5a]2037    {
2038      J_A=VdStrictGB(J_A,n,m_a);
2039    }
[7fe9f8b]2040  else
2041    {
[0e8a5a]2042      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a);
2043      setring HomWeyl;
2044      matrix J_A=fetch(B,J_A);
2045      J_A=nHomogenize(J_A);
2046      def ringofSyz=Sres(transpose(J_A),lengthofres);
2047      setring ringofSyz;
2048      matrix J_A=transpose(matrix(RES[2]));
2049      matrix zerom;
2050      J_A=subst(J_A,h,1);
2051      logens=ncols(J_A)+1;
2052      list rofA;
2053      for (i=3; i<=n+si+1; i++)
[7fe9f8b]2054        {
[0e8a5a]2055          if (size(RES)>=i)
[7fe9f8b]2056            {
[0e8a5a]2057              zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
2058              if (RES[i]!=zerom)
2059                {
2060                  rofA[i-2]=matrix(RES[i]);// resolution for C[1]
2061                  if (i==3)
2062                    {
2063                      if (nrows(rofA[i-2])-logens+1!=nrows(J_A))
2064                        {
2065                          nr=nrows(J_A)+logens-1;
2066                          nc=ncols(rofA[i-2]);
2067                          rofA[i-2]=matrix(rofA[i-2],nr,nc);
2068                        }
2069                    }
2070                  if (i!=3)
2071                    {
2072                      if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3]))
2073                        {
2074                          nr=nrows(rofA[i-3])+logens-1;
2075                          nc=ncols(rofA[i-2]);
2076                          rofA[i-2]=matrix(rofA[i-2],nr,nc);
2077                        }
2078                    }
2079                  i1=intvec(logens..nrows(rofA[i-2]));
2080                  i2=intvec(1..ncols(rofA[i-2]));
2081                  rofA[i-2]=transpose(submat(rofA[i-2],i1,i2));
2082                  logens=logens+ncols(rofA[i-2]);
2083                  rofA[i-2]=subst(rofA[i-2],h,1);
2084                }
2085              else
2086                {
2087                  rofA[i-2]=list();
2088                }
[7fe9f8b]2089            }
[0e8a5a]2090          else
[7fe9f8b]2091            {
[0e8a5a]2092              rofA[i-2]=list();
[7fe9f8b]2093            }
2094        }
[0e8a5a]2095      if(size(rofA[1])==0)
[7fe9f8b]2096        {
[0e8a5a]2097          omitemptylist=1;
[7fe9f8b]2098        }
[0e8a5a]2099      setring B;
2100      J_A=fetch(ringofSyz,J_A);
2101      if (omitemptylist!=1)
2102        {
2103          list rofA=fetch(ringofSyz,rofA);
2104        }
2105      omitemptylist=0;
2106      kill HomWeyl;
2107      kill ringofSyz;
[7fe9f8b]2108    }
2109  J_AC=transpose(storePi[1]);
2110  for (i=2; i<= size(storePi); i++)
[0e8a5a]2111    {
2112      J_AC=concat(J_AC, transpose(storePi[i]));
2113    }
[7fe9f8b]2114  J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC));
2115  list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a));
2116  Vdstrict[7]=list(m_b);
[0e8a5a]2117  Vdstrict[8]=list(v);
[7fe9f8b]2118  if(Syzstring=="Sres")
2119    {
[0e8a5a]2120      Vdstrict[9]=rofA;
2121      if(size(#)==0)
2122        {
2123          Vdstrict[10]=rofC;
2124        }
[7fe9f8b]2125    }
2126  return (Vdstrict);
2127}
2128
2129////////////////////////////////////////////////////////////////////////////////////
2130
[0e8a5a]2131static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL)
[7fe9f8b]2132{
[0e8a5a]2133  /* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact
2134     pieces.
2135     We asume that we are given two sequences of the form coker(L[i][1])->
2136     coker(L[i][3])->coker(L[i][5]) with differentials L[i][2] and L[i][4] such
2137     that L[1][3]=L[2][1].We are going to transform them to V_d-strict sequences
2138     J_D->J_A->J_F and J_A->J_B->J_C*/
[7fe9f8b]2139  int omitemptylist;
2140  int lengthofres=sizeL+d-1;
2141  int logens;
2142  def B=basering;
2143  matrix J_F=L[1][5];
2144  matrix J_D=L[1][1];
2145  matrix f_FA=L[1][4];
[0e8a5a]2146  /*We find new presentations coker(J_DF) and coker(J_DFC)  for L[1][4]=L[2][1]
2147     and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */
[7fe9f8b]2148  matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA)));
2149  matrix J_DF=divdr(f_DFA,L[1][3]);//coker(J_DF) is isomorphic to coker(L[2][1]);
2150  matrix J_C=L[2][5];
2151  matrix f_CB=L[2][4];
2152  matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB)));
2153  matrix J_DFC=divdr(f_DFCB,L[2][3]);//coker(J_DFC) are coker(L[2][3)]) isomorphic
[0e8a5a]2154  /* find a shift vector on the range of J_F such that the first sequence is
2155     exact*/
[7fe9f8b]2156  matrix P=matrixLift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C);
2157  list storePi;
2158  matrix Pi[1][ncols(J_DFC)];
2159  int i; int j;
2160  for (i=1; i<=nrows(J_C); i++)
2161    {
[0e8a5a]2162      for (j=1; j<=nrows(J_DFC);j++)
2163        {
2164          Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC)));
2165        }
2166      storePi[i]=Pi;
2167      Pi=0;
[7fe9f8b]2168    }
2169  intvec m_a;
2170  list findMin;
2171  list noMin;
2172  int comMin;
2173  int nr,nc;
2174  intvec i1,i2;
2175  for (i=1; i<=ncols(J_DF); i++)
2176    {
[0e8a5a]2177      for (j=1; j<=size(storePi);j++)
2178        {
2179          if (storePi[j][1,i]!=0)
2180            {
2181              comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v);
2182              comMin=comMin-VdDeg(storePi[j][1,i],d,intvec(0));
2183              findMin[size(findMin)+1]=comMin;
2184            }
2185        }
2186      if (size(findMin)!=0)
2187        {
2188          m_a[i]=Min(findMin);// shift vector for L[2][1]
2189          findMin=list();
2190          noMin[i]=0;
2191        }
2192      else
2193        {
2194          noMin[i]=1;
2195        }
[7fe9f8b]2196    }
[0e8a5a]2197  if (size(m_a) < ncols(J_DF))
[7fe9f8b]2198    {
[0e8a5a]2199      m_a[ncols(J_DF)]=0;
[7fe9f8b]2200    }
2201  intvec m_f=m_a[ncols(J_D)+1..size(m_a)];
[0e8a5a]2202  /* Computation of a V_d-strict Groebner basis of J_F=L[1][5]:
2203     if Syzstring=="Vdres" this is done using the method of weighted homogenization
2204     (Prop. 3.9 [OT])
2205     else we use the homogenized Weyl algerba for Groebner basis computations
2206     (Prop 9.9 [OT]), in this case we already compute resolutions
2207     (Thm. 9.10 in [OT]) to omit extra Groebner basis  computations later on*/
[7fe9f8b]2208  if (Syzstring=="Vdres")
[0e8a5a]2209    {
2210      J_F=VdStrictGB(J_F,d,m_f);
2211    }
[7fe9f8b]2212  else
2213    {
[0e8a5a]2214      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f);
2215      setring HomWeyl;
2216      matrix J_F=fetch(B,J_F);
2217      J_F=nHomogenize(J_F);
2218      def ringofSyz=Sres(transpose(J_F),lengthofres);
2219      setring ringofSyz;
2220      matrix J_F=transpose(matrix(RES[2]));
2221      J_F=subst(J_F,h,1);
2222      logens=ncols(J_F)+1;
2223      list rofF;
2224      for (i=3; i<=d+sizeL+1; i++)
[7fe9f8b]2225        {
[0e8a5a]2226          if (size(RES)>=i)
[7fe9f8b]2227            {
[0e8a5a]2228              if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
2229                {
2230                  rofF[i-2]=(matrix(RES[i]));// resolution for J_F
2231                  if (i==3)
2232                    {
2233                      if (nrows(rofF[i-2])-logens+1!=nrows(J_F))
2234                        {
2235                          nr=nrows(J_F)+logens-1;
2236                          nc=ncols(rofF[i-2]);
2237                          rofF[i-2]=matrix(rofF[i-2],nr,nc);
2238                        }
2239                    }
2240                  if (i!=3)
2241                    {
2242                      if (nrows(rofF[i-2])-logens+1!=nrows(rofF[i-3]))
2243                        {
2244                          nr=nrows(rofF[i-3])+logens-1;
2245                          rofF[i-2]=matrix(rofF[i-2],nr,ncols(rofF[i-2]));
2246                        }
2247                    }
2248                  i1=intvec(logens..nrows(rofF[i-2]));
2249                  i2=intvec(1..ncols(rofF[i-2]));
2250                  rofF[i-2]=transpose(submat(rofF[i-2],i1,i2));
2251                  logens=logens+ncols(rofF[i-2]);
2252                  rofF[i-2]=subst(rofF[i-2],h,1);
2253                }
2254              else
2255                {
2256                  rofF[i-2]=list();
2257                }
[7fe9f8b]2258            }
[0e8a5a]2259          else
[7fe9f8b]2260            {
[0e8a5a]2261              rofF[i-2]=list();
[7fe9f8b]2262            }
2263        }
[0e8a5a]2264      if(size(rofF[1])==0)
[7fe9f8b]2265        {
[0e8a5a]2266          omitemptylist=1;
[7fe9f8b]2267        }
[0e8a5a]2268      setring B;
2269      J_F=fetch(ringofSyz,J_F);
2270      if (omitemptylist!=1)
2271        {
2272          list rofF=fetch(ringofSyz,rofF);
2273        }
2274      omitemptylist=0;
2275      kill HomWeyl;
2276      kill ringofSyz;
[7fe9f8b]2277    }
2278  /*find shift vectors on the range of J_D*/
2279  P=matrixLift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F);
2280  list storePinew;
2281  matrix Pidf[1][ncols(J_DF)];
2282  for (i=1; i<=nrows(J_F); i++)
2283    {
[0e8a5a]2284      for (j=1; j<=nrows(J_DF);j++)
2285        {
2286          Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF)));
2287        }
2288      storePinew[i]=Pidf;
2289      Pidf=0;
[7fe9f8b]2290    }
2291  intvec m_d;
2292  for (i=1; i<=ncols(J_D); i++)
[0e8a5a]2293    {
2294      for (j=1; j<=size(storePinew);j++)
2295        {
2296          if (storePinew[j][1,i]!=0)
2297            {
2298              comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f);
2299              comMin=comMin-VdDeg(storePinew[j][1,i],d,intvec(0));
2300              findMin[size(findMin)+1]=comMin;
2301            }
2302        }
2303      if (size(findMin)!=0)
2304        {
2305          if (noMin[i]==0)
2306            {
2307              m_d[i]=Min(insert(findMin,m_a[i]));
2308              m_a[i]=m_d[i];
2309            }
2310          else
2311            {
2312              m_d[i]=Min(findMin);
2313              m_a[i]=m_d[i];
2314            }
2315        }
[7fe9f8b]2316      else
[0e8a5a]2317        {
2318          m_d[i]=m_a[i];
2319        }
2320      findMin=list();
[7fe9f8b]2321    }
[0e8a5a]2322  /* compute a V_d-strict Groebner basis (and resolution of J_D if
2323     Syzstring!='Vdres') for J_D*/
2324  if (Syzstring=="Vdres")
[7fe9f8b]2325    {
[0e8a5a]2326      J_D=VdStrictGB(J_D,d,m_d);
[7fe9f8b]2327    }
2328  else
2329    {
[0e8a5a]2330      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d);
2331      setring HomWeyl;
2332      matrix J_D=fetch(B,J_D);
2333      J_D=nHomogenize(J_D);
2334      def ringofSyz=Sres(transpose(J_D),lengthofres);
2335      setring ringofSyz;
2336      matrix J_D=transpose(matrix(RES[2]));
2337      J_D=subst(J_D,h,1);
2338      logens=ncols(J_D)+1;
2339      list rofD;
2340      for (i=3; i<=d+sizeL+1; i++)
[7fe9f8b]2341        {
[0e8a5a]2342          if (size(RES)>=i)
[7fe9f8b]2343            {
[0e8a5a]2344              if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
2345                {
2346                  rofD[i-2]=(matrix(RES[i]));// resolution for J_D
2347                  if (i==3)
2348                    {
2349                      if (nrows(rofD[i-2])-logens+1!=nrows(J_D))
2350                        {
2351                          nr=nrows(J_D)+logens-1;
2352                          rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
2353                        }
2354                    }
2355                  if (i!=3)
2356                    {
2357                      if (nrows(rofD[i-2])-logens+1!=nrows(rofD[i-3]))
2358                        {
2359                          nr=nrows(rofD[i-3])+logens-1;
2360                          rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
2361                        }
2362                    }
2363                  i1=intvec(logens..nrows(rofD[i-2]));
2364                  i2=intvec(1..ncols(rofD[i-2]));
2365                  rofD[i-2]=transpose(submat(rofD[i-2],i1,i2));
2366                  logens=logens+ncols(rofD[i-2]);
2367                  rofD[i-2]=subst(rofD[i-2],h,1);
2368                }
2369              else
2370                {
2371                  rofD[i-2]=list();
2372                }
[7fe9f8b]2373            }
[0e8a5a]2374          else
[7fe9f8b]2375            {
[0e8a5a]2376              rofD[i-2]=list();
[7fe9f8b]2377            }
2378        }
[0e8a5a]2379      if(size(rofD[1])==0)
[7fe9f8b]2380        {
[0e8a5a]2381          omitemptylist=1;
[7fe9f8b]2382        }
[0e8a5a]2383      setring B;
2384      J_D=fetch(ringofSyz,J_D);
2385      if (omitemptylist!=1)
2386        {
2387          list rofD=fetch(ringofSyz,rofD);
2388        }
2389      omitemptylist=0;
2390      kill HomWeyl;
2391      kill ringofSyz;
[7fe9f8b]2392    }
[0e8a5a]2393  /* compute new matrices for J_A and J_B  such that their rows form a V_d-strict
2394     Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and
2395     nrows(J_B)=nrows(J_A)+nrows(J_C)*/
[7fe9f8b]2396  J_DF=transpose(storePinew[1]);
2397  for (i=2; i<=nrows(J_F); i++)
[0e8a5a]2398    {
2399      J_DF=concat(J_DF,transpose(storePinew[i]));
2400    }
[7fe9f8b]2401  J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF));
2402  J_DFC=transpose(storePi[1]);
2403  for (i=2; i<=nrows(J_C); i++)
[0e8a5a]2404    {
2405      J_DFC=concat(J_DFC,transpose(storePi[i]));
2406    }
[7fe9f8b]2407  J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC));
2408  intvec m_b=m_a,v;
2409  matrix zero[ncols(J_D)][ncols(J_F)];
2410  matrix g_DA=concat(unitmat(ncols(J_D)),zero);
2411  matrix g_AF=transpose(concat(transpose(zero),unitmat(ncols(J_F))));
2412  matrix zero1[ncols(J_DF)][ncols(J_C)];
2413  matrix g_AB=concat(unitmat(ncols(J_DF)),zero1);
2414  matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C))));
[0e8a5a]2415  list out;
[7fe9f8b]2416  out[1]=list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F));
2417  out[1]=out[1]+list(list(m_d),list(m_a),list(m_f));
2418  out[2]=list(list(J_DF),list(g_AB),list(J_DFC),list(g_BC),list(J_C));
2419  out[2]=out[2]+list(list(m_a),list(m_b),list(v));
2420  if (Syzstring=="Sres")
[0e8a5a]2421    {
2422      out[3]=rofD;
2423      out[4]=rofF;
2424    }
[7fe9f8b]2425  return(out);
2426}
2427
2428////////////////////////////////////////////////////////////////////////////////////
2429
2430static proc VdStrictDoubleComplexes(list L,int d,string Syzstring)
2431{
[0e8a5a]2432  /* We compute  V_d-strict resolutions over the V_d-strict short exact pieces from
2433     the procedure shortExactPiecesToVdStrict.
2434     We use Algorithms 3.14 and 3.15 in [W2]*/
[7fe9f8b]2435  int i,k,c,j,l,totaldeg,comparedegs,SBcom,verk;
[0e8a5a]2436  intvec fordegs;
[7fe9f8b]2437  intvec n_b,i1,i2;
2438  matrix rem,forML,subm,zerom,unitm,subm2;
2439  matrix J_B;
2440  list store;
2441  int t=size(L)+d;
2442  int vd1,vd2,nr,nc;
2443  def B=basering;
2444  int n=nvars(B) div 2;
2445  intvec v;
2446  list forhW;
2447  if (Syzstring=="Sres")
2448    {
[0e8a5a]2449    /*we already computed some of the resolutions in the procedure
2450      shortExactPiecesToVdStrict*/
2451      matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew;
2452      for (k=1; k<=(size(L)+d-1); k++)
2453        {
2454          L[1][1][1][k+1]=list();
2455          L[1][1][2][k+1]=list();
2456          L[1][1][6][k+1]=list();
2457        }
2458      L[1][1][6][size(L)+d+1]=list();
2459      matrix mem;
2460      for (i=2; i<=d+size(L)+1; i++)
2461        {;
2462          v=0;
2463          if(size(L[1][1][3][i-1])!=0)
2464            {
2465              if(i!=d+size(L)+1)
2466                {
2467                  /*horizontal differential*/
2468                  L[1][1][4][i-1]=unitmat(nrows(L[1][1][3][i-1]));
2469                }
2470              for (j=1; j<=nrows(L[1][1][3][i-1]); j++)
2471                {
2472                  mem=submat(L[1][1][3][i-1],j,intvec(1..ncols(L[1][1][3][i-1])));
2473                  v[j]=VdDeg(mem,d,L[1][1][7][i-1]);
2474                }
2475              L[1][1][7][i]=v;//new shift vector
2476              L[1][1][8][i]=v;
2477              L[1][2][6][i]=v;
2478            }
2479          else
2480            {
2481              if (i!=d+size(L)+1)
2482                {
2483                  L[1][1][4][i-1]=list();
2484                }
2485              L[1][1][7][i]=list();
2486              L[1][1][8][i]=list();
2487              L[1][2][6][i]=list();
2488            }
2489        }
2490      if (size(L[1][1][3][d+size(L)])!=0)
2491        {
2492          /*horizontal differential*/
2493          L[1][1][4][d+size(L)]=unitmat(nrows(L[1][1][3][d+size(L)]));
2494        }
2495      else
2496        {
2497          L[1][1][4][d+size(L)]=list();
2498        }
2499      for (k=1; k<size(L); k++)
2500        {
2501          /* We build a V_d-strict resolution for coker(L[k][2][1][1])->
2502             coker(L[k][2][3][1])->coker(L[k][2][5][1]) using the resolution
2503             obtained for coker(L[k][1][3][1]).
2504             L[k][2][i][j] will be the jth module in the resolution of L[k][2][i][1]
2505             for i=1,3,5.
2506             L[k][2][i+5][j] will be the jth  shift vector in the resolution of
2507             L[k][2][i][1](this holds also for the case Syzstring=="Vdres")*/
2508          for (i=2; i<=d+size(L); i++)
2509            {
2510              v=0;
2511              if (size(L[k][2][5][i-1])!=0)
2512                {
2513                  for (j=1; j<=nrows(L[k][2][5][i-1]); j++)
2514                    {
2515                      i1=intvec(1..ncols(L[k][2][5][i-1]));
2516                      mem=submat(L[k][2][5][i-1],j,i1);
2517                      v[j]=VdDeg(mem,d,L[k][2][8][i-1]);
2518                    }
2519                  /*next shift vector in th resolution of coker(L[k][2][5][1])*/
2520                  L[k][2][8][i]=v;
2521                }
[7fe9f8b]2522              else
2523                {
[0e8a5a]2524                  L[k][2][8][i]=list();
[7fe9f8b]2525                }
[0e8a5a]2526              /* we build step by step a resolution for coker(L[k][2][5][1]) using
2527                 the resolutions of coker(L[k][2][1][1]) and coker(L[k][2][5][1])*/
2528              if (size(L[k][2][5][i])!=0)
[7fe9f8b]2529                {
[0e8a5a]2530                  if (size(L[k][2][1][i])!=0 or size(L[k][2][1][i-1])!=0)
2531                    {
2532                      L[k][2][3][i]=transpose(syz(transpose(L[k][2][3][i-1])));
2533                      nr= nrows(L[k][2][1][i-1]);
2534                      nc=ncols(L[k][2][5][i]);
2535                      Pold=matrixLift(L[k][2][3][i]*prodr(nr,nc), L[k][2][5][i]);
2536                      matrix Pi[1][ncols(L[k][2][3][i])];
2537                       for (l=1; l<=nrows(L[k][2][5][i]); l++)
2538                        {
2539                          for (j=1; j<=nrows(L[k][2][3][i]); j++)
2540                            {
2541                              i2=intvec(1..ncols(L[k][2][3][i]));
2542                              Pi=Pi+Pold[l,j]*submat(L[k][2][3][i],j,i2);
2543                            }
2544                          if (l==1)
2545                            {
2546                              Picombined=transpose(Pi);
2547                            }
2548                          else
2549                            {
2550                              Picombined=concat(Picombined,transpose(Pi));
2551                            }
2552                          Pi=0;
2553                        }
2554                       kill Pi;
2555                       Picombined=transpose(Picombined);
2556                       if (size(L[k][2][1][i])!=0)
2557                        {
2558                          if (i==2)
2559                            {
2560                              containsndeg=(0:ncols(L[k][2][1][1]));
2561                            }
2562                          containsndeg=nDeg(L[k][2][1][i-1],containsndeg);
2563                          forhW=list(L[k][2][6][i],containsndeg);
2564                          def HomWeyl=makeHomogenizedWeyl(n,forhW);
2565                          setring HomWeyl;
2566                          list L=fetch(B,L);
2567                          matrix M=L[k][2][1][i];
2568                          module Mmod;
2569                          list forM=nHomogenize(M,containsndeg,1);
2570                          M=forM[1];
2571                          totaldeg=forM[2];
2572                          kill forM;
2573                          matrix Maorig=fetch(B,Picombined);
2574                          matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M)));
2575                          matrix mem,subm,zerom;
2576                          matrix Pinew;
2577                          M=transpose(M);
2578                          SBcom=0;
2579                          for (l=1; l<=nrows(Ma); l++)
2580                            {
2581                              zerom=matrix(0,1,(ncols(Maorig)-ncols(Ma)));
2582                              i1=(ncols(Ma)+1..ncols(Maorig));
2583                              if (submat(Maorig,l,i1)==zerom)
2584                                {
2585                                  for (cc=1; cc<=ncols(Ma); cc++)
2586                                    {
2587                                      Maorig[l,cc]=0;
2588                                    }
2589                                }
2590                              i2=(ncols(Ma)+1..ncols(Maorig));
2591                              i1=(1..ncols(Ma));
2592                              if (VdDeg(submat(Maorig,l,i1),d,L[k][2][6][i])>
2593                                  VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]) and
2594                                  submat(Maorig,l,i1)!=matrix(0,1,ncols(Ma)))
2595                                {
2596                                  /*V_d-Grad is to big--> we make it smaller using
2597                                    Vdnormal form computations*/
2598                                  if (SBcom==0)
2599                                  {
2600                                    Mmod=slimgb(M);
2601                                    M=Mmod;
2602                                    SBcom=1;
2603                                  }
2604                                  //print("Reduzierung des V_d-Grades(Stelle1)");
2605                                  i2=(ncols(Ma)+1..ncols(Maorig));
2606                                  vd1=VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]);
2607                                  mem=submat(Ma,l,(1..ncols(Ma)));
2608                                  mem=nHomogenize(mem,containsndeg);
2609                                  mem=h^totaldeg*mem;
2610                                  mem=transpose(mem);
2611                                  mem=reduce(mem,Mod);//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2612                                  matrix jt=transpose(subst(mem,h,1));
2613                                  setring B;
2614                                  matrix jt=fetch(HomWeyl,jt);
2615                                  matrix need=fetch(HomWeyl,Maorig);
2616                                  need=submat(need,l,(1..ncols(need)));
2617                                  i1=L[k][2][6][i];
2618                                  i2=L[k][2][8][i];
2619                                  jt=VdNormalForm(need,L[k][2][1][i],d,i1,i2);
2620                                  setring HomWeyl;
2621                                  mem=fetch(B,jt);
2622                                  mem=transpose(mem);
2623                                  if (l==1)
2624                                    {
2625                                      Pinew=mem;
2626                                    }
2627                                  else
2628                                    {
2629                                      Pinew=concat(Pinew,mem);
2630                                    }
2631                                  vd2=VdDeg(transpose(mem),d,L[k][2][6][i]);
2632                                  if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem)))
2633                                    {//should not happen!!
2634                                      //print("Reduzierung fehlgeschlagen!!(Stelle1)");
2635                                    }
2636                                }
2637                              else
2638                                {
2639                                  if (l==1)
2640                                    {
2641                                      Pinew=transpose(submat(Ma,l,(1..ncols(Ma))));
2642                                    }
2643                                  else
2644                                    {
2645                                      subm=transpose(submat(Ma,l,(1..ncols(Ma))));
2646                                      Pinew=concat(Pinew,subm);
2647                                    }
2648                                }
2649                            }
2650                          Pinew=subst(Pinew,h,1);
2651                          Pinew=transpose(Pinew);
2652                          setring B;
2653                          Pinew=fetch(HomWeyl,Pinew);
2654                          kill HomWeyl;
2655                          L[k][2][3][i]=concat(Pinew,L[k][2][5][i]);
2656                          subm=transpose(L[k][2][3][i]);
2657                          subm=concat(transpose(L[k][2][1][i]),subm);
2658                          L[k][2][3][i]=transpose(subm);
2659                        }
2660                      else
2661                        {
2662                          L[k][2][3][i]=Picombined;
2663                        }
2664                      L[k+1][1][1][i]=L[k][2][5][i];
2665                      nr=nrows(L[k][2][1][i-1]);
2666                      nc=ncols(L[k][2][5][i]);
2667                      L[k][2][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
2668                      L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nc);
2669                      v=L[k][2][6][i],L[k][2][8][i];
2670                      L[k][2][7][i]=v;
2671                      L[k+1][1][6][i]=L[k][2][8][i];
2672                    }
[7fe9f8b]2673                  else
[0e8a5a]2674                    {
2675                      L[k][2][3][i]=L[k][2][5][i];
2676                      L[k][2][2][i]=list();
2677                      L[k][2][7][i]=L[k][2][8][i];
2678                      L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1]));
2679                      L[k+1][1][6][i]=L[k][2][8][i];
2680                      L[k+1][1][1][i]=L[k][2][5][i];
2681                    }
[7fe9f8b]2682                }
[0e8a5a]2683              else
[7fe9f8b]2684                {
[0e8a5a]2685                  if (size(L[k][2][1][i])!=0)
2686                    {
2687                      if (size(L[k][2][5][i-1])!=0)
2688                        {
2689                          nr=nrows(L[k][2][5][i-1]);
2690                          L[k][2][3][i]=concat(L[k][2][1][i],matrix(0,1,nr));
2691                          v=L[k][2][6][i],L[k][2][8][i];
2692                          L[k][2][7][i]=v;
2693                          nc=nrows(L[k][2][1][i-1]);
2694                          L[k][2][2][i]=concat(unitmat(nc),matrix(0,nc,nr));
2695                          L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nr);
2696                        }
2697                      else
2698                        {
2699                          L[k][2][3][i]=L[k][2][1][i];
2700                          L[k][2][7][i]=L[k][2][6][i];
2701                          L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1]));
2702                          L[k][2][4][i]=list();
2703                        }
2704                      L[k+1][1][1][i]=L[k][2][5][i];
2705                      L[k+1][1][6][i]=L[k][2][8][i];
2706                    }
[7fe9f8b]2707                  else
[0e8a5a]2708                    {
2709                      L[k][2][3][i]=list();
2710                      if (size(L[k][2][6][i])!=0)
2711                        {
2712                          if (size(L[k][2][8][i])!=0)
2713                            {
2714                              v=L[k][2][6][i],L[k][2][8][i];
2715                              L[k][2][7][i]=v;
2716                              nr=nrows(L[k][2][1][i-1]);
2717                              nc=nrows(L[k][2][5][i-1]);
2718                              L[k][2][2][i]=concat(unitmat(nc),matrix(0,nr,nc));
2719                              L[k][2][4][i]=prodr(nr,nrows(L[k][2][5][i-1]));
2720                            }
2721                          else
2722                            {
2723                              L[k][2][7][i]=L[k][2][6][i];
2724                              L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1]));
2725                              L[k][2][4][i]=list();
2726                            }
2727                        }
2728                      else
2729                        {
2730                          if (size(L[k][2][8][i])!=0)
2731                            {
2732                              L[k][2][7][i]=L[k][2][8][i];
2733                              L[k][2][2][i]=list();
2734                              L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1]));
2735                            }
2736                          else
2737                            {
2738                              L[k][2][7][i]=list();
2739                              L[k][2][2][i]=list();
2740                              L[k][2][4][i]=list();
2741                            }
2742                        }
2743                      L[k+1][1][1][i]=L[k][2][5][i];
2744                      L[k+1][1][6][i]=L[k][2][8][i];
2745                    }
[7fe9f8b]2746                }
2747            }
[0e8a5a]2748          i=d+size(L)+1;
2749          v=0;
2750          if (size(L[k][2][5][i-1])!=0)
[7fe9f8b]2751            {
[0e8a5a]2752              for (j=1; j<=nrows(L[k][2][5][i-1]); j++)
2753                {
2754                  mem=submat(L[k][2][5][i-1],j,intvec(1..ncols(L[k][2][5][i-1])));
2755                  v[j]=VdDeg(mem,d,L[k][2][8][i-1]);
2756                }
2757              L[k][2][8][i]=v;
2758              if (size(L[k][2][6][i])!=0)
2759                {
2760                  v=L[k][2][6][i],L[k][2][8][i];
2761                  L[k][2][7][i]=v;
2762                }
[7fe9f8b]2763              else
[0e8a5a]2764                {
2765                  L[k][2][7][i]=L[k][2][8][i];
2766                }
[7fe9f8b]2767            }
[0e8a5a]2768          else
[7fe9f8b]2769            {
[0e8a5a]2770              L[k][2][8][i]=list();
2771              L[k][2][7][i]=L[k][2][6][i];
2772            }
2773          L[k+1][1][6][i]=L[k][2][8][i];
2774          /* now we build V_d-strict resolutions for the sequences
2775             coker(L[k+1][1][1][1])->coker(L[k+1][1][3][1])->coker(L[k+1][1][5][i])
2776             using the resolutions  for coker(L[k][2][5][1]) we just obtained
2777             (works exactly the same as above)*/
2778          for (i=2; i<=d+size(L); i++)
2779            {
2780              v=0;
2781              if (size(L[k+1][1][5][i-1])!=0)
2782                {
2783                  for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++)
2784                    {
2785                      i1=intvec(1..ncols(L[k+1][1][5][i-1]));
2786                      mem=submat(L[k+1][1][5][i-1],j,i1);
2787                      v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]);
2788                    }
2789                  L[k+1][1][8][i]=v;
2790                }
[7fe9f8b]2791              else
2792                {
[0e8a5a]2793                  L[k+1][1][8][i]=list();
[7fe9f8b]2794                }
[0e8a5a]2795              if (size(L[k+1][1][5][i])!=0)
[7fe9f8b]2796                {
[0e8a5a]2797                  if (size(L[k+1][1][1][i])!=0 or size(L[k+1][1][1][i-1])!=0)
2798                    {
2799                      L[k+1][1][3][i]=transpose(syz(transpose(L[k+1][1][3][i-1])));
2800                      nr=nrows(L[k+1][1][1][i-1]);
2801                      nc=ncols(L[k+1][1][5][i]);
2802                      Pold=matrixLift(L[k+1][1][3][i]*prodr(nr,nc),L[k+1][1][5][i]);
2803                      matrix Pi[1][ncols(L[k+1][1][3][i])];
2804                      for (l=1; l<=nrows(L[k+1][1][5][i]); l++)
2805                        {
2806                          for (j=1; j<=nrows(L[k+1][1][3][i]); j++)
2807                            {
2808                              i2=intvec(1..ncols(L[k+1][1][3][i]));
2809                              Pi=Pi+Pold[l,j]*submat(L[k+1][1][3][i],j,i2);
2810                            }
2811                          if (l==1)
2812                            {
2813                              Picombined=transpose(Pi);
2814                            }
2815                          else
2816                            {
2817                              Picombined=concat(Picombined,transpose(Pi));
2818                            }
2819                          Pi=0;
2820                        }
2821                      kill Pi;
2822                      Picombined=transpose(Picombined);
2823                      if(size(L[k+1][1][1][i])!=0)
2824                        {
2825                          if (i==2)
2826                            {
2827                              containsndeg=(0:ncols(L[k+1][1][1][i-1]));
2828                            }
2829                          containsndeg=nDeg(L[k+1][1][1][i-1],containsndeg);
2830                          forhW=list(L[k+1][1][6][i], containsndeg);
2831                          def HomWeyl=makeHomogenizedWeyl(n,forhW);
2832                          setring HomWeyl;
2833                          list L=fetch(B,L);
2834                          matrix M=L[k+1][1][1][i];
2835                          module Mmod;
2836                          list forM=nHomogenize(M,containsndeg,1);
2837                          M=forM[1];
2838                          totaldeg=forM[2];
2839                          kill forM;
2840                          matrix Maorig=fetch(B,Picombined);
2841                          matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M)));
2842                          Ma=nHomogenize(Ma,containsndeg);
2843                          matrix mem,subm,zerom,subm2;
2844                          matrix Pinew;
2845                          M=transpose(M);
2846                          SBcom=0;
2847                          for (l=1; l<=nrows(Ma); l++)
2848                            {
2849                              i2=(ncols(Ma)+1..ncols(Maorig));
2850                              nc=ncols(Maorig)-ncols(Ma);
2851                              if (submat(Maorig,l,i2)==matrix(0,1,nc))
2852                                {
2853                                  for (cc=1; cc<=ncols(Ma); cc++)
2854                                    {
2855                                      Maorig[l,cc]=0;
2856                                    }
2857                                }
2858                              i1=(1..ncols(Ma));
2859                              i2=L[k+1][1][8][i];
2860                              subm=submat(Maorig,l,i1);
2861                              subm2=submat(Maorig,l,(ncols(Ma)+1..ncols(Maorig)));
2862                              if (VdDeg(subm,d,L[k+1][1][6][i])>VdDeg(subm2,d,i2)
2863                                  and subm!=matrix(0,1,ncols(Ma)))
2864                                {
2865                                  //print("Reduzierung des Vd-Grades (Stelle2)");
2866                                  if (SBcom==0)
2867                                    {
2868                                      Mmod=slimgb(M);
2869                                      M=Mmod;
2870                                      SBcom=1;
2871                                    }
2872                                  vd1=VdDeg(subm2,d,L[k+1][1][8][i]);
2873                                  mem=submat(Ma,l,(1..ncols(Ma)));
2874                                  mem=nHomogenize(mem,containsndeg);
2875                                  mem=h^totaldeg*mem;
2876                                  mem=transpose(mem);
2877                                  mem=reduce(mem,Mmod);
2878                                  if (l==1)
2879                                    {
2880                                      Pinew=mem;
2881                                    }
2882                                  else
2883                                    {
2884                                      Pinew=concat(Pinew,mem);
2885                                    }
2886                                  vd2=VdDeg(transpose(mem),d,L[k+1][1][6][i]);
2887                                  if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem)))
2888                                    {//should not happen
2889                                      //print("Reduzierung fehlgeschlagen!!!!(Stelle2)");
2890                                    }
2891                                }
2892                              else
2893                                {
2894                                  if (l==1)
2895                                    {
2896                                      Pinew=transpose(submat(Ma,l,(1..ncols(Ma))));
2897                                    }
2898                                  else
2899                                    {
2900                                      subm=transpose(submat(Ma,l,(1..ncols(Ma))));
2901                                      Pinew=concat(Pinew,subm);
2902                                    }
2903                                }
2904                            }
2905                          Pinew=subst(Pinew,h,1);
2906                          Pinew=transpose(Pinew);
2907                          setring B;
2908                          Pinew=fetch(HomWeyl,Pinew);
2909                          kill HomWeyl;
2910                          L[k+1][1][3][i]=concat(Pinew,L[k+1][1][5][i]);
2911                          subm=transpose(L[k+1][1][1][i]);
2912                          subm2=transpose(L[k+1][1][3][i]);
2913                          L[k+1][1][3][i]=transpose(concat(subm,subm2));
2914                        }
2915                      else
2916                        {
2917                          L[k+1][1][3][i]=Picombined;
2918                        }
2919                      L[k+1][2][1][i]=L[k+1][1][3][i];
2920                      nr=nrows(L[k+1][1][1][i-1]);
2921                      nc=ncols(L[k+1][1][5][i]);
2922                      L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
2923                      L[k+1][1][4][i]=prodr(nr,nc);
2924                      v=L[k+1][1][6][i],L[k+1][1][8][i];
2925                      L[k+1][1][7][i]=v;
2926                      L[k+1][2][6][i]=L[k+1][1][7][i];
2927                    }
[7fe9f8b]2928                  else
[0e8a5a]2929                    {
2930                      L[k+1][1][3][i]=L[k+1][1][5][i];
2931                      L[k+1][1][2][i]=list();
2932                      L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1]));
2933                      L[k+1][1][7][i]=L[k+1][1][8][i];
2934                      L[k+1][2][6][i]=L[k+1][1][7][i];
2935                      L[k+1][2][1][i]=L[k+1][1][3][i];
2936                    }
[7fe9f8b]2937                }
[0e8a5a]2938              else
[7fe9f8b]2939                {
[0e8a5a]2940                  if (size(L[k+1][1][1][i])!=0)
2941                    {
2942                      if (size(L[k+1][1][5][i-1])!=0)
2943                        {
2944                          zerom=matrix(0,1,nrows(L[k+1][1][5][i-1]));
2945                          L[k+1][1][3][i]=concat(L[k+1][1][1][i],zerom);
2946                          v=L[k+1][1][6][i],L[k+1][1][8][i];
2947                          L[k+1][1][7][i]=v;
2948                          nr=nrows(L[k+1][1][1][i-1]);
2949                          nc=nrows(L[k+1][1][5][i-1]);
2950                          L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
2951                          L[k+1][1][4][i]=prodr(nr,nc);
2952                        }
2953                      else
2954                        {
2955                          L[k+1][1][3][i]=L[k+1][1][1][i];
2956                          L[k+1][1][7][i]=L[k+1][1][6][i];
2957                          L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1]));
2958                          L[k+1][1][4][i]=list();
2959                        }
2960                      L[k+1][2][1][i]=L[k+1][1][3][i];
2961                      L[k+1][2][6][i]=L[k+1][1][7][i];
2962                    }
[7fe9f8b]2963                  else
[0e8a5a]2964                    {
2965                      L[k+1][1][3][i]=list();
2966                      if (size(L[k+1][1][6][i])!=0)
2967                        {
2968                          if (size(L[k+1][1][8][i])!=0)
2969                            {
2970                              v=L[k+1][1][6][i],L[k+1][1][8][i];
2971                              L[k+1][1][7][i]=v;
2972                              nr=nrows(L[k+1][1][1][i-1]);
2973                              nc=nrows(L[k+1][1][5][i-1]);
2974                              L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
2975                              L[k+1][1][4][i]=prodr(nr,nrows(L[k+1][1][5][i-1]));
2976                            }
2977                          else
2978                            {
2979                              L[k+1][1][7][i]=L[k+1][1][6][i];
2980                              L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1]));
2981                              L[k+1][1][4][i]=list();
2982                            }
2983                        }
2984                      else
2985                        {
2986                          if (size(L[k+1][1][8][i])!=0)
2987                            {
2988                              L[k+1][1][7][i]=L[k+1][1][8][i];
2989                              L[k+1][1][2][i]=list();
2990                              L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1]));
2991                            }
2992                          else
2993                            {
2994                              L[k+1][1][7][i]=list();
2995                              L[k+1][1][2][i]=list();
2996                              L[k+1][1][4][i]=list();
2997                            }
2998                        }
2999
3000                      L[k+1][2][1][i]=L[k+1][1][3][i];
3001                      L[k+1][2][6][i]=L[k+1][1][7][i];
3002                    }
[7fe9f8b]3003                }
[0e8a5a]3004            }
3005          i=size(L)+d+1;
3006          v=0;
3007          if (size(L[k+1][1][5][i-1])!=0)
[7fe9f8b]3008            {
[0e8a5a]3009              for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++)
3010                {
3011                  i1=intvec(1..ncols(L[k+1][1][5][i-1]));
3012                  mem=submat(L[k+1][1][5][i-1],j,i1);
3013                  v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]);
3014                }
3015              L[k+1][1][8][i]=v;
3016              if (size(L[k+1][1][6][i])!=0)
3017                {
3018                  v=L[k+1][1][6][i],L[k+1][1][8][i];
3019                  L[k+1][1][7][i]=v;
3020                }
[7fe9f8b]3021              else
[0e8a5a]3022                {
3023                  L[k+1][1][7][i]=L[k+1][1][8][i];
3024                }
[7fe9f8b]3025            }
[0e8a5a]3026          else
[7fe9f8b]3027            {
[0e8a5a]3028              L[k+1][1][8][i]=list();
3029              L[k+1][1][7][i]=L[k+1][1][8][i];
[7fe9f8b]3030            }
[0e8a5a]3031          L[k+1][2][6][i]=L[k+1][1][7][i];
[7fe9f8b]3032        }
[0e8a5a]3033      for (k=1; k<=(size(L)+d); k++)
[7fe9f8b]3034        {
[0e8a5a]3035          L[size(L)][2][5][k]=list();
3036          L[size(L)][2][4][k]=list();
3037          L[size(L)][2][8][k]=list();
3038          L[size(L)][2][3][k]=L[size(L)][2][1][k];
3039          L[size(L)][2][7][k]=L[size(L)][2][6][k];
[7fe9f8b]3040        }
[0e8a5a]3041      L[size(L)][2][7][size(L)+d+1]=L[size(L)][2][6][size(L)+d+1];
3042      L[size(L)][2][8][size(L)+d+1]=list();
3043      /* building the resolution of the last short exact piece*/
3044      for (i=2; i<=d+size(L); i++)
[7fe9f8b]3045        {
[0e8a5a]3046          v=0;
3047          if(size(L[size(L)][2][1][i-1])!=0)
3048            {
3049              L[size(L)][2][2][i]=unitmat(nrows(L[size(L)][2][1][i-1]));
3050            }
3051          else
3052            {
3053              L[size(L)][2][2][i-1]=list();
3054            }
[7fe9f8b]3055        }
[0e8a5a]3056      return(L);
[7fe9f8b]3057    }
3058  /*case Syzstring=="Vdres"*/
3059  list forVd;
3060  for (k=1; k<=(size(L)+d); k++)//?????
3061    {
[0e8a5a]3062      /* we compute a V_d-strict resolution for the first short exact piece*/
3063      L[1][1][1][k+1]=list();
3064      L[1][1][2][k+1]=list();
3065      L[1][1][6][k+1]=list();
3066      if (size(L[1][1][3][k])!=0)
3067        {
3068          for (i=1; i<=nrows(L[1][1][3][k]); i++)
3069            {
3070              rem=submat(L[1][1][3][k],i,(1..ncols(L[1][1][3][k])));
3071              n_b[i]=VdDeg(rem,d,L[1][1][7][k]);
3072            }
3073          J_B=transpose(syz(transpose(L[1][1][3][k])));
3074          L[1][1][7][k+1]=n_b;
3075          L[1][1][8][k+1]=n_b;
3076          L[1][1][4][k+1]=unitmat(nrows(L[1][1][3][k]));
3077          if (J_B!=matrix(0,nrows(J_B),ncols(J_B)))
3078            {
3079              J_B=VdStrictGB(J_B,d,n_b);
3080              L[1][1][3][k+1]=J_B;
3081              L[1][1][5][k+1]=J_B;
3082            }
3083          else
3084            {
3085              L[1][1][3][k+1]=list();
3086              L[1][1][5][k+1]=list();
3087            }
3088          n_b=0;
3089        }
[7fe9f8b]3090      else
3091        {
[0e8a5a]3092          L[1][1][3][k+1]=list();
3093          L[1][1][5][k+1]=list();
3094          L[1][1][7][k+1]=list();
3095          L[1][1][8][k+1]=list();
3096          L[1][1][4][k+1]=list();
[7fe9f8b]3097        }
[0e8a5a]3098      /* we compute step by step V_d-strict resolutions over
3099         coker(L[i][2][1][1])->coker(L[i][2][3][1])->coker(L[i][2][1][5])
3100         and coker(L[i+1][1][1][1])->coker(L[i+1][1][3][1])->coker(L[i+1][1][1][5])
3101         using the already computed resolutions for coker(L[i][2][1][1])=
3102         coker(L[i][1][3][1]) and coker(L[i+1][1][1][1])=coker(L[i][2][5][1])*/
3103      for (i=1; i<size(L); i++)
[7fe9f8b]3104        {
[0e8a5a]3105          forVd[1]=L[i][2][1][k];
3106          forVd[2]=L[i][2][2][k];
3107          forVd[3]=L[i][2][3][k];
3108          forVd[4]=L[i][2][4][k];
3109          forVd[5]=L[i][2][5][k];
3110          forVd[6]=L[i][2][6][k];
3111          forVd[7]=L[i][2][7][k];
3112          forVd[8]=L[i][2][8][k];
3113          store=toVdStrict2x3Complex(forVd,d,L[i][1][3][k+1],L[i][1][7][k+1]);
3114          for (j=1; j<=8; j++)
3115            {
3116              L[i][2][j][k+1]=store[j];
3117            }
3118          forVd[1]=L[i+1][1][1][k];
3119          forVd[2]=L[i+1][1][2][k];
3120          forVd[3]=L[i+1][1][3][k];
3121          forVd[4]=L[i+1][1][4][k];
3122          forVd[5]=L[i+1][1][5][k];
3123          forVd[6]=L[i+1][1][6][k];
3124          forVd[7]=L[i+1][1][7][k];
3125          forVd[8]=L[i+1][1][8][k];
3126          store=toVdStrict2x3Complex(forVd,d,L[i][2][5][k+1],L[i][2][8][k+1]);
3127          for (j=1; j<=8; j++)
3128            {
3129              L[i+1][1][j][k+1]=store[j];
3130            }
[7fe9f8b]3131        }
[0e8a5a]3132      if (size(L[size(L)][1][7][k+1])!=0)
[7fe9f8b]3133        {
[0e8a5a]3134          L[size(L)][2][4][k+1]=list();
3135          L[size(L)][2][5][k+1]=list();
3136          L[size(L)][2][6][k+1]=L[size(L)][1][7][k+1];
3137          L[size(L)][2][7][k+1]=L[size(L)][1][7][k+1];
3138          L[size(L)][2][8][k+1]=list();
3139          L[size(L)][2][2][k+1]=unitmat(size(L[size(L)][1][7][k+1]));
3140          if (size(L[size(L)][1][3][k+1])!=0)
3141            {
3142              L[size(L)][2][1][k+1]=L[size(L)][1][3][k+1];
3143              L[size(L)][2][3][k+1]=L[size(L)][1][3][k+1];
3144            }
3145          else
3146            {
3147              L[size(L)][2][1][k+1]=list();
3148              L[size(L)][2][3][k+1]=list();
3149            }
[7fe9f8b]3150        }
3151      else
3152        {
[0e8a5a]3153          for (j=1; j<=8; j++)
3154            {
3155              L[size(L)][2][j][k+1]=list();
3156            }
[7fe9f8b]3157        }
[0e8a5a]3158    }
3159  k=t;
3160  intvec n_c;
3161  intvec vn_b;
3162  list N_b;
3163  int n;
3164  /*computation of the shift vectors*/
3165  for (i=1; i<=size(L); i++)
3166    {
3167      for (n=1; n<=2; n++)
[7fe9f8b]3168        {
[0e8a5a]3169          if (i==1 and n==1)
3170            {
3171              L[i][n][6][k+1]=list();
3172            }
3173          else
3174            {
3175              if (n==1)
3176                {
3177                  L[i][1][6][k+1]=L[i-1][2][8][k+1];
3178                }
3179              else
3180                {
3181                  L[i][2][6][k+1]=L[i][1][7][k+1];
3182                }
3183            }
3184          N_b[1]=L[i][n][6][k+1];
3185          if (size(L[i][n][5][k])!=0)
3186            {
3187              for (j=1; j<=nrows(L[i][n][5][k]); j++)
3188                {
3189                  rem=submat(L[i][n][5][k],j,(1..ncols(L[i][n][5][k])));
3190                  n_c[j]=VdDeg(rem,d,L[i][n][8][k]);
3191                }
3192              L[i][n][8][k+1]=n_c;
3193            }
3194          else
3195            {
3196              L[i][n][8][k+1]=list();
3197            }
3198          N_b[2]=L[i][n][8][k+1];
3199          n_c=0;
3200          if (size(N_b[1])!=0)
3201            {
3202              vn_b=N_b[1];
3203              if (size(N_b[2])!=0)
3204                {
3205                  vn_b=vn_b,N_b[2];
3206                }
3207              L[i][n][7][k+1]=vn_b;
3208            }
3209          else
3210            {
3211              if (size(N_b[2])!=0)
3212                {
3213                  L[i][n][7][k+1]=N_b[2];
3214                }
3215              else
3216                {
3217                  L[i][n][7][k+1]=list();
3218                }
3219            }
[7fe9f8b]3220        }
3221    }
3222  return(L);
3223}
3224
[0e8a5a]3225////////////////////////////////////////////////////////////////////////////////////
[7fe9f8b]3226
3227static proc toVdStrict2x3Complex(list L,int d,list #)
3228{
[0e8a5a]3229  /* We build a one-step free resolution over a V_d-strict short exact piece
3230     (Algorithm 3.14 in [W2]).
3231     This procedure is called from the procedure VdStrictDoubleComplexes
3232     if Syzstring=='Vdres'*/
3233  matrix rem;
[7fe9f8b]3234  int i,j,cc;
[0e8a5a]3235  int nr;
[7fe9f8b]3236  list J_A=list(list());
3237  list J_B=list(list());
3238  list J_C=list(list());
3239  list g_AB=list(list());
3240  list g_BC=list(list());
3241  list n_a=list(list());
3242  list n_b=list(list());
3243  list n_c=list(list());
3244  intvec n_b1;
3245  matrix fromnf;
3246  intvec i1,i2;
3247  /* compute a one step V_d-strict resolution for L[5]*/
3248  if (size(L[5])!=0)
3249    {
[0e8a5a]3250      intvec n_c1;
3251      for (i=1; i<=nrows(L[5]); i++)
3252        {
3253          rem=submat(L[5],i,intvec(1..ncols(L[5])));
3254          n_c1[i]=VdDeg(rem,d, L[8]);//new shift vector
3255        }
3256      n_c[1]=n_c1;
3257      J_C[1]=transpose(syz(transpose(L[5])));
3258      if (J_C[1]!=matrix(0,nrows(J_C[1]),ncols(J_C[1])))
3259        {
3260          J_C[1]=VdStrictGB(J_C[1],d,n_c1);
3261          if (size(#[2])!=0)// new shift vector for the resolution of L[1]
3262            {
3263              n_a[1]=#[2];
3264              n_b1=n_a[1],n_c[1];
3265              n_b[1]=n_b1;
3266              matrix zero[nrows(L[1])][nrows(L[5])];
3267              g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));
3268              if (size(#[1])!=0)
3269                {
3270                  J_A=#[1];// one step V_d-strict resolution for L[1]
3271                  /* use resolutions of L[1] and L[5] to build a resolution for
3272                     L[3]*/
3273                  J_B[1]=transpose(matrix(syz(transpose(L[3]))));
3274                  matrix P=matrixLift(J_B[1]*prodr(nrows(L[1]),nrows(L[5])),J_C[1]);
3275                  matrix Pi[1][ncols(J_B[1])];
3276                  matrix Picombined;
3277                  for (i=1; i<=nrows(J_C[1]); i++)
3278                    {
3279                      for (j=1; j<=nrows(J_B[1]);j++)
3280                        {
3281                          Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1])));
3282                        }
3283                      if (i==1)
3284                        {
3285                          Picombined=transpose(Pi);
3286                        }
3287                      else
3288                        {
3289                          Picombined=concat(Picombined,transpose(Pi));
3290                        }
3291                      Pi=0;
3292                    }
3293                  Picombined=transpose(Picombined);
3294                  fromnf=VdNormalForm(Picombined,J_A[1],d,n_a[1],n_c[1]);
3295                  i1=intvec(1..nrows(Picombined));
3296                  i2=intvec((ncols(J_A[1])+1)..ncols(Picombined));
3297                  Picombined=concat(fromnf,submat(Picombined,i1,i2));
3298                  J_B[1]=transpose(matrix(J_A[1],nrows(J_A[1]),ncols(J_B[1])));
3299                  J_B[1]=transpose(concat(J_B[1],transpose(Picombined)));
3300                  g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
3301                }
3302              else//L[1] is already a resolution
3303                {
3304                  //compute a resolution for L[3]
3305                  J_B=transpose(matrix(syz(transpose(L[3]))));
3306                  matrix P=matrixLift(J_B[1]*prodr(nrows(L[1]),nrows(L[5])),J_C[1]);
3307                  matrix Pi[1][ncols(J_B[1])];
3308                  matrix Picombined;
3309                  for (i=1; i<=nrows(J_C[1]); i++)
3310                    {
3311                      for (j=1; j<=nrows(J_B[1]);j++)
3312                        {
3313                          Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1])));
3314                        }
3315                      if (i==1)
3316                        {
3317                          Picombined=transpose(Pi);
3318                        }
3319                      else
3320                        {
3321                          Picombined=concat(Picombined,transpose(Pi));
3322                        }
3323                      Pi=0;
3324                    }
3325                  Picombined=transpose(Picombined);
3326                  J_B[1]=Picombined;
3327                  g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
3328                }
3329            }
3330          else
3331            {
3332              n_b=n_c[1];
3333              J_B[1]=J_C[1];
3334              g_BC=unitmat(ncols(J_C[1]));
3335            }
3336        }
[7fe9f8b]3337      else
[0e8a5a]3338        {
3339          J_C=list(list());// L[5] is already a resolution
3340          if (size(#[2])!=0)
3341            {
3342              matrix zero[nrows(L[1])][nrows(L[5])];
3343              g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
3344              n_a[1]=#[2];
3345              n_b1=n_a[1],n_c[1];
3346              n_b[1]=n_b1;
3347              g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));
3348              if (size(#[1])!=0)
3349                {
3350                  J_A=#[1];
3351                  /*resolution of L[3]*/
3352                  nr=nrows(J_A[1]);
3353                  J_B=concat(J_A[1],matrix(0,nr,nrows(L[3])-nrows(L[1])));
3354                }
3355            }
3356          else
3357            {
3358              n_b=n_c[1];
3359              g_BC=unitmat(ncols(L[5]));
3360            }
3361        }
[7fe9f8b]3362    }
3363  else// L[5]=list();
3364    {
[0e8a5a]3365      if (size(#[2])!=0)
3366        {
3367          n_a[1]=#[2];
3368          n_b=n_a[1];
3369          g_AB=unitmat(size(n_b[1]));
3370          if (size(#[1])!=0)
3371            {
3372              J_A=#[1];
3373              J_B[1]=J_A[1];// resolution of L[3] equals that of L[1]
3374            }
3375        }
[7fe9f8b]3376    }
3377  list out=(J_A[1],g_AB[1],J_B[1],g_BC[1],J_C[1],n_a[1],n_b[1],n_c[1]);
3378  return (out);
3379}
3380
3381////////////////////////////////////////////////////////////////////////////////////
3382
3383static proc assemblingDoubleComplexes(list L)
3384{
[0e8a5a]3385  /* The input is the output of VdStrictDoubleComplexes, we assemble the
3386     resolutions of the L[i][2][3][1] to obtain a V_d-strict free Cartan-Eilenberg
3387     resolution with modules P^i_j (1<=i<=size(L), j>=0) for the seqeunce
3388     coker(L[1][2][3][1])->...->coker(L[size(L)][2][3][1])*/
[7fe9f8b]3389  list out;
3390  int i,j,k,l,oldj,newj,nr,nc;
3391  for (i=1; i<=size(L); i++)
[0e8a5a]3392    {
3393      out[i]=list(list());
3394      out[i][1][1]=ncols(L[i][2][3][1]);//rank of module P^i_0
3395      if (size(L[i][2][5][1])!=0)
3396        {
3397          /*horizontal differential P^i_0->P^(i+1)_0*/
3398          nc=ncols(L[i][2][5][1]);
3399          out[i][1][4]=prodr(ncols(L[i][2][3][1])-ncols(L[i][2][5][1]),nc);
3400        }
[7fe9f8b]3401      else
[0e8a5a]3402        {
3403          /*horizontal differential P^i_0->0*/
3404          out[i][1][4]=matrix(0,ncols(L[i][2][3][1]),1);
3405        }
3406      oldj=newj;
3407      for (j=1; j<=size(L[i][2][3]);j++)
3408        {
3409          out[i][j][2]=L[i][2][7][j];//shift vector of P^i_{j-1}
3410          if (size(L[i][2][3][j])==0)
3411            {
3412              newj =j;
3413              break;
3414            }
3415          out[i][j+1]=list();
3416          out[i][j+1][1]=nrows(L[i][2][3][j]);//rank of the module P^i_j
3417          out[i][j+1][3]=L[i][2][3][j];//vertical differential P^i_j->P^(i+1)_j
3418          if (size(L[i][2][5][j])!=0)
3419            {
3420              //horizonal differential P^i_j->P^(i-1)_j
3421              nr=nrows(L[i][2][3][j])-nrows(L[i][2][5][j]);
3422              out[i][j+1][4]=(-1)^j*prodr(nr,nrows(L[i][2][5][j]));
3423            }
3424          else
3425            {
3426              /*horizontal differential P^i_j->P^(i-1)_j*/
3427              out[i][j+1][4]=matrix(0,nrows(L[i][2][3][j]),1);
3428            }
3429          if(j==size(L[i][2][3]))
3430            {
3431              out[i][j+1][2]=L[i][2][7][j+1];//shift vector of P^i_j
3432              newj=j+1;
3433            }
3434        }
3435      if (i>1)
3436        {
3437
3438          for (k=1; k<=Min(list(oldj,newj)); k++)
3439            {
3440              /*horizonal differential P^(i-1)_(k-1)->P^i_(k-1)*/
3441              nr=nrows(out[i-1][k][4]);
3442              out[i-1][k][4]=matrix(out[i-1][k][4],nr,out[i][k][1]);
3443            }
3444          for (k=newj+1; k<=oldj; k++)
3445            {
3446              /*no differential needed*/
3447              out[i-1][k]=delete(out[i-1][k],4);
3448            }
3449        }
[7fe9f8b]3450    }
3451  return (out);
3452}
3453
[0e8a5a]3454////////////////////////////////////////////////////////////////////////////////////
[7fe9f8b]3455
3456static proc totalComplex(list L);
3457{
3458  /* Input is the output of assemblingDoubleComplexes.
[0e8a5a]3459     We obtain a complex C^1[m^1]->...->C^(r)[m^r]  with differentials d^i and
3460     shift  vectors m^i (where C^r is placed in degree size(L)-1).
3461     This complex is dercribed in the list out as follows:
3462     rank(C^i)=out[3*i-2]; m_i=out[3*i-1] and d^i=out[3*i]*/
[7fe9f8b]3463  list out;intvec rem1;intvec v; list remsize; int emp;
3464  int i; int j; int c; int d; matrix M; int k; int l;
3465  int n=nvars(basering) div 2;
3466  list K;
3467  for (i=1; i<=n+1; i++)
3468    {
[0e8a5a]3469      K[i]=list();
[7fe9f8b]3470    }
[0e8a5a]3471  L=K+L;
3472  for (i=1; i<=size(L); i++)
[7fe9f8b]3473    {
[0e8a5a]3474      emp=0;
3475      if (size(L[i])!=0)
[7fe9f8b]3476        {
[0e8a5a]3477          out[3*i-2]=L[i][1][1];
3478          v=L[i][1][1];
3479          rem1=L[i][1][2];
[7fe9f8b]3480          emp=1;
3481        }
[0e8a5a]3482      else
[7fe9f8b]3483        {
[0e8a5a]3484          out[3*i-2]=0;
3485          v=0;
[7fe9f8b]3486        }
[0e8a5a]3487      for (j=i+1; j<=size(L); j++)
3488        {
3489          if (size(L[j])>=j-i+1)
3490            {
3491              out[3*i-2]=out[3*i-2]+L[j][j-i+1][1];
3492              if (emp==0)
3493                {
3494                  rem1=L[j][j-i+1][2];
3495                  emp=1;
3496                }
3497              else
3498                {
3499                  rem1=rem1,L[j][j-i+1][2];
3500                }
3501              v[size(v)+1]=L[j][j-i+1][1];
3502            }
3503          else
3504            {
3505              v[size(v)+1]=0;
3506            }
3507        }
3508      out[3*i-1]=rem1;
3509      v[size(v)+1]=0;
3510      remsize[i]=v;
[7fe9f8b]3511    }
3512  int o1;
3513  int o2;
3514  for (i=1; i<=size(L)-1; i++)
3515    {
[0e8a5a]3516      o1=1;
3517      o2=1;
3518      if (size(out[3*i-2])!=0)
3519        {
3520          o1=out[3*i-2];
3521        }
3522      if (size(out[3*i+1])!=0)
[7fe9f8b]3523        {
[0e8a5a]3524          o2=out[3*i+1];
[7fe9f8b]3525        }
[0e8a5a]3526      M=matrix(0,o1,o2);
3527      if (size(L[i])!=0)
[7fe9f8b]3528        {
[0e8a5a]3529          if (size(L[i][1][4])!=0)
[7fe9f8b]3530            {
[0e8a5a]3531              M=matrix(L[i][1][4],o1,o2);
[7fe9f8b]3532            }
[0e8a5a]3533        }
3534      c=remsize[i][1];
3535      for (j=i+1; j<=size(L); j++)
3536        {
3537          if (remsize[i][j-i+1]!=0)
3538            {
3539              for (k=c+1; k<=c+remsize[i][j-i+1]; k++)
3540                {
3541                  for (l=d+1; l<=d+remsize[i+1][j-i];l++)
3542                    {
3543                      M[k,l]=L[j][j-i+1][3][(k-c),(l-d)];
3544                    }
3545                }
3546              d=d+remsize[i+1][j-i];
3547              if (remsize[i+1][j-i+1]!=0)
3548                {
3549                  for (k=c+1; k<=c+remsize[i][j-i+1]; k++)
3550                    {
3551                      for (l=d+1; l<=d+remsize[i+1][j-i+1];l++)
3552                        {
3553                          M[k,l]=L[j][j-i+1][4][k-c,l-d];
3554                        }
3555                    }
3556                  c=c+remsize[i][j-i+1];
3557                }
3558            }
3559          else
3560            {
3561              d=d+remsize[i+1][j-i];
3562            }
3563        }
3564      out[3*i]=M;
3565      d=0; c=0;
[7fe9f8b]3566    }
3567  out[3*size(L)]=matrix(0,out[3*size(L)-2],1);
3568  return (out);
3569
3570}
3571
[0e8a5a]3572////////////////////////////////////////////////////////////////////////////////////
[7fe9f8b]3573//COMPUTATION OF THE BLOBAL B-FUNCTION
[0e8a5a]3574////////////////////////////////////////////////////////////////////////////////////
[7fe9f8b]3575
3576static proc globalBFun(list L,list #)
3577{
[0e8a5a]3578  /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]),
3579    where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an
3580    intvec of size n_i.
3581    We compute bounds for the minimal and maximal integer roots of the b-functions
3582    of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def.
3583    6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal
3584    intersection (cf. Remark 6.1.7 in [R] 2012).
3585    This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this
3586    proc is only called from the procedure deRhamCohomology and the input comes
3587    originally from the procedure toVdstrictFreeComplex*/
[7fe9f8b]3588  if (size(#)==0)//# may contain the Syzstring
[0e8a5a]3589    {
3590      string Syzstring="Sres";
3591    }
[7fe9f8b]3592  else
[0e8a5a]3593    {
3594      string Syzstring=#[1];
3595    }
[7fe9f8b]3596  int i,j;
3597  def W=basering;
3598  int n=nvars(W) div 2;
3599  list G0;
3600  ideal I;
3601  for (j=1; j<=size(L); j++)
3602    {
[0e8a5a]3603      G0[j]=list();
3604      for (i=1; i<=ncols(L[j][1]); i++)
3605        {
3606          G0[j][i]=I;
3607        }
[7fe9f8b]3608    }
3609  list out;
3610  ideal I; poly f;
3611  intvec i1;
3612  for (j=1; j<=size(L); j++)
3613    {
[0e8a5a]3614      /*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict
3615        Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1]
3616        is already a V_d-strict Groebner basis, because it was obtained by the
3617        procedure toVdStrictFreeComplex*/
3618      if (L[j][2]!=intvec(0:size(L[j][2])) or Syzstring=="noCE")
3619              {
3620          if (Syzstring=="Vdres")
3621            {
3622              L[j][1]=VdStrictGB(L[j][1],n);
3623            }
3624          else
3625            {
3626              def HomWeyl=makeHomogenizedWeyl(n);
3627              setring HomWeyl;
3628              list L=fetch(W,L);
3629              L[j][1]=nHomogenize(L[j][1]);
3630              L[j][1]=transpose(matrix(slimgb(transpose(L[j][1]))));
3631              L[j][1]=subst(L[j][1],h,1);
3632              setring W;
3633              L=fetch(HomWeyl,L);
3634              kill HomWeyl;
3635            }
3636        }
3637      for (i=1; i<=ncols(L[j][1]); i++)
3638        {
3639          G0[j][i]=I;
3640        }
3641      for (i=1; i<=nrows(L[j][1]); i++)
3642        {
3643          /*computes the terms of maximal V_d-degree of the biggest non-zero
3644            component of submat(L[j][1],i,(1..ncols(L[j][1])))*/
3645          i1=(1..ncols(L[j][1]));
3646          out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);
3647          // f=L[j][1][i,out[2]];
3648          G0[j][out[2]]=G0[j][out[2]],out[1];
3649          G0[j][out[2]]=compress(G0[j][out[2]]);
3650        }
[7fe9f8b]3651    }
3652  list save;
3653  int l;
3654  list weights;
[0e8a5a]3655  /*bFctIdealModified computes the intersection of G0[j][i] and
3656    x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/
[7fe9f8b]3657  for (j=1; j<=size(G0); j++)
3658    {
[0e8a5a]3659      for (i=1; i<=size(G0[j]); i++)
3660        {
3661          G0[j][i]=bFctIdealModified(G0[j][i]);
3662        }
3663      for (i=1; i<=size(G0[j]); i++)
[7fe9f8b]3664        {
[0e8a5a]3665          weights=list();
3666          if (size(G0[j][i])!=0)
3667            {
3668              for (l=i; l<=size(G0[j]); l++)
3669                {
3670                  weights[size(weights)+1]=L[j][2][l];
3671                }
3672              G0[j][i]=list(G0[j][i][1]+Min(weights),G0[j][i][2]+Max(weights));
3673            }
[7fe9f8b]3674        }
3675    }
3676  list allmin;
3677  list allmax;
3678  for (j=1; j<=size(G0); j++)
3679    {
[0e8a5a]3680      for (i=1; i<=size(G0[j]); i++)
3681        {
3682          if (size(G0[j][i])!=0)
3683            {
3684              allmin[size(allmin)+1]=G0[j][i][1];
3685              allmax[size(allmax)+1]=G0[j][i][2];
3686            }
3687        }
[7fe9f8b]3688    }
3689  list minmax=list(Min(allmin),Max(allmax));
3690  return(minmax);
3691}
3692
[0e8a5a]3693////////////////////////////////////////////////////////////////////////////////////
[7fe9f8b]3694
3695static proc exactGlobalBFun(list L,list #)
3696{
[0e8a5a]3697  /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]),
3698    where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an
3699    intvec of size n_i.
3700    We compute bounds for the minimal and maximal integer roots of the b-functions
3701    of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def.
3702    6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal
3703    intersection (cf. Remark 6.1.7 in [R] 2012).
3704    This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this
3705    proc is only called from the procedure deRhamCohomology and the input comes
3706    originally from the procedure toVdstrictFreeComplex*/
[7fe9f8b]3707  if (size(#)==0)//# may contain the Syzstring
[0e8a5a]3708    {
3709      string Syzstring="Sres";
3710    }
[7fe9f8b]3711  else
[0e8a5a]3712    {
3713      string Syzstring=#[1];
3714    }
[7fe9f8b]3715  int i,j,k;
3716  def W=basering;
3717  int n=nvars(W) div 2;
3718  list G0;
3719  ideal I;
3720  for (j=1; j<=size(L); j++)
3721    {
[0e8a5a]3722      G0[j]=list();
3723      for (i=1; i<=ncols(L[j][1]); i++)
3724        {
3725          G0[j][i]=I;
3726        }
[7fe9f8b]3727    }
3728  list out;
3729  matrix M;
3730  ideal I; poly f;
3731  intvec i1;
3732  for (j=1; j<=size(L); j++)
[0e8a5a]3733    {
3734      M=L[j][1];
3735      /*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict
3736        Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1]
3737        is already a V_d-strict Groebner basis, because it was obtained by the
3738        procedure toVdStrictFreeComplex*/
3739      for (k=1; k<=ncols(L[j][1]); k++)
[7fe9f8b]3740        {
[0e8a5a]3741          L[j][1]=permcol(M,1,k);
3742          if (Syzstring=="Vdres")
3743            {
3744              L[j][1]=VdStrictGB(L[j][1],n);
3745            }
3746          else
3747            {
3748              def HomWeyl=makeHomogenizedWeyl(n);
3749              setring HomWeyl;
3750              list L=fetch(W,L);
3751              L[j][1]=nHomogenize(L[j][1]);
3752              L[j][1]=transpose(matrix(slimgb(transpose(L[j][1]))));
3753              L[j][1]=subst(L[j][1],h,1);
3754              setring W;
3755              L=fetch(HomWeyl,L);
3756              kill HomWeyl;
3757            }
3758          for (i=1; i<=nrows(L[j][1]); i++)
3759            {
3760              /*computes the terms of maximal V_d-degree of the biggest non-zero
3761                component of submat(L[j][1],i,(1..ncols(L[j][1])))*/
3762              i1=(1..ncols(L[j][1]));
3763              out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);
3764              if (out[2]==1)
3765                {
3766                  G0[j][k]=G0[j][k],out[1];
3767                  G0[j][k]=compress(G0[j][k]);
3768                }
3769            }
[7fe9f8b]3770        }
3771    }
3772  list save;
3773  int l;
3774  list weights;
[0e8a5a]3775  /*bFctIdealModified computes the intersection of G0[j][i] and
3776    x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/
3777  for (j=1; j<=size(G0); j++)
3778    {
3779      for (i=1; i<=size(G0[j]); i++)
3780        {
3781          G0[j][i]=bFctIdealModified(G0[j][i]);
3782        }
3783      for (i=1; i<=size(G0[j]); i++)
3784        {
3785          if (size(G0[j][i])!=0)
3786            {
3787              G0[j][i]=list(G0[j][i][1]+L[j][2][i],G0[j][i][2]+L[j][2][i]);
3788            }
3789        }
3790    }
3791  list allmin;
3792  list allmax;
[7fe9f8b]3793  for (j=1; j<=size(G0); j++)
3794    {
[0e8a5a]3795      for (i=1; i<=size(G0[j]); i++)
3796        {
3797          if (size(G0[j][i])!=0)
3798            {
3799              allmin[size(allmin)+1]=G0[j][i][1];
3800              allmax[size(allmax)+1]=G0[j][i][2];
3801            }
3802        }
3803    }
3804  list minmax=list(Min(allmin),Max(allmax));
3805  return(minmax);
3806}
3807
3808////////////////////////////////////////////////////////////////////////////////////
3809
3810////////////////////////////////////////////////////////////////////////////////////
3811
3812static proc exactGlobalBFunIntegration(list L,list #)
3813{
3814  /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]),
3815    where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an
3816    intvec of size n_i.
3817    We compute bounds for the minimal and maximal integer roots of the b-functions
3818    of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def.
3819    6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal
3820    intersection (cf. Remark 6.1.7 in [R] 2012).
3821    This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this
3822    proc is only called from the procedure deRhamCohomology and the input comes
3823    originally from the procedure toVdstrictFreeComplex*/
3824  string Syzstring="Sres";
3825  int i,j,k;
3826  def W=basering;
3827  int n=nvars(W) div 2;
3828//   def C=makeConverseWeyl(n);
3829//   setring C;
3830//   ideal Jn=x(1);
3831//   for (i=2; i<=nvars(basering) div 2; i++)
3832//     {
3833//       Jn=Jn,var(nvars(basering) div 2 + i);
3834//     }
3835//   for (i=1; i<=nvars(basering) div 2; i++)
3836//     {
3837//       Jn=Jn,var(i);
3838//     }
3839//   map transtc=W,Jn;
3840//   list L=transtc(L);
3841  list G0;
3842  ideal I;
3843  for (j=1; j<=size(L); j++)
[7fe9f8b]3844    {
[0e8a5a]3845      G0[j]=list();
3846      for (i=1; i<=ncols(L[j][1]); i++)
3847        {
3848          G0[j][i]=I;
3849        }
3850    }
3851  list out;
3852  matrix M;
3853  ideal I;
3854  poly f;
3855  intvec i1;
3856  for (j=1; j<=size(L); j++)
3857    {
3858      M=L[j][1];
3859      /*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict
3860        Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1]
3861        is already a V_d-strict Groebner basis, because it was obtained by the
3862        procedure toVdStrictFreeComplex*/
3863      for (k=1; k<=ncols(L[j][1]); k++)
3864        {
3865          L[j][1]=permcol(M,1,k);
3866          def HomWeyl=makeHomogenizedWeylTilde(n);
3867          setring HomWeyl;
3868          list L=fetch(W,L);
3869          L[j][1]=nHomogenize(L[j][1]);
3870          L[j][1]=transpose(matrix(slimgb(transpose(L[j][1]))));
3871          L[j][1]=subst(L[j][1],h,1);
3872          setring W;
3873          L=fetch(HomWeyl,L);
3874          kill HomWeyl;
3875          for (i=1; i<=nrows(L[j][1]); i++)
3876            {
3877              /*computes the terms of maximal V_d-degree of the biggest non-zero
3878                component of submat(L[j][1],i,(1..ncols(L[j][1])))*/
3879              i1=(1..ncols(L[j][1]));
[08fa62]3880              out=VdDegTilde(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);//hier koennte es evtl noch einen Fehler geben!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[0e8a5a]3881              f=L[j][1][i,out[2]];
3882              if (out[2]==1)
3883                {
3884                  G0[j][k]=G0[j][k],out[1];
3885                  G0[j][k]=compress(G0[j][k]);
3886                }
3887            }
3888        }
3889    }
3890  list save;
3891  int l;
3892  list weights;
3893  /*bFctIdealModified computes the intersection of G0[j][i] and
3894    x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/
3895  for (j=1; j<=size(G0); j++)
3896    {
3897      for (i=1; i<=size(G0[j]); i++)
3898        {
3899          G0[j][i]=bFctIdealModified(G0[j][i],1);
3900        }
3901      for (i=1; i<=size(G0[j]); i++)
3902        {
3903          if (size(G0[j][i])!=0)
3904            {
3905              G0[j][i]=list(G0[j][i][1]+L[j][2][i],G0[j][i][2]+L[j][2][i]);
3906            }
3907        }
[7fe9f8b]3908    }
3909  list allmin;
3910  list allmax;
3911  for (j=1; j<=size(G0); j++)
3912    {
[0e8a5a]3913      for (i=1; i<=size(G0[j]); i++)
3914        {
3915          if (size(G0[j][i])!=0)
3916            {
3917              allmin[size(allmin)+1]=G0[j][i][1];
3918              allmax[size(allmax)+1]=G0[j][i][2];
3919            }
3920        }
[7fe9f8b]3921    }
3922  list minmax=list(Min(allmin),Max(allmax));
3923  return(minmax);
3924}
3925
3926////////////////////////////////////////////////////////////////////////////////////
3927
[0e8a5a]3928static proc bFctIdealModified (ideal I, list #)
[7fe9f8b]3929{/*modified version of the procedure bfunIdeal from bfun.lib*/
[0e8a5a]3930  int tilde;
3931  if (size(#)!=0)
3932    {
3933      tilde=#[1];
3934    }
[7fe9f8b]3935  def B= basering;
3936  int n = nvars(B) div 2;
3937  intvec w=(1:n);
[0e8a5a]3938  //  if (tilde==0)
3939  // {
3940      I= initialIdealW(I,-w,w);
3941      //    }
3942//   else
3943//     {
3944//       I= initialIdealW(I,w,-w);
3945//     }
[7fe9f8b]3946  poly s; int i;
[0e8a5a]3947  if (tilde==0)
3948    {
3949      for (i=1; i<=n; i++)
3950        {
3951          s=s+x(i)*D(i);
3952        }
3953    }
3954  else
3955    {
3956      for (i=1; i<=n; i++)
3957        {
3958          s=s-D(i)*x(i);
3959        }
3960    }
[7fe9f8b]3961  /*pIntersect computes the intersection on s and I*/
3962  vector b = pIntersect(s,I);
3963  list RL = ringlist(B); RL = RL[1..4];
3964  RL[2] = list(safeVarName("s"));
3965  RL[3] = list(list("dp",intvec(1)),list("C",intvec(0)));
3966  def @S = ring(RL); setring @S;
3967  vector b = imap(B,b);
3968  poly bs = vec2poly(b);
3969  ring r=0,s,dp;
3970  poly bs=imap(@S,bs);
3971  /*find minimal and maximal integer root*/
3972  ideal allfac=factorize(bs,1);
3973  list allfacs;
3974  for (i=1; i<=ncols(allfac); i++)
[0e8a5a]3975    {
3976      allfacs[i]=allfac[i];
3977    }
[7fe9f8b]3978  number testzero;
3979  list zeros;
3980  for (i=1; i<=size(allfacs); i++)
3981    {
[0e8a5a]3982      if (deg(allfacs[i])==1)
3983        {
3984          testzero=number(subst(allfacs[i],s,0))/leadcoef(allfacs[i]);
3985          if (testzero-int(testzero)==0)
3986            {
3987              zeros[size(zeros)+1]=int(-1)*int(testzero);
3988            }
3989        }
[7fe9f8b]3990    }
3991  if (size(zeros)!=0)
[0e8a5a]3992    {
3993      list minmax=(Min(zeros),Max(zeros));
3994    }
[7fe9f8b]3995  else
[0e8a5a]3996    {
3997      list minmax=list();
3998    }
[7fe9f8b]3999  setring B;
4000  return(minmax);
4001}
4002
4003////////////////////////////////////////////////////////////////////////////////////
4004
4005static proc safeVarName (string s)
4006{/* from the library "bfun.lib"*/
4007  string S = "," + charstr(basering) + "," + varstr(basering) + ",";
4008  s = "," + s + ",";
4009  while (find(S,s) <> 0)
4010  {
4011    s[1] = "@";
4012    s = "," + s;
4013  }
4014  s = s[2..size(s)-1];
4015  return(s)
4016}
4017
4018////////////////////////////////////////////////////////////////////////////////////
4019
4020static proc globalBFunOT(list L,list #)
4021{
[0e8a5a]4022  /*this proc is currently not used since globalBFun computes the same output and is
4023    faster, however globalBFun works only for non-zero b-functions!*/
4024  /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]),
4025    where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an
4026    intvec of size n_i.
4027    We compute bounds for the minimal and maximal integer roots of the b-functions
4028    of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def.
4029    6.1.1 in [R]) using Algorithm 6.1.6 in [R].*/
[7fe9f8b]4030  if (size(#)==0)
[0e8a5a]4031    {
4032      string Syzstring="Sres";
4033    }
[7fe9f8b]4034  else
[0e8a5a]4035    {
4036      string Syzstring=#[1];
4037    }
[7fe9f8b]4038  int i; int j;
4039  def W=basering;
4040  int n=nvars(W) div 2;
4041  list G0;
4042  ideal I;
4043  intvec i1;
4044  for (j=1; j<=size(L); j++)
4045    {
[0e8a5a]4046      G0[j]=list();
4047      for (i=1; i<=ncols(L[j][1]); i++)
4048        {
4049          G0[j][i]=I;
4050        }
[7fe9f8b]4051    }
4052  list out;
4053  for (j=1; j<=size(L); j++)
4054    {
[0e8a5a]4055      if (L[j][2]!=intvec(0:size(L[j][2])))
4056              {
4057          if (Syzstring=="Vdres")
4058            {
4059              L[j][1]=VdStrictGB(L[j][1],n);
4060            }
4061          else
4062            {
4063              def HomWeyl=makeHomogenizedWeyl(n);
4064              setring HomWeyl;
4065              list L=fetch(W,L);
4066              L[j][1]=nHomogenize(L[j][1]);
4067              L[j][1]=transpose(matrix(slimgb(transpose(L[j][1]))));
4068              L[j][1]=subst(L[j][1],h,1);
4069              setring W;
4070              L=fetch(HomWeyl,L);
4071              kill HomWeyl;
4072            }
4073        }
4074      for (i=1; i<=nrows(L[j][1]); i++)
4075        {
4076          i1=(1..ncols(L[j][1]));
4077          out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);
4078          G0[j][out[2]][size(G0[j][out[2]])+1]=(out[1]);
4079        }
4080    }
[7fe9f8b]4081  list Data=ringlist(W);
4082  for (i=1; i<=n; i++)
[0e8a5a]4083    {
4084      Data[2][2*n+i]=Data[2][i];
4085      Data[2][3*n+i]=Data[2][n+i];
4086      Data[2][i]="v("+string(i)+")";
4087      Data[2][n+i]="w("+string(i)+")";
4088    }
[7fe9f8b]4089  Data[3][1][1]="M";
4090  intvec mord=(0:16*n^2);
4091  mord[1..2*n]=(1:2*n);
4092  mord[6*n+1..8*n]=(1:2*n);
4093  for (i=0; i<=2*n-2; i++)
[0e8a5a]4094    {
4095      mord[(3+i)*4*n-i]=-1;
4096      mord[(2*n+2+i)*4*n-2*n-i]=-1;
4097    }
[7fe9f8b]4098  Data[3][1][2]=mord;
4099  matrix Ones=UpOneMatrix(4*n);
4100  Data[5]=Ones;
4101  matrix con[2*n][2*n];
4102  Data[6]=transpose(concat(con,transpose(concat(con,Data[6]))));
4103  def Wuv=ring(Data);
4104  setring Wuv;
4105  list G0=imap(W,G0); list G3; poly lterm;intvec lexp;
[0e8a5a]4106  list G1,G2,LL;
4107  intvec e,f;
4108  int  kapp,k,l;
4109  poly h;
[7fe9f8b]4110  ideal I;
4111  for (l=1; l<=size(G0); l++)
[0e8a5a]4112    {
4113      G1[l]=list();  G2[l]=list(); G3[l]=list();
4114      for (i=1; i<=size(G0[l]); i++)
4115        {
4116          for (j=1; j<=ncols(G0[l][i]);j++)
4117            {
4118                    G0[l][i][j]=mHom(G0[l][i][j]);
4119            }
4120          for (j=1; j<=nvars(Wuv) div 4; j++)
4121            {
4122              G0[l][i][size(G0[l][i])+1]=1-v(j)*w(j);
4123            }
4124          G1[l][i]=slimgb(G0[l][i]);
4125          G2[l][i]=I;
4126          G3[l][i]=list();
4127          for (j=1; j<=ncols(G1[l][i]); j++)
4128            {
4129              e=leadexp(G1[l][i][j]);
4130              f=e[1..2*n];
4131              if (f==intvec(0:(2*n)))
4132                {
4133                  for (k=1; k<=n; k++)
4134                    {
4135                      kapp=-e[2*n+k]+e[3*n+k];
4136                      if (kapp>0)
4137                        {
4138                          G1[l][i][j]=(x(k)^kapp)*G1[l][i][j];
4139                        }
4140                      if (kapp<0)
4141                        {
4142                          G1[l][i][j]=(D(k)^(-kapp))*G1[l][i][j];
4143                        }
4144                    }
4145                  G2[l][i][size(G2[l][i])+1]=G1[l][i][j];
4146                  G3[l][i][size(G3[l][i])+1]=list();
4147                  while (G1[l][i][j]!=0)
4148                    {
4149                      lterm=lead(G1[l][i][j]);
4150                      G1[l][i][j]=G1[l][i][j]-lterm;
4151                      lexp=leadexp(lterm);
4152                      lexp=lexp[2*n+1..3*n];
4153                      LL=list(lexp,leadcoef(lterm));
4154                      G3[l][i][size(G3[l][i])][size(G3[l][i][size(G3[l][i])])+1]=LL;
4155                    }
4156                }
4157            }
4158        }
[7fe9f8b]4159    }
4160  ring r=0,(s(1..n)),dp;
4161  ideal I;
4162  map G3forr=Wuv,I;
4163  list G3=G3forr(G3);
4164  poly fs,gs;
4165  int a;
4166  list G4;
4167  for (l=1; l<=size(G3); l++)
4168    {
[0e8a5a]4169      G4[l]=list();
4170      for (i=1; i<=size(G3[l]);i++)
[7fe9f8b]4171        {
[0e8a5a]4172          G4[l][i]=I;
4173
4174          for (j=1; j<=size(G3[l][i]); j++)
[7fe9f8b]4175            {
[0e8a5a]4176              fs=0;
4177              for (k=1; k<=size(G3[l][i][j]); k++)
4178                {
4179                  gs=1;
4180                  for (a=1; a<=n; a++)
4181                    {
4182                      if (G3[l][i][j][k][1][a]!=0)
4183                        {
4184                          gs=gs*permuteVar(list(G3[l][i][j][k][1][a]),a);
4185                        }
4186                    }
4187                  gs=gs*G3[l][i][j][k][2];
4188                  fs=fs+gs;
4189                }
4190              G4[l][i]=G4[l][i],fs;
[7fe9f8b]4191            }
4192        }
4193    }
4194  if (n==1)
4195    {
[0e8a5a]4196      ring rnew=0,t,dp;
[7fe9f8b]4197    }
[0e8a5a]4198  else
[7fe9f8b]4199    {
[0e8a5a]4200      ring rnew=0,(t,s(2..n)),dp;
[7fe9f8b]4201    }
[0e8a5a]4202  ideal Iformap;
4203  Iformap[1]=t;
4204   poly forel=1;
4205   for (i=2; i<=n; i++)
4206     {
4207       Iformap[1]=Iformap[1]-s(i);
4208       Iformap[i]=s(i);
4209       forel=forel*s(i);
4210     }
4211   map rtornew=r,Iformap;
4212   list G4=rtornew(G4);
4213   list getintvecs=fetch(W,L);
4214   ideal J;
4215   option(redSB);
4216   for (l=1; l<=size(G4); l++)
4217     {
4218       J=1;
4219       for (i=1; i<=size(G4[l]); i++)
4220         {
4221           G4[l][i]=eliminate(G4[l][i],forel);
4222           J=intersect(J,G4[l][i]);
4223         }
4224       G4[l]=poly(std(J)[1]);
4225     }
4226   list minmax;
4227   list mini=list();
4228   list maxi=list();
4229   list L=fetch(W,L);
4230   for (i=1; i<=size(G4); i++)
4231     {
[76d26c]4232       minmax[i]=minIntRootD(G4[i],1);
[0e8a5a]4233       if (size(minmax[i])!=0)
4234         {
4235           mini=insert(mini,minmax[i][1]+Min(L[i][2]));
4236           maxi=insert(maxi,minmax[i][2]+Max(L[i][2]));
4237         }
4238     }
4239   mini=Min(mini);
4240   maxi=Max(maxi);
4241   minmax=list(mini[1],maxi[1]);
4242   option(none);
[7fe9f8b]4243  return(minmax);
4244}
4245
4246////////////////////////////////////////////////////////////////////////////////////
4247//COMPUTATION OF THE COHOMOLOGY
4248////////////////////////////////////////////////////////////////////////////////////
4249
4250static proc findCohomology(list L,int le)
4251{
[0e8a5a]4252/*computes the cohomology of the complex (D^i,d^i) given by D^i=C^L[2*i-1] and
[7fe9f8b]4253  d^i=L[2*i]*/
4254  def R=basering;
4255  ring r=0,(x),dp;
4256  list L=imap(R,L);
4257  list out;
4258  int i, ker, im;
4259  matrix S;
4260  option(returnSB);
4261  option(redSB);
4262  for (i=2; i<=size(L); i=i+2)
4263    {
[0e8a5a]4264      if (L[i-1]==0)
4265        {
4266          out[i div 2]=0;
4267          im=0;
4268        }
4269      else
4270        {
4271          S=matrix(syz(transpose(L[i])));
4272          if (S!=matrix(0,nrows(S),ncols(S)))
4273            {
4274              ker=ncols(S);
4275              out[i div 2]=ker-im;
4276              im=L[i-1]-ker;
4277            }
4278          else
4279            {
[08fa62]4280              out[i div 2]=0;////achtung geaendert??????????????????????????????????????????????????!!!!!!!!!!!!!!!!!!!!!!!!!war mal out[i-1]
[0e8a5a]4281              im=L[i-1];
4282            }
4283        }
4284    }
4285  option(none);
4286  while (size(out)>le)
4287    {
4288      out=delete(out,1);
[7fe9f8b]4289    }
[0e8a5a]4290  setring R;
4291  return(out);
4292}
4293
4294////////////////////////////////////////////////////////////////////////////////////
4295
4296
4297static proc findCohomologyDiffForms(list L,int le)
4298{
4299  /*computes the cohomology of the complex (D^i,d^i) given by D^i=C^L[2*i-1] and
4300    d^i=L[2*i]*/
4301  def R=basering;
4302  list outdiffforms=list(var(1));
4303  ring r=0,(x),dp;
4304  list L=imap(R,L);
4305  list out;
4306  list outdiffforms;
4307  int i, ker, im, j;
4308  matrix S;
4309  matrix concreteimage=matrix(0);
4310  module concreteimagemod=concreteimage;
4311  option(returnSB);
4312  option(redSB);
4313  matrix redS;
4314  for (i=2; i<=size(L); i=i+2)
[7fe9f8b]4315    {
[0e8a5a]4316      if (L[i-1]==0)
4317        {
4318          out[i div 2]=0;
4319          im=0;
4320          concreteimage=matrix(0);
4321          concreteimagemod=concreteimage;
4322          outdiffforms[i div 2]=list();
4323        }
[7fe9f8b]4324      else
[0e8a5a]4325        {
4326          S=matrix(transpose(syz(transpose(L[i]))));
4327          if (S!=matrix(0,nrows(S),ncols(S)))
4328            {
4329              ker=nrows(S);
4330              out[i div 2]=ker-im;
4331              if(out[i div 2]==0)
4332                {
4333                  outdiffforms[i div 2]=list();
4334                }
4335              else
4336                {
4337                  outdiffforms[i div 2]=list();
4338                  if (concreteimage==matrix(0))
4339                    {
4340                      for (j=1; j<=nrows(S); j++)
4341                        {
4342                          outdiffforms[ i div 2][j]=submat(S,j,intvec(1..ncols(S)));
4343                        }
4344                    }
4345                  else
4346                    {
4347                      redS=transpose(std(reduce(transpose(S),concreteimagemod)));
4348                      for (j=1; j<=nrows(redS); j++)
4349                        {
4350                          if (submat(redS,j, intvec(1..ncols(redS)))!=matrix(0,1,ncols(redS)))
4351                            {
4352                              outdiffforms[i div 2][size(outdiffforms[i div 2])+1]=submat(redS,j, intvec(1..ncols(redS)));
4353                            }
4354                        }
4355                    }
4356                }
4357              im=L[i-1]-ker;
4358              concreteimagemod=std(transpose(L[i]));
4359              concreteimage=concreteimagemod;
4360              concreteimage=transpose(concreteimage);
4361
4362
4363              //concreteimage=transpose(std(transpose(L[i])));//Achtung:hier wieder das Problem mit no Standard basis!!!!!!!!!!!!!
4364            }
4365          else
4366            {
4367              out[i div 2]=0;
4368              outdiffforms[i div 2]=0;
4369              im=L[i-1];
4370              concreteimagemod=std(transpose(L[i]));
4371              concreteimage=concreteimagemod;
4372              concreteimage=transpose(concreteimage);
4373              //concreteimage=transpose(std(transpose(L[i])));
4374            }
4375        }
[7fe9f8b]4376    }
4377  option(none);
4378  while (size(out)>le)
[0e8a5a]4379    {
4380      out=delete(out,1);
4381      outdiffforms=delete(outdiffforms,1);
4382    }
[7fe9f8b]4383  setring R;
[0e8a5a]4384  outdiffforms=imap(r,outdiffforms);
4385  list outall=list(out,outdiffforms);
4386  option(noredSB);
4387  option(noreturnSB);
4388  return(outall);
[7fe9f8b]4389}
4390
[0e8a5a]4391
4392
[7fe9f8b]4393////////////////////////////////////////////////////////////////////////////////////
4394//AUXILIARY PROCEDURES
4395////////////////////////////////////////////////////////////////////////////////////
4396
[0e8a5a]4397static proc findPreimage(matrix m, matrix n)
[7fe9f8b]4398{
[0e8a5a]4399  def W=basering;//input wird in spaltenform angenommen, output in zeilenform
4400  list rl=ringlist(W);
4401  list rlnew=rl;
4402  rlnew[3][1]=rl[3][2];
4403  rlnew[3][2]=rl[3][1];
4404  def Wnew=ring(rlnew);
4405  setring Wnew;
4406  matrix m=imap(W,m);
4407  matrix n=imap(W,n);
4408  def Opp=opposite(Wnew);
4409  setring Opp;
4410  matrix m=oppose(Wnew,m);
4411  matrix n=oppose(Wnew,n);
4412  option(redSB);
4413  //matrix m=imap(W,m);
4414  //  matrix n=imap(W,n);
4415  int i;
4416  matrix preim;
4417  if (n!=matrix(0,nrows(n),ncols(n)))
4418    {
4419      matrix con=concat(m,n);
4420      matrix s=syz(con);
4421      for (i=1; i<=ncols(s); i++)
4422        {
4423          if (s[nrows(s),i]==1)
4424            {
4425              preim=(-1)*submat(s,1..ncols(m),i);
4426              break;
4427            }
4428        }
4429    }
4430  else
4431    {
4432      matrix s=syz(m);
4433      preim=submat(s,1..ncols(m),1);
4434    }
4435  option(noredSB);
4436  setring Wnew;
4437  matrix preim=oppose(Opp,preim);
4438  setring W;
4439  matrix preim=imap(Wnew,preim);
4440  return(transpose(preim));
4441}
4442
4443////////////////////////////////////////////////////////////////////////////////////
4444
4445static proc divdr(matrix m,matrix n, list #)
4446{
4447  if (n!=matrix(0,nrows(n),ncols(n)))
4448    {
4449      m=transpose(m);
4450      n=transpose(n);
4451      matrix con=concat(m,n);
4452      matrix s=syz(con);
4453      s=submat(s,1..ncols(m),1..ncols(s));
4454      s=transpose(compress(s));
4455    }
4456  else
4457    {
4458      matrix s=transpose(syz(transpose(m)));
4459    }
4460  int i;
4461  matrix g;
4462  matrix sm;
4463  if (size(#)!=0)
4464    {
4465      for (i=1; i<=nrows(s); i++)
4466        {
4467          g=deletecol(transpose(s),i);
4468          sm=transpose(submat(s,i,intvec(1..ncols(s))));
4469          sm=reduce(sm,slimgb(g));
4470          if (sm==matrix(0,nrows(sm),ncols(sm)))
4471            {
4472              s=g;
4473              s=transpose(s);
4474              i=i-1;
4475            }
4476        }
4477    }
[7fe9f8b]4478  return(s);
4479}
4480////////////////////////////////////////////////////////////////////////////////////
4481
4482static proc matrixLift(matrix M,matrix N)
4483{
4484  intvec v=option(get);
4485  option(none);
4486  matrix l=transpose(lift(transpose(M),transpose(N)));
4487  option(set,v);
4488  return(l);
4489}
4490
4491////////////////////////////////////////////////////////////////////////////////////
4492
4493static proc VdStrictGB (matrix M,int d,list #)
4494"USAGE:VdStrictGB(M,d[,v]); M a matrix, d an integer, v an optional intvec
4495ASSUME:-basering is the nth Weyl algebra D_n @*
4496       -1<=d<=n @*
4497       -v (if given) is the shift vector on the range of M (in particular,
4498        size(v)=ncols(M)); otherwise v is assumed to be the zero shift vector
[0e8a5a]4499RETURN:matrix N; the rows of N form a V_d-strict Groebner basis with respect to v
[7fe9f8b]4500       for the module generated by the rows of M
4501"
[0e8a5a]4502{
[7fe9f8b]4503  if (M==matrix(0,nrows(M),ncols(M)))
4504    {
4505      return (matrix(0,1,ncols(M)));
4506    }
4507  intvec op=option(get);
4508  def W =basering;
4509  int ncM=ncols(M);
4510  list Data=ringlist(W);
4511  Data[2]=list("nhv")+Data[2];
4512  Data[3][3]=Data[3][1];
4513  Data[3][1]=list("dp",intvec(1));
4514  matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2]));
4515  Data[5]=re;
4516  int k,l;
4517  Data[6]=transpose(concat(matrix(0,1,1),transpose(concat(matrix(0,1,1),Data[6]))));
4518  def Whom=ring(Data);// D_n[nhv] with the new commuative variable nhv
4519  setring Whom;
4520  matrix Mnew=imap(W,M);
4521  intvec v;
4522  if (size(#)!=0)
4523    {
4524      v=#[1];
4525    }
4526  if (size(v) < ncM)
4527    {
4528      v=v,0:(ncM-size(v));
4529    }
4530  Mnew=homogenize(Mnew, d, v);//homogenization of M with respect to the new variable
4531  Mnew=transpose(Mnew);
4532  Mnew=slimgb(Mnew);// computes a Groebner basis of the homogenzition of M
4533  Mnew=subst(Mnew,nhv,1);// substitution of 1 gives V_d-strict Groebner basis  of M
4534  Mnew=compress(Mnew);
4535  Mnew=transpose(Mnew);
4536  setring W;
4537  M=imap(Whom,Mnew);
4538  option(set,op);
4539  return(M);
4540}
4541
4542////////////////////////////////////////////////////////////////////////////////////
4543
4544static proc VdNormalForm(matrix F,matrix M,int d,intvec v,list #)
[0e8a5a]4545"USAGE:VdNormalForm(F,M,d,v[,w]); F and M matrices, d int, v intvec, w an optional
[7fe9f8b]4546       intvec
4547ASSUME:-basering is the nth Weyl algebra D_n @*
4548       -F a n_1 x n_2-matrix and M a m_1 x m_2-matrix with m_2<=n_2 @*
4549       -d is an integer between 1 and n @*
4550       -v is a shift vector for D_n^(m_2) and hence size(v)=m_2 @*
4551       -w is a shift vector for D_n^(m_1-m_2) and hence size(v)=m_1-m_2 @*
4552RETURN:a n_1 x n_2-matrix N such that:@*
[0e8a5a]4553       -If no optional intvec w is given:(N[i,1],..,N[i,m_2]) is a V_d-strict normal
4554        form of (F[i,1],...,F[i,m_2]) with respect to a V_d-strict Groebner basis of
[7fe9f8b]4555        the rows of M and the shift vector v
[0e8a5a]4556       -If w is given:(N[i,1],..,N[i,m_2]) is chosen such that
[7fe9f8b]4557        Vddeg((N[i,1],...,N[i,m_2])[v])<=Vddeg((F[i,m_2+1],...,F[i,m_1])[v]);
4558       -N[i,j]=F[i,j] for j>m_2
4559"
[0e8a5a]4560{
[7fe9f8b]4561  int SBcom;
4562  def W =basering;
4563  int c=ncols(M);
4564  matrix keepF=F;
4565  if (size(#)!=0)
[0e8a5a]4566    {
4567      intvec w=#[1];
4568    }
[7fe9f8b]4569  F=submat(F,intvec(1..nrows(F)),intvec(1..c));
4570  list Data=ringlist(W);
4571  Data[2]=list("nhv")+Data[2];
4572  Data[3][3]=Data[3][1];
4573  Data[3][1]=list("dp",intvec(1));
4574  matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2]));
4575  Data[5]=re;
4576  int k,l,nr,nc;
4577  matrix rep[size(Data[2])][size(Data[2])];
4578  for (l=size(Data[2])-1;l>=1; l--)
4579    {
[0e8a5a]4580      for (k=l-1; k>=1;k--)
4581        {
4582          rep[k+1,l+1]=Data[6][k,l];
4583        }
[7fe9f8b]4584    }
4585  Data[6]=rep;
4586  def Whom=ring(Data);//new ring D_n[nvh] this new commuative variable nhv
4587  setring Whom;
4588  matrix Mnew=imap(W,M);
4589  list forMnew=homogenize(Mnew,d,v,1);//commputes homogenization of M;
4590  Mnew=forMnew[1];
4591  int rightexp=forMnew[2];
4592  matrix Fnew=imap(W,F);
4593  matrix keepF=imap(W,keepF);
4594  matrix Fb;
4595  int cc;
4596  intvec i1,i2;
4597  matrix zeromat,subm1,subm2,zeromat2;
4598  for (l=1; l<=nrows(Fnew); l++)
4599    {
[0e8a5a]4600      if (size(#)!=0)
4601        {
4602          subm2=submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF)));
4603          zeromat2=matrix(0,1,ncols(subm2));
4604          if (submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF)))==zeromat2)
4605            {
4606              for (cc=1; cc<=ncols(Fnew); c++)
4607                {
4608                  Fnew[l,cc]=0;
4609                }
4610            }
4611          i1=intvec(1..ncols(Fnew));
4612          subm1=submat(Fnew,l,i1);
4613          subm2=submat(keepF,l,(ncols(Fnew)+1)..ncols(keepF));
4614          zeromat=matrix(0,1,ncols(Fnew));
4615          if (VdDegnhv(subm1,d,v)>VdDegnhv(subm2,d,w)
4616              and submat(Fnew,l,intvec(1..ncols(Fnew)))!=zeromat)
4617            {
[08fa62]4618              //print("Reduzierung des V_d-Grades noetig");
[0e8a5a]4619              /*We need to reduce the V_d-degree. First we homogenize the
4620                lth row of Fnew*/
4621              Fb=homogenize(subm1,d,v)*(nhv^rightexp);
4622              if (SBcom==0)
4623                {
4624                  /*computes a V_d-strict standard basis*/
4625                  Mnew=slimgb(transpose(Mnew));//
4626                  SBcom=1;
4627                }
4628              /*computes a V_d-strict normal form for FB*/
4629              Fb=transpose(reduce(transpose(Fb),Mnew));
4630              if (VdDegnhv(Fb,d,v)> VdDegnhv(subm2,d,w)
4631                  and Fb!=matrix(0,nrows(Fb),ncols(Fb)))//should not happen
4632                {
4633                  //print("Reduzierung fehlgeschlagen!!!!!!!!!!!!!!!!");
4634                }
4635            }
4636          else
4637            {
4638              /*condition on V_ddeg already satisfied -> no normal form
4639                computation is needed*/
4640              Fb=submat(Fnew,l,intvec(1..ncols(Fnew)));
4641            }
4642        }
4643      else
4644        {
4645          Fb=homogenize(submat(Fnew,l,intvec(1..ncols(Fnew))),d,v);
4646          if (SBcom==0)
4647            {
4648              Mnew=slimgb(transpose(Mnew));// computes a V_d-strict Groebner basis
4649              SBcom=1;
4650            }
4651          Fb=transpose(reduce(transpose(Fb),Mnew));//normal form
4652        }
4653      for (k=1; k<=ncols(Fnew);k++)
4654        {
4655          Fnew[l,k]=Fb[1,k];
4656        }
[7fe9f8b]4657    }
4658  Fnew=subst(Fnew,nhv,1);//obtain normal form in D_n
4659  setring W;
4660  F=imap(Whom,Fnew);
4661  return(F);
4662}
4663
4664////////////////////////////////////////////////////////////////////////////////////
4665
4666static proc homogenize (matrix M,int d,intvec v,list #)
4667{
4668  /* we compute the F[v]-homogenization of each row of M (cf. Def. 3.4 in [OT])*/
4669  if (M==matrix(0,nrows(M),ncols(M)))
[0e8a5a]4670    {
4671      return(M);
4672    }
[7fe9f8b]4673  int i,l,s, kmin, nhvexp;
[0e8a5a]4674  poly f;
[7fe9f8b]4675  intvec vnm;
4676  list findmin,maxnhv,rempoly,remk,rem1,rem2;
4677  int n=(nvars(basering)-1) div 2;
[0e8a5a]4678  for (int k=1; k<=nrows(M); k++)
4679    {
4680      for (l=1; l<=ncols (M); l++)
4681        {
4682          f=M[k,l];
4683          s=size(f);
4684          for (i=1; i<=s; i++)
4685            {
4686              vnm=leadexp(f);
4687              vnm=vnm[n+2..n+d+1]-vnm[2..d+1];
4688              kmin=sum(vnm)+v[l];
4689              rem1[size(rem1)+1]=lead(f);
4690              rem2[size(rem2)+1]=kmin;
4691              findmin=insert(findmin,kmin);
4692              f=f-lead(f);
4693            }
4694          rempoly[l]=rem1;
4695          remk[l]=rem2;
4696          rem1=list();
4697          rem2=list();
4698        }
4699      if (size(findmin)!=0)
4700        {
4701          kmin=Min(findmin);
4702        }
4703      for (l=1; l<=ncols(M); l++)
4704        {
4705          if (M[k,l]!=0)
4706            {
4707              M[k,l]=0;
4708              for (i=1; i<=size(rempoly[l]);i++)
4709                {
4710                  nhvexp=remk[l][i]-kmin;
4711                  M[k,l]=M[k,l]+nhv^(nhvexp)*rempoly[l][i];
4712                  maxnhv[size(maxnhv)+1]=nhvexp;
4713                }
4714            }
4715        }
4716      rempoly=list();
4717      remk=list();
4718      findmin=list();
4719    }
[7fe9f8b]4720  maxnhv=Max(maxnhv);
4721  nhvexp=maxnhv[1];
4722  if (size(#)!=0)
[0e8a5a]4723    {
4724      return(list(M,nhvexp));//only needed for normal form computations
4725    }
[7fe9f8b]4726  return(M);
4727}
4728
4729////////////////////////////////////////////////////////////////////////////////////
4730
4731static proc soldr (matrix M,matrix N)
[0e8a5a]4732{
4733  /* We compute a ncols(M) x nrows(M)-matrix C such that
4734     C[i,1]M_1+...+C[i,nrows(M)]M_(nrows(M))= e_i mod im(N),
4735     where e_i is the ith basis element on the range of M, M_j denotes the jth row
4736     of M and im(N) is generated by the rows of N */
[7fe9f8b]4737  int n=nrows(M);
4738  int q=ncols(M);
4739  matrix S=concat(transpose(M),transpose(N));
4740  def W=basering;
4741  list Data=ringlist(W);
4742  list Save=Data[3];
4743  Data[3]=list(list("c",0),list("dp",intvec(1..nvars(W))));
4744  def Wmod=ring(Data);
4745  setring Wmod;
4746  matrix Smod=imap(W,S);
4747  matrix E[q][1];
4748  matrix Smod2,Smodnew;
4749  option(returnSB);
4750  int i,j;
4751  for (i=1;i<=q;i++)
[0e8a5a]4752    {
4753      E[i,1]=1;
4754      Smod2=concat(E,Smod);
4755      Smod2=syz(Smod2);
4756      E[i,1]=0;
4757      for (j=1;j<=ncols(Smod2);j++)
4758        {
4759          if (Smod2[1,j]==1)
4760            {
4761              Smodnew=concat(Smodnew,(-1)*(submat(Smod2,intvec(2..n+1),j)));
4762              break;
4763            }
4764        }
[7fe9f8b]4765    }
4766  Smodnew=transpose(submat(Smodnew,intvec(1..n),intvec(2..q+1)));
4767  setring W;
4768  matrix  Snew=imap(Wmod,Smodnew);
4769  option(none);
4770  return (Snew);
4771}
4772
4773////////////////////////////////////////////////////////////////////////////////////
4774
4775static proc prodr (int k,int l)
4776{
4777  if (k==0)
[0e8a5a]4778    {
4779      matrix P=unitmat(l);
4780      return (P);
4781    }
[7fe9f8b]4782  matrix O[l][k];
4783  matrix P=transpose(concat(O,unitmat(l)));
4784  return (P);
4785}
4786
4787////////////////////////////////////////////////////////////////////////////////////
4788
4789static proc VdDeg(matrix M,int d,intvec v,list #)
4790{
4791  /* We assume that the basering it the nth Weyl algebra and that  M is a 1 x r-
[0e8a5a]4792     matrix.
4793     We compute the V_d-deg of M with respect to the shift vector v,
4794     i.e V_ddeg(M)=max (V_ddeg(M_i)+v[i]), where k=V_ddeg(M_i) if k is the minimal
4795     integer, such that M_i can be expressed as a sum of operators
4796     x(1)^(a_1)*...*x(n)^(a_n)*D(1)^(b_1)*...*D(n)^(b_n) with
4797     a_1+..+a_d+k>=b_1+..+b_d*/
[7fe9f8b]4798  int i, j, etoint;
4799  int n=nvars(basering) div 2;
4800  intvec  e;
4801  list findmax;
4802  int c=ncols(M);
4803  poly l;
4804  list positionpoly,positionVd;
4805  for (i=1; i<=c; i++)
4806    {
[0e8a5a]4807      positionpoly[i]=list();
4808      positionVd[i]=list();
4809      while (M[1,i]!=0)
4810        {
4811          l=lead(M[1,i]);
4812          positionpoly[i][size(positionpoly[i])+1]=l;
4813          e=leadexp(l);
4814          e=-e[1..d]+e[n+1..n+d];
4815          e=sum(e)+v[i];
4816          etoint=e[1];
4817          positionVd[i][size(positionVd[i])+1]=etoint;
4818          findmax[size(findmax)+1]=etoint;
4819          M[1,i]=M[1,i]-l;
4820        }
4821    }
4822  if (size(findmax)!=0)
4823    {
4824      int maxVd=Max(findmax);
4825      if (size(#)==0)
4826        {
4827          return (maxVd);
4828        }
4829    }
4830  else // M is 0-modul
4831    {
4832      return(int(0));
4833    }
4834  l=0;
4835  for (i=c; i>=1; i--)
4836    {
4837      for (j=1; j<=size(positionVd[i]); j++)
4838        {
4839          if (positionVd[i][j]==maxVd)
4840            {
4841              l=l+positionpoly[i][j];
4842            }
4843        }
4844      if (l!=0)
4845        {
4846          /*returns the largest component that has maximal V_d-degree
4847            and its terms of maximal V_d-deg (needed for globalBFun)*/
4848          return (list(l,i));
4849        }
4850    }
4851}
4852
4853////////////////////////////////////////////////////////////////////////////////////
4854
4855static proc VdDegTilde(matrix M,int d,intvec v,list #)
4856{
4857  /* We assume that the basering it the nth Weyl algebra and that  M is a 1 x r-
4858     matrix.
4859     We compute the \tilde(V_d)-deg of M with respect to the shift vector v,
4860     i.e \tilde(V_d)deg(M)=max (\tilde(V_d)deg(M_i)+v[i]), where k=\tilde(V_d)deg(M_i) if k is the minimal
4861     integer, such that M_i can be expressed as a sum of operators
4862     x(1)^(a_1)*...*x(n)^(a_n)*D(1)^(b_1)*...*D(n)^(b_n) with
4863     a_1+..+a_d<=b_1+..+b_d+k*/
4864  int i, j, etoint;
4865  int n=nvars(basering) div 2;
4866  intvec  e;
4867  list findmax;
4868  int c=ncols(M);
4869  poly l;
4870  list positionpoly,positionVd;
4871  for (i=1; i<=c; i++)
4872    {
4873      positionpoly[i]=list();
4874      positionVd[i]=list();
4875      while (M[1,i]!=0)
4876        {
4877          l=lead(M[1,i]);
4878          positionpoly[i][size(positionpoly[i])+1]=l;
4879          e=leadexp(l);
4880          e=e[1..d]-e[n+1..n+d];
4881          e=sum(e)+v[i];
4882          etoint=e[1];
4883          positionVd[i][size(positionVd[i])+1]=etoint;
4884          findmax[size(findmax)+1]=etoint;
4885          M[1,i]=M[1,i]-l;
4886        }
[7fe9f8b]4887    }
4888  if (size(findmax)!=0)
4889    {
[0e8a5a]4890      int maxVd=Max(findmax);
4891      if (size(#)==0)
4892        {
4893          return (maxVd);
4894        }
[7fe9f8b]4895    }
4896  else // M is 0-modul
4897    {
[0e8a5a]4898      return(int(0));
[7fe9f8b]4899    }
[0e8a5a]4900  l=0;
4901  for (i=c; i>=1; i--)
[7fe9f8b]4902    {
[0e8a5a]4903      for (j=1; j<=size(positionVd[i]); j++)
4904        {
4905          if (positionVd[i][j]==maxVd)
4906            {
4907              l=l+positionpoly[i][j];
4908            }
4909        }
4910      if (l!=0)
4911        {
4912          /*returns the largest component that has maximal V_d-degree
4913            and its terms of maximal V_d-deg (needed for globalBFun)*/
4914          return (list(l,i));
4915        }
[7fe9f8b]4916    }
4917}
[0e8a5a]4918
[7fe9f8b]4919////////////////////////////////////////////////////////////////////////////////////
4920
4921static proc VdDegnhv(matrix M,int d,intvec v,list #)
4922{
[0e8a5a]4923  /* As the procedure VdDeg, but the basering is the nth Weyl algebra
4924     with a commutative variable nhv*/
[7fe9f8b]4925  int i,j,etoint;
4926  int n=nvars(basering) div 2;
4927  intvec  e;
4928  int etoint;
4929  list findmax;
4930  int c=ncols(M);
4931  poly l;
4932  list positionpoly;
4933  list positionVd;
4934  for (i=1; i<=c; i++)
4935    {
[0e8a5a]4936      positionpoly[i]=list();
4937      positionVd[i]=list();
4938      while (M[1,i]!=0)
4939        {
4940          l=lead(M[1,i]);
4941          positionpoly[i][size(positionpoly[i])+1]=l;
4942          e=leadexp(l);
4943          e=-e[2..d+1]+e[n+2..n+d+1];
4944          e=sum(e)+v[i];
4945          etoint=e[1];
4946          positionVd[i][size(positionVd[i])+1]=etoint;
4947          findmax[size(findmax)+1]=etoint;
4948          M[1,i]=M[1,i]-l;
4949        }
[7fe9f8b]4950    }
4951  if (size(findmax)!=0)
4952    {
[0e8a5a]4953      int maxVd=Max(findmax);
4954      if (size(#)==0)
4955        {
4956          return (maxVd);
4957        }
[7fe9f8b]4958    }
4959  else // M is 0-modul
[0e8a5a]4960    {
4961      return(int(0));
4962    }
[7fe9f8b]4963}
4964
4965////////////////////////////////////////////////////////////////////////////////////
4966
4967static proc deletecol(matrix M,int l)
4968{
[0e8a5a]4969  if (ncols(M)==1)
4970    {
4971      return(M);
4972    }
[7fe9f8b]4973  int s=ncols(M);
4974  if (l==1)
[0e8a5a]4975    {
4976      M=submat(M,(1..nrows(M)),(2..ncols(M)));
4977      return(M);
4978    }
[7fe9f8b]4979  if (l==s)
[0e8a5a]4980    {
4981      M=submat(M,(1..nrows(M)),(1..(ncols(M)-1)));
4982      return(M);
4983    }
[7fe9f8b]4984  intvec v=(1..(l-1)),((l+1)..s);
4985  M=submat(M,(1..nrows(M)),v);
4986  return(M);
4987}
4988
4989////////////////////////////////////////////////////////////////////////////////////
4990
4991static proc mHom(poly f)
4992{/*for globalBFunOT*/
4993  poly g;
4994  poly l;
4995  poly add;
4996  intvec e;
4997  list minint;
4998  list remf;
4999  int i;
5000  int j;
5001  int n=nvars(basering) div 4;
5002  if (f==0)
[0e8a5a]5003    {
5004      return(f);
5005    }
[7fe9f8b]5006  while (f!=0)
[0e8a5a]5007    {
5008      l=lead(f);
5009      e=leadexp(l);
5010      remf[size(remf)+1]=list();
5011      remf[size(remf)][1]=l;
5012      for (i=1; i<=n; i++)
5013        {
5014          remf[size(remf)][i+1]=-e[2*n+i]+e[3*n+i];
5015          if (size(minint)<i)
5016            {
5017              minint[i]=list();
5018            }
5019          minint[i][size(minint[i])+1]=-e[2*n+i]+e[3*n+i];
5020        }
5021      f=f-l;
5022    }
[7fe9f8b]5023  for (i=1; i<=n; i++)
[0e8a5a]5024    {
5025      minint[i]=Min(minint[i]);
5026    }
[7fe9f8b]5027  for (i=1; i<=size(remf); i++)
5028    {
[0e8a5a]5029      add=remf[i][1];
5030      for (j=1; j<=n; j++)
5031        {
5032          add=v(j)^(remf[i][j+1]-minint[j])*add;
5033        }
5034      g=g+add;
[7fe9f8b]5035    }
5036  return (g);
5037}
5038
5039////////////////////////////////////////////////////////////////////////////////////
5040
5041static proc permuteVar(list L,int n)
5042{/*for globalBFunOT*/
5043  if (typeof(L[1])=="intvec")
[0e8a5a]5044    {
5045      intvec v=L[1];
5046    }
[7fe9f8b]5047  else
[0e8a5a]5048    {
5049      intvec v=(1:L[1]),(0:L[1]);
5050    }
[7fe9f8b]5051  int i;int k; int indi=0;
5052  int j;
5053  int s=size(v);
5054  poly e;
5055  intvec fore;
5056  for (i=2; i<=size(v); i=i+2)
[0e8a5a]5057    {
[7fe9f8b]5058
[0e8a5a]5059      if (v[i]!=0)
5060        {
5061           j=i+1;
5062          while (v[j]!=0)
5063            {
5064              j=j+1;
5065            }
5066          v[i]=0;
5067          v[j]=1;
5068          fore=0;
5069          indi=0;
5070          for (k=1; k<=size(v); k++)
5071            {
5072              if (k!=i and k!=j)
5073                {
5074                  if (indi==0)
5075                    {
5076                      indi=1;
5077                      fore[1]=v[k];
5078                    }
5079                  else
5080                    {
5081                      fore[size(fore)+1]=v[k];
5082                    }
5083                }
5084            }
5085          e=e-(j-i)*permutevar(list(fore),n);
[7fe9f8b]5086        }
5087    }
5088  e=e+s(n)^(size(v) div 2);
5089  return (e);
5090}
5091
5092////////////////////////////////////////////////////////////////////////////////////
5093
5094static proc makeHomogenizedWeyl(int n,list #)
5095{
5096  /*modified version of the procedure makeWeyl() from the library nctools.lib*/
5097  /*Creates the nth homogenized Weyl algebra with variables x(1),..,x(n),D(1),..,
[0e8a5a]5098    D(n) and homogenization variable h, i.e. it holds x(i)*D(i)=D(i)*x(1)+h^2.
5099    If # contains on intvec v, we assign weight v[i] to the ith module component.*/
[7fe9f8b]5100  if (n<1)
[0e8a5a]5101    {
5102      print*("Incorrect input");
5103      return();
5104    }
[7fe9f8b]5105  if (n ==1)
[0e8a5a]5106    {
5107      ring @rr = 0,(x(1),D(1),h),dp;
5108    }
[7fe9f8b]5109  else
[0e8a5a]5110    {
5111      ring @rr = 0,(x(1..n),D(1..n),h),dp;
5112    }
[7fe9f8b]5113  setring @rr;
[0e8a5a]5114  int i=0;
[7fe9f8b]5115  if (size(#)==0)
[0e8a5a]5116    {
5117      def @rrr = homogenizedWeyl(i);
5118    }
[7fe9f8b]5119  else
[0e8a5a]5120    {
5121      def @rrr=homogenizedWeyl(i,#);
5122    }
5123  return(@rrr);
5124}
5125
5126////////////////////////////////////////////////////////////////////////////////////
5127
5128static proc makeHomogenizedWeylTilde(int n,list #)
5129{
5130  /*modified version of the procedure makeWeyl() from the library nctools.lib*/
5131  /*Creates the nth homogenized Weyl algebra with variables x(1),..,x(n),D(1),..,
5132    D(n) and homogenization variable h, i.e. it holds x(i)*D(i)=D(i)*x(1)+h^2.
5133    If # contains on intvec v, we assign weight v[i] to the ith module component.*/
5134  if (n<1)
5135    {
5136      print*("Incorrect input");
5137      return();
5138    }
5139  if (n ==1)
5140    {
5141      ring @rr = 0,(x(1),D(1),h),dp;
5142    }
5143  else
5144    {
5145      ring @rr = 0,(x(1..n),D(1..n),h),dp;
5146    }
5147  setring @rr;
5148  int i=1;
5149  if (size(#)==0)
5150    {
5151      def @rrr = homogenizedWeyl(i);
5152    }
5153  else
5154    {
5155      def @rrr=homogenizedWeyl(i,#);
5156    }
5157  return(@rrr);
5158}
5159
5160////////////////////////////////////////////////////////////////////////////////////
5161
5162static proc makeConverseHomogenizedWeylTilde(int n,list #)
5163{
5164  /*modified version of the procedure makeWeyl() from the library nctools.lib*/
5165  /*Creates the nth homogenized Weyl algebra with variables x(1),..,x(n),D(1),..,
5166    D(n) and homogenization variable h, i.e. it holds x(i)*D(i)=D(i)*x(1)+h^2.
5167    If # contains on intvec v, we assign weight v[i] to the ith module component.*/
5168  if (n<1)
5169    {
5170      print*("Incorrect input");
5171      return();
5172    }
5173  if (n ==1)
5174    {
5175      ring @rr = 0,(D(1),x(1),h),dp;
5176    }
5177  else
5178    {
5179      ring @rr = 0,(D(1..n),x(1..n),h),dp;
5180    }
5181  setring @rr;
5182  int i=1;
5183  if (size(#)==0)
5184    {
5185      def @rrr = converseHomogenizedWeyl(i);
5186    }
5187  else
5188    {
5189      def @rrr=converseHomogenizedWeyl(i,#);
5190    }
[7fe9f8b]5191  return(@rrr);
5192}
5193
5194////////////////////////////////////////////////////////////////////////////////////
5195
[0e8a5a]5196static proc converseHomogenizedWeyl (int tilde,list #)
[7fe9f8b]5197{
5198  /*modified version of the procedure Weyl() from the library nctools.lib*/
[0e8a5a]5199  /*Creates a homogenized Weyl algebra structure on the basering. We assume
5200    n=nvars(basering) is odd. The first (n-1)/2 variables will be treated as the
5201    x(i), the next (n-1)/2 as the corresponding differentials D(i) and the last as
5202    the homogenization variable h, i.e. it holds x(i)*D(i)=D(i)*x(1)+h^2.
5203    If # contains on intvec v, we assign weight v[i] to the ith module component.*/
[7fe9f8b]5204  string rname=nameof(basering);
5205  if ( rname == "basering") // i.e. no ring has been set yet
[0e8a5a]5206    {
5207      "You have to call the procedure from the ring";
5208      return();
5209    }
[7fe9f8b]5210  int nv = nvars(basering);
5211  int N = (nv-1) div 2;
5212  if (((nv-1) % 2) != 0)
[0e8a5a]5213    {
5214      "Cannot create homogenized Weyl structure for an even number of generators";
5215      return();
5216    }
[7fe9f8b]5217  matrix @D[nv][nv];
5218  int i;
5219  for ( i=1; i<=N; i++ )
[0e8a5a]5220    {
5221      @D[i,N+i]=-h^2;
5222    }
[7fe9f8b]5223  def @R = nc_algebra(1,@D);
5224  setring @R;
5225  list RL=ringlist(@R);
5226  intvec v;
5227  /*we need this ordering for Groebner basis computations*/
[0e8a5a]5228  if (tilde==0)
5229    {
5230      for (i=1; i<=N; i++)
5231        {
5232          v[i]=-1;
5233          v[N+i]=1;
5234        }
5235    }
5236  else
5237    {
5238      for (i=1; i<=N; i++)
5239        {
5240          v[i]=1;
5241          v[N+i]=-1;
5242        }
5243    }
5244  v[nv]=0;
5245  /* we assign weights to module components*/
5246  if (size(#)!=0)
5247    {
5248      if (typeof(#[1])=="intvec")
5249        {
5250          intvec m=#[1];
5251          for (i=1; i<=size(m); i++)
5252            {
5253              v[size(v)+1]=m[i];//assigns weight m[i] to the ith module component
5254            }
5255          RL[3]=insert(RL[3],list("am",v));
5256        }
5257      else
5258        {
5259          RL[3]=insert(RL[3],list("a",v));
5260        }
5261    }
5262  else
5263    {
5264      RL[3]=insert(RL[3],list("a",v));
5265    }
5266  intvec w=(1:nv);
5267  if (size(#)>=2)
5268    {
5269      if (typeof(#[2])=="intvec")
5270        {
5271          intvec n=#[2];
5272          for (i=1; i<=size(n); i++)
5273            {
5274              w[size(w)+1]=n[i];
5275            }
5276          RL[3]=insert(RL[3],list("am",w));
5277        }
5278      else
5279        {
5280          RL[3]=insert(RL[3],list("a",w));
5281        }
5282    }
5283  else
5284    {
5285      RL[3]=insert(RL[3],list("a",w));
5286    }
5287  /*this ordering is needed for globalBFun and globalBFunOT*/
5288  list saveord=RL[3][3];
5289  RL[3][3]=RL[3][4];
5290  RL[3][4]=saveord;
5291  intvec notforh=(1:(size(RL[3][4][2])-1));
5292  RL[3][4][2]=notforh;
5293  RL[3][5]=list("dp",1);
5294  def @@R=ring(RL);
5295  return(@@R);
5296}
5297///////////////////////////////////////////////////////////////////////////////////
5298
5299static proc homogenizedWeyl (int tilde,list #)
5300{
5301  /*modified version of the procedure Weyl() from the library nctools.lib*/
5302  /*Creates a homogenized Weyl algebra structure on the basering. We assume
5303    n=nvars(basering) is odd. The first (n-1)/2 variables will be treated as the
5304    x(i), the next (n-1)/2 as the corresponding differentials D(i) and the last as
5305    the homogenization variable h, i.e. it holds x(i)*D(i)=D(i)*x(1)+h^2.
5306    If # contains on intvec v, we assign weight v[i] to the ith module component.*/
5307  string rname=nameof(basering);
5308  if ( rname == "basering") // i.e. no ring has been set yet
5309    {
5310      "You have to call the procedure from the ring";
5311      return();
5312    }
5313  int nv = nvars(basering);
5314  int N = (nv-1) div 2;
5315  if (((nv-1) % 2) != 0)
5316    {
5317      "Cannot create homogenized Weyl structure for an even number of generators";
5318      return();
5319    }
5320  matrix @D[nv][nv];
5321  int i;
5322  for ( i=1; i<=N; i++ )
5323    {
5324      @D[i,N+i]=h^2;
5325    }
5326  def @R = nc_algebra(1,@D);
5327  setring @R;
5328  list RL=ringlist(@R);
5329  intvec v;
5330  /*we need this ordering for Groebner basis computations*/
5331  if (tilde==0)
5332    {
5333      for (i=1; i<=N; i++)
5334        {
5335          v[i]=-1;
5336          v[N+i]=1;
5337        }
5338    }
5339  else
5340    {
5341      for (i=1; i<=N; i++)
5342        {
5343          v[i]=1;
5344          v[N+i]=-1;
5345        }
5346    }
[7fe9f8b]5347  v[nv]=0;
5348  /* we assign weights to module components*/
5349  if (size(#)!=0)
5350    {
[0e8a5a]5351      if (typeof(#[1])=="intvec")
5352        {
5353          intvec m=#[1];
5354          for (i=1; i<=size(m); i++)
5355            {
5356              v[size(v)+1]=m[i];//assigns weight m[i] to the ith module component
5357            }
5358          RL[3]=insert(RL[3],list("am",v));
5359        }
5360      else
5361        {
5362          RL[3]=insert(RL[3],list("a",v));
5363        }
[7fe9f8b]5364    }
5365  else
[0e8a5a]5366    {
5367      RL[3]=insert(RL[3],list("a",v));
5368    }
[7fe9f8b]5369  intvec w=(1:nv);
5370  if (size(#)>=2)
5371    {
[0e8a5a]5372      if (typeof(#[2])=="intvec")
5373        {
5374          intvec n=#[2];
5375          for (i=1; i<=size(n); i++)
5376            {
5377              w[size(w)+1]=n[i];
5378            }
5379          RL[3]=insert(RL[3],list("am",w));
5380        }
5381      else
5382        {
5383          RL[3]=insert(RL[3],list("a",w));
5384        }
[7fe9f8b]5385    }
5386  else
[0e8a5a]5387    {
5388      RL[3]=insert(RL[3],list("a",w));
5389    }
[7fe9f8b]5390  /*this ordering is needed for globalBFun and globalBFunOT*/
5391  list saveord=RL[3][3];
5392  RL[3][3]=RL[3][4];
5393  RL[3][4]=saveord;
5394  intvec notforh=(1:(size(RL[3][4][2])-1));
5395  RL[3][4][2]=notforh;
5396  RL[3][5]=list("dp",1);
5397  def @@R=ring(RL);
5398  return(@@R);
5399}
5400
5401////////////////////////////////////////////////////////////////////////////////////
5402
5403static proc nHomogenize (matrix M,list #)
5404{
5405  /* # may contain an intvec v, if no intvec is given, we assume that v=(0:ncols(M))
[0e8a5a]5406     We compute the h[v]-homogenization of the rows of M as in Definition 9.2 [OT]*/
[7fe9f8b]5407  int l; poly f; int s; int i; intvec vnm;int kmin; list findmax;
5408  int n=(nvars(basering)-1) div 2;
5409  list rempoly;
5410  list remk;
5411  list rem1;
5412  list rem2;
5413  list maxhexp;
5414  int hexp;
5415  intvec v=(0:ncols(M));
5416  if (size(#)!=0)
5417    {
[0e8a5a]5418      if (typeof(#[1])=="intvec")
5419        {
5420          v=#[1];
5421        }
[7fe9f8b]5422    }
5423  if (size(v)<ncols(M))
5424    {
[0e8a5a]5425      for (i=size(v)+1; i<=ncols(M); i++)
5426        {
5427          v[i]=0;
5428        }
[7fe9f8b]5429    }
[0e8a5a]5430  for (int k=1; k<=nrows(M); k++)
[7fe9f8b]5431    {
[0e8a5a]5432      for (l=1; l<=ncols (M); l++)
5433        {
5434          f=M[k,l];
5435          s=size(f);
5436          for (i=1; i<=s; i++)
5437            {
5438              vnm=leadexp(f);
5439              kmin=sum(vnm)+v[l];
5440              rem1[size(rem1)+1]=lead(f);
5441              rem2[size(rem2)+1]=kmin;
5442              findmax=insert(findmax,kmin);
5443              f=f-lead(f);
5444            }
5445          rempoly[l]=rem1;
5446          remk[l]=rem2;
5447          rem1=list();
5448          rem2=list();
5449        }
5450      if (size(findmax)!=0)
5451        {
5452          kmin=Max(findmax);
5453        }
5454      else
[7fe9f8b]5455        {
[0e8a5a]5456          kmin=0;
5457        }
5458      for (l=1; l<=ncols(M); l++)
5459        {
5460          if (M[k,l]!=0)
5461            {
5462              M[k,l]=0;
5463              for (i=1; i<=size(rempoly[l]);i++)
5464                {
5465                  hexp=kmin-remk[l][i];
5466                  maxhexp[size(maxhexp)+1]=hexp;
5467                  M[k,l]=M[k,l]+h^hexp*rempoly[l][i];
5468                }
5469            }
[7fe9f8b]5470        }
[0e8a5a]5471      rempoly=list();
5472      remk=list();
5473      findmax=list();
[7fe9f8b]5474    }
5475  if (size(maxhexp)!=0)
[0e8a5a]5476    {
5477      maxhexp=Max(maxhexp);
5478      hexp=maxhexp[1];
5479    }
[7fe9f8b]5480  else
[0e8a5a]5481    {
5482      hexp=0;
5483    }
[7fe9f8b]5484  if (size(#)>1)
[0e8a5a]5485    {
5486      list forreturn=M,hexp;
[7fe9f8b]5487
[0e8a5a]5488      return(forreturn);
5489    }
[7fe9f8b]5490  return(M);
5491}
5492
5493////////////////////////////////////////////////////////////////////////////////////
5494
5495static proc max(int i,int j)
5496{
5497  if(i>j){return(i);}
5498  return(j);
5499}
5500
5501////////////////////////////////////////////////////////////////////////////////////
5502
5503static proc nDeg (matrix M,intvec m)
[0e8a5a]5504{/*we compute an intvec n such that n[i]=max(deg(M[i,j])+m[j]|M[i,j]!=0) (where deg
5505   stands for the total degree) if (M[i,j]!=0 for some j) and n[i]=0 else*/
[7fe9f8b]5506  int i; int j;
5507  intvec n;
5508  list L;
5509  for (i=1; i<=nrows(M); i++)
5510    {
[0e8a5a]5511      L=list();
5512      for (j=1; j<=ncols(M); j++)
5513        {
5514          if (M[i,j]!=0)
5515            {
5516              L=insert(L,deg(M[i,j])+m[j]);
5517            }
5518        }
5519      if (size(L)==0)
5520        {
5521          n[i]=0;
5522        }
5523      else
5524        {
5525          n[i]=Max(L);
5526        }
[7fe9f8b]5527    }
5528  return(n);
5529}
5530
5531////////////////////////////////////////////////////////////////////////////////////
5532
[76d26c]5533static proc minIntRootD(list L,list #)
5534"USAGE:minIntRootD(L [,M]); L list, M optinonal list
[7fe9f8b]5535ASSUME:L a list of univariate polynomials with rational coefficients @*
[0e8a5a]5536       the variable of the polynomial is s if size(#)==0 (needed for proc
[7fe9f8b]5537       MVComplex) and t else (needed for globalBFun)
[0e8a5a]5538RETURN:-if size(#)==0: int i, where i is an integer root of one of the polynomials
[7fe9f8b]5539        and it is minimal with respect to that property@*
[0e8a5a]5540       -if size(#)!=0: list L=(i,j), where i is as above and j is an integer root
5541        of one of the polynomials and is maximal with respect to that property (if
[7fe9f8b]5542        an integer root exists) or L=list() else
5543"
5544{
5545  def B=basering;
5546  if (size(#)==0)
[0e8a5a]5547    {
5548      ring rnew=0,s,dp;
5549    }
[7fe9f8b]5550  else
[0e8a5a]5551    {
5552      ring rnew=0,t,dp;
5553    }
[7fe9f8b]5554  list L=imap(B,L);
5555
5556  int i;
5557  int j;
5558  number isint;
5559  list possmin;
5560  ideal allfac;
5561  list allfacs;
5562  for (i=1; i<=size(L); i++)
5563    {
[0e8a5a]5564      allfac=factorize(L[i],1);
5565      for (j=1; j<=ncols(allfac); j++)
5566        {
5567          allfacs[j]=allfac[j];
5568        }
5569      for (j=1; j<=size(allfacs); j++)
[7fe9f8b]5570        {
[0e8a5a]5571          if (deg(allfacs[j])==1)
5572            {
5573              isint=number(subst(allfacs[j],var(1),0)/leadcoef(allfacs[j]));
5574              if (isint-int(isint)==0)
5575                {
5576                  possmin[size(possmin)+1]=int(isint);
5577                }
5578            }
[7fe9f8b]5579        }
[0e8a5a]5580      allfacs=list();
[7fe9f8b]5581    }
5582  int zerolist;
5583  if (size(possmin)!=0)
[0e8a5a]5584    {
5585      int miniroot=(-1)*Max(possmin);
5586      int maxiroot=(-1)*Min(possmin);
5587    }
[7fe9f8b]5588  else
[0e8a5a]5589    {
5590      zerolist=1;
5591    }
[7fe9f8b]5592  setring B;
5593  if (size(#)==0)
[0e8a5a]5594    {
5595      return(miniroot);
5596    }
5597  else
5598    {
5599      if (zerolist==0)
5600        {
5601          return(list(miniroot,maxiroot));
5602        }
5603      else
5604        {
5605          return(list());
5606        }
5607    }
5608}
5609
5610////////////////////////////////////////////////////////////////////////////////////
5611
5612proc converseWeyl(list #)
5613{
5614  string rname=nameof(basering);
5615  int @chr = 0;
5616  int nv = nvars(basering);
5617  int N = nv div 2;
5618  matrix @D[nv][nv];
5619  int i;
5620  for ( i=1; i<=N; i++)
[7fe9f8b]5621  {
[0e8a5a]5622      @D[i,N+i]=-1;
[7fe9f8b]5623  }
[0e8a5a]5624  def @R = nc_algebra(1,@D);
5625  return(@R);
5626}
5627
5628////////////////////////////////////////////////////////////////////////////////////
5629
5630proc makeConverseWeyl(int n, list #)
5631{
5632  if (n==1)
5633    {
5634      ring @rr = 0,(D(1),x(1)),dp;
5635    }
[7fe9f8b]5636  else
5637  {
[0e8a5a]5638    ring @rr = 0,(D(1..n),x(1..n)),dp;
5639  }
5640  setring @rr;
5641  def @rrr = converseWeyl();
5642  return(@rrr);
5643}
5644
5645////////////////////////////////////////////////////////////////////////////////////
5646
5647proc makeOmega(int n)
5648{
5649  def R=basering;
5650  int i;
5651  int j,k,l;
5652  list omega;
5653  omega[1]=list(list(list()));
5654  omega[2]=list();
5655  for (i=1; i<=n; i++)
[7fe9f8b]5656    {
[0e8a5a]5657      omega[2][i]=list(i);
[7fe9f8b]5658    }
[0e8a5a]5659  for (i=2; i<=n; i++)
[7fe9f8b]5660    {
[0e8a5a]5661      omega[i+1]=list();
5662      for (j=1; j<=size(omega[i]); j++)
5663        {
5664          if (omega[i][j][size(omega[i][j])]<n)
5665            {
5666              for (k=omega[i][j][size(omega[i][j])]+1; k<=n; k++)
5667                {
5668                  omega[i+1][size(omega[i+1])+1]=omega[i][j];
5669                  omega[i+1][size(omega[i+1])][size( omega[i+1][size(omega[i+1])])+1]=k;
5670                }
5671            }
5672        }
[7fe9f8b]5673    }
[0e8a5a]5674  list omegamaps;
5675  matrix om;
5676  list lms;
5677  omegamaps[1]=matrix(0,n,1);
5678  for (i=1; i<=n; i++)
5679    {
5680      omegamaps[1][i,1]=var(n+i);
5681    }
5682  for (i=2; i<=n; i++)
5683    {
5684      om=matrix(0,size(omega[i+1]),size(omega[i]));
5685      for (k=1; k<=size(omega[i]); k++)
5686        {
5687          for (l=1; l<=size(omega[i+1]); l++)
5688            {
5689              lms=LMSubset(omega[i][k],omega[i+1][l],1);
5690              om[l,k]=lms[2]*var(n+lms[1]);
5691            }
5692        }
5693      omegamaps[i]=om;
5694    }
5695  omegamaps[n+1]=matrix(0,1,1);
5696  list allomega;
5697  for (i=1; i<=n+1; i++)
5698    {
5699      allomega[2*i]=omega[n+2-i];
5700      allomega[2*i-1]=omegamaps[n+2-i];
5701    }
5702  return(allomega);
5703}
5704
5705////////////////////////////////////////////////////////////////////////////////////
5706
5707static proc makeDoubleComplex(list L, list M, list Q, list G)
5708{
5709  list doublecomplex;
5710  int i,j,k,l;
5711  int s1;
5712  int s2;
5713  int c;
5714  int d;
5715  list gens=list();
5716  for (i=1; i<=size(L) div 2; i++)
5717    {
5718      doublecomplex[i]=list();
5719      for (j=1; j<=size(M) div 2; j++)
5720        {
5721          doublecomplex[i][j]=list();
5722          doublecomplex[i][j]=list(M[2*j]+list(L[2*i-1]));
5723          gens=list();
5724          doublecomplex[i][j][6]=G[i];
5725          if (size(Q[i])!=0)
5726            {
5727              doublecomplex[i][j][4]=tensor(unitmat(size(M[2*j])),Q[i]);
5728              for (c=1; c<=size(M[2*j]); c++)
5729                {
5730                  for (d=1; d<=ncols(Q[i]); d++)
5731                    {
5732                      gens[size(gens)+1]=list(M[2*j][c],d);
5733                    }
5734                }
5735              doublecomplex[i][j][5]=gens;
5736            }
5737          else
5738            {
5739              doublecomplex[i][j][4]=list();
5740              doublecomplex[i][j][5]=list();
5741            }
5742          if (size(Q[i])!=0)
5743            {
5744              if (Q[i]==matrix(0,nrows(Q[i]),ncols(Q[i])))
5745                {
5746                  doublecomplex[i][j][4]=list();
5747                }
5748            }
5749          if (j!=1)
5750            {
5751              s1=(size(doublecomplex[i][j-1][1])-1)*doublecomplex[i][j-1][1][size(doublecomplex[i][j-1][1])];
5752              s2=(size(doublecomplex[i][j][1])-1)*doublecomplex[i][j][1][size(doublecomplex[i][j][1])];
5753              if (s1==0 or s2==0)
5754                {
5755                  doublecomplex[i][j-1][3]=list();
5756                }
5757              else
5758                {
5759                  doublecomplex[i][j-1][3]=tensor(M[2*j-1],unitmat(L[2*i-1]));
5760                }
5761
5762            }
5763          if (j==size(M) div 2)
5764            {
5765              doublecomplex[i][j][3]=list();
5766            }
5767          if (i!=1)
5768            {
5769              s1=(size(doublecomplex[i-1][j][1])-1)*doublecomplex[i-1][j][1][size(doublecomplex[i-1][j][1])];
5770              s2=(size(doublecomplex[i][j][1])-1)*doublecomplex[i][j][1][size(doublecomplex[i][j][1])];
5771              if (s1==0 or s2==0)
5772                {
5773                  doublecomplex[i-1][j][2]=list();
5774                }
5775              else
5776                {
5777                  doublecomplex[i-1][j][2]=tensor(unitmat(size(M[2*j])),L[2*(i-1)]);
5778                }
5779            }
5780          if (i==size(L) div 2)
5781            {
5782              doublecomplex[i][j][2]=list();
5783            }
5784        }
5785    }
5786  return(doublecomplex);
[7fe9f8b]5787}
5788
[0e8a5a]5789////////////////////////////////////////////////////////////////////////////////////
5790
5791static proc transferDiffforms(matrix m, list L)
5792{
5793  int i;
5794  list transfered;
5795  if (size(L[4])==0)
5796    {
5797      return(list());
5798    }
5799  if (size(L[5])==0)
5800    {
5801      return(list());
5802    }
5803  m=m*L[4];
5804  list transferedm=list();
5805  int si=L[5][size(L[5])][2];//Anzahl der direkten Summanden in \oplus R_F_I
5806  matrix fortrans=matrix(0,1,si);
5807  list omegagen=list();
5808  list save=list();
5809  int t;
5810  int c;
5811  int j;
5812  list converteddiff;
5813  vector w;
5814  poly p=1;
5815  for (i=1; i<=ncols(m); i++)
5816    {
5817      if (m[1,i]!=0)
5818        {
5819          if (size(omegagen)==0)
5820            {
5821              omegagen=L[5][i][1];
5822              fortrans[1,L[5][i][2]]= fortrans[1,L[5][i][2]]+m[1,i];
5823            }
5824          else
5825            {
5826              t=0;
5827              for (j=1; j<=size(omegagen);j++)
5828                {
5829                  if (size(omegagen[j])!=0)
5830                    {
5831                      if (omegagen[j]!=L[5][i][1][j])
5832                        {
5833                          t=1;
5834                        }
5835                    }
5836                }
5837              if (t==0)
5838                {
5839                  fortrans[1,L[5][i][2]]= fortrans[1,L[5][i][2]]+m[1,i];
5840                }
5841              else
5842                {
5843                  converteddiff=list();
5844                  for (j=1; j<=ncols(fortrans); j++)
5845                    {
5846                      if (fortrans[1,j]!=0)
5847                        {
5848                          w=[p,L[6][j]];
5849                          converteddiff[j]=dmodActionRat(fortrans[1,j],w);
5850                        }
5851                      else
5852                        {
5853                          converteddiff[j]=0;
5854                        }
5855
5856                    }
5857                  save[size(save)+1]=list(converteddiff,omegagen);
5858                  omegagen=L[5][i][1];
5859                  fortrans=matrix(0,1,si);
5860                  fortrans[1,L[5][i][2]]= fortrans[1,L[5][i][2]]+m[1,i];
5861                }
5862            }
5863        }
5864    }
5865  if (fortrans==matrix(0,1,si))
5866    {
5867      return(list());
5868    }
5869  converteddiff=list();
5870  for (j=1; j<=ncols(fortrans); j++)
5871    {
5872      if (fortrans[1,j]!=0)
5873        {
5874          w=[p,L[6][j]];
5875          converteddiff[j]=dmodActionRat(fortrans[1,j],w);
5876        }
5877      else
5878        {
5879          converteddiff[j]=0;
5880        }
5881    }
5882  save[size(save)+1]=list(converteddiff,omegagen);
5883  return(save);
5884}
[7fe9f8b]5885
5886////////////////////////////////////////////////////////////////////////////////////
5887////////////////////////////////////////////////////////////////////////////////////
5888////////////////////////////////////////////////////////////////////////////////////
5889/*
5890////////////////////////////////////////////////////////////////////////////////////
5891FURTHER EXAMPLES FOR TESTING THE PROCEDURES
5892////////////////////////////////////////////////////////////////////////////////////
5893LIB "derham.lib";
5894
5895//----------------------------------------
5896//EXAMPLE 1
5897//----------------------------------------
5898ring r=0,(x,y),dp;
5899poly f=y2-x3-2x+3;
5900list L=deRhamCohomology(f);
5901L;
5902kill r;
5903
5904//----------------------------------------
5905//EXAMPLE 2
5906//----------------------------------------
5907ring r=0,(x,y),dp;
5908poly f=y2-x3-x;
5909list L=deRhamCohomology(f);
5910L;
5911kill r;
5912
5913//----------------------------------------
5914//EXAMPLE 3
5915//----------------------------------------
5916ring r=0,(x,y),dp;
5917list C=list(x2-1,(x+1)*y,y*(y2+2x+1));
5918list L=deRhamCohomology(C);
5919L;
5920kill r;
5921
5922//----------------------------------------
5923//EXAMPLE 4
5924//----------------------------------------
5925ring r=0,(x,y,z),dp;
5926list C=list(x*(x-1),y,z*(z-1),z*(x-1));
5927list L=deRhamCohomology(C);
5928L;
5929kill r;
5930
5931//----------------------------------------
5932//EXAMPLE 5
5933//----------------------------------------
5934ring r=0,(x,y,z),dp;
5935list C=list(x*y,y*z);
5936list L=deRhamCohomology(C,"Vdres");
5937L;
5938kill r;
5939
5940//----------------------------------------
5941//EXAMPLE 6
5942//----------------------------------------
5943ring r=0,(x,y,z,u),dp;
5944list C=list(x,y,z,u);
5945list L=deRhamCohomology(C);
5946L;
5947kill r;
5948
5949//----------------------------------------
5950//EXAMPLE 7
5951//----------------------------------------
5952ring r=0,(x,y,z),dp;
5953poly f=x3+y3+z3;
5954list L=deRhamCohomology(f);
5955L;
5956kill r;
5957
5958//----------------------------------------
5959//EXAMPLE 8
5960//----------------------------------------
5961ring r=0,(x,y,z),dp;
5962poly f=x2+y2+z2;
5963list L=deRhamCohomology(f,"Vdres");
5964L;
5965kill r;
5966
5967//----------------------------------------
5968//EXAMPLE 9
5969//----------------------------------------
5970ring r=0,(x,y,z,u),dp;
5971list C=list(x2+y2+z2,u);
5972list L=deRhamCohomology(C);
5973L;
5974kill r;
5975
5976
5977//----------------------------------------
5978//EXAMPLE 10
5979//----------------------------------------
5980ring r=0,(x,y,z),dp;
5981list C=list((x*(y-1),y2-1));
5982list L=deRhamCohomology(C);
5983L;
5984kill r;
5985
5986
5987*/
Note: See TracBrowser for help on using the repository browser.