Changeset 0e8a5a in git


Ignore:
Timestamp:
Sep 26, 2013, 2:16:29 PM (10 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
12d61707b29797f7bf3e42e5b9f0a01e937da42b
Parents:
89b2d289ea9d528c84c54a5e203ee55bb49e600c
Message:
new version of derham.lib
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/derham.lib

    r89b2d28 r0e8a5a  
    1 /////////////////////////////////////////////////////////////////////////////////
    2 version="version derham.lib 4.0.0.0 Jun_2013 "; // $Id$
     1///////////////////////////////////////////////////////////////////////////////
     2version="";
    33category="Noncommutative";
    44info="
    55LIBRARY:  derham.lib      Computation of deRham cohomology
    6 
    76AUTHORS:  Cornelia Rottner, rottner@mathematik.uni-kl.de
    8 
    9 OVERVIEW:
    10 A library for computing the de Rham cohomology of complements of complex affine
     7OVERVIEW:
     8PROCEDURES:
     9
     10";
     11
     12
     13LIB "nctools.lib";
     14LIB "matrix.lib";
     15LIB "qhmoduli.lib";
     16LIB "general.lib";
     17LIB "dmod.lib";
     18LIB "bfun.lib";
     19LIB "dmodapp.lib";
     20LIB "poly.lib";
     21
     22/////////////////////////////////////////////////////////////////////////////
     23
     24static proc divdr(matrix m, matrix n)
     25{
     26  m=transpose(m);
     27  n=transpose(n);
     28  matrix con=concat(m,n);
     29  matrix s=syz(con);
     30  s=submat(s,1..ncols(m),1..ncols(s));
     31  s=transpose(compress(s));
     32  return(s);
     33}
     34
     35/////////////////////////////////////////////////////////////////////////////
     36
     37
     38static proc matrixlift(matrix M, matrix N)
     39{
     40  // option(noreturnSB);
     41  matrix l=transpose(lift(transpose(M),transpose(N)));
     42  return(l);
     43}
     44
     45///////////////////////////////////////////////////////////////////////////////
     46
     47proc shortexactpieces(list #)
     48{
     49  matrix  Bnew= divdr(#[2],#[3]);
     50  matrix Bold=Bnew;
     51  matrix Z=divdr(Bnew,#[1]);
     52  list bzh; list zcb;
     53  bzh=list(list(),list(),Z,unitmat(ncols(Z)),Z);
     54  zcb=(Z, Bnew, #[1], unitmat(ncols(#[1])), Bnew);
     55  list sep;
     56  sep[1]=(list(bzh,zcb));
     57  int i;
     58  list out;
     59  for (i=3; i<=size(#)-2; i=i+2)
     60    {
     61      out=bzhzcb(Bold, #[i-1] , #[i], #[i+1], #[i+2]);
     62      sep[size(sep)+1]=out[1];
     63      Bold=out[2];
     64    }
     65  bzh=(divdr(#[size(#)-2], #[size(#)-1]),#[size(#)-2], #[size(#)-1],unitmat(ncols(#[size(#)-1])),transpose(concat(transpose(#[size(#)-2]),transpose(#[size(#)-1]))));
     66  zcb=(#[size(#)-1], unitmat(ncols(#[size(#)-1])), #[size(#)-1],list(),list());
     67  sep[size(sep)+1]=list(bzh,zcb);
     68  return(sep);
     69}
     70
     71////////////////////////////////////////////////////////////////////////////////////////
     72
     73static proc bzhzcb (matrix Bold, matrix f0, matrix C1, matrix f1,matrix C2)
     74{
     75  matrix Bnew=divdr(f1,C2);
     76  matrix Z= divdr(Bnew,C1);
     77  matrix lift1= matrixlift(Bnew,f0);
     78  list bzh=(Bold, lift1, Z, unitmat(ncols(Z)), transpose(concat(transpose(lift1),transpose(Z))));
     79  list zcb=(Z, Bnew, C1, unitmat(ncols(C1)),Bnew);
     80  list out=(list(bzh, zcb), Bnew);
     81  return(out);
     82}
     83
     84//////////////////////////////////////////////////////////////////////////////////////
     85
     86proc VdstrictGB (matrix M, int d ,list #);
     87"USAGE:VdstrictGB(M,d[,v]); M a matrix, d an integer, v an optional intvec(shift vector)
     88RETURN:matrix M; the rows of M foem a Vd-strict Groebner basis for imM
     89ASSUME:1<=d<=nvars(basering)/2; size(v)=ncols(M)
     90"
     91{
     92  if (M==matrix(0,nrows(M),ncols(M)))
     93    {
     94      return (matrix(0,1,ncols(M)));
     95    }
     96  def W =basering;
     97  int ncM=ncols(M);
     98  list Data=ringlist(W);
     99  Data[2]=list("nhv")+Data[2];
     100  Data[3][3]=Data[3][1];
     101  Data[3][1]=Data[3][2];
     102  Data[3][2]=list("dp",intvec(1));
     103  matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2]));
     104  Data[5]=re;
     105  int k; int l;
     106  Data[6]=transpose(concat(matrix(0,1,1),transpose(concat(matrix(0,1,1),Data[6]))));
     107  def Whom=ring(Data);
     108  setring Whom;
     109  matrix Mnew=imap(W,M);
     110  intvec v;
     111  if (size(#)!=0)
     112    {
     113      v=#[1];
     114    }
     115  if (size(v) < ncM)
     116    {
     117      v=v,0:(ncM-size(v));
     118    }
     119  Mnew=homogenize(Mnew, d, v);
     120  Mnew=transpose(Mnew);
     121  Mnew=std(Mnew);
     122  Mnew=subst(Mnew,nhv,1);
     123  Mnew=transpose(Mnew);
     124  setring W;
     125  M=imap(Whom,Mnew);
     126  return(M);
     127}
     128
     129////////////////////////////////////////////////////////////////////////////////////
     130
     131static proc Vdnormalform(matrix F, matrix M, int d, intvec v)
     132{
     133  def W =basering;
     134  int c=ncols(M);
     135  F=submat(F,intvec(1..nrows(F)),intvec(1..c));
     136  list Data=ringlist(W);
     137  Data[2]=list("nhv")+Data[2];
     138  Data[3][3]=Data[3][1];
     139  Data[3][1]=Data[3][2];
     140  Data[3][2]=list("dp",intvec(1));
     141  matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2]));
     142  Data[5]=re;
     143  int k;
     144  int l;
     145  matrix rep[size(Data[2])][size(Data[2])];
     146  for (l=size(Data[2])-1;l>=1; l--)
     147    {
     148      for (k=l-1; k>=1;k--)
     149        {
     150          rep[k+1,l+1]=Data[6][k,l];
     151        }
     152    }
     153  Data[6]=rep;
     154  def Whom=ring(Data);
     155  setring Whom;
     156  matrix Mnew=imap(W,M);
     157  Mnew=(homogenize(Mnew, d, v));//doppelte Berechung unnötig->muss noch geändert werden!!!
     158  matrix Fnew=imap(W,F);
     159  matrix Fb;
     160  for (l=1; l<=nrows(Fnew); l++)
     161    {
     162      Fb=homogenize(submat(Fnew,l,intvec(1..ncols(Fnew))),d,v);
     163      Fb=transpose(reduce(transpose(Fb),std(transpose(Mnew))));// doppelte Berechnung unnötig, unterdrückt aber Fehler meldung
     164      for (k=1; k<=ncols(Fnew);k++)
     165        {
     166          Fnew[l,k]=Fb[1,k];
     167        }
     168    }
     169  Fnew=subst(Fnew,nhv,1);
     170  setring W;
     171  F=imap(Whom,Fnew);
     172  return(F);
     173}
     174
     175
     176///////////////////////////////////////////////////////////////////////////////
     177
     178static proc homogenize (matrix M, int d, intvec v)
     179{
     180  int l; poly f; int s; int i; intvec vnm;int kmin; list findmin;
     181  int n=(nvars(basering)-1) div 2;
     182  list rempoly;
     183  list remk;
     184  list rem1;
     185  list rem2;
     186  for (int k=1; k<=nrows(M); k++) //man könnte auch paralell immer weiter homogenisieren, d.h. immer ein enues Minimum finden und das dann machen
     187    {
     188      for (l=1; l<=ncols (M); l++)
     189        {
     190          f=M[k,l];
     191          s=size(f);
     192          for (i=1; i<=s; i++)
     193            {
     194              vnm=leadexp(f);
     195              vnm=vnm[n+2..n+d+1]-vnm[2..d+1];
     196              kmin=sum(vnm)+v[l];
     197              rem1[size(rem1)+1]=lead(f);
     198              rem2[size(rem2)+1]=kmin;
     199              findmin=insert(findmin,kmin);
     200              f=f-lead(f);
     201            }
     202          rempoly[l]=rem1;
     203          remk[l]=rem2;
     204          rem1=list();
     205          rem2=list();
     206        }
     207      if (size(findmin)!=0)
     208        {
     209          kmin=Min(findmin);
     210        }
     211      for (l=1; l<=ncols(M); l++)
     212        {
     213          if (M[k,l]!=0)
     214            {
     215              M[k,l]=0;
     216              for (i=1; i<=size(rempoly[l]);i++)
     217                {
     218                  M[k,l]=M[k,l]+nhv^(remk[l][i]-kmin)*rempoly[l][i];
     219                }
     220            }
     221        }
     222      rempoly=list();
     223      remk=list();
     224      findmin=list();
     225    }
     226  return(M);
     227}
     228
     229
     230//////////////////////////////////////////////////////////////////////////////////////
     231
     232static proc soldr (matrix M, matrix N)
     233{
     234  int n=nrows(M);
     235  int q=ncols(M);
     236  matrix S=concat(transpose(M),transpose(N));
     237  def W=basering;
     238  list Data=ringlist(W);
     239  list Save=Data[3];
     240  Data[3]=list(list("c",0),list("dp",intvec(1..nvars(W))));
     241  def Wmod=ring(Data);
     242  setring Wmod;
     243  matrix Smod=imap(W,S);
     244  matrix E[q][1];
     245  matrix Smod2;
     246  matrix Smodnew;
     247  option(returnSB);
     248  int i; int j;
     249  for (i=1;i<=q;i++)
     250    {
     251      E[i,1]=1;
     252      Smod2=concat(E,Smod);
     253      print (Smod2);
     254      Smod2=syz(Smod2);
     255      E[i,1]=0;
     256      for (j=1;j<=ncols(Smod2);j++)
     257        {
     258          if (Smod2[1,j]==1)
     259            {
     260              Smodnew=concat(Smodnew,(-1)*(submat(Smod2,intvec(2..n+1),j)));
     261              break;
     262            }
     263        }
     264    }
     265  Smodnew=transpose(submat(Smodnew,intvec(1..n),intvec(2..q+1)));
     266  setring W;
     267  matrix  Snew=imap(Wmod,Smodnew);
     268  return (Snew);
     269}
     270
     271
     272/////////////////////////////////////////////////////////////////////////////
     273
     274proc toVdstrictsequence (list C,int n, intvec v)
     275
     276{
     277  matrix J_C=VdstrictGB(C[5],n,list(v));
     278  matrix J_A=C[1];
     279  matrix f_CB=C[4];
     280  matrix f_ACB=transpose(concat(transpose(C[2]),transpose(f_CB)));
     281  matrix J_AC=divdr(f_ACB,C[3]);
     282  matrix P=matrixlift(J_AC * prodr(ncols(C[1]),ncols(C[5])) ,J_C);
     283  list storePi;
     284  matrix Pi[1][ncols(J_AC)];
     285  int i;int j;
     286  for (i=1; i<=nrows(J_C); i++)
     287    {
     288      for (j=1; j<=nrows(J_AC);j++)
     289        {
     290          Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC)));
     291        }
     292      storePi[i]=Pi;
     293      Pi=0;
     294    }
     295  intvec m_a;
     296  list findMin;
     297  int comMin;
     298  for (i=1; i<=ncols(J_A); i++)
     299    {
     300      for (j=1; j<=size(storePi);j++)
     301        {
     302          if (storePi[j][1,i]!=0)
     303            {
     304              comMin=Vddeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v)-Vddeg(storePi[j][1,i],n,intvec(0));
     305              findMin[size(findMin)+1]=comMin;
     306            }
     307        }
     308      if (size(findMin)!=0)
     309        {
     310          m_a[i]=Min(findMin);
     311          findMin=list();
     312        }
     313      else
     314        {
     315          m_a[i]=0;
     316        }
     317    }
     318  matrix zero[ncols(J_A)][ncols(J_C)];
     319
     320  matrix g_AB=concat(unitmat(ncols(J_A)),zero);
     321  matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C)))));
     322  intvec m_b=m_a,v;
     323
     324  J_A=VdstrictGB(J_A,n,m_a);
     325  J_AC=transpose(storePi[1]);
     326  for (i=2; i<= size(storePi); i++)
     327    {
     328      J_AC=concat(J_AC, transpose(storePi[i]));
     329    }
     330  J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC));
     331
     332  list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a),list(m_b),list(v));
     333  return (Vdstrict);
     334}
     335
     336
     337/////////////////////////////////////////////////////////////////////////
     338
     339static proc prodr (int k, int l)
     340{
     341  if (k==0)
     342    {
     343      matrix P=unitmat(l);
     344      return (P);
     345    }
     346  matrix O[l][k];
     347  matrix P=transpose(concat(O,unitmat(l)));
     348  return (P);
     349}
     350
     351/////////////////////////////////////////////////////////////////////////
     352
     353proc Vddeg(matrix M, int d, intvec v, list #)//Aternative: in WHom Leadmonom ausrechnen!!
     354  "USAGE: Vddeg(M,d,v); M 1xr-matrix, d int, v intvec of size r
     355RETURN:int; the Vd-degree of M
     356"
     357{
     358  int i;int j;
     359  int n=nvars(basering) div 2;
     360  intvec  e;
     361  int etoint;
     362  list findmax;
     363  int c=ncols(M);
     364  poly l;
     365  list positionpoly;
     366  list positionVd;
     367  for (i=1; i<=c; i++)
     368    {
     369      positionpoly[i]=list();
     370      positionVd[i]=list();
     371      while (M[1,i]!=0)
     372        {
     373          l=lead(M[1,i]);
     374          positionpoly[i][size(positionpoly[i])+1]=l;
     375          e=leadexp(l);
     376          e=-e[1..d]+e[n+1..n+d];
     377          e=sum(e)+v[i];
     378          etoint=e[1];
     379          positionVd[i][size(positionVd[i])+1]=etoint;
     380          findmax[size(findmax)+1]=etoint;
     381          M[1,i]=M[1,i]-l;
     382        }
     383    }
     384  if (size(findmax)!=0)
     385    {
     386      int maxVd=Max(findmax);
     387      if (size(#)==0)
     388        {
     389          return (maxVd);
     390        }
     391    }
     392  else // M is 0-modul
     393    {
     394      return(int(0));
     395    }
     396  l=0;
     397  for (i=c; i>=1; i--)
     398    {
     399      for (j=1; j<=size(positionVd[i]); j++)
     400        {
     401          if (positionVd[i][j]==maxVd)
     402            {
     403              l=l+positionpoly[i][j];
     404            }
     405        }
     406      if (l!=0)
     407        {
     408          return (list(l,i));
     409        }
     410    }
     411
     412}
     413
     414
     415///////////////////////////////////////////////////////////////////////////////
     416
     417proc toVdstrictsequences (list L,int d, intvec v)
     418  "USAGE: toVdstrictsequences(L,d,v); L list, d int, v intvec, L contains two lists of short exact sequences(D,f_DA,A,f_AF,F) and (A,f_AB,B,f_BC,C), v is a shift vector on the range of C
     419RETURN: list of two lists; each lists contains Vd-strict exact sequences with corresponding shift vectors
     420"
     421{
     422  matrix J_F=L[1][5];
     423  matrix J_D=L[1][1];
     424  matrix f_FA=L[1][4];
     425  matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA)));
     426  matrix J_DF=divdr(f_DFA,L[1][3]);
     427  matrix J_C=L[2][5];
     428  matrix f_CB=L[2][4];
     429  matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB)));
     430  matrix J_DFC=divdr(f_DFCB,L[2][3]);
     431  matrix P=matrixlift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C);
     432  list storePi;
     433  matrix Pi[1][ncols(J_DFC)];
     434  int i; int j;
     435  for (i=1; i<=nrows(J_C); i++)
     436    {
     437
     438      for (j=1; j<=nrows(J_DFC);j++)
     439        {
     440          Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC)));
     441        }
     442      storePi[i]=Pi;
     443      Pi=0;
     444    }
     445  intvec m_a;
     446  list findMin;
     447  list noMin;
     448  int comMin;
     449  for (i=1; i<=ncols(J_DF); i++)
     450    {
     451      for (j=1; j<=size(storePi);j++)
     452        {
     453          if (storePi[j][1,i]!=0)
     454            {
     455              comMin=Vddeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v)-Vddeg(storePi[j][1,i],d,intvec(0));
     456              findMin[size(findMin)+1]=comMin;
     457            }
     458        }
     459      if (size(findMin)!=0)
     460        {
     461          m_a[i]=Min(findMin);
     462          findMin=list();
     463          noMin[i]=0;
     464        }
     465      else
     466        {
     467          noMin[i]=1;
     468        }
     469    }
     470  if (size(m_a) < ncols(J_DF))
     471    {
     472      m_a[ncols(J_DF)]=0;
     473    }
     474  intvec m_f=m_a[ncols(J_D)+1..size(m_a)];
     475  J_F=VdstrictGB(J_F,d,m_f);
     476  P=matrixlift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F);// selbe Prinzip wie oben--> evtl auslagern
     477  list storePinew;
     478  matrix Pidf[1][ncols(J_DF)];
     479  for (i=1; i<=nrows(J_F); i++)
     480    {
     481      for (j=1; j<=nrows(J_DF);j++)
     482        {
     483          Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF)));
     484        }
     485      storePinew[i]=Pidf;
     486      Pidf=0;
     487    }
     488  intvec m_d;
     489  for (i=1; i<=ncols(J_D); i++)
     490    {
     491      for (j=1; j<=size(storePinew);j++)
     492        {
     493          if (storePinew[j][1,i]!=0)
     494            {
     495              comMin=Vddeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f)-Vddeg(storePinew[j][1,i],d,intvec(0));
     496              findMin[size(findMin)+1]=comMin;
     497            }
     498        }
     499      if (size(findMin)!=0)
     500        {
     501          if (noMin[i]==0)
     502            {
     503              m_d[i]=Min(insert(findMin,m_a[i]));
     504              m_a[i]=m_d[i];
     505            }
     506          else
     507            {
     508              m_d[i]=Min(findMin);
     509              m_a[i]=m_d[i];
     510            }
     511        }
     512      else
     513        {
     514          m_d[i]=m_a[i];
     515        }
     516      findMin=list();
     517    }
     518  J_D=VdstrictGB(J_D,d,m_d);
     519  J_DF=transpose(storePinew[1]);
     520  for (i=2; i<=nrows(J_F); i++)
     521    {
     522      J_DF=concat(J_DF,transpose(storePinew[i]));
     523    }
     524  J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF));
     525  J_DFC=transpose(storePi[1]);
     526  for (i=2; i<=nrows(J_C); i++)
     527    {
     528      J_DFC=concat(J_DFC,transpose(storePi[i]));
     529    }
     530  J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC));
     531  intvec m_b=m_a,v;
     532  matrix zero[ncols(J_D)][ncols(J_F)];
     533  matrix g_DA=concat(unitmat(ncols(J_D)),zero);
     534  matrix g_AF=transpose(concat(transpose(zero),unitmat(ncols(J_F))));
     535  matrix zero1[ncols(J_DF)][ncols(J_C)];
     536  matrix g_AB=concat(unitmat(ncols(J_DF)),zero1);
     537  matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C))));
     538  list out=(list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F),list(m_d),list(m_a),list(m_f)),list(list(J_DF),list(g_AB),list(J_DFC),list(g_BC),list(J_C),list(m_a),list(m_b),list(v)));
     539  return(out),
     540    }
     541
     542///////////////////////////////////////////////////////////////////////////////////////////
     543
     544proc shortexactpiecestoVdstrict(list C, int d,list #)
     545
     546{
     547
     548  int s =size(C);
     549  if (size(#)==0)
     550    {
     551      intvec v=0:ncols(C[s][1][5]);
     552    }
     553  else
     554    {
     555      intvec v=#[1];
     556    }
     557  list out;
     558  out[s]=list(toVdstrictsequence(C[s][1],d,v));
     559  out[s][2]=list(list(out[s][1][3][1]),list(unitmat(ncols(out[s][1][3][1]))),list(out[s][1][3][1]),list(list()),list(list()));
     560  out[s][2][6]=list(out[s][1][7][1]);
     561  out[s][2][7]=list(out[s][2][6][1]);
     562  out[s][2][8]=list(list());
     563  int i;
     564  for (i=s-1; i>=2; i--)
     565    {
     566      C[i][2][5]=out[i+1][1][1][1];
     567      out[i]=toVdstrictsequences(C[i],d,out[i+1][1][6][1]);
     568    }
     569  out[1]=list(list());
     570  out[1][2]=toVdstrictsequence(C[1][2],d,out[2][1][6][1]);
     571  out[1][1][3]=list(out[1][2][1][1]);
     572  out[1][1][5]=list(out[1][2][1][1]);
     573  out[1][1][4]=list(unitmat(ncols(out[1][1][3][1])));
     574  out[1][1][7]=list(out[1][2][6][1]);
     575  out[1][1][8]=list(out[1][2][6][1]);
     576  out[1][1][1]=list(list());
     577  out[1][1][2]=list(list());
     578  out[1][1][6]=list(list());
     579  list Hi;
     580  for (i=1; i<=size(out); i++)
     581    {
     582      Hi[i]=list(out[i][1][5][1],out[i][1][8][1]);
     583    }
     584  list outall;
     585  outall[1]=out;
     586  print (out);
     587  outall[2]=Hi;
     588  return(outall);
     589
     590
     591}
     592
     593///////////////////////////////////////////////////////////////////////////////////////////
     594
     595proc toVdstrict2x3complex(list L, int d, list #)
     596{
     597  matrix rem; int i; int j;
     598  list J_A=list(list());
     599  list J_B=list(list());
     600  list J_C=list(list());
     601  list g_AB=list(list());
     602  list g_BC=list(list());
     603  list n_a=list(list());
     604  list n_b=list(list());
     605  list n_c=list(list());
     606  intvec n_b1;
     607  if (size(L[5])!=0)
     608    {
     609      intvec n_c1;
     610      for (i=1; i<=nrows(L[5]); i++)
     611        {
     612          rem=submat(L[5],i,intvec(1..ncols(L[5])));
     613          n_c1[i]=Vddeg(rem,d, L[8]);
     614        }
     615      n_c[1]=n_c1;
     616      J_C[1]=transpose(syz(transpose(L[5])));
     617      if (J_C[1]!=matrix(0,nrows(J_C[1]),ncols(J_C[1])))
     618        {
     619          J_C[1]=VdstrictGB(J_C[1],d,n_c1);
     620          if (size(#[2])!=0)
     621            {
     622              n_a[1]=#[2];
     623              n_b1=n_a[1],n_c[1];
     624              n_b[1]=n_b1;
     625              matrix zero[nrows(L[1])][nrows(L[5])];
     626              g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));
     627              if (size(#[1])!=0)
     628                {
     629                  J_A=#[1];
     630                  J_B=transpose(matrix(syz(transpose(L[3]))));
     631                  matrix P=matrixlift(J_B[1] * prodr(nrows(L[1]),nrows(L[5])) ,J_C[1]);
     632
     633                  matrix Pi[1][ncols(J_B[1])];
     634                  matrix Picombined;
     635                  for (i=1; i<=nrows(J_C[1]); i++)
     636                    {
     637                      for (j=1; j<=nrows(J_B[1]);j++)
     638                        {
     639                          Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1])));
     640
     641                        }
     642                      if (i==1)
     643                        {
     644                          Picombined=transpose(Pi);
     645                        }
     646                      else
     647                        {
     648                          Picombined=concat(Picombined,transpose(Pi));
     649                        }
     650                      Pi=0;
     651                    }
     652                  Picombined=transpose(Picombined);
     653                  Picombined=concat(Vdnormalform(Picombined,J_A[1],d,n_a[1]),submat(Picombined,intvec(1..nrows(Picombined)),intvec((ncols(J_A[1])+1)..ncols(Picombined))));
     654                  J_B[1]=transpose(concat(transpose(matrix(J_A[1],nrows(J_A[1]),ncols(J_B[1]))),transpose(Picombined)));
     655                  g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
     656                }
     657              else
     658                {
     659                  J_B[1]=concat(matrix(0,nrows(J_C[1]),nrows(L[3])-nrows(L[5])),J_C[1]);
     660                  g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
     661                }
     662            }
     663          else
     664            {
     665              n_b=n_c[1];
     666              J_B[1]=J_C[1];
     667              g_BC=unitmat(ncols(J_C[1]));
     668            }
     669        }
     670      else
     671        {
     672          J_C=list(list());
     673          if (size(#[2])!=0)
     674            {
     675              matrix zero[nrows(L[1])][nrows(L[5])];
     676              g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
     677              n_a[1]=#[2];
     678              n_b1=n_a[1],n_c[1];
     679              n_b[1]=n_b1;
     680              g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));;
     681
     682              if (size(#[1])!=0)
     683                {
     684                  J_A=#[1];
     685                  J_B=concat(J_A[1],matrix(0,nrows(J_A[1]),nrows(L[3])-nrows(L[1])));
     686                }
     687            }
     688          else
     689            {
     690              n_b=n_c[1];
     691              g_BC=unitmat(ncols(L[5]));
     692            }
     693
     694        }
     695    }
     696  else
     697    {
     698      if (size(#[2])!=0)
     699        {
     700          n_a[1]=#[2];
     701          n_b=n_a[1];
     702          g_AB=unitmat(size(n_b[1]));
     703          if (size(#[1])!=0)
     704            {
     705              J_A=#[1];
     706              J_B[1]=J_A[1];
     707            }
     708        }
     709    }
     710  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]);
     711  return (out);
     712}
     713
     714
     715//////////////////////////////////////////////////////////////////////////
     716
     717proc Vdstrictdoublecompexes(list L, int d)
     718{
     719  int i; int k; int c; int j;
     720  intvec n_b;
     721  matrix rem;
     722  matrix J_B;
     723  list store;
     724  int t=size(L)+nvars(basering) div 2-2;
     725  for (k=1; k<=(size(L)+nvars(basering) div 2-3); k++)//
     726    {
     727      L[1][1][1][k+1]=list();
     728      L[1][1][2][k+1]=list();
     729      L[1][1][6][k+1]=list();
     730      if (size(L[1][1][3][k])!=0)
     731        {
     732          for (i=1; i<=nrows(L[1][1][3][k]); i++)
     733            {
     734              rem=submat(L[1][1][3][k],i,(1..ncols(L[1][1][3][k])));
     735              n_b[i]=Vddeg(rem,d,L[1][1][7][k]);
     736            }
     737          J_B=transpose(syz(transpose(L[1][1][3][k])));
     738          L[1][1][7][k+1]=n_b;
     739          L[1][1][8][k+1]=n_b;
     740          L[1][1][4][k+1]=unitmat(nrows(L[1][1][3][k]));
     741          if (J_B!=matrix(0,nrows(J_B),ncols(J_B)))
     742            {
     743              J_B=VdstrictGB(J_B,d,n_b);
     744              L[1][1][3][k+1]=J_B;
     745              L[1][1][5][k+1]=J_B;
     746            }
     747          else
     748            {
     749              L[1][1][3][k+1]=list();
     750              L[1][1][5][k+1]=list();
     751            }
     752          n_b=0;
     753        }
     754      else
     755        {
     756          L[1][1][3][k+1]=list();
     757          L[1][1][5][k+1]=list();
     758          L[1][1][7][k+1]=list();
     759          L[1][1][8][k+1]=list();
     760          L[1][1][4][k+1]=list();
     761        }
     762      for (i=1; i<size(L); i++)
     763        {
     764          store=toVdstrict2x3complex(list(L[i][2][1][k],L[i][2][2][k],L[i][2][3][k],L[i][2][4][k],L[i][2][5][k],L[i][2][6][k],L[i][2][7][k],L[i][2][8][k]),d,L[i][1][3][k+1],L[i][1][7][k+1]);
     765          for (j=1; j<=8; j++)
     766            {
     767              L[i][2][j][k+1]=store[j];
     768            }
     769
     770          store=toVdstrict2x3complex(list(L[i+1][1][1][k],L[i+1][1][2][k],L[i+1][1][3][k],L[i+1][1][4][k],L[i+1][1][5][k],L[i+1][1][6][k],L[i+1][1][7][k],L[i+1][1][8][k]),d,L[i][2][5][k+1],L[i][2][8][k+1]);
     771
     772          for (j=1; j<=8; j++)
     773            {
     774              L[i+1][1][j][k+1]=store[j];
     775            }
     776        }
     777      if (size(L[size(L)][1][7][k+1])!=0)
     778        {
     779          L[size(L)][2][4][k+1]=list();
     780          L[size(L)][2][5][k+1]=list();
     781          L[size(L)][2][6][k+1]=L[size(L)][1][7][k+1];
     782          L[size(L)][2][7][k+1]=L[size(L)][1][7][k+1];
     783          L[size(L)][2][8][k+1]=list();
     784          L[size(L)][2][2][k+1]=unitmat(size(L[size(L)][1][7][k+1]));
     785
     786          if (size(L[size(L)][1][3][k+1])!=0)
     787            {
     788              L[size(L)][2][1][k+1]=L[size(L)][1][3][k+1];
     789              L[size(L)][2][3][k+1]=L[size(L)][1][3][k+1];
     790            }
     791          else
     792            {
     793              L[size(L)][2][1][k+1]=list();
     794              L[size(L)][2][3][k+1]=list();
     795            }
     796        }
     797      else
     798        {
     799          for (j=1; j<=8; j++)
     800            {
     801              L[size(L)][2][j][k+1]=list();
     802            }
     803        }
     804    }
     805
     806
     807  k=t;
     808  intvec n_c;
     809  intvec vn_b;
     810  list N_b;
     811  int n;
     812  for (i=1; i<=size(L); i++)
     813    {
     814      for (n=1; n<=2; n++)
     815        {
     816          if (i==1 and n==1)
     817            {
     818              L[i][n][6][k+1]=list();
     819            }
     820          else
     821            {
     822              if (n==1)
     823                {
     824                  L[i][1][6][k+1]=L[i-1][2][8][k+1];
     825                }
     826              else
     827                {
     828                  L[i][2][6][k+1]=L[i][1][7][k+1];
     829                }
     830            }
     831          N_b[1]=L[i][n][6][k+1];
     832          if (size(L[i][n][5][k])!=0)
     833            {
     834              for (j=1; j<=nrows(L[i][n][5][k]); j++)
     835                {
     836                  rem=submat(L[i][n][5][k],j,(1..ncols(L[i][n][5][k])));
     837                  n_c[j]=Vddeg(rem,d,L[i][n][8][k]);
     838                }
     839              L[i][n][8][k+1]=n_c;
     840            }
     841          else
     842            {
     843              L[i][n][8][k+1]=list();
     844            }
     845          N_b[2]=L[i][n][8][k+1];
     846          n_c=0;
     847          if (size(N_b[1])!=0)
     848            {
     849              vn_b=N_b[1];
     850              if (size(N_b[2])!=0)
     851                {
     852                  vn_b=vn_b,N_b[2];
     853                }
     854              L[i][n][7][k+1]=vn_b;
     855            }
     856          else
     857            {
     858              if (size(N_b[2])!=0)
     859                {
     860                  L[i][n][7][k+1]=N_b[2];
     861                }
     862              else
     863                {
     864                  L[i][n][7][k+1]=list();
     865                }
     866            }
     867
     868        }
     869    }
     870  return(L);
     871}
     872
     873////////////////////////////////////////////////////////////////////////////
     874
     875proc assemblingdoublecomplexes(list L)
     876{
     877  list out;
     878  int i; int j;int k;int l; int oldj; int newj;
     879  for (i=1; i<=size(L); i++)
     880    {
     881      out[i]=list(list());
     882      out[i][1][1]=ncols(L[i][2][3][1]);
     883      if (size(L[i][2][5][1])!=0)
     884        {
     885          out[i][1][4]=prodr(ncols(L[i][2][3][1])-ncols(L[i][2][5][1]),ncols(L[i][2][5][1]));
     886        }
     887      else
     888        {
     889          out[i][1][4]=matrix(0,ncols(L[i][2][3][1]),1);
     890        }
     891
     892      oldj=newj;
     893      for (j=1; j<=size(L[i][2][3]);j++)
     894        {
     895          out[i][j][2]=L[i][2][7][j];
     896          if (size(L[i][2][3][j])==0)
     897            {
     898              newj =j;
     899              break;
     900            }
     901          out[i][j+1]=list();
     902          out[i][j+1][1]=nrows(L[i][2][3][j]);
     903          out[i][j+1][3]=L[i][2][3][j];
     904          if (size(L[i][2][5][j])!=0)
     905            {
     906              out[i][j+1][4]=(-1)^j*prodr(nrows(L[i][2][3][j])-nrows(L[i][2][5][j]),nrows(L[i][2][5][j]));
     907            }
     908          else
     909            {
     910              out[i][j+1][4]=matrix(0,nrows(L[i][2][3][j]),1);
     911            }
     912          if(j==size(L[i][2][3]))
     913            {
     914              out[i][j+1][2]=L[i][2][7][j+1];
     915              newj=j+1;
     916            }
     917        }
     918      if (i>1)
     919        {
     920          for (k=1; k<=Min(list(oldj,newj)); k++)
     921            {
     922              out[i-1][k][4]=matrix(out[i-1][k][4],nrows(out[i-1][k][4]),out[i][k][1]);
     923            }
     924          for (k=newj+1; k<=oldj; k++)
     925            {
     926              out[i-1][k]=delete(out[i-1][k],4);
     927            }
     928        }
     929    }
     930  return (out);
     931}
     932
     933//////////////////////////////////////////////////////////////////////////////
     934
     935proc totalcomplex(list L);
     936{
     937  list out;intvec rem1;intvec v; list remsize; int emp;
     938  int i; int j; int c; int d; matrix M; int k; int l;
     939  int n=nvars(basering) div 2;
     940  list K;
     941  for (i=1; i<=n; i++)
     942    {
     943      K[i]=list();
     944    }
     945  L=K+L;
     946  for (i=1; i<=size(L); i++)
     947    {
     948      emp=0;
     949      if (size(L[i])!=0)
     950        {
     951          out[3*i-2]=L[i][1][1];
     952          v=L[i][1][1];
     953          rem1=L[i][1][2];
     954          emp=1;
     955        }
     956      else
     957        {
     958          out[3*i-2]=0;
     959          v=0;
     960        }
     961
     962      for (j=i+1; j<=size(L); j++)
     963        {
     964          if (size(L[j])>=j-i+1)
     965            {
     966              out[3*i-2]=out[3*i-2]+L[j][j-i+1][1];
     967              if (emp==0)
     968                {
     969                  rem1=L[j][j-i+1][2];
     970                  emp=1;
     971                }
     972              else
     973                {
     974                  rem1=rem1,L[j][j-i+1][2];
     975                }
     976              v[size(v)+1]=L[j][j-i+1][1];
     977            }
     978          else
     979            {
     980              v[size(v)+1]=0;
     981            }
     982        }
     983      out[3*i-1]=rem1;
     984      v[size(v)+1]=0;
     985      remsize[i]=v;
     986    }
     987  int o1;
     988  int o2;
     989  for (i=1; i<=size(L)-1; i++)
     990    {
     991      o1=1;
     992      o2=1;
     993      if (size(out[3*i-2])!=0)
     994        {
     995          o1=out[3*i-2];
     996        }
     997      if (size(out[3*i+1])!=0)
     998        {
     999          o2=out[3*i+1];
     1000        }
     1001      M=matrix(0,o1,o2);
     1002      if (size(L[i])!=0)
     1003        {
     1004          if (size(L[i][1][4])!=0)
     1005            {
     1006              M=matrix(L[i][1][4],o1,o2);
     1007            }
     1008        }
     1009      c=remsize[i][1];
     1010      // d=remsize[i+1][1];
     1011      for (j=i+1; j<=size(L); j++)
     1012        {
     1013          if (remsize[i][j-i+1]!=0)
     1014            {
     1015              for (k=c+1; k<=c+remsize[i][j-i+1]; k++)
     1016                {
     1017                  for (l=d+1; l<=d+remsize[i+1][j-i];l++)
     1018                    {
     1019                      M[k,l]=L[j][j-i+1][3][(k-c),(l-d)];
     1020                    }
     1021                }
     1022              d=d+remsize[i+1][j-i];
     1023              if (remsize[i+1][j-i+1]!=0)
     1024                {
     1025                  for (k=c+1; k<=c+remsize[i][j-i+1]; k++)
     1026                    {
     1027                      for (l=d+1; l<=d+remsize[i+1][j-i+1];l++)
     1028                        {
     1029                          M[k,l]=L[j][j-i+1][4][k-c,l-d];
     1030                        }
     1031                    }
     1032                  c=c+remsize[i][j-i+1];
     1033                }
     1034            }
     1035          else
     1036            {
     1037              d=d+remsize[i+1][j-i];
     1038            }
     1039        }
     1040      out[3*i]=M;
     1041      d=0; c=0;
     1042    }
     1043  out[3*size(L)]=matrix(0,out[3*size(L)-2],1);
     1044  return (out);
     1045}
     1046
     1047/////////////////////////////////////////////////////////////////////////////////////
     1048
     1049proc toVdstrictfreecomplex(list L,list #)
     1050{
     1051  def B=basering;
     1052  int n=nvars(B) div 2+2;
     1053  int d=nvars(B) div 2;
     1054  intvec v;
     1055  list out;list outall;
     1056  int i;int j;
     1057  matrix mem;
     1058  int k;
     1059  if (size(#)!=0)
     1060    {
     1061      for (i=1; i<=size(#); i++)
     1062        {
     1063          if (typeof(#[i])==intvec)
     1064            {
     1065              v=#[i];
     1066            }
     1067          if (typeof(#[i])==int)
     1068            {
     1069              d=#[i];
     1070            }
     1071        }
     1072    }
     1073  if (size(L)==2)
     1074    {
     1075      v=(0:ncols(L[1]));
     1076      out[3*n-1]=v;
     1077      out[3*n-2]=ncols(L[1]);
     1078      out[3*n]=L[2];
     1079      out[3*n-3]=VdstrictGB(L[1],d,v);
     1080      for (i=n-1; i>=1; i--)
     1081        {
     1082          out[3*i-2]=nrows(out[3*i]);
     1083          v=0;
     1084          for (j=1; j<=out[3*i-2]; j++)
     1085            {
     1086              mem=submat(out[3*i],j,intvec(1..ncols(out[3*i])));
     1087              v[j]=Vddeg(mem,d, out[3*i+2]);
     1088            }
     1089          out[3*i-1]=v;
     1090          if (i!=1)
     1091            {
     1092              out[3*i-3]=transpose(syz(transpose(out[3*i])));
     1093              if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3])))
     1094                {
     1095                  out[3*i-3]=VdstrictGB(out[3*i-3],d,out[3*i-1]);
     1096                }
     1097              else
     1098                {
     1099                  out[3*i-3]=matrix(0,1,ncols(out[3*i-3]));
     1100                  out[3*i-4]=intvec(0);
     1101                  out[3*i-5]=int(0);
     1102                  for (j=i-2; j>=1; j--)
     1103                    {
     1104                      out[3*j]=matrix(0,1,1);
     1105                      out[3*j-1]=intvec(0);
     1106                      out[3*j-2]=int(0);
     1107                    }
     1108                  break;
     1109                }
     1110            }
     1111        }
     1112      outall[1]=out;
     1113      outall[2]=list(list(out[3*n-3],out[3*n-1]));
     1114      return(outall);
     1115    }
     1116  out=shortexactpieces(L);
     1117  list rem;
     1118  if (v!=intvec(0:size(v)))
     1119    {
     1120      rem=shortexactpiecestoVdstrict(out,d,v);
     1121    }
     1122  else
     1123    {
     1124      rem=shortexactpiecestoVdstrict(out,d);
     1125    }
     1126  out=Vdstrictdoublecompexes(rem[1],d);
     1127  out=assemblingdoublecomplexes(out);
     1128  out=totalcomplex(out);
     1129  outall[1]=out;
     1130  outall[2]=rem[2];
     1131  return (outall);
     1132}
     1133
     1134////////////////////////////////////////////////////////////////////////////////
     1135
     1136proc derhamcohomology(list L)
     1137{
     1138  def R=basering;
     1139  int n=nvars(R);int le=2*size(L)+n-1;
     1140  def W=makeWeyl(n);
     1141  setring W;
     1142  list man=ringlist(W);
     1143  if (n==1)
     1144    {
     1145      man[2][1]="x(1)";
     1146      man[2][2]="D(1)";
     1147      def Wi=ring(man);
     1148      setring Wi;
     1149      kill W;
     1150      def W=Wi;
     1151      setring W;
     1152      list man=ringlist(W);
     1153    }
     1154  man[2][size(man[2])+1]="s";;
     1155  man[3][3]=man[3][2];
     1156  man[3][2]=list("dp",intvec(1));
     1157  matrix N=UpOneMatrix(size(man[2]));
     1158  man[5]=N;
     1159  matrix M[1][1];
     1160  man[6]=transpose(concat(transpose(concat(man[6],M)),M));
     1161  def Ws=ring(man); setring R;  int r=size(L); int i;  int j;int k; int l; int count;  list Fi; list subsets; list maxnum; list bernsteinpolys; list annideals; list minint; list diffmaps;
     1162  for (i=1; i<=r; i++)
     1163    {
     1164      Fi[i]=list(); bernsteinpolys[i]=list(); annideals[i]=list(); subsets[i]=list();
     1165      maxnum[i]=list();
     1166      Fi[1][i]=L[i];
     1167      maxnum[1][i]=i;
     1168      subsets[1][i]=intvec(i);
     1169    }
     1170  intvec v;
     1171  for (i=2; i<=r; i++)
     1172    {
     1173      count=1;
     1174      for (j=1; j<=size(Fi[i-1]);j++)
     1175        {
     1176          for (k=maxnum[i-1][j]+1; k<=r; k++)
     1177            {
     1178              maxnum[i][count]=k;
     1179              v=subsets[i-1][j],k;
     1180              subsets[i][count]=v;
     1181              Fi[i][count]=lcm(Fi[i-1][j],L[k]);/////////
     1182              count=count+1;
     1183            }
     1184        }
     1185    }
     1186  for (i=1; i<=r; i++)
     1187    {
     1188      for (j=1; j<=size(Fi[i]); j++)
     1189        {
     1190          bernsteinpolys[i][j]=bfct(Fi[i][j])[1];
     1191          for (k=1; k<=ncols(bernsteinpolys[i][j]); k++)
     1192            {
     1193              if (isInt(number(bernsteinpolys[i][j][k]))==1)
     1194                {
     1195                  minint[size(minint)+1]=int(bernsteinpolys[i][j][k]);
     1196                }
     1197            }
     1198          def D=Sannfs(Fi[i][j]);
     1199          setring Ws;
     1200          annideals[i][j]=fetch(D,LD);
     1201          kill D;
     1202          setring R;
     1203        }
     1204    }
     1205  int m=Min(minint);
     1206  list zw;
     1207  for (i=1; i<r; i++)
     1208    {
     1209      diffmaps[i]=matrix(0,size(subsets[i]),size(subsets[i+1]));
     1210      for (j=1; j<=size(subsets[i]); j++)
     1211        {
     1212          for (k=1; k<=size(subsets[i+1]); k++)
     1213            {
     1214              zw=mysubset(subsets[i][j],subsets[i+1][k]);
     1215              diffmaps[i][j,k]=zw[2]*(L[zw[1]]/gcd(L[zw[1]],Fi[i][j]))^(-m);
     1216            }
     1217        }
     1218    }
     1219  diffmaps[r]=matrix(0,1,1);
     1220  setring Ws;
     1221  for (i=1; i<=r; i++)
     1222    {
     1223      for (j=1; j<=size(annideals[i]); j++)
     1224        {
     1225          annideals[i][j]=subst(annideals[i][j],s,m);
     1226        }
     1227    }
     1228  setring W;
     1229  list annideals=imap(Ws,annideals);
     1230  list diffmaps=fetch(R,diffmaps);
     1231  list fortoVdstrict;
     1232  ideal IFourier=var(n+1);
     1233  for (i=2;i<=n;i++)
     1234    {
     1235      IFourier=IFourier,var(n+i);
     1236    }
     1237  for (i=1; i<=n;i++)
     1238    {
     1239      IFourier=IFourier,-var(i);
     1240    }
     1241  map cFourier=W,IFourier;
     1242  matrix sup;
     1243  for (i=1; i<=r; i++)
     1244    {
     1245      sup=matrix(annideals[i][1]);
     1246      fortoVdstrict[2*i-1]=transpose(cFourier(sup));
     1247      for (j=2; j<=size(annideals[i]); j++)
     1248        {
     1249          sup=matrix(annideals[i][j]);
     1250          fortoVdstrict[2*i-1]=dsum(fortoVdstrict[2*i-1],transpose(cFourier(sup)));
     1251        }
     1252      sup=diffmaps[i];
     1253      fortoVdstrict[2*i]=cFourier(sup);
     1254    }
     1255  list rem=toVdstrictfreecomplex(fortoVdstrict);
     1256  list newcomplex=rem[1];
     1257  list minmaxk=globalbfun(rem[2]);
     1258  if (size(minmaxk)==0)
     1259    {
     1260      return (0);
     1261    }
     1262  list truncatedcomplex; list shorten; list  upto;
     1263  for (i=1; i<=size(newcomplex) div 3; i++)
     1264    {
     1265      shorten[3*i-1]=list();
     1266      for (j=1; j<=size(newcomplex[3*i-1]); j++)
     1267        {
     1268          shorten[3*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1,minmaxk[2]-newcomplex[3*i-1][j]+1);
     1269          upto[size(upto)+1]=shorten[3*i-1][j][2];
     1270          if (shorten[3*i-1][j][2]<=0)
     1271            {
     1272              shorten[3*i-1][j]=list();
     1273            }
     1274          else
     1275            {
     1276              if (shorten[3*i-1][j][1]<=0)
     1277                {
     1278                  shorten[3*i-1][j][1]=1;
     1279                }
     1280            }
     1281        }
     1282    }
     1283  int iupto=Max(upto);
     1284  if (iupto<=0)
     1285    {
     1286      /////die Kohomologie ist dann überall 0, muss noch entsprechend ausgegeben werden
     1287    }
     1288  list allpolys;
     1289  allpolys[1]=list(1);
     1290  list minvar;
     1291  minvar[1]=list(1);
     1292  for (i=1; i<=iupto-1; i++)
     1293    {
     1294      allpolys[i+1]=list();
     1295      minvar[i+1]=list();
     1296      for (k=1; k<=size(allpolys[i]); k++)
     1297        {
     1298          for (j=minvar[i][k]; j<=nvars(W) div 2; j++)
     1299            {
     1300              allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j);
     1301              minvar[i+1][size(minvar[i+1])+1]=j;
     1302            }
     1303        }
     1304    }
     1305  list keepformatrix;list sizetruncom;int stc;list fortrun;
     1306  for (i=1; i<=size(newcomplex) div 3; i++)
     1307    {
     1308      truncatedcomplex[2*i-1]=list();
     1309      sizetruncom[2*i-1]=list();
     1310      sizetruncom[2*i]=list();
     1311      truncatedcomplex[2*i]=newcomplex[3*i];
     1312      v=0;count=0;
     1313      sizetruncom[2*i][1]=0;
     1314      for (j=1; j<=newcomplex[3*i-2]; j++)
     1315        {
     1316          if (size(shorten[3*i-1][j])!=0)
     1317            {
     1318              fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]);
     1319              truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1];
     1320              count=count+fortrun[2];
     1321              sizetruncom[2*i-1][size(sizetruncom[2*i-1])+1]=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1);
     1322              sizetruncom[2*i][size(sizetruncom[2*i])+1]=count;
     1323              if (v!=0)
     1324                {
     1325                  v[size(v)+1]=j;
     1326                }
     1327              else
     1328                {
     1329                  v[1]=j;
     1330                }
     1331            }
     1332        }
     1333
     1334      if (v!=0)
     1335        {
     1336          truncatedcomplex[2*i]=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i]));
     1337          if (i!=1)
     1338            {
     1339              truncatedcomplex[2*(i-1)]=submat(truncatedcomplex[2*(i-1)],1..nrows(truncatedcomplex[2*(i-1)]),v);
     1340            }
     1341        }
     1342      else
     1343        {
     1344          truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i]));
     1345          if (i!=1)
     1346            {
     1347              truncatedcomplex[2*(i-1)]=matrix(0,nrows(truncatedcomplex[2*(i-1)]),1);
     1348            }
     1349        }
     1350    }
     1351  int b;int d;poly form;poly lform; poly nform;int ideg;int kplus; int lplus;
     1352  for (i=1; i<size(truncatedcomplex) div 2; i++)
     1353    {
     1354      M=matrix(0,max(1,sizetruncom[2*i][size(sizetruncom[2*i])]),sizetruncom[2*i+2][size(sizetruncom[2*i+2])]);
     1355      for (k=1; k<=size(truncatedcomplex[2*i-1]);k++)
     1356        {
     1357          for (l=1; l<=size(truncatedcomplex[2*(i+1)-1]); l++)
     1358            {
     1359              if (size(sizetruncom[2*i])!=1)//?
     1360                {
     1361                  for (j=1; j<=size(truncatedcomplex[2*i-1][k]); j++)
     1362                    {
     1363                      for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++)
     1364                        {
     1365                          form=truncatedcomplex[2*i-1][k][j][b][1]*truncatedcomplex[2*i][k,l];
     1366                          while (form!=0)
     1367                            {
     1368                              lform=lead(form);
     1369                              v=leadexp(lform);
     1370                              v=v[1..n];
     1371                              if (v==(0:n))
     1372                                {
     1373                                  ideg=deg(lform)-sizetruncom[2*(i+1)-1][l][1];
     1374                                  if (ideg>=0)
     1375                                    {
     1376                                      for (d=1; d<=size(truncatedcomplex[2*(i+1)-1][l][ideg+1]);d++)
     1377                                        {
     1378                                          if (leadmonom(lform)==truncatedcomplex[2*(i+1)-1][l][ideg+1][d][1])
     1379                                            {
     1380                                              M[sizetruncom[2*i][k]+truncatedcomplex[2*i-1][k][j][b][2],sizetruncom[2*(i+1)][l]+truncatedcomplex[2*(i+1)-1][l][ideg+1][d][2]]=leadcoef(lform);
     1381                                              break;
     1382                                            }
     1383                                        }
     1384                                    }
     1385                                }
     1386                              form=form-lform;
     1387                            }
     1388                        }
     1389                    }
     1390                }
     1391            }
     1392        }
     1393      truncatedcomplex[2*i]=M;
     1394      truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])];
     1395    }
     1396  truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])];
     1397  if (truncatedcomplex[2*i-1]!=0)
     1398    {
     1399      truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1);
     1400    }
     1401  setring R;
     1402  list truncatedcomplex=imap(W,truncatedcomplex);
     1403  list derhamhom=findhomology(truncatedcomplex,le);
     1404  return (derhamhom);
     1405}
     1406
     1407///////////////////////////////////
     1408static proc sublist(list L, int m, int n)
     1409{
     1410  list out;
     1411  int i; int j;
     1412  int count;
     1413  for (i=m; i<=n; i++)
     1414    {
     1415      out[size(out)+1]=list();
     1416      for (j=1; j<=size(L[i]); j++)
     1417        {
     1418          count=count+1;
     1419          out[size(out)][j]=list(L[i][j],count);
     1420        }
     1421    }
     1422  list o=list(out,count);
     1423  return(o);
     1424}
     1425
     1426//////////////////////////////////////////////////////////////////////////
     1427static proc mysubset(intvec L, intvec M)
     1428{
     1429  int i;
     1430  int j=1;
     1431  list position=(M[size(M)],(-1)^(size(L)));
     1432  for (i=1; i<=size(L); i++)
     1433    {
     1434      if (L[i]!=M[j])
     1435        {
     1436          if (L[i]!=M[j+1] or j!=i)
     1437            {
     1438              return (L[i],0);
     1439            }
     1440          else
     1441            {
     1442              position=(M[i],(-1)^(i-1));
     1443              j=j+i;
     1444            }
     1445        }
     1446      j=j+1;
     1447    }
     1448  return (position);
     1449}
     1450
     1451
     1452
     1453////////////////////////////////////////////////////////////////////////////
     1454
     1455proc globalbfun(list L)
     1456{
     1457  int i; int j;
     1458  def W=basering;
     1459  int n=nvars(W) div 2;
     1460  list G0;
     1461  ideal I;
     1462  for (j=1; j<=size(L); j++)
     1463    {
     1464      G0[j]=list();
     1465      for (i=1; i<=ncols(L[j][1]); i++)
     1466        {
     1467          G0[j][i]=I;
     1468        }
     1469    }
     1470  list out;
     1471  for (j=1; j<=size(L); j++)
     1472    {
     1473      for (i=1; i<=nrows(L[j][1]); i++)
     1474        {
     1475          out=Vddeg(submat(L[j][1],i,(1..ncols(L[j][1]))),n,L[j][2],1);
     1476          G0[j][out[2]][size(G0[j][out[2]])+1]=(out[1]);
     1477        }
     1478    }
     1479  list Data=ringlist(W);
     1480  for (i=1; i<=n; i++)
     1481    {
     1482      Data[2][2*n+i]=Data[2][i];
     1483      Data[2][3*n+i]=Data[2][n+i];
     1484      Data[2][i]="v("+string(i)+")";
     1485      Data[2][n+i]="w("+string(i)+")";
     1486    }
     1487  Data[3][1][1]="M";
     1488  intvec mord=(0:16*n^2);
     1489  mord[1..2*n]=(1:2*n);
     1490  mord[6*n+1..8*n]=(1:2*n);
     1491  for (i=0; i<=2*n-2; i++)
     1492    {
     1493      mord[(3+i)*4*n-i]=-1;
     1494      mord[(2*n+2+i)*4*n-2*n-i]=-1;
     1495    }
     1496  Data[3][1][2]=mord;//ordering mh?????????
     1497  matrix Ones=UpOneMatrix(4*n);
     1498  Data[5]=Ones;
     1499  matrix con[2*n][2*n];
     1500  Data[6]=transpose(concat(con,transpose(concat(con,Data[6]))));
     1501
     1502  def Wuv=ring(Data);
     1503  setring Wuv;
     1504  list G0=imap(W,G0); list G3; poly lterm;intvec lexp;
     1505  list G1;  list G2; intvec e; intvec f; int  kapp; int k; int l; poly h; ideal I;
     1506  for (l=1; l<=size(G0); l++)
     1507    {
     1508      G1[l]=list();  G2[l]=list(); G3[l]=list();
     1509      for (i=1; i<=size(G0[l]); i++)
     1510        {
     1511          for (j=1; j<=ncols(G0[l][i]);j++)
     1512            {
     1513              G0[l][i][j]=mhom(G0[l][i][j]);
     1514            }
     1515          for (j=1; j<=nvars(Wuv) div 4; j++)
     1516            {
     1517              G0[l][i][size(G0[l][i])+1]=1-v(j)*w(j);
     1518            }
     1519          G1[l][i]=std(G0[l][i]);
     1520          G2[l][i]=I;
     1521          G3[l][i]=list();
     1522          for (j=1; j<=ncols(G1[l][i]); j++)
     1523            {
     1524              e=leadexp(G1[l][i][j]);
     1525              f=e[1..2*n];
     1526              if (f==intvec(0:(2*n)))
     1527                {
     1528                  for (k=1; k<=n; k++)
     1529                    {
     1530                      kapp=-e[2*n+k]+e[3*n+k];
     1531                      if (kapp>0)
     1532                        {
     1533                          G1[l][i][j]=(x(k)^kapp)*G1[l][i][j];
     1534                        }
     1535                      if (kapp<0)
     1536                        {
     1537                          G1[l][i][j]=(D(k)^(-kapp))*G1[l][i][j];
     1538                        }
     1539                    }
     1540                  G2[l][i][size(G2[l][i])+1]=G1[l][i][j];
     1541                  G3[l][i][size(G3[l][i])+1]=list();
     1542                  while (G1[l][i][j]!=0)
     1543                    {
     1544                      lterm=lead(G1[l][i][j]);
     1545                      G1[l][i][j]=G1[l][i][j]-lterm;
     1546                      lexp=leadexp(lterm);
     1547                      lexp=lexp[2*n+1..3*n];
     1548                      G3[l][i][size(G3[l][i])][size(G3[l][i][size(G3[l][i])])+1]=list(lexp,leadcoef(lterm));
     1549                    }
     1550
     1551                }
     1552            }
     1553        }
     1554    }
     1555  ring r=0,(s(1..n)),dp;
     1556  ideal I;
     1557  map G3forr=Wuv,I;
     1558  list G3=G3forr(G3);
     1559  poly fs;
     1560  poly gs;
     1561  int a;
     1562  list G4;
     1563  for (l=1; l<=size(G3); l++)
     1564    {
     1565      G4[l]=list();
     1566      for (i=1; i<=size(G3[l]);i++)
     1567        {
     1568          G4[l][i]=I;
     1569          for (j=1; j<=size(G3[l][i]); j++)
     1570            {
     1571              fs=0;
     1572              for (k=1; k<=size(G3[l][i][j]); k++)
     1573                {
     1574                  gs=1;
     1575                  for (a=1; a<=n; a++)
     1576                    {
     1577                      if (G3[l][i][j][k][1][a]!=0)
     1578                        {
     1579                          gs=gs*permutevar(list(G3[l][i][j][k][1][a]),a);
     1580                        }
     1581                    }
     1582                  gs=gs*G3[l][i][j][k][2];
     1583                  fs=fs+gs;
     1584                }
     1585              G4[l][i]=G4[l][i],fs;
     1586            }
     1587        }
     1588    }
     1589  if (n==1)
     1590    {
     1591      ring rnew=0,t,dp;
     1592    }
     1593  else
     1594    {
     1595      ring rnew=0,(t,s(2..n)),dp;
     1596    }
     1597  ideal Iformap;
     1598  Iformap[1]=t;
     1599   poly forel=1;
     1600   for (i=2; i<=n; i++)
     1601     {
     1602       Iformap[1]=Iformap[1]-s(i);
     1603       Iformap[i]=s(i);
     1604       forel=forel*s(i);
     1605     }
     1606   map rtornew=r,Iformap;
     1607   list G4=rtornew(G4);
     1608   list getintvecs=fetch(W,L);
     1609   ideal J;
     1610   option(redSB);
     1611   for (l=1; l<=size(G4); l++)
     1612     {
     1613       J=1;
     1614       for (i=1; i<=size(G4[l]); i++)
     1615         {
     1616           G4[l][i]=eliminate(G4[l][i],forel);
     1617           G4[l][i]=subst(G4[l][i],t,t-getintvecs[l][2][i]);
     1618           J=intersect(J,G4[l][i]);
     1619         }
     1620       G4[l]=poly(std(J)[1]);
     1621     }
     1622   list minmax=minmaxintroot(G4);//besser factorize nehmen
     1623   // Fall: keine Nullstelle muss noch weiter beruecksichtigt werden
     1624  return(minmax);
     1625}
     1626
     1627
     1628
     1629
     1630//////////////////////////////////////////////////////////////////////////
     1631
     1632proc minmaxintroot(list L);
     1633{
     1634  int i; int j; int k; int l; int sa; int s; number d; poly f; poly rest; list a0; list possk; list alldiv; intvec e;
     1635  possk[1]=list();
     1636  for (i=1; i<=size(L); i++)
     1637    {
     1638      d=1;
     1639      f=L[i];
     1640      while (f!=0)
     1641        {
     1642          rest=lead(f);
     1643          d=d*denominator(leadcoef(rest));
     1644          f=f-rest;
     1645        }
     1646      e=leadexp(rest);
     1647      if (e[1]!=0)
     1648        {
     1649          rest=rest/(t^(e[1]));
     1650          possk[1][size(possk[1])+1]=i;
     1651        }
     1652      a0[i]=int(absValue(d*rest));
     1653    }
     1654  int m=Max(a0);
     1655  for (i=2; i<=m+1; i++)
     1656    {
     1657      possk[i]=list();
     1658    }
     1659  list allprimefac;
     1660  for (i=1; i<=size(L); i++)
     1661    {
     1662      allprimefac=primefactors(a0[i]);
     1663      alldiv=1;
     1664      possk[2][size(possk[2])+1]=i;
     1665
     1666      for (j=1; j<=size(allprimefac[1]); j++)
     1667        {
     1668          s=size(alldiv);
     1669          for (k=1; k<=s; k++)
     1670            {
     1671              for (l=1; l<=allprimefac[2][j]; l++)
     1672                {
     1673                  alldiv[size(alldiv)+1]=alldiv[k]*allprimefac[1][j]^l;
     1674                  possk[alldiv[size(alldiv)]+1][size(possk[alldiv[size(alldiv)]+1])+1]=i;
     1675                }
     1676            }
     1677        }
     1678    }
     1679  int mink;
     1680  int maxk;
     1681  int indi;
     1682  for (i=m+1; i>=1; i--)
     1683    {
     1684      if (size(possk[i])!=0)
     1685        {
     1686          for (j=1; j<=size(possk[i]); j++)
     1687            {
     1688              if (subst(L[possk[i][j]],t,(i-1))==0)
     1689                {
     1690                  maxk=i-1;
     1691                  indi=1;
     1692                  break;
     1693                }
     1694            }
     1695        }
     1696      if (maxk!=0)
     1697        {
     1698          break;
     1699        }
     1700    }
     1701  int indi2;
     1702  for (i=m+1; i>=1; i--)
     1703    {
     1704      if (size(possk[i])!=0)
     1705        {
     1706          for (j=1; j<=size(possk[i]); j++)
     1707            {
     1708              if (subst(L[possk[i][j]],t,-(i-1))==0)
     1709                {
     1710                  mink=-i+1;
     1711                  indi2=1;
     1712                  break;
     1713                }
     1714            }
     1715        }
     1716      if (mink!=0)
     1717        {
     1718          break;
     1719        }
     1720    }
     1721  list mima=mink,maxk;
     1722  if (indi==0)
     1723    {
     1724      if (indi2==0)
     1725        {
     1726          mima=list();//es gibt keine ganzzahlige NS
     1727        }
     1728      else
     1729        {
     1730          mima[2]=mima[1];
     1731        }
     1732    }
     1733  else
     1734    {
     1735      if (indi2==0)
     1736        {
     1737          mima[1]=mima[2];
     1738        }
     1739    }
     1740  return (mima);
     1741}
     1742
     1743///////////////////////////////////////////////////////
     1744
     1745
     1746proc findhomology(list L, int le)
     1747{
     1748  int li;
     1749  matrix M; matrix N;
     1750  matrix N1;
     1751  matrix lift1;
     1752  list out;
     1753  int i;
     1754  option (redSB);
     1755  for (i=2; i<=size(L); i=i+2)
     1756    {
     1757      if (L[i-1]==0)
     1758        {
     1759          li=0;
     1760          out[i div 2]=0;
     1761        }
     1762      else
     1763        {
     1764
     1765          if (li==0)
     1766            {
     1767
     1768
     1769              li=L[i-1];
     1770              N1=transpose(syz(transpose(L[i])));
     1771              out[i div 2]=matrix(transpose(syz(transpose(N1))));
     1772              out[i div 2]=transpose(matrix(std(transpose(out[i div 2]))));
     1773
     1774            }
     1775
     1776          else
     1777            {
     1778
     1779
     1780              li=L[i-1];
     1781              N1=transpose(syz(transpose(L[i])));
     1782              N=transpose(syz(transpose(N1)));
     1783              lift1=matrixlift(N1,L[i-2]);
     1784              out[i div 2]=transpose(concat(transpose(lift1),transpose(N)));
     1785              out[i div 2]=transpose(matrix(std(transpose(out[i div 2]))));
     1786            }
     1787        }
     1788      if (out[i div 2]!=matrix(0,1,ncols(out[i div 2])))
     1789        {
     1790          out[i div 2]=ncols(out[i div 2])-nrows(out[i div 2]);
     1791        }
     1792      else
     1793        {
     1794          out[i div 2]=ncols(out[i div 2]);
     1795        }
     1796    }
     1797  if (size(out)>le)
     1798    {
     1799      out=delete(out,1);
     1800    }
     1801  return(out);
     1802}
     1803
     1804
     1805
     1806
     1807/////////////////////////////////////////////////////////////////////
     1808
     1809static proc mhom(poly f)
     1810{
     1811  poly g;
     1812  poly l;
     1813  poly add;
     1814  intvec e;
     1815  list minint;
     1816  list remf;
     1817  int i;
     1818  int j;
     1819  int n=nvars(basering) div 4;
     1820  if (f==0)
     1821    {
     1822      return(f);
     1823    }
     1824  while (f!=0)
     1825    {
     1826      l=lead(f);
     1827      e=leadexp(l);
     1828      remf[size(remf)+1]=list();
     1829      remf[size(remf)][1]=l;
     1830      for (i=1; i<=n; i++)
     1831        {
     1832          remf[size(remf)][i+1]=-e[2*n+i]+e[3*n+i];
     1833          if (size(minint)<i)
     1834            {
     1835              minint[i]=list();
     1836            }
     1837          minint[i][size(minint[i])+1]=-e[2*n+i]+e[3*n+i];
     1838        }
     1839      f=f-l;
     1840    }
     1841  for (i=1; i<=n; i++)
     1842    {
     1843      minint[i]=Min(minint[i]);
     1844    }
     1845  for (i=1; i<=size(remf); i++)
     1846    {
     1847      add=remf[i][1];
     1848      for (j=1; j<=n; j++)
     1849        {
     1850          add=v(j)^(remf[i][j+1]-minint[j])*add;
     1851        }
     1852      g=g+add;
     1853    }
     1854  return (g);
     1855}
     1856
     1857
     1858
     1859
     1860//////////////////////////////////////////////////////////////////////////
     1861
     1862static proc permutevar(list L,int n)
     1863{
     1864  if (typeof(L[1])=="intvec")
     1865    {
     1866      intvec v=L[1];
     1867    }
     1868  else
     1869    {
     1870      intvec v=(1:L[1]),(0:L[1]);
     1871    }
     1872  int i;int k; int indi=0;
     1873  int j;
     1874  int s=size(v);
     1875  poly e;
     1876  intvec fore;
     1877  for (i=2; i<=size(v); i=i+2)
     1878    {
     1879      if (v[i]!=0)
     1880        {
     1881          j=i+1;
     1882          while (v[j]!=0)
     1883            {
     1884              j=j+1;
     1885            }
     1886          v[i]=0;
     1887          v[j]=1;
     1888          fore=0;
     1889          indi=0;
     1890          for (k=1; k<=size(v); k++)
     1891            {
     1892              if (k!=i and k!=j)
     1893                {
     1894                  if (indi==0)
     1895                    {
     1896                      indi=1;
     1897                      fore[1]=v[k];
     1898                    }
     1899                  else
     1900                    {
     1901                      fore[size(fore)+1]=v[k];
     1902                    }
     1903                }
     1904            }
     1905          e=e-(j-i)*permutevar(list(fore),n);
     1906        }
     1907    }
     1908  e=e+s(n)^(size(v) div 2);
     1909  return (e);
     1910}
     1911
     1912///////////////////////////////////////////////////////////////////////////////
     1913static proc max(int i,int j)
     1914{
     1915  if(i>j){return(i);}
     1916  return(j);
     1917}
     1918
     1919////////////////////////////////////////////////////////////////////////////////////
     1920version="$Id$";
     1921category="Noncommutative";
     1922info="
     1923LIBRARY:  derham.lib      Computation of deRham cohomology
     1924
     1925AUTHORS:  Cornelia Rottner, rottner@mathematik.uni-kl.de
     1926
     1927OVERVIEW:
     1928A library for computing the de Rham cohomology of complements of complex affine
    111929varieties.
    121930
    131931
    14 REFERENCES: 
    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
     1932REFERENCES:
     1933[OT] Oaku, T.; Takayama, N.: Algorithms of D-modules - restriction, tensor product,
     1934     localzation, and local cohomology groups}, J. Pure Appl. Algebra 156, 267-308
    171935     (2001)
    181936[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,
     1937[W1] Walther, U.: Algorithmic computation of local cohomology modules and the local
     1938     cohomological dimension of algebraic varieties}, J. Pure Appl. Algebra 139,
    211939     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)
     1940[W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of
     1941     Complex Affine Varieties}, J. Symbolic Computation 29, 796-839 (2000)
     1942[W3] Walther, U.: Computing the cup product structure for complements of complex
     1943     affine varieties, J. Pure Appl. Algebra 164, 247-273 (2001)
    241944
    251945
    … …  
    391959LIB "poly.lib";
    401960LIB "schreyer.lib";
     1961LIB "dmodloc.lib";
     1962
    411963
    421964////////////////////////////////////////////////////////////////////////////////////
    431965
    441966proc deRhamCohomology(list L,list #)
    45 "USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices 
     1967"USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices
    461968        optional list consisting of one up to three strings @*
    471969        The optional strings may be one of the strings@*
    48         -'Vdres': compute the Cartan-Eilenberg resolutions via V_d-homogenization
     1970        -'noCE': compute quasi-isomorphic complexes without using Cartan-Eilenberg
     1971         resolutionsq@*
     1972        -'Vdres': compute quasi-isomorphic complexes using Cartan-Eilenberg
     1973         resolutions; the CE resolutions are computed via V__d-homogenization
    491974         and without using Schreyer's method @*
    50         -'Sres': compute the Cartan-Eilenberg resolutions in the homogenized Weyl
    51          algebra using Schreyer's method@*
     1975        -'Sres': compute quasi-isomorphic complexes using Cartan-Eilenberg
     1976         resolutions in the homogenized Weyl algebra via Schreyer's method@*
    521977        one of the strings@*
    53         -'iterativeloc': compute localizations by factorizing the polynomials and 
     1978        -'iterativeloc': compute localizations by factorizing the polynomials and
    541979         sucessive localization of the factors @*
    551980        -'no iterativeloc': compute localizations by directly localizing the
    … …  
    601985        -'exactroots' computes the minimal and maximal integer root of the global
    611986         b-function
    62         The default is 'Sres', 'iterativeloc' and 'onlybounds'.
     1987        The default is 'noCE', 'iterativeloc' and 'onlybounds'.
    631988ASSUME: -The basering must be a polynomial ring over the field of rational numbers@*
    64 RETURN: list, where the ith entry is the (i-1)st de Rham cohomology group of the 
    65         complement of the complex affine variety given by the polynomials in L   
     1989RETURN: list, where the ith entry is the (i-1)st de Rham cohomology group of the
     1990        complement of the complex affine variety given by the polynomials in L
    661991EXAMPLE:example deRhamCohomology; shows an example
    671992"
    … …  
    762001  int n=nvars(R);
    772002  int le=size(L)+n;
    78   string Syzstring="Sres";
     2003  string Syzstring="noCE";
    792004  int onlybounds=1;
     2005  int diffforms;
    802006  for (i=1; i<=size(#); i++)
    81   {
    82     if (#[i]=="Vdres")
    83     {
    84       Syzstring="Vdres";
    85     }
    86     if (#[i]=="noiterativeloc")
    87     {
    88       recursiveloc=0;
    89     }
    90     if (#[i]=="exactroots")
    91     {
    92       onlybounds=0;
    93     }
    94   }
     2007    {
     2008      if (#[i]=="Sres")
     2009        {
     2010          Syzstring="Sres";
     2011        }
     2012      if (#[i]=="Vdres")
     2013        {
     2014          Syzstring="Vdres";
     2015        }
     2016      if (#[i]=="noiterativeloc")
     2017        {
     2018          recursiveloc=0;
     2019        }
     2020      if (#[i]=="exactroots")
     2021        {
     2022          onlybounds=0;
     2023        }
     2024      if (#[i]=="diffforms")
     2025        {
     2026          diffforms=1;
     2027        }
     2028    }
    952029  for (i=1; i<=size(L); i++)
    96   {
    97     if (L[i]==0)
    98     {
    99       L=delete(L,i);
    100       i=i-1;
    101     }
    102   }
     2030    {
     2031      if (L[i]==0)
     2032        {
     2033          L=delete(L,i);
     2034          i=i-1;
     2035        }
     2036    }
    1032037  if (size(L)==0)
    104   {
    105     return (list(0));
    106   }
     2038    {
     2039      return (list(0));//////////////////////////////////////////////////////////////////stimmt das jetzt?!??????????????????????????????????
     2040    }
    1072041  for (i=1; i<= size(L); i++)
    108   {
    109     if (leadcoef(L[i])-L[i]==0)
    110     {
    111       return(list(1));   
    112     }
    113   }
     2042    {
     2043      if (leadcoef(L[i])-L[i]==0)
     2044        {
     2045          return(list(1));    ///////////////////////////////////////////////////////////////stimmt das jetzt?!????????????????????????????????????
     2046        }
     2047    }
    1142048  if (size(L)==0)
    115   {
    116     /*the complement of the variety given by the input is the whole space*/
    117     return(list(1));
    118   }
     2049    {
     2050      /*the complement of the variety given by the input is the whole space*/
     2051      return(list(1));
     2052    }
    1192053  for (i=1; i<=size(L); i++)
    120   {
    121     if (typeof(L[i])!="poly")
    122     {
    123       print("The input list must consist of polynomials");
    124       retrun();
    125     }
    126   }
     2054    {
     2055      if (typeof(L[i])!="poly")
     2056        {
     2057          print("The input list must consist of polynomials");
     2058          return();
     2059        }
     2060    }
     2061  if (size(L)==1 and Syzstring=="noCE")
     2062    {
     2063      Syzstring="Sres";
     2064    }
    1272065  /* 1st step: compute the Mayer-Vietoris Complex and its Fourier transform*/
    1282066  def W=MVComplex(L,recursiveloc);//new ring that contains the MV complex
    1292067  setring W;
    1302068  list fortoVdstrict=MV;
    131   ideal IFourier=var(n+1);
    132   for (i=2;i<=n;i++)
    133   {
    134     IFourier=IFourier,var(n+i);
    135   }
    136   for (i=1; i<=n;i++)
    137   {
    138     IFourier=IFourier,-var(i);
    139   }
    140   map cFourier=W,IFourier;
    141   matrix sup;
    142   for (i=1; i<=size(MV); i++)
    143   {
    144     sup=fortoVdstrict[i];
    145     /*takes the Fourier transform of the MV complex*/
    146     fortoVdstrict[i]=cFourier(sup);
    147   }
    148   /* 2nd step: Compute a V_d-strict free complex that is quasi-isomorphic to the
    149   complex fortoVdstrict
    150   The 1st entry of the list rem will be the quasi-isomorphic complex, the 2nd
    151   entry contains the cohomology modules and is needed for the computation of the
    152   global b-function*/
    153   list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring);
     2069  if (diffforms==0)
     2070    {
     2071      ideal IFourier=var(n+1);
     2072      for (i=2;i<=n;i++)
     2073        {
     2074          IFourier=IFourier,var(n+i);
     2075        }
     2076      for (i=1; i<=n;i++)
     2077        {
     2078          IFourier=IFourier,-var(i);
     2079        }
     2080      map cFourier=W,IFourier;
     2081      matrix sup;
     2082      for (i=1; i<=size(MV); i++)
     2083        {
     2084          sup=fortoVdstrict[i];
     2085          /*takes the Fourier transform of the MV complex*/
     2086          fortoVdstrict[i]=cFourier(sup);
     2087        }
     2088    }
     2089  /* 2nd step: Compute a V_d-strict free complex that is quasi-isomorphic to the
     2090     complex fortoVdstrict
     2091     The 1st entry of the list rem will be the quasi-isomorphic complex, the 2nd
     2092     entry contains the cohomology modules and is needed for the computation of the
     2093     global b-function*/
     2094  if (Syzstring=="noCE")
     2095    {
     2096      list rem=quasiisomorphicVdComplex(fortoVdstrict,diffforms);
     2097      list quasiiso=rem[3];
     2098    }
     2099  else
     2100    {
     2101      list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring,diffforms);
     2102      if (diffforms==1)
     2103        {
     2104          list quasiiso=list(matrix(1,1,1));
     2105        }
     2106    }
    1542107  list newcomplex=rem[1];
    155   ////////////////////////////////////////////////////////////////////////////////////
    156   /* 3rd step: Compute the  bounds for the minimal and maximal integer root of the 
    157   global b-function of newcomplex(i.e. compute the lcm of the b-functions of its
    158   cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal
    159   integer root.
    160 
    161   If we compute only the bounds, we omit additional Groebner basis computations.
    162   However this leads to a higher-dimensional truncated complex.
    163 
    164   Note that the  cohomology modules are already contained in rem[2].
    165   minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/
    166   if (onlybounds==1)
    167   {
    168     list minmaxk=globalBFun(rem[2],Syzstring);
    169   }
     2108////////////////////////////////////////////////////////////////////////////////////
     2109  /* 3rd step: Compute the  bounds for the minimal and maximal integer root of the
     2110     global b-function of newcomplex(i.e. compute the lcm of the b-functions of its
     2111     cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal
     2112     integer root.
     2113
     2114     If we compute only the bounds, we omit additional Groebner basis computations.
     2115     However this leads to a higher-dimensional truncated complex.
     2116
     2117     Note that the  cohomology modules are already contained in rem[2].
     2118     minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/
     2119  if (diffforms==1)
     2120    {
     2121      list minmaxk=exactGlobalBFunIntegration(rem[2]);
     2122    }
    1702123  else
    171   {
    172     list minmaxk=exactGlobalBFun(rem[2],Syzstring);
    173   }
     2124    {
     2125      if (onlybounds==1)
     2126        {
     2127          list minmaxk=globalBFun(rem[2],Syzstring);
     2128        }
     2129      else
     2130        {
     2131          list minmaxk=exactGlobalBFun(rem[2],Syzstring);
     2132        }
     2133    }
    1742134  if (size(minmaxk)==0)
    175   {
    176     return (0);
    177   }
    178   /*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where
    179   C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i-2],
    180   m^i=newcomplex[3*i-1], d^i=newcomplex[3*i]), using Thm 5.7 in [W1]:
    181   The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set
    182   (0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and
    183   if it is placed in component k it holds that
    184   minmaxk[1]-m^i[k]<=deg(P_(i_j))<=minmaxk[2]-m^i[k]*/
     2135    {
     2136      return (0);
     2137    }
     2138  ///////////////////////////////////////////////////////////////////////////Bis hierhin angepasst
     2139  /*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where
     2140    C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i-2],
     2141    m^i=newcomplex[3*i-1], d^i=newcomplex[3*i]), using Thm 5.7 in [W1]:
     2142    The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set
     2143    (0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and
     2144    if it is placed in component k it holds that
     2145    minmaxk[1]-m^i[k]<=deg(P_(i_j))<=minmaxk[2]-m^i[k]*/
    1852146  int k,l;
    1862147  list truncatedcomplex,shorten,upto;
    1872148  for (i=1; i<=size(newcomplex) div 3; i++)
    188   {
    189     shorten[3*i-1]=list();
    190     for (j=1; j<=size(newcomplex[3*i-1]); j++)
    191     {
    192       /*shorten[3*i-1][j][k]=minmaxk[k]-m^i[j]+1 (for k=1,2) if this value is
    193       positive otherwise we will set it to be list();
    194       we added +1, because we will use a list, where we put in position l
    195       polys of degree l+1*/
    196       shorten[3*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1);
    197       shorten[3*i-1][j][2]=minmaxk[2]-newcomplex[3*i-1][j]+1;
    198       upto[size(upto)+1]=shorten[3*i-1][j][2];
    199       if (shorten[3*i-1][j][2]<=0)
    200       {
    201         shorten[3*i-1][j]=list();
    202       }
    203       else
    204       {
    205         if (shorten[3*i-1][j][1]<=0)
    206         {
    207           shorten[3*i-1][j][1]=1;
    208         }
    209       }
    210     }
    211   }
     2149    {
     2150      shorten[3*i-1]=list();
     2151      for (j=1; j<=size(newcomplex[3*i-1]); j++)
     2152        {
     2153          /*shorten[3*i-1][j][k]=minmaxk[k]-m^i[j]+1 (for k=1,2) if this value is
     2154            positive otherwise we will set it to be list();
     2155.-            we added +1, because we will use a list, where we put in position l
     2156            polys of degree l+1*/
     2157          shorten[3*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1);
     2158          if (diffforms==1)
     2159            {
     2160              shorten[3*i-1][j][1]=1;
     2161            }
     2162          shorten[3*i-1][j][2]=minmaxk[2]-newcomplex[3*i-1][j]+1;
     2163          upto[size(upto)+1]=shorten[3*i-1][j][2];
     2164          if (shorten[3*i-1][j][2]<=0)
     2165            {
     2166              shorten[3*i-1][j]=list();
     2167            }
     2168          else
     2169            {
     2170              if (shorten[3*i-1][j][1]<=0)
     2171                {
     2172                  shorten[3*i-1][j][1]=1;
     2173                }
     2174            }
     2175        }
     2176    }
    2122177  int iupto=Max(upto);//maximal degree +1 of the polynomials we have to consider
    2132178  if (iupto<=0)
    214   {
    215     return(list(0));
    216   }
     2179    {
     2180      return(list(0));
     2181    }
    2172182  list allpolys;
    2182183  /*allpolys[i] will consist list of all monomials in D(1),...,D(n) of degree i-1*/
    2192184  allpolys[1]=list(1);
    2202185  list minvar;
     2186  list keepv;
    2212187  minvar[1]=list(1);
    2222188  for (i=1; i<=iupto-1; i++)
    223   {
    224     allpolys[i+1]=list();
    225     minvar[i+1]=list();
    226     for (k=1; k<=size(allpolys[i]); k++)
    227     {
    228       for (j=minvar[i][k]; j<=nvars(W) div 2; j++)
    229       {
    230         allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j);
    231         minvar[i+1][size(minvar[i+1])+1]=j;
    232       }
    233     }
    234   }
     2189    {
     2190      allpolys[i+1]=list();
     2191      minvar[i+1]=list();
     2192      for (k=1; k<=size(allpolys[i]); k++)
     2193        {
     2194          for (j=minvar[i][k]; j<=nvars(W) div 2; j++)
     2195            {
     2196              if (diffforms==0)
     2197                {
     2198                  allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j);
     2199                }
     2200              else
     2201                {
     2202                  allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*x(j);
     2203                }
     2204              minvar[i+1][size(minvar[i+1])+1]=j;
     2205            }
     2206        }
     2207    }
    2352208  list keepformatrix,sizetruncom,fortrun,fst;
    2362209  int count,stc;
    2372210  intvec v,forin;
    2382211  matrix subm;
     2212  list keepcount;
     2213  list passendespoly;
    2392214  /*now we compute the truncation*/
    2402215  for (i=1; i<=size(newcomplex) div 3; i++)
    241   {
    242     /*truncatedcomplex[2*i-1] will contain all the generators for the truncation
    243     of D_n/(x(1),..,x(n))\otimes C[i]*/
    244     truncatedcomplex[2*i-1]=list();
    245     sizetruncom[2*i-1]=list();
    246     sizetruncom[2*i]=list();
    247     /*truncatedcomplex[2*i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i])
    248     ->trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/
    249     truncatedcomplex[2*i]=newcomplex[3*i];
    250     v=0;count=0;
    251     sizetruncom[2*i][1]=0;
    252     for (j=1; j<=newcomplex[3*i-2]; j++)
    253     {
    254       if (size(shorten[3*i-1][j])!=0)
    255       {
    256         fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]);
    257         truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1];
    258         count=count+fortrun[2];
    259         fst=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1);
    260         sizetruncom[2*i-1][size(sizetruncom[2*i-1])+1]=fst;
    261         sizetruncom[2*i][size(sizetruncom[2*i])+1]=count;
    262         if (v!=0)
    263         {
    264           v[size(v)+1]=j;
    265         }
    266         else
    267         {
    268           v[1]=j;
    269         }
    270       }
    271     }
    272     if (v!=0)
    273     {
    274       subm=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i]));
    275       truncatedcomplex[2*i]=subm;
    276       if (i!=1)
    277       {
    278         i1=1..nrows(truncatedcomplex[2*(i-1)]);
    279         subm=submat(truncatedcomplex[2*(i-1)],i1,v);
    280         truncatedcomplex[2*(i-1)]=subm;
    281       }
    282     }
    283     else
    284     {
    285       truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i]));
    286       if (i!=1)
    287       {
    288         nr=nrows(truncatedcomplex[2*(i-1)]);
    289         truncatedcomplex[2*(i-1)]=matrix(0,nr,1);
    290       }
    291     }         
    292   }
     2216    {
     2217      /*truncatedcomplex[2*i-1] will contain all the generators for the truncation
     2218        of D_n/(x(1),..,x(n))\otimes C[i]*/
     2219      truncatedcomplex[2*i-1]=list();
     2220      sizetruncom[2*i-1]=list();
     2221      sizetruncom[2*i]=list();
     2222      passendespoly[i]=list();
     2223      /*truncatedcomplex[2*i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i])
     2224        ->trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/
     2225      truncatedcomplex[2*i]=newcomplex[3*i];
     2226      v=0;count=0;
     2227      sizetruncom[2*i][1]=0;
     2228      for (j=1; j<=newcomplex[3*i-2]; j++)
     2229        {
     2230          if (size(shorten[3*i-1][j])!=0)
     2231            {
     2232              fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]);
     2233              truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1];
     2234              for (k=1; k<=size(fortrun[1]); k++)
     2235                {
     2236                  for (l=1; l<=size(fortrun[1][k]); l++)
     2237                    {
     2238                      passendespoly[i][size(passendespoly[i])+1]=list(fortrun[1][k][l][1],j);
     2239                    }
     2240                }
     2241              count=count+fortrun[2];
     2242              fst=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1);
     2243              sizetruncom[2*i-1][size(sizetruncom[2*i-1])+1]=fst;
     2244              sizetruncom[2*i][size(sizetruncom[2*i])+1]=count;
     2245              if (v!=0)
     2246                {
     2247                  v[size(v)+1]=j;
     2248                }
     2249              else
     2250                {
     2251                  v[1]=j;
     2252                }
     2253            }
     2254        }
     2255      if (v!=0)
     2256        {
     2257          keepv[i]=v;
     2258          subm=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i]));
     2259          truncatedcomplex[2*i]=subm;
     2260          if (i!=1)
     2261            {
     2262              i1=1..nrows(truncatedcomplex[2*(i-1)]);
     2263              subm=submat(truncatedcomplex[2*(i-1)],i1,v);
     2264              truncatedcomplex[2*(i-1)]=subm;
     2265            }
     2266        }
     2267      else
     2268        {
     2269          keepv[i]=list();
     2270          truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i]));
     2271          if (i!=1)
     2272            {
     2273              nr=nrows(truncatedcomplex[2*(i-1)]);
     2274              truncatedcomplex[2*(i-1)]=matrix(0,nr,1);
     2275            }
     2276        }
     2277    }
     2278  list keeptruncatedcomplex=truncatedcomplex;
    2932279  matrix M;
    2942280  int st,pi,pj;
    2952281  poly ptc;
    2962282  int b,d,ideg,kplus,lplus;
     2283  int z;
    2972284  poly form,lform,nform;
    2982285  /*computation of the maps*/
     2286  if (diffforms==1)
     2287    {
     2288      def ConvWeyl=makeConverseWeyl(nvars(basering) div 2);
     2289      setring ConvWeyl;
     2290      poly form,lform,nform;
     2291      poly ptc;
     2292      list truncatedcomplex;
     2293      matrix M;
     2294      ideal I=x(1);
     2295      for (i=2; i<=nvars(basering) div 2; i++)
     2296        {
     2297          I=I,var(nvars(basering) div 2 + i);
     2298        }
     2299      for (i=1; i<=nvars(basering) div 2; i++)
     2300        {
     2301          I=I,var(i);
     2302        }
     2303      map transtc=W,I;
     2304      truncatedcomplex=transtc(truncatedcomplex);
     2305    }
    2992306  for (i=1; i<size(truncatedcomplex) div 2; i++)
    300   {
    301     nr=max(1,sizetruncom[2*i][size(sizetruncom[2*i])]);
    302     nc=max(1,sizetruncom[2*i+2][size(sizetruncom[2*i+2])]);
    303     M=matrix(0,nr,nc);   
    304     for (k=1; k<=size(truncatedcomplex[2*i-1]);k++)
    305     {
    306       for (l=1; l<=size(truncatedcomplex[2*(i+1)-1]); l++)
    307       {
    308         if (size(sizetruncom[2*i])!=1)
    309         {
    310           for (j=1; j<=size(truncatedcomplex[2*i-1][k]); j++)     
    311           {
    312             for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++)
    313             {
    314               form=truncatedcomplex[2*i-1][k][j][b][1];
    315               form=form*truncatedcomplex[2*i][k,l];
    316               while (form!=0)
    317               {                 
    318                 lform=lead(form);
    319                 v=leadexp(lform);
    320                 v=v[1..n];                   
    321                 if (v==(0:n))
    322                 {                         
    323                   ideg=deg(lform)-sizetruncom[2*(i+1)-1][l][1];                     
    324                   if (ideg>=0)
    325                   {         
    326                     nr=ideg+1;
    327                     st=size(truncatedcomplex[2*(i+1)-1][l][nr]);
    328                     for (d=1; d<=st;d++)
     2307    {
     2308      nr=max(1,sizetruncom[2*i][size(sizetruncom[2*i])]);
     2309      nc=max(1,sizetruncom[2*i+2][size(sizetruncom[2*i+2])]);
     2310      M=matrix(0,nr,nc);
     2311      for (k=1; k<=size(truncatedcomplex[2*i-1]);k++)
     2312        {
     2313          for (l=1; l<=size(truncatedcomplex[2*(i+1)-1]); l++)
     2314            {
     2315              if (size(sizetruncom[2*i])!=1)
     2316                {
     2317                  for (j=1; j<=size(truncatedcomplex[2*i-1][k]); j++)
    3292318                    {
    330                       nc=2*(i+1)-1;
    331                       ptc=truncatedcomplex[nc][l][ideg+1][d][1];
    332                       if (leadmonom(lform)==ptc)
    333                       {
    334                         nr=2*i-1;
    335                         pi=truncatedcomplex[nr][k][j][b][2];
    336                         pi=pi+sizetruncom[2*i][k];
    337                         nc=2*(i+1)-1;
    338                         nr=ideg+1;
    339                         pj=truncatedcomplex[nc][l][nr][d][2];
    340                         pj=pj+sizetruncom[2*(i+1)][l];
    341                         M[pi,pj]=leadcoef(lform);
    342                         break;
    343                       }
     2319                      for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++)
     2320                        {
     2321                          form=truncatedcomplex[2*i-1][k][j][b][1];
     2322                          form=form*truncatedcomplex[2*i][k,l];
     2323
     2324
     2325                          for (z=1; z<=nvars(basering) div 2; z++)//neu
     2326                            {//
     2327                              form=subst(form,var(z),0);//
     2328                            }//
     2329
     2330                          while (form!=0)
     2331                            {
     2332                              lform=lead(form);
     2333                              v=leadexp(lform);
     2334                              v=v[1..n];
     2335                              // if (v==(0:n))
     2336                              //{
     2337                                  ideg=deg(lform)-sizetruncom[2*(i+1)-1][l][1];
     2338                                  if (ideg>=0)
     2339                                    {
     2340                                      nr=ideg+1;
     2341                                      st=size(truncatedcomplex[2*(i+1)-1][l][nr]);
     2342                                      for (d=1; d<=st;d++)
     2343                                        {
     2344                                          nc=2*(i+1)-1;
     2345                                          ptc=truncatedcomplex[nc][l][ideg+1][d][1];
     2346                                          if (leadmonom(lform)==ptc)
     2347                                            {
     2348                                              nr=2*i-1;
     2349                                              pi=truncatedcomplex[nr][k][j][b][2];
     2350                                              pi=pi+sizetruncom[2*i][k];
     2351                                              nc=2*(i+1)-1;
     2352                                              nr=ideg+1;
     2353                                              pj=truncatedcomplex[nc][l][nr][d][2];
     2354                                              pj=pj+sizetruncom[2*(i+1)][l];
     2355                                              M[pi,pj]=leadcoef(lform);
     2356                                              break;
     2357                                            }
     2358                                        }
     2359                                    }
     2360                                  //        }
     2361
     2362                              form=form-lform;
     2363                            }
     2364                        }
    3442365                    }
    345                   }
    346                 }                 
    347                 form=form-lform;
    348               }
    349             }
    350           }
    351         }
    352       }
    353     }   
    354     truncatedcomplex[2*i]=M;
    355     truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])];
    356   }
     2366                }
     2367            }
     2368        }
     2369      truncatedcomplex[2*i]=M;
     2370      truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])];
     2371    }
    3572372  truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])];
    3582373  if (truncatedcomplex[2*i-1]!=0)
    359   {
    360     truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1);
     2374    {
     2375      truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1);
     2376    }
     2377  if (diffforms==1)
     2378    {
     2379      setring W;
     2380    truncatedcomplex=imap(ConvWeyl,truncatedcomplex);
    3612381  }
    3622382  setring R;
    363   list truncatedcomplex=imap(W,truncatedcomplex);
    364   /*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex,
     2383 list truncatedcomplex=imap(W,truncatedcomplex);
     2384/*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex,
    3652385  i.e. D^i=C^truncatedcomplex[2*i-1] and d^i=truncatedcomplex[2*i]*/
    366   list derhamhom=findCohomology(truncatedcomplex,le);
    367   option(set,saveoptions);
    368   return (derhamhom);
     2386 if (diffforms==0)
     2387   {
     2388     list derhamhom=findCohomology(truncatedcomplex,le);
     2389     option(set,saveoptions);
     2390     return (derhamhom);
     2391   }
     2392 list outall=findCohomologyDiffForms(truncatedcomplex,le);
     2393 setring W;
     2394 list dimanddiff=imap(R,outall);
     2395 list alldiffforms=dimanddiff[2];
     2396 while(size(alldiffforms)<size(passendespoly))
     2397   {
     2398     passendespoly=delete(passendespoly,1);
     2399   }
     2400 list newdiffforms;
     2401 matrix Diff;
     2402 for (i=1; i<=size(alldiffforms); i++)
     2403   {
     2404     newdiffforms[i]=list();
     2405     for (j=1; j<=size(alldiffforms[i]); j++)
     2406       {
     2407         Diff=matrix(0,1,newcomplex[3*(i+size(newcomplex) div 3 - size(alldiffforms))-2]);
     2408         for (k=1; k<=ncols(alldiffforms[i][j]); k++)
     2409           {
     2410             if (alldiffforms[i][j][1,k]!=0)
     2411               {
     2412                 Diff[1,passendespoly[i][k][2]]=Diff[1,passendespoly[i][k][2]]+alldiffforms[i][j][1,k]*passendespoly[i][k][1];
     2413               }
     2414           }
     2415         newdiffforms[i][j]=Diff;
     2416       }
     2417   }
     2418 list omegacomplex=makeOmega(nvars(W) div 2);
     2419 list newcomplexmod;
     2420 for (i=1; i<=size(newcomplex) div 3; i++)
     2421   {
     2422     newcomplexmod[2*i-1]=newcomplex[3*i-2];
     2423     newcomplexmod[2*i]=newcomplex[3*i];
     2424   }
     2425 while (size(dimanddiff[1])<size(newcomplexmod) div 2)
     2426   {
     2427     newcomplexmod=delete(newcomplexmod,1);
     2428     newcomplexmod=delete(newcomplexmod,1);
     2429   }
     2430 while (size(dimanddiff[1])<size(quasiiso))
     2431   {
     2432     quasiiso=delete(quasiiso,1);
     2433   }
     2434 while (size(dimanddiff[1])>size(generators))
     2435   {
     2436     generators=insert(generators,list());
     2437   }
     2438 while (size(dimanddiff[1])>size(quasiiso))
     2439   {
     2440     quasiiso=insert(quasiiso,list());
     2441   }
     2442 int keepsign;
     2443 list derhamdiff;
     2444 list doublecom=makeDoubleComplex(newcomplexmod,omegacomplex,quasiiso,generators);
     2445 matrix diffform;
     2446 int stopping;
     2447 int p;
     2448 matrix convert;
     2449 list interim;
     2450 list correspondingposition;
     2451 list allforms=list();
     2452 for (i=1; i<=size(newdiffforms); i++)
     2453   {
     2454     derhamdiff[i]=list();
     2455     allforms[i]=list();
     2456     for (j=1; j<=size(newdiffforms[i]); j++)
     2457       {
     2458         allforms[i][j]=list();
     2459         keepsign=1;
     2460         derhamdiff[i][j]=0;
     2461         diffform=newdiffforms[i][j];//Zeilenform
     2462         correspondingposition=doublecom[i][1];//needed fpr transformation process
     2463         interim=transferDiffforms(diffform,correspondingposition);
     2464         if (size(interim)!=0)
     2465           {
     2466             allforms[i][j][size(allforms[i][j])+1]=interim;
     2467           }
     2468         stopping=0;
     2469         p=1;
     2470         for (k=i; k<=size(newdiffforms); k++)
     2471           {
     2472             keepsign=(-1)*keepsign;
     2473             if (stopping==0)
     2474               {
     2475                 if (size(doublecom[k][p][2])==0)
     2476                   {
     2477                     stopping=1;
     2478                   }
     2479                 else
     2480                   {
     2481                     if (size(doublecom[k+1][p][3])!=0)
     2482                       {
     2483                         diffform=diffform*doublecom[k][p][2];//Spaltenform
     2484                         if (diffform!=matrix(0,nrows(diffform),ncols(diffform)))
     2485                           {
     2486                              diffform=findPreimage(doublecom[k+1][p][3],transpose(diffform));//Zeilenform
     2487                             correspondingposition=doublecom[k+1][p+1];//needed for transformation process
     2488                             interim=transferDiffforms(keepsign*diffform,correspondingposition);
     2489                             if (size(interim)!=0)
     2490                               {
     2491                                 allforms[i][j][size(allforms[i][j])+1]=interim;
     2492                               }
     2493                             p=p+1;
     2494                           }
     2495                         else
     2496                           {
     2497                             stopping=1;
     2498                           }
     2499                       }
     2500                     else
     2501                       {
     2502                         stopping=1;
     2503                       }
     2504                   }
     2505               }
     2506           }
     2507       }
     2508   }
     2509 setring R;
     2510 list allforms=fetch(W,allforms);
     2511 option(set,saveoptions);
     2512 return (allforms);
    3692513}
    3702514
    … …  
    3852529        -L is a list of non-constant polynomials
    3862530RETURN: ring W: the nth Weyl algebra @*
    387         W contains a list MV, which represents the Mayer-Vietrois complex (C^i,d^i) of the 
     2531        W contains a list MV, which represents the Mayer-Vietrois complex (C^i,d^i) of the
    3882532        polynomials contained in L as follows:@*
    3892533        the C^i are given  by D_n^ncols(C[2*i-1])/im(C[2*i-1]) and the differentials
    … …  
    3932537{
    3942538  /* We follow algorithm 3.2.5 in [R],if #!=0 we use also  Remark 3.2.6 in [R] for
    395   an additional iterative localization*/
     2539     an additional iterative localization*/
    3962540  def R=basering;
    3972541  int i;
    3982542  int iterative=1;
    3992543  if (size(#)!=0)
    400   {
    401     iterative=#[1];
    402   }
     2544    {
     2545      iterative=#[1];
     2546    }
    4032547  for (i=1; i<=size(L); i++)
    404   {
    405     if (L[i]==0)
    406     {
    407       print("localization with respect to 0 not possible");
    408       return();
    409     }
    410     if (leadcoef(L[i])-L[i]==0)
    411     {
    412       print("polynomials must be non-constant");
    413       return();
    414     }
    415   }
     2548    {
     2549      if (L[i]==0)
     2550        {
     2551          print("localization with respect to 0 not possible");
     2552          return();
     2553        }
     2554      if (leadcoef(L[i])-L[i]==0)
     2555        {
     2556          print("polynomials must be non-constant");
     2557          return();
     2558        }
     2559    }
    4162560  if (iterative==1)
    417   {
    418     /*compute the localizations by factorizing the polynomials and iterative
    419     localization of the factors */
    420     for (i=1; i<=size(L); i++)
    421     {
    422       L[i]=factorize(L[i],1);
    423     }
    424   }
     2561    {
     2562      /*compute the localizations by factorizing the polynomials and iterative
     2563        localization of the factors */
     2564      for (i=1; i<=size(L); i++)
     2565        {
     2566          L[i]=factorize(L[i],1);
     2567        }
     2568    }
    4252569  int r=size(L);
    4262570  int n=nvars(basering);
    … …  
    4282572  /*construct the ring Ws*/
    4292573  def W=makeWeyl(n);
    430   setring W; 
     2574  setring W;
    4312575  list man=ringlist(W);
    4322576  if (n==1)
    433   {
    434     man[2][1]="x(1)";
    435     man[2][2]="D(1)";
    436     def Wi=ring(man);
    437     setring Wi;
    438     kill W;
    439     def W=Wi;
    440     setring W;
    441     list man=ringlist(W);
    442   }
     2577    {
     2578      man[2][1]="x(1)";
     2579      man[2][2]="D(1)";
     2580      def Wi=ring(man);
     2581      setring Wi;
     2582      kill W;
     2583      def W=Wi;
     2584      setring W;
     2585      list man=ringlist(W);
     2586    }
    4432587  man[2][size(man[2])+1]="s";;
    4442588  man[3][3]=man[3][2];
    … …  
    4482592  matrix M[1][1];
    4492593  man[6]=transpose(concat(transpose(concat(man[6],M)),M));
    450   def Ws=ring(man); 
     2594  def Ws=ring(man);
    4512595  setring Ws;
    452   int j,k,l,c; 
    453   list L=fetch(R,L); 
    454   list Cech; 
     2596  int j,k,l,c;
     2597  list L=fetch(R,L);
     2598  list Cech;
    4552599  ideal J=var(1+n);
    4562600  for (i=2; i<=n; i++)
    457   {
    458     J=J,var(i+n);
    459   }
     2601    {
     2602      J=J,var(i+n);
     2603    }
    4602604  Cech[1]=list(J);
    4612605  list Theta, remminroots;
    … …  
    4682612  int rmr;
    4692613  if (iterative==0)
    470   {/*computation of the modules of the MV complex*/
    471     for (i=1; i<=r; i++)
    472     {
    473       findminintroot=list();
    474       Cech[i+1]=list();
    475       Theta[i+1]=list();
    476       k1=1;
    477       for (j=1; j<=i; j++)
    478       {
    479         k1[size(k1)+1]=size(Theta[j+1]);
    480         for (k=1; k<=k1[j]; k++)
    481         {
    482           Theta[j+1][size(Theta[j+1])+1]=list(Theta[j][k][1]+list(i));
    483           Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i];
    484           /*We compute the s-parametric annihilator J(s)  and the b-function
    485           of the polynomial L[i] and Cech[i][k] to localize the module 
    486           D_n/(D(1),...,D(n))[L[i]^(-1)]\otimes D_n^c/im(Cech[i][k]),
    487           where c=ncols(Cech[i][k]) and the im(Cech[i][k]) is generated by
    488           the rows of the matrix.
    489           If we plug the minimal integer root r(or a smaller integer
    490           value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic to
    491           the above localization*/
    492           rem=SannfsIBM(L[i],Cech[j][k]);
    493           Cech[j+1][size(Cech[j+1])+1]=rem[1];
    494           findminintroot[size(findminintroot)+1]=rem[2];
    495         }
    496       }
    497       /* we compute the minimal root of all b-functions of L[i] computed above,
    498       because we want to plug in the same root r in all s-parametric
    499       annihilators we computed for L[i]  ->this will ensure  we can compute
    500       the maps of the MV complex*/
    501       minroot=minIntRoot0(findminintroot);
    502       for (j=1; j<=i; j++)
    503       {
    504         for (k=1; k<=k1[j]; k++)
    505         {
    506           sk=size(Cech[j+1])+1-k;
    507           Cech[j+1][size(Cech[j+1])+1-k]=subst(Cech[j+1][sk],s,minroot);
    508         }
    509       }
    510       remminroots[i]=minroot;
    511     }
    512     Cech=delete(Cech,1);
    513     Theta=delete(Theta,1);
    514     list zw;
    515     poly reme;
    516     /*computation of the maps of the MV complex*/
    517     for (i=1; i<r; i++)
    518     {
    519       diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
    520       for (j=1; j<=size(Cech[i]); j++)
    521       {
    522         for (k=1; k<=size(Cech[i+1]); k++)
    523         {
    524           zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
    525           if (zw[2]!=0)
    526           {
    527             rmr=-remminroots[zw[1]];
    528             reme=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
    529             zw[2]=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
    530             diffmaps[i][j,k]=zw[2];
    531           }
    532         }
    533       }
    534     }
    535     diffmaps[r]=matrix(0,1,1);
    536   }
     2614    {/*computation of the modules of the MV complex*/
     2615      for (i=1; i<=r; i++)
     2616        {
     2617          findminintroot=list();
     2618          Cech[i+1]=list();
     2619          Theta[i+1]=list();
     2620          k1=1;
     2621          for (j=1; j<=i; j++)
     2622            {
     2623              k1[size(k1)+1]=size(Theta[j+1]);
     2624              for (k=1; k<=k1[j]; k++)
     2625                {
     2626                  Theta[j+1][size(Theta[j+1])+1]=list(Theta[j][k][1]+list(i));
     2627                  Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i];
     2628                  /*We compute the s-parametric annihilator J(s)  and the b-function
     2629                    of the polynomial L[i] and Cech[i][k] to localize the module
     2630                    D_n/(D(1),...,D(n))[L[i]^(-1)]\otimes D_n^c/im(Cech[i][k]),
     2631                    where c=ncols(Cech[i][k]) and the im(Cech[i][k]) is generated by
     2632                    the rows of the matrix.
     2633                    If we plug the minimal integer root r(or a smaller integer
     2634                    value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic to
     2635                    the above localization*/
     2636                  rem=SannfsIBM(L[i],Cech[j][k]);
     2637                  Cech[j+1][size(Cech[j+1])+1]=rem[1];
     2638                  findminintroot[size(findminintroot)+1]=rem[2];
     2639                }
     2640            }
     2641          /* we compute the minimal root of all b-functions of L[i] computed above,
     2642             because we want to plug in the same root r in all s-parametric
     2643             annihilators we computed for L[i]  ->this will ensure  we can compute
     2644             the maps of the MV complex*/
     2645          minroot=minIntRoot(findminintroot);
     2646          for (j=1; j<=i; j++)
     2647            {
     2648              for (k=1; k<=k1[j]; k++)
     2649                {
     2650                  sk=size(Cech[j+1])+1-k;
     2651                  Cech[j+1][size(Cech[j+1])+1-k]=subst(Cech[j+1][sk],s,minroot);
     2652                }
     2653            }
     2654          remminroots[i]=minroot;
     2655        }
     2656      Cech=delete(Cech,1);
     2657      Theta=delete(Theta,1);
     2658      list zw;
     2659      poly reme;
     2660      /*computation of the maps of the MV complex*/
     2661      for (i=1; i<r; i++)
     2662        {
     2663          diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
     2664          for (j=1; j<=size(Cech[i]); j++)
     2665            {
     2666              for (k=1; k<=size(Cech[i+1]); k++)
     2667                {
     2668                  zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
     2669                  if (zw[2]!=0)
     2670                    {
     2671                      rmr=-remminroots[zw[1]];
     2672                      reme=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
     2673                      zw[2]=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
     2674                      diffmaps[i][j,k]=zw[2];
     2675                    }
     2676                }
     2677            }
     2678        }
     2679      diffmaps[r]=matrix(0,1,1);
     2680    }
     2681  list generators;
    5372682  if (iterative==1)
    538   {
    539     for (i=1; i<=r;i++)
    540     {
    541       Cech[i+1]=list();
    542       Theta[i+1]=list();
    543       k1=1;
    544       for (c=1; c<=size(L[i]); c++)
    545       {
    546         findminintroot=list();
    547         for (j=1; j<=i; j++)
    548         {
    549           if (c==1)
    550           {
    551             k1[size(k1)+1]=size(Theta[j+1]);
    552           }
    553           for (k=1; k<=k1[j]; k++)
    554           {
    555             /*We compute the s-parametric annihilator J(s)  und the b-
    556             function of the polynomial L[i][c] and Cech[i][k] to
    557             localize the module D_n/(D(1),...,D(n))[L[i][c]^(-1)]\otimes
    558             D_n^c/im(Cech[i][k]), where c=ncols(Cech[i][k]).
    559             If we plug the minimal integer root r(or a smaller integer
    560             value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic
    561             to the above localization*/
    562             if (c==1)
    563             {
    564               rmr=size(Theta[j+1])+1;
    565               Theta[j+1][rmr]=list(Theta[j][k][1]+list(i));
    566               Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i][c];
    567               rem=SannfsIBM(L[i][c],Cech[j][k]);
    568               Cech[j+1][size(Cech[j+1])+1]=rem[1];
    569               findminintroot[size(findminintroot)+1]=rem[2];
    570             }
    571             else
    572             {
    573               st=size(Theta[j+1])-k1[j]+k;
    574               Theta[j+1][st][2]=Theta[j+1][st][2]*L[i][c];
    575               rem=SannfsIBM(L[i][c],Cech[j+1][size(Cech[j+1])-k1[j]+k]);
    576               Cech[j+1][size(Cech[j+1])-k1[j]+k]=rem[1];
    577               findminintroot[size(findminintroot)+1]=rem[2];
    578             }
    579           }
    580         }
    581         /* we compute the minimal root of all b-functions of L[i][c]
    582         computed above,because we want to plug in the same root r in all
    583         s-parametric annihilators we computed for L[i]  ->this will
    584         ensure  we can compute the maps of the MV complex*/
    585         minroot=minIntRoot0(findminintroot);
    586         for (j=1; j<=i; j++)
    587         {
    588           for (k=1; k<=k1[j]; k++)
    589           {
    590             st=size(Cech[j+1])+1-k;
    591             Cech[j+1][st]=subst(Cech[j+1][st],s,minroot);
    592           }
    593         }
    594         if (c==1)
    595         {
    596           remminroots[i]=list();
    597         }
    598         remminroots[i][c]=minroot;
    599       }
    600     }
    601     Cech=delete(Cech,1);
    602     Theta=delete(Theta,1);
    603     list zw;
    604     poly reme;
    605     /*maps of the MV Complex*/
    606     for (i=1; i<r; i++)
    607     {
    608       diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
    609       for (j=1; j<=size(Cech[i]); j++)
    610       {
    611         for (k=1; k<=size(Cech[i+1]); k++)
    612         {
    613           zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
    614           if (zw[2]!=0)
    615           {
    616             reme=1;
    617             for (c=1; c<=size(L[zw[1]]);c++)
    618             {
    619               reme=reme*L[zw[1]][c]^(-remminroots[zw[1]][c]);
    620             }
    621             diffmaps[i][j,k]=zw[2]*reme;
    622           }
    623         }
    624       }
    625     }   
    626     diffmaps[r]=matrix(0,1,1);
    627   }
     2683    {
     2684      for (i=1; i<=r;i++)
     2685        {
     2686          generators[i]=list();////////////////////////////////////////////////////////////////////
     2687          Cech[i+1]=list();
     2688          Theta[i+1]=list();
     2689          k1=1;
     2690          for (c=1; c<=size(L[i]); c++)
     2691            {
     2692              findminintroot=list();
     2693              for (j=1; j<=i; j++)
     2694                {
     2695                  if (c==1)
     2696                    {
     2697                      k1[size(k1)+1]=size(Theta[j+1]);
     2698                    }
     2699                  for (k=1; k<=k1[j]; k++)
     2700                    {
     2701                      /*We compute the s-parametric annihilator J(s)  und the b-
     2702                        function of the polynomial L[i][c] and Cech[i][k] to
     2703                        localize the module D_n/(D(1),...,D(n))[L[i][c]^(-1)]\otimes
     2704                        D_n^c/im(Cech[i][k]), where c=ncols(Cech[i][k]).
     2705                        If we plug the minimal integer root r(or a smaller integer
     2706                        value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic
     2707                        to the above localization*/
     2708                      if (c==1)
     2709                        {
     2710                          rmr=size(Theta[j+1])+1;
     2711                          Theta[j+1][rmr]=list(Theta[j][k][1]+list(i));
     2712                          Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i][c];
     2713                          rem=SannfsIBM(L[i][c],Cech[j][k]);
     2714                          Cech[j+1][size(Cech[j+1])+1]=rem[1];
     2715                          findminintroot[size(findminintroot)+1]=rem[2];
     2716                        }
     2717                      else
     2718                        {
     2719                          st=size(Theta[j+1])-k1[j]+k;
     2720                          Theta[j+1][st][2]=Theta[j+1][st][2]*L[i][c];
     2721                          rem=SannfsIBM(L[i][c],Cech[j+1][size(Cech[j+1])-k1[j]+k]);
     2722                          Cech[j+1][size(Cech[j+1])-k1[j]+k]=rem[1];
     2723                          findminintroot[size(findminintroot)+1]=rem[2];
     2724                        }
     2725                    }
     2726                }
     2727                /* we compute the minimal root of all b-functions of L[i][c]
     2728                   computed above,because we want to plug in the same root r in all
     2729                   s-parametric annihilators we computed for L[i]  ->this will
     2730                   ensure  we can compute the maps of the MV complex*/
     2731              minroot=minIntRoot(findminintroot);
     2732              for (j=1; j<=i; j++)
     2733                {
     2734                  for (k=1; k<=k1[j]; k++)
     2735                    {
     2736                      st=size(Cech[j+1])+1-k;
     2737                      Cech[j+1][st]=subst(Cech[j+1][st],s,minroot);
     2738                    }
     2739                }
     2740              if (c==1)
     2741                {
     2742                  remminroots[i]=list();
     2743                }
     2744              remminroots[i][c]=minroot;
     2745            }
     2746        }
     2747      Cech=delete(Cech,1);
     2748      Theta=delete(Theta,1);
     2749      list zw;
     2750      poly reme;
     2751      /*maps of the MV Complex*/
     2752      for (i=1; i<r; i++)
     2753        {
     2754          diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
     2755          for (j=1; j<=size(Cech[i]); j++)
     2756            {
     2757              for (k=1; k<=size(Cech[i+1]); k++)
     2758                {
     2759                  zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
     2760                  if (zw[2]!=0)
     2761                    {
     2762                      reme=1;
     2763                      for (c=1; c<=size(L[zw[1]]);c++)
     2764                        {
     2765                          reme=reme*L[zw[1]][c]^(-remminroots[zw[1]][c]);
     2766                        }
     2767                      diffmaps[i][j,k]=zw[2]*reme;
     2768                    }
     2769                }
     2770            }
     2771        }
     2772      diffmaps[r]=matrix(0,1,1);
     2773      for (i=1; i<=r; i++)
     2774        {
     2775          for (j=1; j<=size(Theta[i]); j++)
     2776            {
     2777              generators[i][j]=1;
     2778              for (c=1; c<=size(Theta[i][j][1]); c++)
     2779                {
     2780                  for (k=1; k<=size(L[Theta[i][j][1][c]]); k++)
     2781                    {
     2782                      generators[i][j]=generators[i][j]*L[Theta[i][j][1][c]][k]^((-1)*remminroots[Theta[i][j][1][c]][k]);
     2783                    }
     2784                }
     2785            }
     2786        }
     2787    }
    6282788  setring W;
    6292789  /*map the modules and maps to the Weyl algebra*/
    6302790  list diffmaps=imap(Ws,diffmaps);
    6312791  list Cechmodules=imap(Ws,Cech);
     2792  if (iterative==1)
     2793    {
     2794      list Theta=imap(Ws,Theta);
     2795      list generators=imap(Ws,generators);
     2796    }
    6322797  list Cech;
    6332798  matrix sup;
    6342799  for (i=1; i<=r; i++)
    635   {
    636     sup=transpose(matrix(Cechmodules[i][1]));
    637     Cech[2*i-1]=sup;
    638     for (j=2; j<=size(Cechmodules[i]); j++)
    639     {
    640       sup=transpose(matrix(Cechmodules[i][j]));
    641       Cech[2*i-1]=dsum(Cech[2*i-1],sup);
    642     }
    643     sup=matrix(diffmaps[i]);
    644     Cech[2*i]=sup;
    645   }
     2800    {
     2801      sup=transpose(matrix(Cechmodules[i][1]));
     2802      Cech[2*i-1]=sup;
     2803      for (j=2; j<=size(Cechmodules[i]); j++)
     2804        {
     2805          sup=transpose(matrix(Cechmodules[i][j]));
     2806          Cech[2*i-1]=dsum(Cech[2*i-1],sup);
     2807        }
     2808      sup=matrix(diffmaps[i]);
     2809      Cech[2*i]=sup;
     2810    }
    6462811  list MV=Cech;
     2812  if (iterative==1)
     2813    {
     2814      export Theta;
     2815      export generators;
     2816    }
    6472817  export MV;
     2818
    6482819  return (W);
    6492820}
    6502821
    6512822example
    652 { "EXAMPLE:"; 
     2823{ "EXAMPLE:";
    6532824  ring r = 0,(x,y,z),dp;
    6542825  list L=xy,xz;
    … …  
    6622833static proc SannfsIBM(poly F,ideal myJ)
    6632834"USAGE: SannfsIBM(f,J), F poly, J ideal
    664 ASSUME: basering is D_n[s], where D_n is the Weyl algebra and s and extra 
     2835ASSUME: basering is D_n[s], where D_n is the Weyl algebra and s and extra
    6652836        commutative variable@*
    666         f is a polynomial in the variables x(1),...,x(n) with rational coefficients 
     2837        f is a polynomial in the variables x(1),...,x(n) with rational coefficients
    6672838        @*
    668         J is holonomic and f-saturated     
     2839        J is holonomic and f-saturated
    6692840RETURN  AlList of the form (K,g), where K is an ideal and g a univariant polynomial
    6702841        in  the variable s. K is the s-parametric annihilator of F and J and g is
    … …  
    6732844{
    6742845  /*modified version of the procedure SannfsBM from the library dmod.lib: SannfsBM
    675   computes the s-parametric annihilator for J=(x_1,...,x_n)*/
    676   /* We use Algorithm 3.1.12 in[R] to compute the s-parametric 
    677       annihilator. Then we use the s-parametric annihilator to compute the b-function
    678       via Algorithm 4.7 in [W1].*/
     2846    computes the s-parametric annihilator for J=(x_1,...,x_n)*/
     2847  /* We use Algorithm 3.1.12 in[R] to compute the s-parametric
     2848     annihilator. Then we use the s-parametric annihilator to compute the b-function
     2849     via Algorithm 4.7 in [W1].*/
    6792850  /* We assume that the basering the the nth Weyl algebra D_n. We create the ring
    680       D_n[s,t], where t*s=s*t-t*/
     2851     D_n[s,t], where t*s=s*t-t*/
    6812852  def save = basering;
    6822853  int N = nvars(basering)-1;
    … …  
    6882859  list tmp;
    6892860  intvec iv;
    690   L[1] = RL[1]; 
    691   L[4] = RL[4]; 
     2861  L[1] = RL[1];
     2862  L[4] = RL[4];
    6922863  list Name  = RL[2];
    6932864  Name=delete(Name,size(Name));
    … …  
    6962867  RName[2] = "s";
    6972868  list DName;
    698   for(i=1;i<=N div 2;i++)
     2869 for(i=1;i<=N div 2;i++)
    6992870  {
    700     DName[i] = var(N div 2+i); 
     2871    DName[i] = var(N div 2+i);
    7012872    Name=delete(Name,N div 2+1);
    7022873  }
    … …  
    7062877  L[2]   = NName;
    7072878  kill NName;
    708   tmp[1]  = "lp"; 
     2879  tmp[1]  = "lp";
    7092880  iv      = 1,1;
    710   tmp[2]  = iv; 
     2881  tmp[2]  = iv;
    7112882  Lord[1] = tmp;
    712   tmp[1]  = "dp"; 
     2883  tmp[1]  = "dp";
    7132884  s       = "iv=";
    7142885  for(i=1;i<=Nnew;i++)
    … …  
    7422913  ideal myJ=imap(save,myJ);
    7432914  for (i=1; i<=N div 2; i++)
    744   {
    745     myJ=subst(myJ,D(i),D(i)+diff(F,x(i))*t);
    746   }
     2915    {
     2916      myJ=subst(myJ,D(i),D(i)+diff(F,x(i))*t);
     2917    }
    7472918  ideal I = t*F+s;
    7482919  I=I,myJ;//the s-parametric annihilator in D_n[s,t]
    … …  
    7632934
    7642935////////////////////////////////////////////////////////////////////////////////////
    765 //COMPUTATION OF QA UASI-ISOMORPHIC V_D-STRICT FREE COMPLEX
     2936//COMPUTATION OF A QUASI-ISOMORPHIC V_D-STRICT FREE COMPLEX
    7662937////////////////////////////////////////////////////////////////////////////////////
    7672938
     2939static proc quasiisomorphicVdComplex(list L,list #)
     2940"USAGE: quasiisomorphicVdComplex(L[,df]); L a list of the form (M_1,f_1,...,M_s,f_s),
     2941        where the M_i and f_i are matrices
     2942ASSUME: Basering is the Weyl algebra D_n @*
     2943        (M_1,f_1,...,M_s,f_s) represents a complex 0->D_n^(r_1)/im(M_1)->
     2944        D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i)
     2945        is generated by the rows of M_i. In particular it hold:@*
     2946        - The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @*
     2947        -the image of M_1*f_i is contained in the image of M_(i+1) @*
     2948        d is an integer between 1 and n. If no value for d is given, it is assumed
     2949        to be n @*
     2950        df is an optional int, if df equals 1 a \tilde(V_d)-strict complex
     2951        will be computed (instead of a V_d-strict one) (for a definition see [W3])
     2952RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @*
     2953        L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,i_s,g_s,m_s) such that:@*
     2954        -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs
     2955         of size i_j@*
     2956        -D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0  is a V_d-strict complex
     2957         with differentials m_i that is quasi-isomorphic to the complex given by L@*
     2958        L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and
     2959        the n_i are shift vectors such that:@*
     2960        -coker(H_i) is the ith cohomology group of the complex given by L_1@*
     2961        -the n_i are the shift vectors of the coker(H_i)
     2962THEORY: We follow Proposition 3.2 and Corollary 3.3 in [W3]
     2963"
     2964{
     2965  int tilde;
     2966  if (size(#)!=0)
     2967    {
     2968      tilde=#[1];
     2969    }
     2970  def B=basering;
     2971  int n=nvars(B) div 2 + 1;//+1 müsste stimmen! bitte kontrollieren!
     2972  int d=nvars(B) div 2;
     2973  int r=size(L) div 2;
     2974  int lonc=n+r;
     2975  int Kiold=0;
     2976  matrix kerold;
     2977  // matrix kernew=out[r][2][2];
     2978  matrix kernew=diag(1,ncols(L[size(L)-1]));
     2979  module mL;
     2980  int i;
     2981  int k;
     2982  matrix testm;
     2983  int Kinew=nrows(kernew);
     2984  int Jiold=0;
     2985  int Jinew=0;
     2986  matrix Niold;
     2987  matrix Ninew;
     2988  list newcomplex;
     2989  int Aiold=Kinew;
     2990  matrix savediv;
     2991  newcomplex[3*lonc-2]=Kinew;
     2992  newcomplex[3*lonc-1]=intvec(0:Kinew);
     2993  newcomplex[3*lonc]=matrix(0,Kinew,1);
     2994  list quasiiso;
     2995  quasiiso[lonc]=diag(1,Kinew);
     2996  matrix invimage;
     2997  matrix keralpha;
     2998  intvec v;
     2999  int j;
     3000  matrix sc;
     3001  matrix fnc;
     3002  int indk;
     3003  int indj;
     3004  int Aiold;
     3005  list saveres;
     3006  matrix Liplus;
     3007  for (i=r-1; i>=0; i--)
     3008    {
     3009      indk=0;
     3010      indj=0;
     3011      Kiold=Kinew;
     3012      kerold=kernew;
     3013      if (i!=0)
     3014        {
     3015          // kernew=divdr(L[2*i],L[2*i+1],1);
     3016          kernew=divdr(L[2*i],L[2*i+1]);
     3017          mL=slimgb(transpose(L[2*i-1]));
     3018          for (k=1; k<=nrows(kernew); k++)
     3019            {
     3020              testm=reduce(transpose(submat(kernew,k,intvec(1..ncols(kernew)))),mL);
     3021              if (testm==matrix(0,nrows(testm),ncols(testm)))
     3022                {
     3023                  kernew=transpose(deletecol(transpose(kernew),k));
     3024                  k=k-1;
     3025                }
     3026            }
     3027          Kinew=nrows(kernew);
     3028          if (kernew==matrix(0,nrows(kernew),ncols(kernew)))
     3029            {
     3030              Kinew=0;
     3031              indk=1;
     3032            }
     3033        }
     3034      else
     3035        {
     3036          Kinew=0;
     3037          indk=1;
     3038        }
     3039      Jiold=Jinew;
     3040      Niold=Ninew;
     3041      keralpha=transpose(syz(transpose(newcomplex[3*(i+n)+3])));
     3042      if (i!=0)
     3043        {
     3044          invimage=divdr(quasiiso[n+i+1],transpose(concat(transpose(L[2*i]),transpose(L[2*i+1]))));
     3045          Ninew=vdStrictIntersect(keralpha,invimage,newcomplex[3*(n+i+1)-1],tilde);//////////////
     3046        }
     3047      else
     3048        {
     3049          invimage=divdr(quasiiso[n+i+1],L[2*i+1]);
     3050          saveres=vdStrictIntersectPlus(keralpha,invimage,newcomplex[3*(n+i+1)-1],tilde);////////////////////////
     3051
     3052          ///////////////////BIS HIERHIN VERALLGEMEINERT////////////////////////////////////////////////////////////////////
     3053
     3054
     3055          Ninew=saveres[1];
     3056        }
     3057      Jinew=nrows(Ninew);
     3058      if (Ninew==matrix(0,nrows(Ninew),ncols(Ninew)))
     3059        {
     3060          Jinew=0;
     3061          indk=1;
     3062        }
     3063      newcomplex[3*(n+i)-2]=Kinew+Jinew;
     3064      v=0;
     3065      if (indk==0)
     3066        {
     3067          v=(0:Kinew);
     3068          if (indj==0)
     3069            {
     3070              fnc=transpose(concat(transpose(matrix(0,Kinew,Kiold+Jiold)),transpose(Ninew)));
     3071            }
     3072          else
     3073            {
     3074              fnc=matrix(0,Kinew,Kiold+Jiold);
     3075            }
     3076        }
     3077      else
     3078        {
     3079          if (indj==0)
     3080            {
     3081              fnc=Ninew;
     3082            }
     3083          else
     3084            {
     3085              fnc=matrix(0,1,Kiold+Jiold);
     3086              newcomplex[3*(n+i)-2]=1;
     3087            }
     3088        }
     3089      Aiold=Jinew+Kinew;
     3090      if (Aiold==0)
     3091        {
     3092          Aiold=1;
     3093        }
     3094      newcomplex[3*(n+i)]=fnc;
     3095      for (j=1; j<=Jinew; j++)
     3096        {
     3097          if (tilde==0)
     3098            {
     3099              v[Kinew+j]=VdDeg(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]);
     3100            }
     3101          else
     3102            {
     3103              v[Kinew+j]=VdDegTilde(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]);
     3104            }
     3105        }
     3106      newcomplex[3*(n+i)-1]=v;
     3107      if (i==0)
     3108        {
     3109          quasiiso[n+i]=matrix(0,Jinew,1);
     3110        }
     3111      else
     3112        {
     3113          if (indj==0)
     3114            {
     3115              sc=submat(fnc,intvec(Kinew+1..nrows(fnc)),intvec(1..ncols(fnc)))*quasiiso[n+i+1];
     3116              Liplus=transpose(concat(transpose(L[2*i]),transpose(L[2*i+1])));
     3117              sc=matrixLift(Liplus,sc);//stimmt das jetzt
     3118              sc=submat(sc,intvec(1..nrows(sc)),intvec(1..nrows(L[2*i])));
     3119              if (indk==0)
     3120                {
     3121                  //pi=kernew
     3122                  quasiiso[n+i]=transpose(concat(transpose(kernew),transpose(sc)));
     3123                }
     3124              else
     3125                {
     3126                  quasiiso[n+i]=sc;
     3127                }
     3128            }
     3129          else
     3130            {
     3131              if (indk==0)
     3132                {
     3133                  quasiiso[n+i]=kernew;
     3134                }
     3135              else
     3136                {
     3137                  quasiiso[n+i]=matrix(0,1,ncols(kernew));
     3138                }
     3139            }
     3140        }
     3141    }
     3142  for (i=1; i<=n-1; i++)
     3143    {
     3144      quasiiso[n-i]=list();
     3145      if (size(saveres[2][i])!=0)
     3146        {
     3147          newcomplex[3*(n-i)]=saveres[2][i];
     3148          newcomplex[3*(n-i)-2]=nrows(saveres[2][i]);
     3149          v=0;
     3150          for (j=1; j<=newcomplex[3*(n-i)-2]; j++)
     3151            {
     3152              if (tilde==0)
     3153                {
     3154                  v[j]=VdDeg(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(n-i)+2]);
     3155                }
     3156              else
     3157                {
     3158                  v[j]=VdDegTilde(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(n-i)+2]);
     3159                }
     3160            }
     3161          newcomplex[3*(n-i)-1]=v;
     3162        }
     3163      else
     3164        {
     3165          newcomplex[3*(n-i)]=matrix(0,1,1);
     3166          if (newcomplex[3*(n-i)+1]!=0)
     3167            {
     3168              newcomplex[3*(n-i)]=matrix(0,1,newcomplex[3*(n-i)+1]);
     3169            }
     3170          newcomplex[3*(n-i)-2]=int(0);
     3171          newcomplex[3*(n-i)-1]=intvec(0);
     3172        }
     3173    }
     3174  list result;
     3175  result[1]=newcomplex;
     3176  result[2]=list();
     3177  list forsep;
     3178  for (i=1; i<=size(L) div 2+1; i++)
     3179    {
     3180      forsep[2*i]=newcomplex[3*(n+i-1)];
     3181      forsep[2*i-1]=matrix(0,1,nrows(forsep[2*i]));
     3182    }
     3183  forsep=shortExactPieces(forsep);
     3184  list listofHis;
     3185  matrix forVd;
     3186  for (i=1; i<=size(L) div 2; i++)
     3187    {
     3188      v=0;
     3189      listofHis[i]=list(forsep[i+1][1][5]);
     3190      forVd=forsep[i+1][2][2];
     3191      for (j=1; j<=nrows(forVd); j++)
     3192        {
     3193          if (tilde==0)
     3194            {
     3195              v[j]=VdDeg(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)-1]);
     3196            }
     3197          else
     3198            {
     3199              v[j]=VdDegTilde(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)-1]);
     3200            }
     3201        }
     3202      listofHis[i][2]=v;
     3203    }
     3204  result[2]=listofHis;
     3205  result[3]=quasiiso;
     3206  return(result);
     3207}
     3208
     3209////////////////////////////////////////////////////////////////////////////////////
     3210
     3211static proc vdStrictIntersect(matrix M, matrix N, intvec v, int tilde)
     3212{
     3213  def B=basering;
     3214  option(returnSB);//                    alternative:erst intersect und dann SB-Berechung mit slimgb
     3215  if (tilde==0)
     3216    {
     3217      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v);
     3218    }
     3219  else
     3220    {
     3221      def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v);
     3222    }
     3223  setring HomWeyl;
     3224  matrix M=fetch(B,M);
     3225  matrix N=fetch(B,N);
     3226  M=nHomogenize(M);
     3227  N=nHomogenize(N);
     3228  matrix vdintersection=transpose(intersect(transpose(M),transpose(N)));
     3229  vdintersection=subst(vdintersection,h,1);
     3230  setring B;
     3231  matrix vdintersection=fetch(HomWeyl,vdintersection);
     3232  option(noreturnSB);
     3233  return(vdintersection);
     3234}
     3235
     3236////////////////////////////////////////////////////////////////////////////////////
     3237
     3238static proc vdStrictIntersectPlus(matrix M, matrix N, intvec v, int tilde)
     3239{
     3240  def B=basering;
     3241  int n=nvars(B) div 2;
     3242  matrix vdint=transpose(intersect(transpose(M),transpose(N)));
     3243  if (tilde==0)
     3244    {
     3245      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v);
     3246    }
     3247  else
     3248    {
     3249      def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v);
     3250    }
     3251  setring HomWeyl;
     3252  matrix vdint=fetch(B,vdint);
     3253  matrix N=fetch(B,N);
     3254  vdint=nHomogenize(vdint);
     3255  intvec i1;
     3256  intvec i2;
     3257  int i;
     3258  int nr;
     3259  int nc;
     3260  def ringofSyz=Sres(transpose(vdint),n);////////////////////////////////////////////////////////////////
     3261  setring ringofSyz;
     3262  matrix vdint=transpose(matrix(RES[2]));
     3263  vdint=subst(vdint,h,1);
     3264  int logens=ncols(vdint)+1;
     3265  int omitemptylist;
     3266  matrix zerom;
     3267  list rofA;
     3268  for (i=3; i<=n+3; i++)////////////////////////////////////////////////////////////////////////////n und si müssen noch definiert werden
     3269    {
     3270      if (size(RES)>=i)
     3271        {
     3272          zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
     3273          if (RES[i]!=zerom)
     3274            {
     3275              rofA[i-2]=(matrix(RES[i]));
     3276              if (i==3)
     3277                {
     3278                  if (nrows(rofA[i-2])-logens+1!=nrows(vdint))
     3279                    {
     3280                      //build the resolution
     3281                      nr=nrows(vdint)+logens-1;
     3282                      nc=ncols(rofA[i-2]);
     3283                      rofA[i-2]=matrix(rofA[i-2],nr,nc);
     3284                    }
     3285
     3286                }
     3287              if (i!=3)
     3288                {
     3289                  if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3]))
     3290                    {
     3291                      nr=nrows(rofA[i-3])+logens-1;
     3292                      nc=ncols(rofA[i-2]);
     3293                      rofA[i-2]=matrix(rofA[i-2],nr,nc);
     3294                    }
     3295                }
     3296              i1=intvec(logens..nrows(rofA[i-2]));
     3297              i2=intvec(1..ncols(rofA[i-2]));
     3298              rofA[i-2]=transpose(submat(rofA[i-2],i1,i2));
     3299              logens=logens+ncols(rofA[i-2]);
     3300              rofA[i-2]=subst(rofA[i-2],h,1);
     3301            }
     3302          else
     3303            {
     3304              rofA[i-2]=list();
     3305            }
     3306        }
     3307      else
     3308        {
     3309          rofA[i-2]=list();
     3310        }
     3311    }
     3312  if(size(rofA[1])==0)
     3313    {
     3314      omitemptylist=1;
     3315    }
     3316  setring B;
     3317  vdint=fetch(ringofSyz,vdint);
     3318  if (omitemptylist!=1)
     3319    {
     3320      list rofA=fetch(ringofSyz,rofA);
     3321    }
     3322  kill HomWeyl;
     3323  kill ringofSyz;
     3324  return(list(vdint,rofA));
     3325}
     3326
     3327////////////////////////////////////////////////////////////////////////////////////
     3328
    7683329static proc toVdStrictFreeComplex(list L,string Syzstring,list #)
    769 "USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form 
    770         (M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a 
     3330"USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form
     3331        (M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a
    7713332        string, d an optional integer
    7723333ASSUME: Basering is the Weyl algebra D_n @*
    7733334        (M_1,f_1,...,M_s,f_s) represents a complex 0->D_n^(r_1)/im(M_1)->
    774         D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i) 
     3335        D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i)
    7753336        is generated by the rows of M_i. In particular it hold:@*
    7763337        - The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @*
    7773338        -the image of M_1*f_i is contained in the image of M_(i+1) @*
    778         d is an integer between 1 and n. If no value for d is given, it is assumed
    779         to be n @*
     3339        d is an optional integer which indices in the case size(L)=2, whether a
     3340        V_d-strict or \tilde(V_d)-strict will be computed@*
    7803341        Syzstring is either: @*
    781         -'Sres' (computes the resolutions and Groebner bases in the homogenized 
    782          Weyl algebra using Schreyer's method)@* 
     3342        -'Sres' (computes the resolutions and Groebner bases in the homogenized
     3343         Weyl algebra using Schreyer's method)@*
    7833344        or @*
    784         -'Vdres' (computes the resolutions via V_d-homogenization and without 
     3345        -'Vdres' (computes the resolutions via V_d-homogenization and without
    7853346         Schreyer's method)@*
    7863347RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @*
    7873348        L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,i_s,g_s,m_s) such that:@*
    788         -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs 
     3349        -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs
    7893350         of size i_j@*
    7903351        -D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0  is a V_d-strict complex
    7913352         with differentials m_i that is quasi-isomorphic to the complex given by L@*
    792         L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and 
     3353        L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and
    7933354        the n_i are shift vectors such that:@*
    794         -coker(H_i) is the ith cohomology group of the complex given by L_1@* 
     3355        -coker(H_i) is the ith cohomology group of the complex given by L_1@*
    7953356        -the n_i are the shift vectors of the coker(H_i)
    7963357THEORY: We follow Algorithm 3.8 in [W2]
    7973358"
    7983359{
    799   def B=basering; 
     3360  def B=basering;
    8003361  int n=nvars(B) div 2+2;
    8013362  int d=nvars(B) div 2;
    … …  
    8053366  matrix mem;
    8063367  intvec i1,i2;
     3368  int tilde;
    8073369  if (size(#)!=0)
    808   {
    809     for (i=1; i<=size(#); i++)
    810     {
    811       if (typeof(#[i])=="int")
    812       {
    813         if (#[1]>=1 and #[1]<=n)
    814         {
    815           d=#[i];
    816         }
    817       }
    818     }
    819   }
     3370    {
     3371      for (i=1; i<=size(#); i++)
     3372        {
     3373          if (typeof(#[i])=="int")
     3374            {
     3375              tilde=#[i];
     3376            }
     3377        }
     3378    }
    8203379  /* If size(L)=2, our complex consists for only one non-trivial module.
    821   Therefore, we just have to compute a V_d-strict resolution of this module.*/
    822   if (size(L)==2)
    823   {
    824     v=(0:ncols(L[1]));
    825     out[3*n-1]=v;
    826     out[3*n-2]=ncols(L[1]);
    827     out[3*n]=L[2];
    828     if (Syzstring=="Vdres")
    829     {
    830       /*if Syzstring="Vdres", we compute a V_d-strict Groebner basis of L[1]
    831       using F-homogenization (Prop. 3.9 in [OT]); then we compute the syzygies
    832       and make them V_d-strict using Prop  3.9[OT] and so on*/
    833       out[3*n-3]=VdStrictGB(L[1],d,v);   
    834       for (i=n-1; i>=1; i--)
    835       {
    836         out[3*i-2]=nrows(out[3*i]);
    837         v=0;
    838         for (j=1; j<=out[3*i-2]; j++)
    839         {
    840           mem=submat(out[3*i],j,intvec(1..ncols(out[3*i])));
    841           v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector
    842         }
    843         out[3*i-1]=v;
    844         if (i!=1)
    845         {
    846           /*next step in the resolution*/
    847           out[3*i-3]=transpose(syz(transpose(out[3*i])));
    848           if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3])))
    849           {
    850             /*makes the resolution V_d-strict*/
    851             out[3*i-3]=VdStrictGB(out[3*i-3],d,out[3*i-1]);
    852           }
    853           else
    854           {
    855             /*resolution is already computed*/
    856             out[3*i-3]=matrix(0,1,ncols(out[3*i-3]));
    857             out[3*i-4]=intvec(0);
    858             out[3*i-5]=int(0);
    859             for (j=i-2; j>=1; j--)
    860             {
    861               out[3*j]=matrix(0,1,1);
    862               out[3*j-1]=intvec(0);
    863               out[3*j-2]=int(0);
    864             }
    865             break;
    866           }
    867         }
    868       }
    869     }
    870     else
    871     {
    872       /*in the case Syzstring!="Vdres" we compute the resolution in the
    873       homogenized Weyl algebra using Thm 9.10 in[OT]*/
    874       def HomWeyl=makeHomogenizedWeyl(d);
    875       setring HomWeyl;
    876       list L=fetch(B,L);
    877       L[1]=nHomogenize(L[1]);
    878       list out=fetch(B,out);
    879       out[3*n-3]=L[1];
    880       /*computes a ring with a list RES; RES is a V_d-strict resolution of
    881       L[1]*/
    882       def ringofSyz=Sres(transpose(L[1]),d);
    883       setring ringofSyz;
    884       int logens=2;
    885       matrix mem;
    886       list out=fetch(HomWeyl,out);   
    887       out[3*n-3]=transpose(matrix(RES[2]));
    888       out[3*n-3]=subst(out[3*n-3],h,1);
    889       for (i=n-1; i>=1; i--)
    890       {
    891         out[3*i-2]=nrows(out[3*i]);
    892         v=0;     
    893         for (j=1; j<=out[3*i-2]; j++)
    894         {
    895           mem=submat(out[3*i],j,intvec(1..ncols(out[3*i])));
    896           v[j]=VdDeg(mem,d, out[3*i+2]);
    897         }
    898         out[3*i-1]=v;//shift vector such that the resolution RES is V_d-strict
    899         if (i!=1)
    900         {
    901           indi=0;
    902           if (size(RES)>=n-i+2)
    903           {
    904             nr=nrows(matrix(RES[n-i+2]));
    905             mem=matrix(0,nr,ncols(matrix(RES[n-i+2])));
    906             if (matrix(RES[n-i+2])!=mem)
    907             {
    908               indi=1;
    909               out[3*i-3]=(matrix(RES[n-i+2]));                   
    910               if (nrows(out[3*i-3])-logens+1!=nrows(out[3*i]))
    911               {
    912                 mem=out[3*i-3];
    913                 out[3*i-3]=matrix(mem,nrows(mem)+logens-1,ncols(mem));
    914               }
    915               mem=out[3*i-3];
    916               i1=intvec(logens..nrows(mem));
    917               mem=submat(mem,i1,intvec(1..ncols(mem)));
    918               out[3*i-3]=transpose(mem);
    919               out[3*i-3]=subst(out[3*i-3],h,1);   
    920               logens=logens+ncols(out[3*i-3]);
    921             }               
    922           }
    923           if(indi==0)
    924           {
    925             out[3*i-3]=matrix(0,1,nrows(out[3*i]));
    926             out[3*i-4]=intvec(0);
    927             out[3*i-5]=int(0);
    928             for (j=i-2; j>=1; j--)
    929             {
    930               out[3*j]=matrix(0,1,1);
    931               out[3*j-1]=intvec(0);
    932               out[3*j-2]=int(0);
    933             }
    934             break;
    935           }
    936         }
    937       }
    938       setring B;
    939       out=fetch(ringofSyz,out);//contains the V_d-strict resolution
    940       kill ringofSyz;
    941     }
    942     outall[1]=out;
    943     outall[2]=list(list(out[3*n-3],out[3*n-1]));
    944     return(outall);
    945   }
    946   /*case size(L)>2: We compute a quasi-isomorphic free complex following Alg 3.8 in
    947   [W2]*/
    948   /* We denote the complex given by L as (C^i,d^i).
    949   We start by computing in the proc shortExaxtPieces representations for the
    950   short exact sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1), where the B^i, Z^i
    951   and H^i are coboundaries, cocycles and cohomology groups, respectively.*/
    952   out=shortExactPieces(L);
     3380     Therefore, we just have to compute a V_d-strict resolution of this module.*/
     3381  if (size(L)==2)
     3382    {
     3383      v=(0:ncols(L[1]));
     3384      out[3*n-1]=v;
     3385      out[3*n-2]=ncols(L[1]);
     3386      out[3*n]=L[2];
     3387      if (Syzstring=="Vdres")
     3388        {
     3389          /*if Syzstring="Vdres", we compute a V_d-strict Groebner basis of L[1]
     3390            using F-homogenization (Prop. 3.9 in [OT]); then we compute the syzygies
     3391            and make them V_d-strict using Prop  3.9[OT] and so on*/
     3392          out[3*n-3]=VdStrictGB(L[1],d,v);
     3393          for (i=n-1; i>=1; i--)
     3394            {
     3395              out[3*i-2]=nrows(out[3*i]);
     3396              v=0;
     3397              for (j=1; j<=out[3*i-2]; j++)
     3398                {
     3399                  mem=submat(out[3*i],j,intvec(1..ncols(out[3*i])));
     3400                  v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector
     3401                }
     3402              out[3*i-1]=v;
     3403              if (i!=1)
     3404                {
     3405                  /*next step in the resolution*/
     3406                  out[3*i-3]=transpose(syz(transpose(out[3*i])));
     3407                  if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3])))
     3408                    {
     3409                      /*makes the resolution V_d-strict*/
     3410                      out[3*i-3]=VdStrictGB(out[3*i-3],d,out[3*i-1]);
     3411                    }
     3412                  else
     3413                    {
     3414                      /*resolution is already computed*/
     3415                      out[3*i-3]=matrix(0,1,ncols(out[3*i-3]));
     3416                      out[3*i-4]=intvec(0);
     3417                      out[3*i-5]=int(0);
     3418                      for (j=i-2; j>=1; j--)
     3419                        {
     3420                          out[3*j]=matrix(0,1,1);
     3421                          out[3*j-1]=intvec(0);
     3422                          out[3*j-2]=int(0);
     3423                        }
     3424                      break;
     3425                    }
     3426                }
     3427            }
     3428        }
     3429      else
     3430        {
     3431          /*in the case Syzstring!="Vdres" we compute the resolution in the
     3432            homogenized Weyl algebra using Thm 9.10 in[OT]*/
     3433          if (tilde==0)
     3434            {
     3435              def HomWeyl=makeHomogenizedWeyl(d);
     3436            }
     3437          else
     3438            {
     3439              def HomWeyl=makeHomogenizedWeylTilde(d);
     3440            }
     3441          setring HomWeyl;
     3442          list L=fetch(B,L);
     3443          L[1]=nHomogenize(L[1]);
     3444          list out=fetch(B,out);
     3445          out[3*n-3]=L[1];
     3446          /*computes a ring with a list RES; RES is a V_d-strict resolution of
     3447            L[1]*/
     3448          def ringofSyz=Sres(transpose(L[1]),d);
     3449          setring ringofSyz;
     3450          int logens=2;
     3451          matrix mem;
     3452          list out=fetch(HomWeyl,out);
     3453          out[3*n-3]=transpose(matrix(RES[2]));
     3454          out[3*n-3]=subst(out[3*n-3],h,1);
     3455          for (i=n-1; i>=1; i--)
     3456            {
     3457              out[3*i-2]=nrows(out[3*i]);
     3458              v=0;
     3459              for (j=1; j<=out[3*i-2]; j++)
     3460                {
     3461                  mem=submat(out[3*i],j,intvec(1..ncols(out[3*i])));
     3462                  if (tilde==0)
     3463                    {
     3464                      v[j]=VdDeg(mem,d, out[3*i+2]);
     3465                    }
     3466                  else
     3467                    {
     3468                      v[j]=VdDegTilde(mem,d, out[3*i+2]);
     3469                    }
     3470                }
     3471              out[3*i-1]=v;//shift vector such that the resolution RES is V_d-strict
     3472              if (i!=1)
     3473                {
     3474                  indi=0;
     3475                  if (size(RES)>=n-i+2)
     3476                    {
     3477                      nr=nrows(matrix(RES[n-i+2]));
     3478                      mem=matrix(0,nr,ncols(matrix(RES[n-i+2])));
     3479                      if (matrix(RES[n-i+2])!=mem)
     3480                        {
     3481                          indi=1;
     3482                          out[3*i-3]=(matrix(RES[n-i+2]));
     3483                          if (nrows(out[3*i-3])-logens+1!=nrows(out[3*i]))
     3484                            {
     3485                              mem=out[3*i-3];
     3486                              out[3*i-3]=matrix(mem,nrows(mem)+logens-1,ncols(mem));
     3487                            }
     3488                          mem=out[3*i-3];
     3489                          i1=intvec(logens..nrows(mem));
     3490                          mem=submat(mem,i1,intvec(1..ncols(mem)));
     3491                          out[3*i-3]=transpose(mem);
     3492                          out[3*i-3]=subst(out[3*i-3],h,1);
     3493                          logens=logens+ncols(out[3*i-3]);
     3494                        }
     3495                    }
     3496                  if(indi==0)
     3497                    {
     3498                      out[3*i-3]=matrix(0,1,nrows(out[3*i]));
     3499                      out[3*i-4]=intvec(0);
     3500                      out[3*i-5]=int(0);
     3501                      for (j=i-2; j>=1; j--)
     3502                        {
     3503                          out[3*j]=matrix(0,1,1);
     3504                          out[3*j-1]=intvec(0);
     3505                          out[3*j-2]=int(0);
     3506                        }
     3507                      break;
     3508                    }
     3509                }
     3510            }
     3511          setring B;
     3512          out=fetch(ringofSyz,out);//contains the V_d-strict resolution
     3513          kill ringofSyz;
     3514        }
     3515      outall[1]=out;
     3516      outall[2]=list(list(out[3*n-3],out[3*n-1]));
     3517      return(outall);
     3518    }
     3519  /*case size(L)>2: We compute a quasi-isomorphic free complex following Alg 3.8 in
     3520    [W2]*/
     3521  /* We denote the complex given by L as (C^i,d^i).
     3522     We start by computing in the proc shortExaxtPieces representations for the
     3523     short exact sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1), where the B^i, Z^i
     3524     and H^i are coboundaries, cocycles and cohomology groups, respectively.*/
     3525  out=shortExactPieces(L);
    9533526  list rem;
    954   /* shortExactpiecesToVdStrict makes the sequences B^i->Z^i->H^i and 
    955   Z^i->C^i->B^(i+1) V_d-strict*/
     3527  /* shortExactpiecesToVdStrict makes the sequences B^i->Z^i->H^i and
     3528     Z^i->C^i->B^(i+1) V_d-strict*/
    9563529  rem=shortExactPiecesToVdStrict(out,d,Syzstring);
    957   /*VdStrictDoubleComplexes computes V_d-strict resolutions over the seqeunces from 
    958   proc shortExactPiecesToVdstrict*/
     3530  /*VdStrictDoubleComplexes computes V_d-strict resolutions over the seqeunces from
     3531    proc shortExactPiecesToVdstrict*/
    9593532  out=VdStrictDoubleComplexes(rem[1],d,Syzstring);
    9603533  for (i=1;i<=size(out); i++)
    961   {
    962     rem[2][i][1]=out[i][1][5][1];
    963     rem[2][i][2]=out[i][1][8][1];
    964   }
    965   /* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences 
    966   Z^i->C^i->B^(i+1)) together to obtain a Cartan-Eilenberg resolution of
    967   (C^i,d^i)*/
     3534    {
     3535      rem[2][i][1]=out[i][1][5][1];
     3536      rem[2][i][2]=out[i][1][8][1];
     3537    }
     3538  /* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences
     3539     Z^i->C^i->B^(i+1)) together to obtain a Cartan-Eilenberg resolution of
     3540     (C^i,d^i)*/
    9683541  out=assemblingDoubleComplexes(out);
    9693542  /*the proc totalComplex takes the total complex of the double complex from the
    970   proc assemblingDoubleComplexes*/
     3543    proc assemblingDoubleComplexes*/
    9713544  out=totalComplex(out);
    9723545  outall[1]=out;
    973   outall[2]=rem[2];//contains the cohomology groups and their shift vectors 
     3546  outall[2]=rem[2];//contains the cohomology groups and their shift vectors
    9743547  return (outall);
    9753548}
    … …  
    9843557  int count;
    9853558  for (i=m; i<=n; i++)
    986   {
    987     out[size(out)+1]=list();
    988     for (j=1; j<=size(L[i]); j++)
    989     {
    990       count=count+1;
    991       out[size(out)][j]=list(L[i][j],count);
    992     }
    993   }
     3559    {
     3560      out[size(out)+1]=list();
     3561      for (j=1; j<=size(L[i]); j++)
     3562        {
     3563          count=count+1;
     3564          out[size(out)][j]=list(L[i][j],count);
     3565        }
     3566    }
    9943567  list o=list(out,count);
    9953568  return(o);
    … …  
    9983571////////////////////////////////////////////////////////////////////////////////////
    9993572
    1000 static proc LMSubset(list L,list M)
     3573static proc LMSubset(list L,list M, list #)
    10013574{
    10023575  int i;
    10033576  int j=1;
    1004   list position=(M[size(M)],(-1)^(size(L)));
     3577  if (size(#)==0)
     3578    {
     3579      list position=(M[size(M)],(-1)^(size(L)));
     3580    }
     3581  else
     3582    {
     3583      list position=(M[size(M)],1);
     3584    }
    10053585  for (i=1; i<=size(L); i++)
    1006   {
    1007     if (L[i]!=M[j])
    1008     {
    1009       if (L[i]!=M[i+1] or j!=i)
    1010       {
    1011         return (L[i],0);
    1012       }
    1013       else
    1014       {
    1015         position=(M[i],(-1)^(i-1));
    1016         j=j+1;
    1017       }
    1018     }
    1019     j=j+1;
    1020 
    1021   }
     3586    {
     3587      if (L[i]!=M[j])
     3588        {
     3589          if (L[i]!=M[i+1] or j!=i)
     3590            {
     3591              return (L[i],0);
     3592            }
     3593          else
     3594            {
     3595              if (size(#)==0)
     3596                {
     3597                  position=(M[i],(-1)^(i-1));
     3598                }
     3599              else
     3600                {
     3601                  position=(M[i],(-1)^(size(L)+1-i));
     3602                }
     3603              j=j+1;
     3604            }
     3605        }
     3606      j=j+1;
     3607
     3608    }
    10223609  return (position);
    10233610}
    … …  
    10283615{
    10293616  /*we follow Section 3.3 in [W2]*/
    1030   /* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex  C=(C^i,d^i) 
    1031   as in the procedure toVdstrictcomplex*/
     3617  /* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex  C=(C^i,d^i)
     3618     as in the procedure toVdstrictcomplex*/
    10323619  matrix  Bnew= divdr(L[2],L[3]);
    10333620  matrix Bold=Bnew;
    … …  
    10383625  list sep;
    10393626  /* the list sep will be of size s such that
    1040   -sep[i]=(sep[i][1],sep[i][2]) is a list of two lists
    1041   -sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)->coker(Z^i)
    1042   ->coker(H^i) represents the short exact seqeuence B^i(C)->Z^i(C)->H^i(C)
    1043   -sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)->coker(C^i)->
    1044   coker(B^(i+1)) represents the short exact seqeuence Z^i(C)->C^i->B^(i+1)(C)*/
     3627     -sep[i]=(sep[i][1],sep[i][2]) is a list of two lists
     3628     -sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)->coker(Z^i)
     3629      ->coker(H^i) represents the short exact seqeuence B^i(C)->Z^i(C)->H^i(C)
     3630     -sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)->coker(C^i)->
     3631      coker(B^(i+1)) represents the short exact seqeuence Z^i(C)->C^i->B^(i+1)(C)*/
    10453632  sep[1]=list(bzh,zcb);
    10463633  int i;
    10473634  list out;
    10483635  for (i=3; i<=size(L)-2; i=i+2)
    1049   {
    1050     /*the proc bzhzcb computes representations for the short exact seqeunces */
    1051     out=bzhzcb(Bold, L[i-1] , L[i], L[i+1], L[i+2]);
    1052     sep[size(sep)+1]=out[1];
    1053     Bold=out[2];
    1054   }
     3636    {
     3637      /*the proc bzhzcb computes representations for the short exact seqeunces */
     3638      out=bzhzcb(Bold, L[i-1] , L[i], L[i+1], L[i+2]);
     3639      sep[size(sep)+1]=out[1];
     3640      Bold=out[2];
     3641    }
    10553642  bzh=(divdr(L[size(L)-2], L[size(L)-1]),L[size(L)-2], L[size(L)-1]);
    10563643  bzh[4]=unitmat(ncols(L[size(L)-1]));
    … …  
    10793666static proc shortExactPiecesToVdStrict(list C,int d,list #)
    10803667{/* We transform the short exact pieces from procedure shortExactPieces to V_d-
    1081 strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in
    1082 [W2].*/
     3668    strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in
     3669    [W2].*/
    10833670  /* If we compute our Groebner bases in the homogenized Weyl algebra, we already
    1084   compute some resolutions it omit additional Groebner basis computations later
    1085   on.*/
     3671     compute some resolutions it omit additional Groebner basis computations later
     3672     on.*/
    10863673  int s =size(C);int i; int j;
    10873674  string Syzstring="Sres";
    1088   intvec v=0:ncols(C[s][1][5]); 
     3675  intvec v=0:ncols(C[s][1][5]);
    10893676  if (size(#)!=0)
    1090   {
    1091     for (i=1; i<=size(#); i++)
    1092     {
    1093       if (typeof(#[i])=="string")
    1094       {
    1095         Syzstring=#[i];
    1096       }
    1097       if (typeof(#[i])=="intvec")
    1098       {
    1099         v=#[i];
    1100       }
    1101     }
    1102   }
     3677    {
     3678      for (i=1; i<=size(#); i++)
     3679        {
     3680          if (typeof(#[i])=="string")
     3681            {
     3682              Syzstring=#[i];
     3683            }
     3684          if (typeof(#[i])=="intvec")
     3685            {
     3686               v=#[i];
     3687            }
     3688        }
     3689    }
    11033690  list out;
    11043691  list forout;
    11053692  if (Syzstring=="Vdres")
    1106   {
    1107     out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s));
    1108   }
     3693    {
     3694      out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s));
     3695    }
    11093696  else
    1110   {
    1111     forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s);
    1112     list resolutionofA=forout[9];
    1113     list resolutionofC=forout[10];
    1114     forout=delete(forout,10);
    1115     forout=delete(forout,9);
    1116     out[s]=list(forout);
    1117     for (i=1; i<=size(resolutionofC); i++)
    1118     {
    1119       out[s][1][5][i+1]=resolutionofC[i];//save the resolutions
    1120       out[s][1][1][i+1]=resolutionofA[i];
    1121     }     
    1122   }
     3697    {
     3698      forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s);
     3699      list resolutionofA=forout[9];
     3700      list resolutionofC=forout[10];
     3701      forout=delete(forout,10);
     3702      forout=delete(forout,9);
     3703      out[s]=list(forout);
     3704      for (i=1; i<=size(resolutionofC); i++)
     3705        {
     3706          out[s][1][5][i+1]=resolutionofC[i];//save the resolutions
     3707          out[s][1][1][i+1]=resolutionofA[i];
     3708        }
     3709    }
    11233710  out[s][2]=list(list(out[s][1][3][1]));
    11243711  out[s][2][2]=list(unitmat(ncols(out[s][1][3][1])));
    … …  
    11323719  list resolutionofF;
    11333720  for (i=s-1; i>=2; i--)
    1134   {
    1135     C[i][2][5]=out[i+1][1][1][1];
    1136     forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s);
    1137     if (Syzstring=="Sres")
    1138     {
    1139       resolutionofD=forout[3];//save the resolutions
    1140       resolutionofF=forout[4];
    1141       forout=delete(forout,4);
    1142       forout=delete(forout,3);
    1143     }
    1144     out[i]=forout; 
    1145     if(Syzstring=="Sres")
    1146     {
    1147       for (j=2; j<=size(out[i+1][1][1]); j++)
    1148       {
    1149         out[i][2][5][j]=out[i+1][1][1][j];
    1150       }       
    1151       for (j=1; j<=size(resolutionofD);j++)
    1152       {
    1153         out[i][1][1][j+1]=resolutionofD[j];
    1154         out[i][1][5][j+1]=resolutionofF[j];
    1155       }                 
    1156     }
    1157   }
     3721    {
     3722      C[i][2][5]=out[i+1][1][1][1];
     3723      forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s);
     3724      if (Syzstring=="Sres")
     3725        {
     3726          resolutionofD=forout[3];//save the resolutions
     3727          resolutionofF=forout[4];
     3728          forout=delete(forout,4);
     3729          forout=delete(forout,3);
     3730        }
     3731      out[i]=forout;
     3732      if(Syzstring=="Sres")
     3733        {
     3734          for (j=2; j<=size(out[i+1][1][1]); j++)
     3735            {
     3736              out[i][2][5][j]=out[i+1][1][1][j];
     3737            }
     3738          for (j=1; j<=size(resolutionofD);j++)
     3739            {
     3740              out[i][1][1][j+1]=resolutionofD[j];
     3741              out[i][1][5][j+1]=resolutionofF[j];
     3742            }
     3743        }
     3744    }
    11583745  out[1]=list(list());//initalize our list
    11593746  C[1][2][5]=out[2][1][1][1];
    11603747  /*Compute the last V_d-strict seqeunce*/
    11613748  if (Syzstring=="Vdres")
    1162   {
    1163     out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
    1164   }
     3749    {
     3750      out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
     3751    }
    11653752  else
    1166   {
    1167     forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
    1168     out[1][2]=delete(forout,9);
    1169     list resolutionofA2=forout[9];
    1170     for (i=1; i<=size(out[2][1][1]); i++)
    1171     {
    1172       /*put the modules for the resolutions in the right spot*/
    1173       out[1][2][5][i]=out[2][1][1][i];
    1174     }
    1175     for (i=1; i<=size(resolutionofA2); i++)
    1176     {
    1177       out[1][2][1][i+1]=resolutionofA2[i];
    1178     }
    1179   }
     3753    {
     3754      forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
     3755      out[1][2]=delete(forout,9);
     3756      list resolutionofA2=forout[9];
     3757      for (i=1; i<=size(out[2][1][1]); i++)
     3758        {
     3759          /*put the modules for the resolutions in the right spot*/
     3760          out[1][2][5][i]=out[2][1][1][i];
     3761        }
     3762      for (i=1; i<=size(resolutionofA2); i++)
     3763        {
     3764          out[1][2][1][i+1]=resolutionofA2[i];
     3765        }
     3766    }
    11803767  out[1][1][3]=list(out[1][2][1][1]);
    11813768  out[1][1][5]=list(out[1][2][1][1]);
    … …  
    11873774  out[1][1][6]=list(list());
    11883775  if (Syzstring=="Sres")
    1189   {
    1190     for (i=1; i<=size(out[1][2][1]); i++)
    1191     {
    1192       out[1][1][3][i]=out[1][2][1][i];
    1193       out[1][1][5][i]=out[1][2][1][i];
    1194     }
    1195   }
     3776    {
     3777      for (i=1; i<=size(out[1][2][1]); i++)
     3778        {
     3779          out[1][1][3][i]=out[1][2][1][i];
     3780          out[1][1][5][i]=out[1][2][1][i];
     3781        }
     3782    }
    11963783  list Hi;
    11973784  for (i=1; i<=size(out); i++)
    1198   {
    1199     Hi[i]=list(out[i][1][5][1],out[i][1][8][1]);
    1200   }
     3785    {
     3786      Hi[i]=list(out[i][1][5][1],out[i][1][8][1]);
     3787    }
    12013788  list outall;
    12023789  outall[1]=out;
    … …  
    12073794////////////////////////////////////////////////////////////////////////////////////
    12083795
    1209 static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #) 
     3796static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #)
    12103797{
    12113798  /*this is the Algorithm 3.11 in [W2]*/
    … …  
    12153802  def B=basering;
    12163803  matrix bi=slimgb(transpose(C[5]));
    1217   /* Computation of a V_d-strict Groebner basis of C[5]: 
    1218   -if Syzstring=="Vdres" this is done using the method of weighted homogenization
    1219   (Prop. 3.9 [OT])
    1220   -else we use the homogenized Weyl algebra for Groebner basis computations
    1221   (Prop 9.9 [OT]),
    1222   in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit
    1223   extra Groebner basis computations later on*/
     3804  /* Computation of a V_d-strict Groebner basis of C[5]:
     3805     -if Syzstring=="Vdres" this is done using the method of weighted homogenization
     3806     (Prop. 3.9 [OT])
     3807     -else we use the homogenized Weyl algebra for Groebner basis computations
     3808     (Prop 9.9 [OT]),
     3809     in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit
     3810     extra Groebner basis computations later on*/
    12243811  int nr,nc;
    12253812  intvec i1,i2;
    12263813  if (Syzstring=="Vdres")
    1227   {
    1228     if(size(#)==0)
    1229     {
    1230       matrix J_C=VdStrictGB(C[5],n,list(v));
    1231     }
    1232     else
    1233     {
    1234       matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
    1235     }
    1236   }
     3814    {
     3815      if(size(#)==0)
     3816        {
     3817          matrix J_C=VdStrictGB(C[5],n,list(v));
     3818        }
     3819      else
     3820        {
     3821          matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
     3822        }
     3823    }
    12373824  else
    1238   {
    1239     if (size(#)==0)
    1240     {
    1241       matrix MC=C[5];
    1242       def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v);
    1243       setring HomWeyl;
    1244       matrix J_C=fetch(B,MC);
    1245       J_C=nHomogenize(J_C);
    1246       /*computation of V_d-strict resolution of C[5]->needed for proc
    1247       VdstrictDoubleComplexes*/
    1248       def ringofSyz=Sres(groebner(transpose(J_C)), lengthofres);
    1249       setring ringofSyz;
    1250       matrix J_C=transpose(matrix(RES[2]));
    1251       J_C=subst(J_C,h,1);
    1252       logens=ncols(J_C)+1;
    1253       matrix zerom;
    1254       list rofC;//will contain resolution of C
    1255       for (i=3; i<=n+si+1; i++)
    1256       {
    1257         if (size(RES)>=i)
    1258         {
    1259           zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
    1260           if (RES[i]!=zerom)
    1261           {
    1262             rofC[i-2]=(matrix(RES[i]));
    1263 
    1264             if (i==3)
    1265             {
    1266               if (nrows(rofC[i-2])-logens+1!=nrows(J_C))
    1267               {
    1268                 //build the resolution
    1269                 nr=nrows(J_C)+logens-1;
    1270                 nc=ncols(rofC[i-2]);
    1271                 rofC[i-2]=matrix(rofC[i-2],nr,nc);
    1272               }
    1273 
    1274             }
    1275             if (i!=3)
    1276             {
    1277               if (nrows(rofC[i-2])-logens+1!=nrows(rofC[i-3]))
    1278               {
    1279                 nr=nrows(rofC[i-3])+logens-1;
    1280                 nc=ncols(rofC[i-2]);
    1281                 rofC[i-2]=matrix(rofC[i-2],nr,nc);
    1282               }
    1283             }
    1284             i1=intvec(logens..nrows(rofC[i-2]));
    1285             i2=intvec(1..ncols(rofC[i-2]));
    1286             rofC[i-2]=transpose(submat(rofC[i-2],i1,i2));
    1287             logens=logens+ncols(rofC[i-2]);
    1288             rofC[i-2]=subst(rofC[i-2],h,1);
    1289           }
    1290           else
    1291           {
    1292             rofC[i-2]=list();
    1293           }
    1294         }
    1295         else
    1296         {
    1297           rofC[i-2]=list();
    1298         }
    1299       }
    1300       if(size(rofC[1])==0)
    1301       {
    1302         omitemptylist=1;
    1303       }
    1304       setring B;
    1305       matrix  J_C=fetch(ringofSyz,J_C);
    1306       if (omitemptylist!=1)
    1307       {
    1308         list rofC=fetch(ringofSyz,rofC);
    1309       }
    1310       omitemptylist=0;
    1311       kill HomWeyl;
    1312       kill ringofSyz;
    1313     }
    1314     else
    1315     {
    1316       matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
    1317     }
    1318   }
     3825    {
     3826      if (size(#)==0)
     3827        {
     3828          matrix MC=C[5];
     3829          def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v);
     3830          setring HomWeyl;
     3831          matrix J_C=fetch(B,MC);
     3832          J_C=nHomogenize(J_C);
     3833          /*computation of V_d-strict resolution of C[5]->needed for proc
     3834            VdstrictDoubleComplexes*/
     3835          def ringofSyz=Sres(transpose(J_C),lengthofres);
     3836          setring ringofSyz;
     3837          matrix J_C=transpose(matrix(RES[2]));
     3838          J_C=subst(J_C,h,1);
     3839          logens=ncols(J_C)+1;
     3840          matrix zerom;
     3841          list rofC;//will contain resolution of C
     3842          for (i=3; i<=n+si+1; i++)
     3843            {
     3844              if (size(RES)>=i)
     3845                {
     3846                  zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
     3847                  if (RES[i]!=zerom)
     3848                    {
     3849                      rofC[i-2]=(matrix(RES[i]));
     3850
     3851                      if (i==3)
     3852                        {
     3853                          if (nrows(rofC[i-2])-logens+1!=nrows(J_C))
     3854                            {
     3855                              //build the resolution
     3856                              nr=nrows(J_C)+logens-1;
     3857                              nc=ncols(rofC[i-2]);
     3858                              rofC[i-2]=matrix(rofC[i-2],nr,nc);
     3859                            }
     3860
     3861                        }
     3862                      if (i!=3)
     3863                        {
     3864                          if (nrows(rofC[i-2])-logens+1!=nrows(rofC[i-3]))
     3865                            {
     3866                              nr=nrows(rofC[i-3])+logens-1;
     3867                              nc=ncols(rofC[i-2]);
     3868                              rofC[i-2]=matrix(rofC[i-2],nr,nc);
     3869                            }
     3870                        }
     3871                      i1=intvec(logens..nrows(rofC[i-2]));
     3872                      i2=intvec(1..ncols(rofC[i-2]));
     3873                      rofC[i-2]=transpose(submat(rofC[i-2],i1,i2));
     3874                      logens=logens+ncols(rofC[i-2]);
     3875                      rofC[i-2]=subst(rofC[i-2],h,1);
     3876                    }
     3877                  else
     3878                    {
     3879                      rofC[i-2]=list();
     3880                    }
     3881                }
     3882              else
     3883                {
     3884                  rofC[i-2]=list();
     3885                }
     3886            }
     3887          if(size(rofC[1])==0)
     3888            {
     3889              omitemptylist=1;
     3890            }
     3891          setring B;
     3892          matrix  J_C=fetch(ringofSyz,J_C);
     3893          if (omitemptylist!=1)
     3894            {
     3895              list rofC=fetch(ringofSyz,rofC);
     3896            }
     3897          omitemptylist=0;
     3898          kill HomWeyl;
     3899          kill ringofSyz;
     3900        }
     3901      else
     3902        {
     3903          matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
     3904        }
     3905    }
    13193906  /* we compute a V_d-strict Groebner basis for C[3]*/
    13203907  matrix J_A=C[1];
    … …  
    13263913  matrix Pi[1][ncols(J_AC)];
    13273914  for (i=1; i<=nrows(J_C); i++)
    1328   {
    1329     for (j=1; j<=nrows(J_AC);j++)
    1330     {
    1331       Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC)));
    1332     }
    1333     storePi[i]=Pi;
    1334     Pi=0;
    1335   }
     3915    {
     3916      for (j=1; j<=nrows(J_AC);j++)
     3917        {
     3918          Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC)));
     3919        }
     3920      storePi[i]=Pi;
     3921      Pi=0;
     3922    }
    13363923  /*we compute the shift vector for C[1]*/
    13373924  intvec m_a;
    13383925  list findMin;
    1339   int comMin; 
     3926  int comMin;
    13403927  for (i=1; i<=ncols(J_A); i++)
    1341   {
    1342     for (j=1; j<=size(storePi);j++)
    1343     {
    1344       if (storePi[j][1,i]!=0)
    1345       {
    1346         comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v);
    1347         comMin=comMin-VdDeg(storePi[j][1,i],n,intvec(0));
    1348         findMin[size(findMin)+1]=comMin;
    1349       }
    1350     }
    1351     if (size(findMin)!=0)
    1352     {
    1353       m_a[i]=Min(findMin);
    1354       findMin=list();
    1355     }
    1356     else
    1357     {
    1358       m_a[i]=0;
    1359     }
    1360   }
     3928    {
     3929      for (j=1; j<=size(storePi);j++)
     3930        {
     3931          if (storePi[j][1,i]!=0)
     3932            {
     3933              comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v);
     3934              comMin=comMin-VdDeg(storePi[j][1,i],n,intvec(0));
     3935              findMin[size(findMin)+1]=comMin;
     3936            }
     3937        }
     3938      if (size(findMin)!=0)
     3939        {
     3940          m_a[i]=Min(findMin);
     3941          findMin=list();
     3942        }
     3943      else
     3944        {
     3945          m_a[i]=0;
     3946        }
     3947    }
    13613948  matrix zero[ncols(J_A)][ncols(J_C)];
    13623949  matrix g_AB=concat(unitmat(ncols(J_A)),zero);
    13633950  matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C)))));
    13643951  intvec m_b=m_a,v;
    1365   /* computation of a V_d-strict Groebner basis of C[1] (and resolution if 
    1366   Syzstring=='Vdres') */
     3952  /* computation of a V_d-strict Groebner basis of C[1] (and resolution if
     3953     Syzstring=='Vdres') */
    13673954  if (Syzstring=="Vdres")
    1368   {
    1369     J_A=VdStrictGB(J_A,n,m_a);
    1370   }
     3955    {
     3956      J_A=VdStrictGB(J_A,n,m_a);
     3957    }
    13713958  else
    1372   {
    1373     def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a);
    1374     setring HomWeyl;
    1375     matrix J_A=fetch(B,J_A);
    1376     J_A=nHomogenize(J_A);
    1377     def ringofSyz=Sres(groebner(transpose(J_A)),lengthofres);
    1378     setring ringofSyz;
    1379     matrix J_A=transpose(matrix(RES[2]));
    1380     matrix zerom;
    1381     J_A=subst(J_A,h,1);
    1382     logens=ncols(J_A)+1;
    1383     list rofA;
    1384     for (i=3; i<=n+si+1; i++)
    1385     {
    1386       if (size(RES)>=i)
    1387       {
    1388         zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
    1389         if (RES[i]!=zerom)
    1390         {
    1391           rofA[i-2]=matrix(RES[i]);// resolution for C[1]
    1392           if (i==3)
    1393           {
    1394             if (nrows(rofA[i-2])-logens+1!=nrows(J_A))
    1395             {
    1396               nr=nrows(J_A)+logens-1;
    1397               nc=ncols(rofA[i-2]);
    1398               rofA[i-2]=matrix(rofA[i-2],nr,nc);
    1399             }
    1400           }
    1401           if (i!=3)
    1402           {
    1403             if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3]))
    1404             {
    1405               nr=nrows(rofA[i-3])+logens-1;
    1406               nc=ncols(rofA[i-2]);
    1407               rofA[i-2]=matrix(rofA[i-2],nr,nc);
    1408             }
    1409           }
    1410           i1=intvec(logens..nrows(rofA[i-2]));
    1411           i2=intvec(1..ncols(rofA[i-2]));
    1412           rofA[i-2]=transpose(submat(rofA[i-2],i1,i2));
    1413           logens=logens+ncols(rofA[i-2]);
    1414           rofA[i-2]=subst(rofA[i-2],h,1);
    1415         }
    1416         else
    1417         {
    1418           rofA[i-2]=list();
    1419         }
    1420       }
    1421       else
    1422       {
    1423         rofA[i-2]=list();
    1424       }
    1425     }
    1426     if(size(rofA[1])==0)
    1427     {
    1428       omitemptylist=1;
    1429     }
    1430     setring B;
    1431     J_A=fetch(ringofSyz,J_A);
    1432     if (omitemptylist!=1)
    1433     {
    1434       list rofA=fetch(ringofSyz,rofA);
    1435     }
    1436     omitemptylist=0;
    1437     kill HomWeyl;
    1438     kill ringofSyz;
    1439   }
     3959    {
     3960      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a);
     3961      setring HomWeyl;
     3962      matrix J_A=fetch(B,J_A);
     3963      J_A=nHomogenize(J_A);
     3964      def ringofSyz=Sres(transpose(J_A),lengthofres);
     3965      setring ringofSyz;
     3966      matrix J_A=transpose(matrix(RES[2]));
     3967      matrix zerom;
     3968      J_A=subst(J_A,h,1);
     3969      logens=ncols(J_A)+1;
     3970      list rofA;
     3971      for (i=3; i<=n+si+1; i++)
     3972        {
     3973          if (size(RES)>=i)
     3974            {
     3975              zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
     3976              if (RES[i]!=zerom)
     3977                {
     3978                  rofA[i-2]=matrix(RES[i]);// resolution for C[1]
     3979                  if (i==3)
     3980                    {
     3981                      if (nrows(rofA[i-2])-logens+1!=nrows(J_A))
     3982                        {
     3983                          nr=nrows(J_A)+logens-1;
     3984                          nc=ncols(rofA[i-2]);
     3985                          rofA[i-2]=matrix(rofA[i-2],nr,nc);
     3986                        }
     3987                    }
     3988                  if (i!=3)
     3989                    {
     3990                      if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3]))
     3991                        {
     3992                          nr=nrows(rofA[i-3])+logens-1;
     3993                          nc=ncols(rofA[i-2]);
     3994                          rofA[i-2]=matrix(rofA[i-2],nr,nc);
     3995                        }
     3996                    }
     3997                  i1=intvec(logens..nrows(rofA[i-2]));
     3998                  i2=intvec(1..ncols(rofA[i-2]));
     3999                  rofA[i-2]=transpose(submat(rofA[i-2],i1,i2));
     4000                  logens=logens+ncols(rofA[i-2]);
     4001                  rofA[i-2]=subst(rofA[i-2],h,1);
     4002                }
     4003              else
     4004                {
     4005                  rofA[i-2]=list();
     4006                }
     4007            }
     4008          else
     4009            {
     4010              rofA[i-2]=list();
     4011            }
     4012        }
     4013      if(size(rofA[1])==0)
     4014        {
     4015          omitemptylist=1;
     4016        }
     4017      setring B;
     4018      J_A=fetch(ringofSyz,J_A);
     4019      if (omitemptylist!=1)
     4020        {
     4021          list rofA=fetch(ringofSyz,rofA);
     4022        }
     4023      omitemptylist=0;
     4024      kill HomWeyl;
     4025      kill ringofSyz;
     4026    }
    14404027  J_AC=transpose(storePi[1]);
    14414028  for (i=2; i<= size(storePi); i++)
    1442   {
    1443     J_AC=concat(J_AC, transpose(storePi[i]));
    1444   }
     4029    {
     4030      J_AC=concat(J_AC, transpose(storePi[i]));
     4031    }
    14454032  J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC));
    14464033  list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a));
    14474034  Vdstrict[7]=list(m_b);
    1448   Vdstrict[8]=list(v); 
     4035  Vdstrict[8]=list(v);
    14494036  if(Syzstring=="Sres")
    1450   {
    1451     Vdstrict[9]=rofA;
    1452     if(size(#)==0)
    1453     {
    1454       Vdstrict[10]=rofC;
    1455     }
    1456   }
     4037    {
     4038      Vdstrict[9]=rofA;
     4039      if(size(#)==0)
     4040        {
     4041          Vdstrict[10]=rofC;
     4042        }
     4043    }
    14574044  return (Vdstrict);
    14584045}
    … …  
    14604047////////////////////////////////////////////////////////////////////////////////////
    14614048
    1462 static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL)     
     4049static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL)
    14634050{
    1464   /* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact 
    1465   pieces.
    1466   We asume that we are given two sequences of the form coker(L[i][1])->
    1467   coker(L[i][3])->coker(L[i][5]) with differentials L[i][2] and L[i][4] such
    1468   that L[1][3]=L[2][1].We are going to transform them to V_d-strict sequences
    1469   J_D->J_A->J_F and J_A->J_B->J_C*/
     4051  /* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact
     4052     pieces.
     4053     We asume that we are given two sequences of the form coker(L[i][1])->
     4054     coker(L[i][3])->coker(L[i][5]) with differentials L[i][2] and L[i][4] such
     4055     that L[1][3]=L[2][1].We are going to transform them to V_d-strict sequences
     4056     J_D->J_A->J_F and J_A->J_B->J_C*/
    14704057  int omitemptylist;
    14714058  int lengthofres=sizeL+d-1;
    … …  
    14754062  matrix J_D=L[1][1];
    14764063  matrix f_FA=L[1][4];
    1477   /*We find new presentations coker(J_DF) and coker(J_DFC)  for L[1][4]=L[2][1] 
    1478   and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */
     4064  /*We find new presentations coker(J_DF) and coker(J_DFC)  for L[1][4]=L[2][1]
     4065     and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */
    14794066  matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA)));
    14804067  matrix J_DF=divdr(f_DFA,L[1][3]);//coker(J_DF) is isomorphic to coker(L[2][1]);
    … …  
    14834070  matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB)));
    14844071  matrix J_DFC=divdr(f_DFCB,L[2][3]);//coker(J_DFC) are coker(L[2][3)]) isomorphic
    1485   /* find a shift vector on the range of J_F such that the first sequence is 
    1486   exact*/
     4072  /* find a shift vector on the range of J_F such that the first sequence is
     4073     exact*/
    14874074  matrix P=matrixLift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C);
    14884075  list storePi;
    … …  
    14904077  int i; int j;
    14914078  for (i=1; i<=nrows(J_C); i++)
    1492   {
    1493     for (j=1; j<=nrows(J_DFC);j++)
    1494     {
    1495       Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC)));
    1496     }
    1497     storePi[i]=Pi;
    1498     Pi=0;
    1499   }
     4079    {
     4080      for (j=1; j<=nrows(J_DFC);j++)
     4081        {
     4082          Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC)));
     4083        }
     4084      storePi[i]=Pi;
     4085      Pi=0;
     4086    }
    15004087  intvec m_a;
    15014088  list findMin;
    … …  
    15054092  intvec i1,i2;
    15064093  for (i=1; i<=ncols(J_DF); i++)
    1507   {
    1508     for (j=1; j<=size(storePi);j++)
    1509     {
    1510       if (storePi[j][1,i]!=0)
    1511       {
    1512         comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v);
    1513         comMin=comMin-VdDeg(storePi[j][1,i],d,intvec(0));
    1514         findMin[size(findMin)+1]=comMin;
    1515       }
    1516     }
    1517     if (size(findMin)!=0)
    1518     {
    1519       m_a[i]=Min(findMin);// shift vector for L[2][1]
    1520       findMin=list();
    1521       noMin[i]=0;
    1522     }
    1523     else
    1524     {
    1525       noMin[i]=1;
    1526     }
    1527   }
     4094    {
     4095      for (j=1; j<=size(storePi);j++)
     4096        {
     4097          if (storePi[j][1,i]!=0)
     4098            {
     4099              comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v);
     4100              comMin=comMin-VdDeg(storePi[j][1,i],d,intvec(0));
     4101              findMin[size(findMin)+1]=comMin;
     4102            }
     4103        }
     4104      if (size(findMin)!=0)
     4105        {
     4106          m_a[i]=Min(findMin);// shift vector for L[2][1]
     4107          findMin=list();
     4108          noMin[i]=0;
     4109        }
     4110      else
     4111        {
     4112          noMin[i]=1;
     4113        }
     4114    }
    15284115  if (size(m_a) < ncols(J_DF))
    1529   {
    1530     m_a[ncols(J_DF)]=0;
    1531   }
     4116    {
     4117      m_a[ncols(J_DF)]=0;
     4118    }
    15324119  intvec m_f=m_a[ncols(J_D)+1..size(m_a)];
    1533   /* Computation of a V_d-strict Groebner basis of J_F=L[1][5]: 
    1534   if Syzstring=="Vdres" this is done using the method of weighted homogenization
    1535   (Prop. 3.9 [OT])
    1536   else we use the homogenized Weyl algerba for Groebner basis computations
    1537   (Prop 9.9 [OT]), in this case we already compute resolutions
    1538   (Thm. 9.10 in [OT]) to omit extra Groebner basis  computations later on*/
     4120  /* Computation of a V_d-strict Groebner basis of J_F=L[1][5]:
     4121     if Syzstring=="Vdres" this is done using the method of weighted homogenization
     4122     (Prop. 3.9 [OT])
     4123     else we use the homogenized Weyl algerba for Groebner basis computations
     4124     (Prop 9.9 [OT]), in this case we already compute resolutions
     4125     (Thm. 9.10 in [OT]) to omit extra Groebner basis  computations later on*/
    15394126  if (Syzstring=="Vdres")
    1540   {
    1541     J_F=VdStrictGB(J_F,d,m_f);
    1542   }
     4127    {
     4128      J_F=VdStrictGB(J_F,d,m_f);
     4129    }
    15434130  else
    1544   {
    1545     def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f);
    1546     setring HomWeyl;
    1547     matrix J_F=fetch(B,J_F);
    1548     J_F=nHomogenize(J_F);
    1549     def ringofSyz=Sres(groebner(transpose(J_F)),lengthofres);   
    1550     setring ringofSyz;
    1551     matrix J_F=transpose(matrix(RES[2]));
    1552     J_F=subst(J_F,h,1);
    1553     logens=ncols(J_F)+1;
    1554     list rofF;
    1555     for (i=3; i<=d+sizeL+1; i++)
    1556     {
    1557       if (size(RES)>=i)
    1558       {
    1559         if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
    1560         {
    1561           rofF[i-2]=(matrix(RES[i]));// resolution for J_F
    1562           if (i==3)
    1563           {
    1564             if (nrows(rofF[i-2])-logens+1!=nrows(J_F))
    1565             {
    1566               nr=nrows(J_F)+logens-1;
    1567               nc=ncols(rofF[i-2]);
    1568               rofF[i-2]=matrix(rofF[i-2],nr,nc);
    1569             }
    1570           }
    1571           if (i!=3)
    1572           {
    1573             if (nrows(rofF[i-2])-logens+1!=nrows(rofF[i-3]))
    1574             {
    1575               nr=nrows(rofF[i-3])+logens-1;
    1576               rofF[i-2]=matrix(rofF[i-2],nr,ncols(rofF[i-2]));
    1577             }
    1578           }
    1579           i1=intvec(logens..nrows(rofF[i-2]));
    1580           i2=intvec(1..ncols(rofF[i-2]));
    1581           rofF[i-2]=transpose(submat(rofF[i-2],i1,i2));
    1582           logens=logens+ncols(rofF[i-2]);
    1583           rofF[i-2]=subst(rofF[i-2],h,1);
    1584         }
    1585         else
    1586         {
    1587           rofF[i-2]=list();
    1588         }
    1589       }
    1590       else
    1591       {
    1592         rofF[i-2]=list();
    1593       }
    1594     }
    1595     if(size(rofF[1])==0)
    1596     {
    1597       omitemptylist=1;
    1598     }
    1599     setring B;
    1600     J_F=fetch(ringofSyz,J_F);
    1601     if (omitemptylist!=1)
    1602     {
    1603       list rofF=fetch(ringofSyz,rofF);
    1604     }
    1605     omitemptylist=0;
    1606     kill HomWeyl;
    1607     kill ringofSyz;
    1608   }
     4131    {
     4132      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f);
     4133      setring HomWeyl;
     4134      matrix J_F=fetch(B,J_F);
     4135      J_F=nHomogenize(J_F);
     4136      def ringofSyz=Sres(transpose(J_F),lengthofres);
     4137      setring ringofSyz;
     4138      matrix J_F=transpose(matrix(RES[2]));
     4139      J_F=subst(J_F,h,1);
     4140      logens=ncols(J_F)+1;
     4141      list rofF;
     4142      for (i=3; i<=d+sizeL+1; i++)
     4143        {
     4144          if (size(RES)>=i)
     4145            {
     4146              if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
     4147                {
     4148                  rofF[i-2]=(matrix(RES[i]));// resolution for J_F
     4149                  if (i==3)
     4150                    {
     4151                      if (nrows(rofF[i-2])-logens+1!=nrows(J_F))
     4152                        {
     4153                          nr=nrows(J_F)+logens-1;
     4154                          nc=ncols(rofF[i-2]);
     4155                          rofF[i-2]=matrix(rofF[i-2],nr,nc);
     4156                        }
     4157                    }
     4158                  if (i!=3)
     4159                    {
     4160                      if (nrows(rofF[i-2])-logens+1!=nrows(rofF[i-3]))
     4161                        {
     4162                          nr=nrows(rofF[i-3])+logens-1;
     4163                          rofF[i-2]=matrix(rofF[i-2],nr,ncols(rofF[i-2]));
     4164                        }
     4165                    }
     4166                  i1=intvec(logens..nrows(rofF[i-2]));
     4167                  i2=intvec(1..ncols(rofF[i-2]));
     4168                  rofF[i-2]=transpose(submat(rofF[i-2],i1,i2));
     4169                  logens=logens+ncols(rofF[i-2]);
     4170                  rofF[i-2]=subst(rofF[i-2],h,1);
     4171                }
     4172              else
     4173                {
     4174                  rofF[i-2]=list();
     4175                }
     4176            }
     4177          else
     4178            {
     4179              rofF[i-2]=list();
     4180            }
     4181        }
     4182      if(size(rofF[1])==0)
     4183        {
     4184          omitemptylist=1;
     4185        }
     4186      setring B;
     4187      J_F=fetch(ringofSyz,J_F);
     4188      if (omitemptylist!=1)
     4189        {
     4190          list rofF=fetch(ringofSyz,rofF);
     4191        }
     4192      omitemptylist=0;
     4193      kill HomWeyl;
     4194      kill ringofSyz;
     4195    }
    16094196  /*find shift vectors on the range of J_D*/
    16104197  P=matrixLift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F);
    … …  
    16124199  matrix Pidf[1][ncols(J_DF)];
    16134200  for (i=1; i<=nrows(J_F); i++)
    1614   {
    1615     for (j=1; j<=nrows(J_DF);j++)
    1616     {
    1617       Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF)));
    1618     }
    1619     storePinew[i]=Pidf;
    1620     Pidf=0;
    1621   }
     4201    {
     4202      for (j=1; j<=nrows(J_DF);j++)
     4203        {
     4204          Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF)));
     4205        }
     4206      storePinew[i]=Pidf;
     4207      Pidf=0;
     4208    }
    16224209  intvec m_d;
    16234210  for (i=1; i<=ncols(J_D); i++)
    1624   {
    1625     for (j=1; j<=size(storePinew);j++)
    1626     {
    1627       if (storePinew[j][1,i]!=0)
    1628       {
    1629         comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f);
    1630         comMin=comMin-VdDeg(storePinew[j][1,i],d,intvec(0));
    1631         findMin[size(findMin)+1]=comMin;
    1632       }
    1633     }
    1634     if (size(findMin)!=0)
    1635     {
    1636       if (noMin[i]==0)
    1637       {
    1638         m_d[i]=Min(insert(findMin,m_a[i]));
    1639         m_a[i]=m_d[i];
    1640       }
     4211    {
     4212      for (j=1; j<=size(storePinew);j++)
     4213        {
     4214          if (storePinew[j][1,i]!=0)
     4215            {
     4216              comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f);
     4217              comMin=comMin-VdDeg(storePinew[j][1,i],d,intvec(0));
     4218              findMin[size(findMin)+1]=comMin;
     4219            }
     4220        }
     4221      if (size(findMin)!=0)
     4222        {
     4223          if (noMin[i]==0)
     4224            {
     4225              m_d[i]=Min(insert(findMin,m_a[i]));
     4226              m_a[i]=m_d[i];
     4227            }
     4228          else
     4229            {
     4230              m_d[i]=Min(findMin);
     4231              m_a[i]=m_d[i];
     4232            }
     4233        }
    16414234      else
    1642       {
    1643         m_d[i]=Min(findMin);
    1644         m_a[i]=m_d[i];
    1645       }
    1646     }
    1647     else
    1648     {
    1649       m_d[i]=m_a[i];
    1650     }
    1651     findMin=list();
    1652   }
    1653   /* compute a V_d-strict Groebner basis (and resolution of J_D if
    1654   Syzstring!='Vdres') for J_D*/
     4235        {
     4236          m_d[i]=m_a[i];
     4237        }
     4238      findMin=list();
     4239    }
     4240  /* compute a V_d-strict Groebner basis (and resolution of J_D if
     4241     Syzstring!='Vdres') for J_D*/
    16554242  if (Syzstring=="Vdres")
    1656   {
    1657     J_D=VdStrictGB(J_D,d,m_d);
    1658   }
     4243    {
     4244      J_D=VdStrictGB(J_D,d,m_d);
     4245    }
    16594246  else
    1660   {   
    1661     def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d);
    1662     setring HomWeyl;
    1663     matrix J_D=fetch(B,J_D);
    1664     J_D=nHomogenize(J_D);
    1665     def ringofSyz=Sres(groebner(transpose(J_D)),lengthofres);
    1666     setring ringofSyz;
    1667     matrix J_D=transpose(matrix(RES[2]));
    1668     J_D=subst(J_D,h,1);
    1669     logens=ncols(J_D)+1;
    1670     list rofD;
    1671     for (i=3; i<=d+sizeL+1; i++)
    1672     {
    1673       if (size(RES)>=i)
    1674       {
    1675         if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
    1676         {
    1677           rofD[i-2]=(matrix(RES[i]));// resolution for J_D
    1678           if (i==3)
    1679           {
    1680             if (nrows(rofD[i-2])-logens+1!=nrows(J_D))
    1681             {
    1682               nr=nrows(J_D)+logens-1;
    1683               rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
    1684             }
    1685           }
    1686           if (i!=3)
    1687           {
    1688             if (nrows(rofD[i-2])-logens+1!=nrows(rofD[i-3]))
    1689             {
    1690               nr=nrows(rofD[i-3])+logens-1;
    1691               rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
    1692             }
    1693           }
    1694           i1=intvec(logens..nrows(rofD[i-2]));
    1695           i2=intvec(1..ncols(rofD[i-2]));
    1696           rofD[i-2]=transpose(submat(rofD[i-2],i1,i2));
    1697           logens=logens+ncols(rofD[i-2]);
    1698           rofD[i-2]=subst(rofD[i-2],h,1);     
    1699         }
    1700         else
    1701         {
    1702           rofD[i-2]=list();
    1703         }
    1704       }
    1705       else
    1706       {
    1707         rofD[i-2]=list();
    1708       }
    1709     }
    1710     if(size(rofD[1])==0)
    1711     {
    1712       omitemptylist=1;
    1713     }
    1714     setring B;
    1715     J_D=fetch(ringofSyz,J_D);
    1716     if (omitemptylist!=1)
    1717     {
    1718       list rofD=fetch(ringofSyz,rofD);
    1719     }
    1720     omitemptylist=0;
    1721     kill HomWeyl;
    1722     kill ringofSyz;
    1723   }
    1724   /* compute new matrices for J_A and J_B  such that their rows form a V_d-strict 
    1725   Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and
    1726   nrows(J_B)=nrows(J_A)+nrows(J_C)*/
     4247    {
     4248      def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d);
     4249      setring HomWeyl;
     4250      matrix J_D=fetch(B,J_D);
     4251      J_D=nHomogenize(J_D);
     4252      def ringofSyz=Sres(transpose(J_D),lengthofres);
     4253      setring ringofSyz;
     4254      matrix J_D=transpose(matrix(RES[2]));
     4255      J_D=subst(J_D,h,1);
     4256      logens=ncols(J_D)+1;
     4257      list rofD;
     4258      for (i=3; i<=d+sizeL+1; i++)
     4259        {
     4260          if (size(RES)>=i)
     4261            {
     4262              if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
     4263                {
     4264                  rofD[i-2]=(matrix(RES[i]));// resolution for J_D
     4265                  if (i==3)
     4266                    {
     4267                      if (nrows(rofD[i-2])-logens+1!=nrows(J_D))
     4268                        {
     4269                          nr=nrows(J_D)+logens-1;
     4270                          rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
     4271                        }
     4272                    }
     4273                  if (i!=3)
     4274                    {
     4275                      if (nrows(rofD[i-2])-logens+1!=nrows(rofD[i-3]))
     4276                        {
     4277                          nr=nrows(rofD[i-3])+logens-1;
     4278                          rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
     4279                        }
     4280                    }
     4281                  i1=intvec(logens..nrows(rofD[i-2]));
     4282                  i2=intvec(1..ncols(rofD[i-2]));
     4283                  rofD[i-2]=transpose(submat(rofD[i-2],i1,i2));
     4284                  logens=logens+ncols(rofD[i-2]);
     4285                  rofD[i-2]=subst(rofD[i-2],h,1);
     4286                }
     4287              else
     4288                {
     4289                  rofD[i-2]=list();
     4290                }
     4291            }
     4292          else
     4293            {
     4294              rofD[i-2]=list();
     4295            }
     4296        }
     4297      if(size(rofD[1])==0)
     4298        {
     4299          omitemptylist=1;
     4300        }
     4301      setring B;
     4302      J_D=fetch(ringofSyz,J_D);
     4303      if (omitemptylist!=1)
     4304        {
     4305          list rofD=fetch(ringofSyz,rofD);
     4306        }
     4307      omitemptylist=0;
     4308      kill HomWeyl;
     4309      kill ringofSyz;
     4310    }
     4311  /* compute new matrices for J_A and J_B  such that their rows form a V_d-strict
     4312     Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and
     4313     nrows(J_B)=nrows(J_A)+nrows(J_C)*/
    17274314  J_DF=transpose(storePinew[1]);
    17284315  for (i=2; i<=nrows(J_F); i++)
    1729   {
    1730     J_DF=concat(J_DF,transpose(storePinew[i]));
    1731   }
     4316    {
     4317      J_DF=concat(J_DF,transpose(storePinew[i]));
     4318    }
    17324319  J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF));
    17334320  J_DFC=transpose(storePi[1]);
    17344321  for (i=2; i<=nrows(J_C); i++)
    1735   {
    1736     J_DFC=concat(J_DFC,transpose(storePi[i]));
    1737   }
     4322    {
     4323      J_DFC=concat(J_DFC,transpose(storePi[i]));
     4324    }
    17384325  J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC));
    17394326  intvec m_b=m_a,v;
    … …  
    17444331  matrix g_AB=concat(unitmat(ncols(J_DF)),zero1);
    17454332  matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C))));
    1746   list out; 
     4333  list out;
    17474334  out[1]=list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F));
    17484335  out[1]=out[1]+list(list(m_d),list(m_a),list(m_f));
    … …  
    17504337  out[2]=out[2]+list(list(m_a),list(m_b),list(v));
    17514338  if (Syzstring=="Sres")
    1752   {
    1753     out[3]=rofD;
    1754     out[4]=rofF;     
    1755   }
     4339    {
     4340      out[3]=rofD;
     4341      out[4]=rofF;
     4342    }
    17564343  return(out);
    17574344}
    … …  
    17614348static proc VdStrictDoubleComplexes(list L,int d,string Syzstring)
    17624349{
    1763   /* We compute  V_d-strict resolutions over the V_d-strict short exact pieces from 
    1764   the procedure shortExactPiecesToVdStrict.
    1765   We use Algorithms 3.14 and 3.15 in [W2]*/
     4350  /* We compute  V_d-strict resolutions over the V_d-strict short exact pieces from
     4351     the procedure shortExactPiecesToVdStrict.
     4352     We use Algorithms 3.14 and 3.15 in [W2]*/
    17664353  int i,k,c,j,l,totaldeg,comparedegs,SBcom,verk;
    1767   intvec fordegs; 
     4354  intvec fordegs;
    17684355  intvec n_b,i1,i2;
    17694356  matrix rem,forML,subm,zerom,unitm,subm2;
    … …  
    17774364  list forhW;
    17784365  if (Syzstring=="Sres")
    1779   {
    1780     /*we already computed some of the resolutions in the procedure
    1781     shortExactPiecesToVdStrict*/
    1782     matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew;
    1783     for (k=1; k<=(size(L)+d-1); k++)
    1784     {
     4366    {
     4367    /*we already computed some of the resolutions in the procedure
     4368      shortExactPiecesToVdStrict*/
     4369      matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew;
     4370      for (k=1; k<=(size(L)+d-1); k++)
     4371        {
     4372          L[1][1][1][k+1]=list();
     4373          L[1][1][2][k+1]=list();
     4374          L[1][1][6][k+1]=list();
     4375        }
     4376      L[1][1][6][size(L)+d+1]=list();
     4377      matrix mem;
     4378      for (i=2; i<=d+size(L)+1; i++)
     4379        {;
     4380          v=0;
     4381          if(size(L[1][1][3][i-1])!=0)
     4382            {
     4383              if(i!=d+size(L)+1)
     4384                {
     4385                  /*horizontal differential*/
     4386                  L[1][1][4][i-1]=unitmat(nrows(L[1][1][3][i-1]));
     4387                }
     4388              for (j=1; j<=nrows(L[1][1][3][i-1]); j++)
     4389                {
     4390                  mem=submat(L[1][1][3][i-1],j,intvec(1..ncols(L[1][1][3][i-1])));
     4391                  v[j]=VdDeg(mem,d,L[1][1][7][i-1]);
     4392                }
     4393              L[1][1][7][i]=v;//new shift vector
     4394              L[1][1][8][i]=v;
     4395              L[1][2][6][i]=v;
     4396            }
     4397          else
     4398            {
     4399              if (i!=d+size(L)+1)
     4400                {
     4401                  L[1][1][4][i-1]=list();
     4402                }
     4403              L[1][1][7][i]=list();
     4404              L[1][1][8][i]=list();
     4405              L[1][2][6][i]=list();
     4406            }
     4407        }
     4408      if (size(L[1][1][3][d+size(L)])!=0)
     4409        {
     4410          /*horizontal differential*/
     4411          L[1][1][4][d+size(L)]=unitmat(nrows(L[1][1][3][d+size(L)]));
     4412        }
     4413      else
     4414        {
     4415          L[1][1][4][d+size(L)]=list();
     4416        }
     4417      for (k=1; k<size(L); k++)
     4418        {
     4419          /* We build a V_d-strict resolution for coker(L[k][2][1][1])->
     4420             coker(L[k][2][3][1])->coker(L[k][2][5][1]) using the resolution
     4421             obtained for coker(L[k][1][3][1]).
     4422             L[k][2][i][j] will be the jth module in the resolution of L[k][2][i][1]
     4423             for i=1,3,5.
     4424             L[k][2][i+5][j] will be the jth  shift vector in the resolution of
     4425             L[k][2][i][1](this holds also for the case Syzstring=="Vdres")*/
     4426          for (i=2; i<=d+size(L); i++)
     4427            {
     4428              v=0;
     4429              if (size(L[k][2][5][i-1])!=0)
     4430                {
     4431                  for (j=1; j<=nrows(L[k][2][5][i-1]); j++)
     4432                    {
     4433                      i1=intvec(1..ncols(L[k][2][5][i-1]));
     4434                      mem=submat(L[k][2][5][i-1],j,i1);
     4435                      v[j]=VdDeg(mem,d,L[k][2][8][i-1]);
     4436                    }
     4437                  /*next shift vector in th resolution of coker(L[k][2][5][1])*/
     4438                  L[k][2][8][i]=v;
     4439                }
     4440              else
     4441                {
     4442                  L[k][2][8][i]=list();
     4443                }
     4444              /* we build step by step a resolution for coker(L[k][2][5][1]) using
     4445                 the resolutions of coker(L[k][2][1][1]) and coker(L[k][2][5][1])*/
     4446              if (size(L[k][2][5][i])!=0)
     4447                {
     4448                  if (size(L[k][2][1][i])!=0 or size(L[k][2][1][i-1])!=0)
     4449                    {
     4450                      L[k][2][3][i]=transpose(syz(transpose(L[k][2][3][i-1])));
     4451                      nr= nrows(L[k][2][1][i-1]);
     4452                      nc=ncols(L[k][2][5][i]);
     4453                      Pold=matrixLift(L[k][2][3][i]*prodr(nr,nc), L[k][2][5][i]);
     4454                      matrix Pi[1][ncols(L[k][2][3][i])];
     4455                       for (l=1; l<=nrows(L[k][2][5][i]); l++)
     4456                        {
     4457                          for (j=1; j<=nrows(L[k][2][3][i]); j++)
     4458                            {
     4459                              i2=intvec(1..ncols(L[k][2][3][i]));
     4460                              Pi=Pi+Pold[l,j]*submat(L[k][2][3][i],j,i2);
     4461                            }
     4462                          if (l==1)
     4463                            {
     4464                              Picombined=transpose(Pi);
     4465                            }
     4466                          else
     4467                            {
     4468                              Picombined=concat(Picombined,transpose(Pi));
     4469                            }
     4470                          Pi=0;
     4471                        }
     4472                       kill Pi;
     4473                       Picombined=transpose(Picombined);
     4474                       if (size(L[k][2][1][i])!=0)
     4475                        {
     4476                          if (i==2)
     4477                            {
     4478                              containsndeg=(0:ncols(L[k][2][1][1]));
     4479                            }
     4480                          containsndeg=nDeg(L[k][2][1][i-1],containsndeg);
     4481                          forhW=list(L[k][2][6][i],containsndeg);
     4482                          def HomWeyl=makeHomogenizedWeyl(n,forhW);
     4483                          setring HomWeyl;
     4484                          list L=fetch(B,L);
     4485                          matrix M=L[k][2][1][i];
     4486                          module Mmod;
     4487                          list forM=nHomogenize(M,containsndeg,1);
     4488                          M=forM[1];
     4489                          totaldeg=forM[2];
     4490                          kill forM;
     4491                          matrix Maorig=fetch(B,Picombined);
     4492                          matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M)));
     4493                          matrix mem,subm,zerom;
     4494                          matrix Pinew;
     4495                          M=transpose(M);
     4496                          SBcom=0;
     4497                          for (l=1; l<=nrows(Ma); l++)
     4498                            {
     4499                              zerom=matrix(0,1,(ncols(Maorig)-ncols(Ma)));
     4500                              i1=(ncols(Ma)+1..ncols(Maorig));
     4501                              if (submat(Maorig,l,i1)==zerom)
     4502                                {
     4503                                  for (cc=1; cc<=ncols(Ma); cc++)
     4504                                    {
     4505                                      Maorig[l,cc]=0;
     4506                                    }
     4507                                }
     4508                              i2=(ncols(Ma)+1..ncols(Maorig));
     4509                              i1=(1..ncols(Ma));
     4510                              if (VdDeg(submat(Maorig,l,i1),d,L[k][2][6][i])>
     4511                                  VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]) and
     4512                                  submat(Maorig,l,i1)!=matrix(0,1,ncols(Ma)))
     4513                                {
     4514                                  /*V_d-Grad is to big--> we make it smaller using
     4515                                    Vdnormal form computations*/
     4516                                  if (SBcom==0)
     4517                                  {
     4518                                    Mmod=slimgb(M);
     4519                                    M=Mmod;
     4520                                    SBcom=1;
     4521                                  }
     4522                                  //print("Reduzierung des V_d-Grades(Stelle1)");
     4523                                  i2=(ncols(Ma)+1..ncols(Maorig));
     4524                                  vd1=VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]);
     4525                                  mem=submat(Ma,l,(1..ncols(Ma)));
     4526                                  mem=nHomogenize(mem,containsndeg);
     4527                                  mem=h^totaldeg*mem;
     4528                                  mem=transpose(mem);
     4529                                  mem=reduce(mem,Mod);//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     4530                                  matrix jt=transpose(subst(mem,h,1));
     4531                                  setring B;
     4532                                  matrix jt=fetch(HomWeyl,jt);
     4533                                  matrix need=fetch(HomWeyl,Maorig);
     4534                                  need=submat(need,l,(1..ncols(need)));
     4535                                  i1=L[k][2][6][i];
     4536                                  i2=L[k][2][8][i];
     4537                                  jt=VdNormalForm(need,L[k][2][1][i],d,i1,i2);
     4538                                  setring HomWeyl;
     4539                                  mem=fetch(B,jt);
     4540                                  mem=transpose(mem);
     4541                                  if (l==1)
     4542                                    {
     4543                                      Pinew=mem;
     4544                                    }
     4545                                  else
     4546                                    {
     4547                                      Pinew=concat(Pinew,mem);
     4548                                    }
     4549                                  vd2=VdDeg(transpose(mem),d,L[k][2][6][i]);
     4550                                  if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem)))
     4551                                    {//should not happen!!
     4552                                      //print("Reduzierung fehlgeschlagen!!(Stelle1)");
     4553                                    }
     4554                                }
     4555                              else
     4556                                {
     4557                                  if (l==1)
     4558                                    {
     4559                                      Pinew=transpose(submat(Ma,l,(1..ncols(Ma))));
     4560                                    }
     4561                                  else
     4562                                    {
     4563                                      subm=transpose(submat(Ma,l,(1..ncols(Ma))));
     4564                                      Pinew=concat(Pinew,subm);
     4565                                    }
     4566                                }
     4567                            }
     4568                          Pinew=subst(Pinew,h,1);
     4569                          Pinew=transpose(Pinew);
     4570                          setring B;
     4571                          Pinew=fetch(HomWeyl,Pinew);
     4572                          kill HomWeyl;
     4573                          L[k][2][3][i]=concat(Pinew,L[k][2][5][i]);
     4574                          subm=transpose(L[k][2][3][i]);
     4575                          subm=concat(transpose(L[k][2][1][i]),subm);
     4576                          L[k][2][3][i]=transpose(subm);
     4577                        }
     4578                      else
     4579                        {
     4580                          L[k][2][3][i]=Picombined;
     4581                        }
     4582                      L[k+1][1][1][i]=L[k][2][5][i];
     4583                      nr=nrows(L[k][2][1][i-1]);
     4584                      nc=ncols(L[k][2][5][i]);
     4585                      L[k][2][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
     4586                      L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nc);
     4587                      v=L[k][2][6][i],L[k][2][8][i];
     4588                      L[k][2][7][i]=v;
     4589                      L[k+1][1][6][i]=L[k][2][8][i];
     4590                    }
     4591                  else
     4592                    {
     4593                      L[k][2][3][i]=L[k][2][5][i];
     4594                      L[k][2][2][i]=list();
     4595                      L[k][2][7][i]=L[k][2][8][i];
     4596                      L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1]));
     4597                      L[k+1][1][6][i]=L[k][2][8][i];
     4598                      L[k+1][1][1][i]=L[k][2][5][i];
     4599                    }
     4600                }
     4601              else
     4602                {
     4603                  if (size(L[k][2][1][i])!=0)
     4604                    {
     4605                      if (size(L[k][2][5][i-1])!=0)
     4606                        {
     4607                          nr=nrows(L[k][2][5][i-1]);
     4608                          L[k][2][3][i]=concat(L[k][2][1][i],matrix(0,1,nr));
     4609                          v=L[k][2][6][i],L[k][2][8][i];
     4610                          L[k][2][7][i]=v;
     4611                          nc=nrows(L[k][2][1][i-1]);
     4612                          L[k][2][2][i]=concat(unitmat(nc),matrix(0,nc,nr));
     4613                          L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nr);
     4614                        }
     4615                      else
     4616                        {
     4617                          L[k][2][3][i]=L[k][2][1][i];
     4618                          L[k][2][7][i]=L[k][2][6][i];
     4619                          L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1]));
     4620                          L[k][2][4][i]=list();
     4621                        }
     4622                      L[k+1][1][1][i]=L[k][2][5][i];
     4623                      L[k+1][1][6][i]=L[k][2][8][i];
     4624                    }
     4625                  else
     4626                    {
     4627                      L[k][2][3][i]=list();
     4628                      if (size(L[k][2][6][i])!=0)
     4629                        {
     4630                          if (size(L[k][2][8][i])!=0)
     4631                            {
     4632                              v=L[k][2][6][i],L[k][2][8][i];
     4633                              L[k][2][7][i]=v;
     4634                              nr=nrows(L[k][2][1][i-1]);
     4635                              nc=nrows(L[k][2][5][i-1]);
     4636                              L[k][2][2][i]=concat(unitmat(nc),matrix(0,nr,nc));
     4637                              L[k][2][4][i]=prodr(nr,nrows(L[k][2][5][i-1]));
     4638                            }
     4639                          else
     4640                            {
     4641                              L[k][2][7][i]=L[k][2][6][i];
     4642                              L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1]));
     4643                              L[k][2][4][i]=list();
     4644                            }
     4645                        }
     4646                      else
     4647                        {
     4648                          if (size(L[k][2][8][i])!=0)
     4649                            {
     4650                              L[k][2][7][i]=L[k][2][8][i];
     4651                              L[k][2][2][i]=list();
     4652                              L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1]));
     4653                            }
     4654                          else
     4655                            {
     4656                              L[k][2][7][i]=list();
     4657                              L[k][2][2][i]=list();
     4658                              L[k][2][4][i]=list();
     4659                            }
     4660                        }
     4661                      L[k+1][1][1][i]=L[k][2][5][i];
     4662                      L[k+1][1][6][i]=L[k][2][8][i];
     4663                    }
     4664                }
     4665            }
     4666          i=d+size(L)+1;
     4667          v=0;
     4668          if (size(L[k][2][5][i-1])!=0)
     4669            {
     4670              for (j=1; j<=nrows(L[k][2][5][i-1]); j++)
     4671                {
     4672                  mem=submat(L[k][2][5][i-1],j,intvec(1..ncols(L[k][2][5][i-1])));
     4673                  v[j]=VdDeg(mem,d,L[k][2][8][i-1]);
     4674                }
     4675              L[k][2][8][i]=v;
     4676              if (size(L[k][2][6][i])!=0)
     4677                {
     4678                  v=L[k][2][6][i],L[k][2][8][i];
     4679                  L[k][2][7][i]=v;
     4680                }
     4681              else
     4682                {
     4683                  L[k][2][7][i]=L[k][2][8][i];
     4684                }
     4685            }
     4686          else
     4687            {
     4688              L[k][2][8][i]=list();
     4689              L[k][2][7][i]=L[k][2][6][i];
     4690            }
     4691          L[k+1][1][6][i]=L[k][2][8][i];
     4692          /* now we build V_d-strict resolutions for the sequences
     4693             coker(L[k+1][1][1][1])->coker(L[k+1][1][3][1])->coker(L[k+1][1][5][i])
     4694             using the resolutions  for coker(L[k][2][5][1]) we just obtained
     4695             (works exactly the same as above)*/
     4696          for (i=2; i<=d+size(L); i++)
     4697            {
     4698              v=0;
     4699              if (size(L[k+1][1][5][i-1])!=0)
     4700                {
     4701                  for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++)
     4702                    {
     4703                      i1=intvec(1..ncols(L[k+1][1][5][i-1]));
     4704                      mem=submat(L[k+1][1][5][i-1],j,i1);
     4705                      v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]);
     4706                    }
     4707                  L[k+1][1][8][i]=v;
     4708                }
     4709              else
     4710                {
     4711                  L[k+1][1][8][i]=list();
     4712                }
     4713              if (size(L[k+1][1][5][i])!=0)
     4714                {
     4715                  if (size(L[k+1][1][1][i])!=0 or size(L[k+1][1][1][i-1])!=0)
     4716                    {
     4717                      L[k+1][1][3][i]=transpose(syz(transpose(L[k+1][1][3][i-1])));
     4718                      nr=nrows(L[k+1][1][1][i-1]);
     4719                      nc=ncols(L[k+1][1][5][i]);
     4720                      Pold=matrixLift(L[k+1][1][3][i]*prodr(nr,nc),L[k+1][1][5][i]);
     4721                      matrix Pi[1][ncols(L[k+1][1][3][i])];
     4722                      for (l=1; l<=nrows(L[k+1][1][5][i]); l++)
     4723                        {
     4724                          for (j=1; j<=nrows(L[k+1][1][3][i]); j++)
     4725                            {
     4726                              i2=intvec(1..ncols(L[k+1][1][3][i]));
     4727                              Pi=Pi+Pold[l,j]*submat(L[k+1][1][3][i],j,i2);
     4728                            }
     4729                          if (l==1)
     4730                            {
     4731                              Picombined=transpose(Pi);
     4732                            }
     4733                          else
     4734                            {
     4735                              Picombined=concat(Picombined,transpose(Pi));
     4736                            }
     4737                          Pi=0;
     4738                        }
     4739                      kill Pi;
     4740                      Picombined=transpose(Picombined);
     4741                      if(size(L[k+1][1][1][i])!=0)
     4742                        {
     4743                          if (i==2)
     4744                            {
     4745                              containsndeg=(0:ncols(L[k+1][1][1][i-1]));
     4746                            }
     4747                          containsndeg=nDeg(L[k+1][1][1][i-1],containsndeg);
     4748                          forhW=list(L[k+1][1][6][i], containsndeg);
     4749                          def HomWeyl=makeHomogenizedWeyl(n,forhW);
     4750                          setring HomWeyl;
     4751                          list L=fetch(B,L);
     4752                          matrix M=L[k+1][1][1][i];
     4753                          module Mmod;
     4754                          list forM=nHomogenize(M,containsndeg,1);
     4755                          M=forM[1];
     4756                          totaldeg=forM[2];
     4757                          kill forM;
     4758                          matrix Maorig=fetch(B,Picombined);
     4759                          matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M)));
     4760                          Ma=nHomogenize(Ma,containsndeg);
     4761                          matrix mem,subm,zerom,subm2;
     4762                          matrix Pinew;
     4763                          M=transpose(M);
     4764                          SBcom=0;
     4765                          for (l=1; l<=nrows(Ma); l++)
     4766                            {
     4767                              i2=(ncols(Ma)+1..ncols(Maorig));
     4768                              nc=ncols(Maorig)-ncols(Ma);
     4769                              if (submat(Maorig,l,i2)==matrix(0,1,nc))
     4770                                {
     4771                                  for (cc=1; cc<=ncols(Ma); cc++)
     4772                                    {
     4773                                      Maorig[l,cc]=0;
     4774                                    }
     4775                                }
     4776                              i1=(1..ncols(Ma));
     4777                              i2=L[k+1][1][8][i];
     4778                              subm=submat(Maorig,l,i1);
     4779                              subm2=submat(Maorig,l,(ncols(Ma)+1..ncols(Maorig)));
     4780                              if (VdDeg(subm,d,L[k+1][1][6][i])>VdDeg(subm2,d,i2)
     4781                                  and subm!=matrix(0,1,ncols(Ma)))
     4782                                {
     4783                                  //print("Reduzierung des Vd-Grades (Stelle2)");
     4784                                  if (SBcom==0)
     4785                                    {
     4786                                      Mmod=slimgb(M);
     4787                                      M=Mmod;
     4788                                      SBcom=1;
     4789                                    }
     4790                                  vd1=VdDeg(subm2,d,L[k+1][1][8][i]);
     4791                                  mem=submat(Ma,l,(1..ncols(Ma)));
     4792                                  mem=nHomogenize(mem,containsndeg);
     4793                                  mem=h^totaldeg*mem;
     4794                                  mem=transpose(mem);
     4795                                  mem=reduce(mem,Mmod);
     4796                                  if (l==1)
     4797                                    {
     4798                                      Pinew=mem;
     4799                                    }
     4800                                  else
     4801                                    {
     4802                                      Pinew=concat(Pinew,mem);
     4803                                    }
     4804                                  vd2=VdDeg(transpose(mem),d,L[k+1][1][6][i]);
     4805                                  if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem)))
     4806                                    {//should not happen
     4807                                      //print("Reduzierung fehlgeschlagen!!!!(Stelle2)");
     4808                                    }
     4809                                }
     4810                              else
     4811                                {
     4812                                  if (l==1)
     4813                                    {
     4814                                      Pinew=transpose(submat(Ma,l,(1..ncols(Ma))));
     4815                                    }
     4816                                  else
     4817                                    {
     4818                                      subm=transpose(submat(Ma,l,(1..ncols(Ma))));
     4819                                      Pinew=concat(Pinew,subm);
     4820                                    }
     4821                                }
     4822                            }
     4823                          Pinew=subst(Pinew,h,1);
     4824                          Pinew=transpose(Pinew);
     4825                          setring B;
     4826                          Pinew=fetch(HomWeyl,Pinew);
     4827                          kill HomWeyl;
     4828                          L[k+1][1][3][i]=concat(Pinew,L[k+1][1][5][i]);
     4829                          subm=transpose(L[k+1][1][1][i]);
     4830                          subm2=transpose(L[k+1][1][3][i]);
     4831                          L[k+1][1][3][i]=transpose(concat(subm,subm2));
     4832                        }
     4833                      else
     4834                        {
     4835                          L[k+1][1][3][i]=Picombined;
     4836                        }
     4837                      L[k+1][2][1][i]=L[k+1][1][3][i];
     4838                      nr=nrows(L[k+1][1][1][i-1]);
     4839                      nc=ncols(L[k+1][1][5][i]);
     4840                      L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
     4841                      L[k+1][1][4][i]=prodr(nr,nc);
     4842                      v=L[k+1][1][6][i],L[k+1][1][8][i];
     4843                      L[k+1][1][7][i]=v;
     4844                      L[k+1][2][6][i]=L[k+1][1][7][i];
     4845                    }
     4846                  else
     4847                    {
     4848                      L[k+1][1][3][i]=L[k+1][1][5][i];
     4849                      L[k+1][1][2][i]=list();
     4850                      L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1]));
     4851                      L[k+1][1][7][i]=L[k+1][1][8][i];
     4852                      L[k+1][2][6][i]=L[k+1][1][7][i];
     4853                      L[k+1][2][1][i]=L[k+1][1][3][i];
     4854                    }
     4855                }
     4856              else
     4857                {
     4858                  if (size(L[k+1][1][1][i])!=0)
     4859                    {
     4860                      if (size(L[k+1][1][5][i-1])!=0)
     4861                        {
     4862                          zerom=matrix(0,1,nrows(L[k+1][1][5][i-1]));
     4863                          L[k+1][1][3][i]=concat(L[k+1][1][1][i],zerom);
     4864                          v=L[k+1][1][6][i],L[k+1][1][8][i];
     4865                          L[k+1][1][7][i]=v;
     4866                          nr=nrows(L[k+1][1][1][i-1]);
     4867                          nc=nrows(L[k+1][1][5][i-1]);
     4868                          L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
     4869                          L[k+1][1][4][i]=prodr(nr,nc);
     4870                        }
     4871                      else
     4872                        {
     4873                          L[k+1][1][3][i]=L[k+1][1][1][i];
     4874                    &