# Changeset 0e8a5a in git

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

### Legend:

Unmodified
 r89b2d28 ///////////////////////////////////////////////////////////////////////////////// version="version derham.lib 4.0.0.0 Jun_2013 "; // $Id$ /////////////////////////////////////////////////////////////////////////////// version=""; category="Noncommutative"; info=" LIBRARY:  derham.lib      Computation of deRham cohomology AUTHORS:  Cornelia Rottner, rottner@mathematik.uni-kl.de OVERVIEW: A library for computing the de Rham cohomology of complements of complex affine OVERVIEW: PROCEDURES: "; LIB "nctools.lib"; LIB "matrix.lib"; LIB "qhmoduli.lib"; LIB "general.lib"; LIB "dmod.lib"; LIB "bfun.lib"; LIB "dmodapp.lib"; LIB "poly.lib"; ///////////////////////////////////////////////////////////////////////////// static proc divdr(matrix m, matrix n) { m=transpose(m); n=transpose(n); matrix con=concat(m,n); matrix s=syz(con); s=submat(s,1..ncols(m),1..ncols(s)); s=transpose(compress(s)); return(s); } ///////////////////////////////////////////////////////////////////////////// static proc matrixlift(matrix M, matrix N) { // option(noreturnSB); matrix l=transpose(lift(transpose(M),transpose(N))); return(l); } /////////////////////////////////////////////////////////////////////////////// proc shortexactpieces(list #) { matrix  Bnew= divdr(#[2],#[3]); matrix Bold=Bnew; matrix Z=divdr(Bnew,#[1]); list bzh; list zcb; bzh=list(list(),list(),Z,unitmat(ncols(Z)),Z); zcb=(Z, Bnew, #[1], unitmat(ncols(#[1])), Bnew); list sep; sep[1]=(list(bzh,zcb)); int i; list out; for (i=3; i<=size(#)-2; i=i+2) { out=bzhzcb(Bold, #[i-1] , #[i], #[i+1], #[i+2]); sep[size(sep)+1]=out[1]; Bold=out[2]; } bzh=(divdr(#[size(#)-2], #[size(#)-1]),#[size(#)-2], #[size(#)-1],unitmat(ncols(#[size(#)-1])),transpose(concat(transpose(#[size(#)-2]),transpose(#[size(#)-1])))); zcb=(#[size(#)-1], unitmat(ncols(#[size(#)-1])), #[size(#)-1],list(),list()); sep[size(sep)+1]=list(bzh,zcb); return(sep); } //////////////////////////////////////////////////////////////////////////////////////// static proc bzhzcb (matrix Bold, matrix f0, matrix C1, matrix f1,matrix C2) { matrix Bnew=divdr(f1,C2); matrix Z= divdr(Bnew,C1); matrix lift1= matrixlift(Bnew,f0); list bzh=(Bold, lift1, Z, unitmat(ncols(Z)), transpose(concat(transpose(lift1),transpose(Z)))); list zcb=(Z, Bnew, C1, unitmat(ncols(C1)),Bnew); list out=(list(bzh, zcb), Bnew); return(out); } ////////////////////////////////////////////////////////////////////////////////////// proc VdstrictGB (matrix M, int d ,list #); "USAGE:VdstrictGB(M,d[,v]); M a matrix, d an integer, v an optional intvec(shift vector) RETURN:matrix M; the rows of M foem a Vd-strict Groebner basis for imM ASSUME:1<=d<=nvars(basering)/2; size(v)=ncols(M) " { if (M==matrix(0,nrows(M),ncols(M))) { return (matrix(0,1,ncols(M))); } def W =basering; int ncM=ncols(M); list Data=ringlist(W); Data[2]=list("nhv")+Data[2]; Data[3][3]=Data[3][1]; Data[3][1]=Data[3][2]; Data[3][2]=list("dp",intvec(1)); matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); Data[5]=re; int k; int l; Data[6]=transpose(concat(matrix(0,1,1),transpose(concat(matrix(0,1,1),Data[6])))); def Whom=ring(Data); setring Whom; matrix Mnew=imap(W,M); intvec v; if (size(#)!=0) { v=#[1]; } if (size(v) < ncM) { v=v,0:(ncM-size(v)); } Mnew=homogenize(Mnew, d, v); Mnew=transpose(Mnew); Mnew=std(Mnew); Mnew=subst(Mnew,nhv,1); Mnew=transpose(Mnew); setring W; M=imap(Whom,Mnew); return(M); } //////////////////////////////////////////////////////////////////////////////////// static proc Vdnormalform(matrix F, matrix M, int d, intvec v) { def W =basering; int c=ncols(M); F=submat(F,intvec(1..nrows(F)),intvec(1..c)); list Data=ringlist(W); Data[2]=list("nhv")+Data[2]; Data[3][3]=Data[3][1]; Data[3][1]=Data[3][2]; Data[3][2]=list("dp",intvec(1)); matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); Data[5]=re; int k; int l; matrix rep[size(Data[2])][size(Data[2])]; for (l=size(Data[2])-1;l>=1; l--) { for (k=l-1; k>=1;k--) { rep[k+1,l+1]=Data[6][k,l]; } } Data[6]=rep; def Whom=ring(Data); setring Whom; matrix Mnew=imap(W,M); Mnew=(homogenize(Mnew, d, v));//doppelte Berechung unnötig->muss noch geändert werden!!! matrix Fnew=imap(W,F); matrix Fb; for (l=1; l<=nrows(Fnew); l++) { Fb=homogenize(submat(Fnew,l,intvec(1..ncols(Fnew))),d,v); Fb=transpose(reduce(transpose(Fb),std(transpose(Mnew))));// doppelte Berechnung unnötig, unterdrückt aber Fehler meldung for (k=1; k<=ncols(Fnew);k++) { Fnew[l,k]=Fb[1,k]; } } Fnew=subst(Fnew,nhv,1); setring W; F=imap(Whom,Fnew); return(F); } /////////////////////////////////////////////////////////////////////////////// static proc homogenize (matrix M, int d, intvec v) { int l; poly f; int s; int i; intvec vnm;int kmin; list findmin; int n=(nvars(basering)-1) div 2; list rempoly; list remk; list rem1; list rem2; 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 { for (l=1; l<=ncols (M); l++) { f=M[k,l]; s=size(f); for (i=1; i<=s; i++) { vnm=leadexp(f); vnm=vnm[n+2..n+d+1]-vnm[2..d+1]; kmin=sum(vnm)+v[l]; rem1[size(rem1)+1]=lead(f); rem2[size(rem2)+1]=kmin; findmin=insert(findmin,kmin); f=f-lead(f); } rempoly[l]=rem1; remk[l]=rem2; rem1=list(); rem2=list(); } if (size(findmin)!=0) { kmin=Min(findmin); } for (l=1; l<=ncols(M); l++) { if (M[k,l]!=0) { M[k,l]=0; for (i=1; i<=size(rempoly[l]);i++) { M[k,l]=M[k,l]+nhv^(remk[l][i]-kmin)*rempoly[l][i]; } } } rempoly=list(); remk=list(); findmin=list(); } return(M); } ////////////////////////////////////////////////////////////////////////////////////// static proc soldr (matrix M, matrix N) { int n=nrows(M); int q=ncols(M); matrix S=concat(transpose(M),transpose(N)); def W=basering; list Data=ringlist(W); list Save=Data[3]; Data[3]=list(list("c",0),list("dp",intvec(1..nvars(W)))); def Wmod=ring(Data); setring Wmod; matrix Smod=imap(W,S); matrix E[q][1]; matrix Smod2; matrix Smodnew; option(returnSB); int i; int j; for (i=1;i<=q;i++) { E[i,1]=1; Smod2=concat(E,Smod); print (Smod2); Smod2=syz(Smod2); E[i,1]=0; for (j=1;j<=ncols(Smod2);j++) { if (Smod2[1,j]==1) { Smodnew=concat(Smodnew,(-1)*(submat(Smod2,intvec(2..n+1),j))); break; } } } Smodnew=transpose(submat(Smodnew,intvec(1..n),intvec(2..q+1))); setring W; matrix  Snew=imap(Wmod,Smodnew); return (Snew); } ///////////////////////////////////////////////////////////////////////////// proc toVdstrictsequence (list C,int n, intvec v) { matrix J_C=VdstrictGB(C[5],n,list(v)); matrix J_A=C[1]; matrix f_CB=C[4]; matrix f_ACB=transpose(concat(transpose(C[2]),transpose(f_CB))); matrix J_AC=divdr(f_ACB,C[3]); matrix P=matrixlift(J_AC * prodr(ncols(C[1]),ncols(C[5])) ,J_C); list storePi; matrix Pi[1][ncols(J_AC)]; int i;int j; for (i=1; i<=nrows(J_C); i++) { for (j=1; j<=nrows(J_AC);j++) { Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC))); } storePi[i]=Pi; Pi=0; } intvec m_a; list findMin; int comMin; for (i=1; i<=ncols(J_A); i++) { for (j=1; j<=size(storePi);j++) { if (storePi[j][1,i]!=0) { comMin=Vddeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v)-Vddeg(storePi[j][1,i],n,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { m_a[i]=Min(findMin); findMin=list(); } else { m_a[i]=0; } } matrix zero[ncols(J_A)][ncols(J_C)]; matrix g_AB=concat(unitmat(ncols(J_A)),zero); matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C))))); intvec m_b=m_a,v; J_A=VdstrictGB(J_A,n,m_a); J_AC=transpose(storePi[1]); for (i=2; i<= size(storePi); i++) { J_AC=concat(J_AC, transpose(storePi[i])); } J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC)); 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)); return (Vdstrict); } ///////////////////////////////////////////////////////////////////////// static proc prodr (int k, int l) { if (k==0) { matrix P=unitmat(l); return (P); } matrix O[l][k]; matrix P=transpose(concat(O,unitmat(l))); return (P); } ///////////////////////////////////////////////////////////////////////// proc Vddeg(matrix M, int d, intvec v, list #)//Aternative: in WHom Leadmonom ausrechnen!! "USAGE: Vddeg(M,d,v); M 1xr-matrix, d int, v intvec of size r RETURN:int; the Vd-degree of M " { int i;int j; int n=nvars(basering) div 2; intvec  e; int etoint; list findmax; int c=ncols(M); poly l; list positionpoly; list positionVd; for (i=1; i<=c; i++) { positionpoly[i]=list(); positionVd[i]=list(); while (M[1,i]!=0) { l=lead(M[1,i]); positionpoly[i][size(positionpoly[i])+1]=l; e=leadexp(l); e=-e[1..d]+e[n+1..n+d]; e=sum(e)+v[i]; etoint=e[1]; positionVd[i][size(positionVd[i])+1]=etoint; findmax[size(findmax)+1]=etoint; M[1,i]=M[1,i]-l; } } if (size(findmax)!=0) { int maxVd=Max(findmax); if (size(#)==0) { return (maxVd); } } else // M is 0-modul { return(int(0)); } l=0; for (i=c; i>=1; i--) { for (j=1; j<=size(positionVd[i]); j++) { if (positionVd[i][j]==maxVd) { l=l+positionpoly[i][j]; } } if (l!=0) { return (list(l,i)); } } } /////////////////////////////////////////////////////////////////////////////// proc toVdstrictsequences (list L,int d, intvec v) "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 RETURN: list of two lists; each lists contains Vd-strict exact sequences with corresponding shift vectors " { matrix J_F=L[1][5]; matrix J_D=L[1][1]; matrix f_FA=L[1][4]; matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA))); matrix J_DF=divdr(f_DFA,L[1][3]); matrix J_C=L[2][5]; matrix f_CB=L[2][4]; matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB))); matrix J_DFC=divdr(f_DFCB,L[2][3]); matrix P=matrixlift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C); list storePi; matrix Pi[1][ncols(J_DFC)]; int i; int j; for (i=1; i<=nrows(J_C); i++) { for (j=1; j<=nrows(J_DFC);j++) { Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC))); } storePi[i]=Pi; Pi=0; } intvec m_a; list findMin; list noMin; int comMin; for (i=1; i<=ncols(J_DF); i++) { for (j=1; j<=size(storePi);j++) { if (storePi[j][1,i]!=0) { comMin=Vddeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v)-Vddeg(storePi[j][1,i],d,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { m_a[i]=Min(findMin); findMin=list(); noMin[i]=0; } else { noMin[i]=1; } } if (size(m_a) < ncols(J_DF)) { m_a[ncols(J_DF)]=0; } intvec m_f=m_a[ncols(J_D)+1..size(m_a)]; J_F=VdstrictGB(J_F,d,m_f); P=matrixlift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F);// selbe Prinzip wie oben--> evtl auslagern list storePinew; matrix Pidf[1][ncols(J_DF)]; for (i=1; i<=nrows(J_F); i++) { for (j=1; j<=nrows(J_DF);j++) { Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF))); } storePinew[i]=Pidf; Pidf=0; } intvec m_d; for (i=1; i<=ncols(J_D); i++) { for (j=1; j<=size(storePinew);j++) { if (storePinew[j][1,i]!=0) { comMin=Vddeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f)-Vddeg(storePinew[j][1,i],d,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { if (noMin[i]==0) { m_d[i]=Min(insert(findMin,m_a[i])); m_a[i]=m_d[i]; } else { m_d[i]=Min(findMin); m_a[i]=m_d[i]; } } else { m_d[i]=m_a[i]; } findMin=list(); } J_D=VdstrictGB(J_D,d,m_d); J_DF=transpose(storePinew[1]); for (i=2; i<=nrows(J_F); i++) { J_DF=concat(J_DF,transpose(storePinew[i])); } J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF)); J_DFC=transpose(storePi[1]); for (i=2; i<=nrows(J_C); i++) { J_DFC=concat(J_DFC,transpose(storePi[i])); } J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC)); intvec m_b=m_a,v; matrix zero[ncols(J_D)][ncols(J_F)]; matrix g_DA=concat(unitmat(ncols(J_D)),zero); matrix g_AF=transpose(concat(transpose(zero),unitmat(ncols(J_F)))); matrix zero1[ncols(J_DF)][ncols(J_C)]; matrix g_AB=concat(unitmat(ncols(J_DF)),zero1); matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C)))); 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))); return(out), } /////////////////////////////////////////////////////////////////////////////////////////// proc shortexactpiecestoVdstrict(list C, int d,list #) { int s =size(C); if (size(#)==0) { intvec v=0:ncols(C[s][1][5]); } else { intvec v=#[1]; } list out; out[s]=list(toVdstrictsequence(C[s][1],d,v)); 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())); out[s][2][6]=list(out[s][1][7][1]); out[s][2][7]=list(out[s][2][6][1]); out[s][2][8]=list(list()); int i; for (i=s-1; i>=2; i--) { C[i][2][5]=out[i+1][1][1][1]; out[i]=toVdstrictsequences(C[i],d,out[i+1][1][6][1]); } out[1]=list(list()); out[1][2]=toVdstrictsequence(C[1][2],d,out[2][1][6][1]); out[1][1][3]=list(out[1][2][1][1]); out[1][1][5]=list(out[1][2][1][1]); out[1][1][4]=list(unitmat(ncols(out[1][1][3][1]))); out[1][1][7]=list(out[1][2][6][1]); out[1][1][8]=list(out[1][2][6][1]); out[1][1][1]=list(list()); out[1][1][2]=list(list()); out[1][1][6]=list(list()); list Hi; for (i=1; i<=size(out); i++) { Hi[i]=list(out[i][1][5][1],out[i][1][8][1]); } list outall; outall[1]=out; print (out); outall[2]=Hi; return(outall); } /////////////////////////////////////////////////////////////////////////////////////////// proc toVdstrict2x3complex(list L, int d, list #) { matrix rem; int i; int j; list J_A=list(list()); list J_B=list(list()); list J_C=list(list()); list g_AB=list(list()); list g_BC=list(list()); list n_a=list(list()); list n_b=list(list()); list n_c=list(list()); intvec n_b1; if (size(L[5])!=0) { intvec n_c1; for (i=1; i<=nrows(L[5]); i++) { rem=submat(L[5],i,intvec(1..ncols(L[5]))); n_c1[i]=Vddeg(rem,d, L[8]); } n_c[1]=n_c1; J_C[1]=transpose(syz(transpose(L[5]))); if (J_C[1]!=matrix(0,nrows(J_C[1]),ncols(J_C[1]))) { J_C[1]=VdstrictGB(J_C[1],d,n_c1); if (size(#[2])!=0) { n_a[1]=#[2]; n_b1=n_a[1],n_c[1]; n_b[1]=n_b1; matrix zero[nrows(L[1])][nrows(L[5])]; g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5]))); if (size(#[1])!=0) { J_A=#[1]; J_B=transpose(matrix(syz(transpose(L[3])))); matrix P=matrixlift(J_B[1] * prodr(nrows(L[1]),nrows(L[5])) ,J_C[1]); matrix Pi[1][ncols(J_B[1])]; matrix Picombined; for (i=1; i<=nrows(J_C[1]); i++) { for (j=1; j<=nrows(J_B[1]);j++) { Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1]))); } if (i==1) { Picombined=transpose(Pi); } else { Picombined=concat(Picombined,transpose(Pi)); } Pi=0; } Picombined=transpose(Picombined); 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)))); J_B[1]=transpose(concat(transpose(matrix(J_A[1],nrows(J_A[1]),ncols(J_B[1]))),transpose(Picombined))); g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); } else { J_B[1]=concat(matrix(0,nrows(J_C[1]),nrows(L[3])-nrows(L[5])),J_C[1]); g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); } } else { n_b=n_c[1]; J_B[1]=J_C[1]; g_BC=unitmat(ncols(J_C[1])); } } else { J_C=list(list()); if (size(#[2])!=0) { matrix zero[nrows(L[1])][nrows(L[5])]; g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); n_a[1]=#[2]; n_b1=n_a[1],n_c[1]; n_b[1]=n_b1; g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));; if (size(#[1])!=0) { J_A=#[1]; J_B=concat(J_A[1],matrix(0,nrows(J_A[1]),nrows(L[3])-nrows(L[1]))); } } else { n_b=n_c[1]; g_BC=unitmat(ncols(L[5])); } } } else { if (size(#[2])!=0) { n_a[1]=#[2]; n_b=n_a[1]; g_AB=unitmat(size(n_b[1])); if (size(#[1])!=0) { J_A=#[1]; J_B[1]=J_A[1]; } } } 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]); return (out); } ////////////////////////////////////////////////////////////////////////// proc Vdstrictdoublecompexes(list L, int d) { int i; int k; int c; int j; intvec n_b; matrix rem; matrix J_B; list store; int t=size(L)+nvars(basering) div 2-2; for (k=1; k<=(size(L)+nvars(basering) div 2-3); k++)// { L[1][1][1][k+1]=list(); L[1][1][2][k+1]=list(); L[1][1][6][k+1]=list(); if (size(L[1][1][3][k])!=0) { for (i=1; i<=nrows(L[1][1][3][k]); i++) { rem=submat(L[1][1][3][k],i,(1..ncols(L[1][1][3][k]))); n_b[i]=Vddeg(rem,d,L[1][1][7][k]); } J_B=transpose(syz(transpose(L[1][1][3][k]))); L[1][1][7][k+1]=n_b; L[1][1][8][k+1]=n_b; L[1][1][4][k+1]=unitmat(nrows(L[1][1][3][k])); if (J_B!=matrix(0,nrows(J_B),ncols(J_B))) { J_B=VdstrictGB(J_B,d,n_b); L[1][1][3][k+1]=J_B; L[1][1][5][k+1]=J_B; } else { L[1][1][3][k+1]=list(); L[1][1][5][k+1]=list(); } n_b=0; } else { L[1][1][3][k+1]=list(); L[1][1][5][k+1]=list(); L[1][1][7][k+1]=list(); L[1][1][8][k+1]=list(); L[1][1][4][k+1]=list(); } for (i=1; i1) { for (k=1; k<=Min(list(oldj,newj)); k++) { out[i-1][k][4]=matrix(out[i-1][k][4],nrows(out[i-1][k][4]),out[i][k][1]); } for (k=newj+1; k<=oldj; k++) { out[i-1][k]=delete(out[i-1][k],4); } } } return (out); } ////////////////////////////////////////////////////////////////////////////// proc totalcomplex(list L); { list out;intvec rem1;intvec v; list remsize; int emp; int i; int j; int c; int d; matrix M; int k; int l; int n=nvars(basering) div 2; list K; for (i=1; i<=n; i++) { K[i]=list(); } L=K+L; for (i=1; i<=size(L); i++) { emp=0; if (size(L[i])!=0) { out[3*i-2]=L[i][1][1]; v=L[i][1][1]; rem1=L[i][1][2]; emp=1; } else { out[3*i-2]=0; v=0; } for (j=i+1; j<=size(L); j++) { if (size(L[j])>=j-i+1) { out[3*i-2]=out[3*i-2]+L[j][j-i+1][1]; if (emp==0) { rem1=L[j][j-i+1][2]; emp=1; } else { rem1=rem1,L[j][j-i+1][2]; } v[size(v)+1]=L[j][j-i+1][1]; } else { v[size(v)+1]=0; } } out[3*i-1]=rem1; v[size(v)+1]=0; remsize[i]=v; } int o1; int o2; for (i=1; i<=size(L)-1; i++) { o1=1; o2=1; if (size(out[3*i-2])!=0) { o1=out[3*i-2]; } if (size(out[3*i+1])!=0) { o2=out[3*i+1]; } M=matrix(0,o1,o2); if (size(L[i])!=0) { if (size(L[i][1][4])!=0) { M=matrix(L[i][1][4],o1,o2); } } c=remsize[i][1]; // d=remsize[i+1][1]; for (j=i+1; j<=size(L); j++) { if (remsize[i][j-i+1]!=0) { for (k=c+1; k<=c+remsize[i][j-i+1]; k++) { for (l=d+1; l<=d+remsize[i+1][j-i];l++) { M[k,l]=L[j][j-i+1][3][(k-c),(l-d)]; } } d=d+remsize[i+1][j-i]; if (remsize[i+1][j-i+1]!=0) { for (k=c+1; k<=c+remsize[i][j-i+1]; k++) { for (l=d+1; l<=d+remsize[i+1][j-i+1];l++) { M[k,l]=L[j][j-i+1][4][k-c,l-d]; } } c=c+remsize[i][j-i+1]; } } else { d=d+remsize[i+1][j-i]; } } out[3*i]=M; d=0; c=0; } out[3*size(L)]=matrix(0,out[3*size(L)-2],1); return (out); } ///////////////////////////////////////////////////////////////////////////////////// proc toVdstrictfreecomplex(list L,list #) { def B=basering; int n=nvars(B) div 2+2; int d=nvars(B) div 2; intvec v; list out;list outall; int i;int j; matrix mem; int k; if (size(#)!=0) { for (i=1; i<=size(#); i++) { if (typeof(#[i])==intvec) { v=#[i]; } if (typeof(#[i])==int) { d=#[i]; } } } if (size(L)==2) { v=(0:ncols(L[1])); out[3*n-1]=v; out[3*n-2]=ncols(L[1]); out[3*n]=L[2]; out[3*n-3]=VdstrictGB(L[1],d,v); for (i=n-1; i>=1; i--) { out[3*i-2]=nrows(out[3*i]); v=0; for (j=1; j<=out[3*i-2]; j++) { mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); v[j]=Vddeg(mem,d, out[3*i+2]); } out[3*i-1]=v; if (i!=1) { out[3*i-3]=transpose(syz(transpose(out[3*i]))); if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3]))) { out[3*i-3]=VdstrictGB(out[3*i-3],d,out[3*i-1]); } else { out[3*i-3]=matrix(0,1,ncols(out[3*i-3])); out[3*i-4]=intvec(0); out[3*i-5]=int(0); for (j=i-2; j>=1; j--) { out[3*j]=matrix(0,1,1); out[3*j-1]=intvec(0); out[3*j-2]=int(0); } break; } } } outall[1]=out; outall[2]=list(list(out[3*n-3],out[3*n-1])); return(outall); } out=shortexactpieces(L); list rem; if (v!=intvec(0:size(v))) { rem=shortexactpiecestoVdstrict(out,d,v); } else { rem=shortexactpiecestoVdstrict(out,d); } out=Vdstrictdoublecompexes(rem[1],d); out=assemblingdoublecomplexes(out); out=totalcomplex(out); outall[1]=out; outall[2]=rem[2]; return (outall); } //////////////////////////////////////////////////////////////////////////////// proc derhamcohomology(list L) { def R=basering; int n=nvars(R);int le=2*size(L)+n-1; def W=makeWeyl(n); setring W; list man=ringlist(W); if (n==1) { man[2][1]="x(1)"; man[2][2]="D(1)"; def Wi=ring(man); setring Wi; kill W; def W=Wi; setring W; list man=ringlist(W); } man[2][size(man[2])+1]="s";; man[3][3]=man[3][2]; man[3][2]=list("dp",intvec(1)); matrix N=UpOneMatrix(size(man[2])); man[5]=N; matrix M[1][1]; man[6]=transpose(concat(transpose(concat(man[6],M)),M)); 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; for (i=1; i<=r; i++) { Fi[i]=list(); bernsteinpolys[i]=list(); annideals[i]=list(); subsets[i]=list(); maxnum[i]=list(); Fi[1][i]=L[i]; maxnum[1][i]=i; subsets[1][i]=intvec(i); } intvec v; for (i=2; i<=r; i++) { count=1; for (j=1; j<=size(Fi[i-1]);j++) { for (k=maxnum[i-1][j]+1; k<=r; k++) { maxnum[i][count]=k; v=subsets[i-1][j],k; subsets[i][count]=v; Fi[i][count]=lcm(Fi[i-1][j],L[k]);///////// count=count+1; } } } for (i=1; i<=r; i++) { for (j=1; j<=size(Fi[i]); j++) { bernsteinpolys[i][j]=bfct(Fi[i][j])[1]; for (k=1; k<=ncols(bernsteinpolys[i][j]); k++) { if (isInt(number(bernsteinpolys[i][j][k]))==1) { minint[size(minint)+1]=int(bernsteinpolys[i][j][k]); } } def D=Sannfs(Fi[i][j]); setring Ws; annideals[i][j]=fetch(D,LD); kill D; setring R; } } int m=Min(minint); list zw; for (i=1; i=0) { for (d=1; d<=size(truncatedcomplex[2*(i+1)-1][l][ideg+1]);d++) { if (leadmonom(lform)==truncatedcomplex[2*(i+1)-1][l][ideg+1][d][1]) { 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); break; } } } } form=form-lform; } } } } } } truncatedcomplex[2*i]=M; truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; } truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; if (truncatedcomplex[2*i-1]!=0) { truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1); } setring R; list truncatedcomplex=imap(W,truncatedcomplex); list derhamhom=findhomology(truncatedcomplex,le); return (derhamhom); } /////////////////////////////////// static proc sublist(list L, int m, int n) { list out; int i; int j; int count; for (i=m; i<=n; i++) { out[size(out)+1]=list(); for (j=1; j<=size(L[i]); j++) { count=count+1; out[size(out)][j]=list(L[i][j],count); } } list o=list(out,count); return(o); } ////////////////////////////////////////////////////////////////////////// static proc mysubset(intvec L, intvec M) { int i; int j=1; list position=(M[size(M)],(-1)^(size(L))); for (i=1; i<=size(L); i++) { if (L[i]!=M[j]) { if (L[i]!=M[j+1] or j!=i) { return (L[i],0); } else { position=(M[i],(-1)^(i-1)); j=j+i; } } j=j+1; } return (position); } //////////////////////////////////////////////////////////////////////////// proc globalbfun(list L) { int i; int j; def W=basering; int n=nvars(W) div 2; list G0; ideal I; for (j=1; j<=size(L); j++) { G0[j]=list(); for (i=1; i<=ncols(L[j][1]); i++) { G0[j][i]=I; } } list out; for (j=1; j<=size(L); j++) { for (i=1; i<=nrows(L[j][1]); i++) { out=Vddeg(submat(L[j][1],i,(1..ncols(L[j][1]))),n,L[j][2],1); G0[j][out[2]][size(G0[j][out[2]])+1]=(out[1]); } } list Data=ringlist(W); for (i=1; i<=n; i++) { Data[2][2*n+i]=Data[2][i]; Data[2][3*n+i]=Data[2][n+i]; Data[2][i]="v("+string(i)+")"; Data[2][n+i]="w("+string(i)+")"; } Data[3][1][1]="M"; intvec mord=(0:16*n^2); mord[1..2*n]=(1:2*n); mord[6*n+1..8*n]=(1:2*n); for (i=0; i<=2*n-2; i++) { mord[(3+i)*4*n-i]=-1; mord[(2*n+2+i)*4*n-2*n-i]=-1; } Data[3][1][2]=mord;//ordering mh????????? matrix Ones=UpOneMatrix(4*n); Data[5]=Ones; matrix con[2*n][2*n]; Data[6]=transpose(concat(con,transpose(concat(con,Data[6])))); def Wuv=ring(Data); setring Wuv; list G0=imap(W,G0); list G3; poly lterm;intvec lexp; list G1;  list G2; intvec e; intvec f; int  kapp; int k; int l; poly h; ideal I; for (l=1; l<=size(G0); l++) { G1[l]=list();  G2[l]=list(); G3[l]=list(); for (i=1; i<=size(G0[l]); i++) { for (j=1; j<=ncols(G0[l][i]);j++) { G0[l][i][j]=mhom(G0[l][i][j]); } for (j=1; j<=nvars(Wuv) div 4; j++) { G0[l][i][size(G0[l][i])+1]=1-v(j)*w(j); } G1[l][i]=std(G0[l][i]); G2[l][i]=I; G3[l][i]=list(); for (j=1; j<=ncols(G1[l][i]); j++) { e=leadexp(G1[l][i][j]); f=e[1..2*n]; if (f==intvec(0:(2*n))) { for (k=1; k<=n; k++) { kapp=-e[2*n+k]+e[3*n+k]; if (kapp>0) { G1[l][i][j]=(x(k)^kapp)*G1[l][i][j]; } if (kapp<0) { G1[l][i][j]=(D(k)^(-kapp))*G1[l][i][j]; } } G2[l][i][size(G2[l][i])+1]=G1[l][i][j]; G3[l][i][size(G3[l][i])+1]=list(); while (G1[l][i][j]!=0) { lterm=lead(G1[l][i][j]); G1[l][i][j]=G1[l][i][j]-lterm; lexp=leadexp(lterm); lexp=lexp[2*n+1..3*n]; G3[l][i][size(G3[l][i])][size(G3[l][i][size(G3[l][i])])+1]=list(lexp,leadcoef(lterm)); } } } } } ring r=0,(s(1..n)),dp; ideal I; map G3forr=Wuv,I; list G3=G3forr(G3); poly fs; poly gs; int a; list G4; for (l=1; l<=size(G3); l++) { G4[l]=list(); for (i=1; i<=size(G3[l]);i++) { G4[l][i]=I; for (j=1; j<=size(G3[l][i]); j++) { fs=0; for (k=1; k<=size(G3[l][i][j]); k++) { gs=1; for (a=1; a<=n; a++) { if (G3[l][i][j][k][1][a]!=0) { gs=gs*permutevar(list(G3[l][i][j][k][1][a]),a); } } gs=gs*G3[l][i][j][k][2]; fs=fs+gs; } G4[l][i]=G4[l][i],fs; } } } if (n==1) { ring rnew=0,t,dp; } else { ring rnew=0,(t,s(2..n)),dp; } ideal Iformap; Iformap[1]=t; poly forel=1; for (i=2; i<=n; i++) { Iformap[1]=Iformap[1]-s(i); Iformap[i]=s(i); forel=forel*s(i); } map rtornew=r,Iformap; list G4=rtornew(G4); list getintvecs=fetch(W,L); ideal J; option(redSB); for (l=1; l<=size(G4); l++) { J=1; for (i=1; i<=size(G4[l]); i++) { G4[l][i]=eliminate(G4[l][i],forel); G4[l][i]=subst(G4[l][i],t,t-getintvecs[l][2][i]); J=intersect(J,G4[l][i]); } G4[l]=poly(std(J)[1]); } list minmax=minmaxintroot(G4);//besser factorize nehmen // Fall: keine Nullstelle muss noch weiter beruecksichtigt werden return(minmax); } ////////////////////////////////////////////////////////////////////////// proc minmaxintroot(list L); { 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; possk[1]=list(); for (i=1; i<=size(L); i++) { d=1; f=L[i]; while (f!=0) { rest=lead(f); d=d*denominator(leadcoef(rest)); f=f-rest; } e=leadexp(rest); if (e[1]!=0) { rest=rest/(t^(e[1])); possk[1][size(possk[1])+1]=i; } a0[i]=int(absValue(d*rest)); } int m=Max(a0); for (i=2; i<=m+1; i++) { possk[i]=list(); } list allprimefac; for (i=1; i<=size(L); i++) { allprimefac=primefactors(a0[i]); alldiv=1; possk[2][size(possk[2])+1]=i; for (j=1; j<=size(allprimefac[1]); j++) { s=size(alldiv); for (k=1; k<=s; k++) { for (l=1; l<=allprimefac[2][j]; l++) { alldiv[size(alldiv)+1]=alldiv[k]*allprimefac[1][j]^l; possk[alldiv[size(alldiv)]+1][size(possk[alldiv[size(alldiv)]+1])+1]=i; } } } } int mink; int maxk; int indi; for (i=m+1; i>=1; i--) { if (size(possk[i])!=0) { for (j=1; j<=size(possk[i]); j++) { if (subst(L[possk[i][j]],t,(i-1))==0) { maxk=i-1; indi=1; break; } } } if (maxk!=0) { break; } } int indi2; for (i=m+1; i>=1; i--) { if (size(possk[i])!=0) { for (j=1; j<=size(possk[i]); j++) { if (subst(L[possk[i][j]],t,-(i-1))==0) { mink=-i+1; indi2=1; break; } } } if (mink!=0) { break; } } list mima=mink,maxk; if (indi==0) { if (indi2==0) { mima=list();//es gibt keine ganzzahlige NS } else { mima[2]=mima[1]; } } else { if (indi2==0) { mima[1]=mima[2]; } } return (mima); } /////////////////////////////////////////////////////// proc findhomology(list L, int le) { int li; matrix M; matrix N; matrix N1; matrix lift1; list out; int i; option (redSB); for (i=2; i<=size(L); i=i+2) { if (L[i-1]==0) { li=0; out[i div 2]=0; } else { if (li==0) { li=L[i-1]; N1=transpose(syz(transpose(L[i]))); out[i div 2]=matrix(transpose(syz(transpose(N1)))); out[i div 2]=transpose(matrix(std(transpose(out[i div 2])))); } else { li=L[i-1]; N1=transpose(syz(transpose(L[i]))); N=transpose(syz(transpose(N1))); lift1=matrixlift(N1,L[i-2]); out[i div 2]=transpose(concat(transpose(lift1),transpose(N))); out[i div 2]=transpose(matrix(std(transpose(out[i div 2])))); } } if (out[i div 2]!=matrix(0,1,ncols(out[i div 2]))) { out[i div 2]=ncols(out[i div 2])-nrows(out[i div 2]); } else { out[i div 2]=ncols(out[i div 2]); } } if (size(out)>le) { out=delete(out,1); } return(out); } ///////////////////////////////////////////////////////////////////// static proc mhom(poly f) { poly g; poly l; poly add; intvec e; list minint; list remf; int i; int j; int n=nvars(basering) div 4; if (f==0) { return(f); } while (f!=0) { l=lead(f); e=leadexp(l); remf[size(remf)+1]=list(); remf[size(remf)][1]=l; for (i=1; i<=n; i++) { remf[size(remf)][i+1]=-e[2*n+i]+e[3*n+i]; if (size(minint)j){return(i);} return(j); } //////////////////////////////////////////////////////////////////////////////////// version="$Id$"; category="Noncommutative"; info=" LIBRARY:  derham.lib      Computation of deRham cohomology AUTHORS:  Cornelia Rottner, rottner@mathematik.uni-kl.de OVERVIEW: A library for computing the de Rham cohomology of complements of complex affine varieties. REFERENCES: [OT] Oaku, T.; Takayama, N.: Algorithms of D-modules - restriction, tensor product, localzation, and local cohomology groups, J. Pure Appl. Algebra 156, 267-308 REFERENCES: [OT] Oaku, T.; Takayama, N.: Algorithms of D-modules - restriction, tensor product, localzation, and local cohomology groups}, J. Pure Appl. Algebra 156, 267-308 (2001) [R]  Rottner, C.: Computing de Rham Cohomology,diploma thesis (2012) [W1] Walther, U.: Algorithmic computation of local cohomology modules and the local cohomological dimension of algebraic varieties, J. Pure Appl. Algebra 139, [W1] Walther, U.: Algorithmic computation of local cohomology modules and the local cohomological dimension of algebraic varieties}, J. Pure Appl. Algebra 139, 303-321 (1999) [W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of Complex Affine Varieties, J. Symbolic Computation 29, 796-839 (2000) [W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of Complex Affine Varieties}, J. Symbolic Computation 29, 796-839 (2000) [W3] Walther, U.: Computing the cup product structure for complements of complex affine varieties, J. Pure Appl. Algebra 164, 247-273 (2001) LIB "poly.lib"; LIB "schreyer.lib"; LIB "dmodloc.lib"; //////////////////////////////////////////////////////////////////////////////////// proc deRhamCohomology(list L,list #) "USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices "USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices optional list consisting of one up to three strings @* The optional strings may be one of the strings@* -'Vdres': compute the Cartan-Eilenberg resolutions via V_d-homogenization -'noCE': compute quasi-isomorphic complexes without using Cartan-Eilenberg resolutionsq@* -'Vdres': compute quasi-isomorphic complexes using Cartan-Eilenberg resolutions; the CE resolutions are computed via V__d-homogenization and without using Schreyer's method @* -'Sres': compute the Cartan-Eilenberg resolutions in the homogenized Weyl algebra using Schreyer's method@* -'Sres': compute quasi-isomorphic complexes using Cartan-Eilenberg resolutions in the homogenized Weyl algebra via Schreyer's method@* one of the strings@* -'iterativeloc': compute localizations by factorizing the polynomials and -'iterativeloc': compute localizations by factorizing the polynomials and sucessive localization of the factors @* -'no iterativeloc': compute localizations by directly localizing the -'exactroots' computes the minimal and maximal integer root of the global b-function The default is 'Sres', 'iterativeloc' and 'onlybounds'. The default is 'noCE', 'iterativeloc' and 'onlybounds'. ASSUME: -The basering must be a polynomial ring over the field of rational numbers@* RETURN: list, where the ith entry is the (i-1)st de Rham cohomology group of the complement of the complex affine variety given by the polynomials in L RETURN: list, where the ith entry is the (i-1)st de Rham cohomology group of the complement of the complex affine variety given by the polynomials in L EXAMPLE:example deRhamCohomology; shows an example " int n=nvars(R); int le=size(L)+n; string Syzstring="Sres"; string Syzstring="noCE"; int onlybounds=1; int diffforms; for (i=1; i<=size(#); i++) { if (#[i]=="Vdres") { Syzstring="Vdres"; } if (#[i]=="noiterativeloc") { recursiveloc=0; } if (#[i]=="exactroots") { onlybounds=0; } } { if (#[i]=="Sres") { Syzstring="Sres"; } if (#[i]=="Vdres") { Syzstring="Vdres"; } if (#[i]=="noiterativeloc") { recursiveloc=0; } if (#[i]=="exactroots") { onlybounds=0; } if (#[i]=="diffforms") { diffforms=1; } } for (i=1; i<=size(L); i++) { if (L[i]==0) { L=delete(L,i); i=i-1; } } { if (L[i]==0) { L=delete(L,i); i=i-1; } } if (size(L)==0) { return (list(0)); } { return (list(0));//////////////////////////////////////////////////////////////////stimmt das jetzt?!?????????????????????????????????? } for (i=1; i<= size(L); i++) { if (leadcoef(L[i])-L[i]==0) { return(list(1)); } } { if (leadcoef(L[i])-L[i]==0) { return(list(1));    ///////////////////////////////////////////////////////////////stimmt das jetzt?!???????????????????????????????????? } } if (size(L)==0) { /*the complement of the variety given by the input is the whole space*/ return(list(1)); } { /*the complement of the variety given by the input is the whole space*/ return(list(1)); } for (i=1; i<=size(L); i++) { if (typeof(L[i])!="poly") { print("The input list must consist of polynomials"); retrun(); } } { if (typeof(L[i])!="poly") { print("The input list must consist of polynomials"); return(); } } if (size(L)==1 and Syzstring=="noCE") { Syzstring="Sres"; } /* 1st step: compute the Mayer-Vietoris Complex and its Fourier transform*/ def W=MVComplex(L,recursiveloc);//new ring that contains the MV complex setring W; list fortoVdstrict=MV; ideal IFourier=var(n+1); for (i=2;i<=n;i++) { IFourier=IFourier,var(n+i); } for (i=1; i<=n;i++) { IFourier=IFourier,-var(i); } map cFourier=W,IFourier; matrix sup; for (i=1; i<=size(MV); i++) { sup=fortoVdstrict[i]; /*takes the Fourier transform of the MV complex*/ fortoVdstrict[i]=cFourier(sup); } /* 2nd step: Compute a V_d-strict free complex that is quasi-isomorphic to the complex fortoVdstrict The 1st entry of the list rem will be the quasi-isomorphic complex, the 2nd entry contains the cohomology modules and is needed for the computation of the global b-function*/ list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring); if (diffforms==0) { ideal IFourier=var(n+1); for (i=2;i<=n;i++) { IFourier=IFourier,var(n+i); } for (i=1; i<=n;i++) { IFourier=IFourier,-var(i); } map cFourier=W,IFourier; matrix sup; for (i=1; i<=size(MV); i++) { sup=fortoVdstrict[i]; /*takes the Fourier transform of the MV complex*/ fortoVdstrict[i]=cFourier(sup); } } /* 2nd step: Compute a V_d-strict free complex that is quasi-isomorphic to the complex fortoVdstrict The 1st entry of the list rem will be the quasi-isomorphic complex, the 2nd entry contains the cohomology modules and is needed for the computation of the global b-function*/ if (Syzstring=="noCE") { list rem=quasiisomorphicVdComplex(fortoVdstrict,diffforms); list quasiiso=rem[3]; } else { list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring,diffforms); if (diffforms==1) { list quasiiso=list(matrix(1,1,1)); } } list newcomplex=rem[1]; //////////////////////////////////////////////////////////////////////////////////// /* 3rd step: Compute the  bounds for the minimal and maximal integer root of the global b-function of newcomplex(i.e. compute the lcm of the b-functions of its cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal integer root. If we compute only the bounds, we omit additional Groebner basis computations. However this leads to a higher-dimensional truncated complex. Note that the  cohomology modules are already contained in rem[2]. minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/ if (onlybounds==1) { list minmaxk=globalBFun(rem[2],Syzstring); } //////////////////////////////////////////////////////////////////////////////////// /* 3rd step: Compute the  bounds for the minimal and maximal integer root of the global b-function of newcomplex(i.e. compute the lcm of the b-functions of its cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal integer root. If we compute only the bounds, we omit additional Groebner basis computations. However this leads to a higher-dimensional truncated complex. Note that the  cohomology modules are already contained in rem[2]. minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/ if (diffforms==1) { list minmaxk=exactGlobalBFunIntegration(rem[2]); } else { list minmaxk=exactGlobalBFun(rem[2],Syzstring); } { if (onlybounds==1) { list minmaxk=globalBFun(rem[2],Syzstring); } else { list minmaxk=exactGlobalBFun(rem[2],Syzstring); } } if (size(minmaxk)==0) { return (0); } /*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i-2], m^i=newcomplex[3*i-1], d^i=newcomplex[3*i]), using Thm 5.7 in [W1]: The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set (0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and if it is placed in component k it holds that minmaxk[1]-m^i[k]<=deg(P_(i_j))<=minmaxk[2]-m^i[k]*/ { return (0); } ///////////////////////////////////////////////////////////////////////////Bis hierhin angepasst /*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i-2], m^i=newcomplex[3*i-1], d^i=newcomplex[3*i]), using Thm 5.7 in [W1]: The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set (0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and if it is placed in component k it holds that minmaxk[1]-m^i[k]<=deg(P_(i_j))<=minmaxk[2]-m^i[k]*/ int k,l; list truncatedcomplex,shorten,upto; for (i=1; i<=size(newcomplex) div 3; i++) { shorten[3*i-1]=list(); for (j=1; j<=size(newcomplex[3*i-1]); j++) { /*shorten[3*i-1][j][k]=minmaxk[k]-m^i[j]+1 (for k=1,2) if this value is positive otherwise we will set it to be list(); we added +1, because we will use a list, where we put in position l polys of degree l+1*/ shorten[3*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1); shorten[3*i-1][j][2]=minmaxk[2]-newcomplex[3*i-1][j]+1; upto[size(upto)+1]=shorten[3*i-1][j][2]; if (shorten[3*i-1][j][2]<=0) { shorten[3*i-1][j]=list(); } else { if (shorten[3*i-1][j][1]<=0) { shorten[3*i-1][j][1]=1; } } } } { shorten[3*i-1]=list(); for (j=1; j<=size(newcomplex[3*i-1]); j++) { /*shorten[3*i-1][j][k]=minmaxk[k]-m^i[j]+1 (for k=1,2) if this value is positive otherwise we will set it to be list(); .-            we added +1, because we will use a list, where we put in position l polys of degree l+1*/ shorten[3*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1); if (diffforms==1) { shorten[3*i-1][j][1]=1; } shorten[3*i-1][j][2]=minmaxk[2]-newcomplex[3*i-1][j]+1; upto[size(upto)+1]=shorten[3*i-1][j][2]; if (shorten[3*i-1][j][2]<=0) { shorten[3*i-1][j]=list(); } else { if (shorten[3*i-1][j][1]<=0) { shorten[3*i-1][j][1]=1; } } } } int iupto=Max(upto);//maximal degree +1 of the polynomials we have to consider if (iupto<=0) { return(list(0)); } { return(list(0)); } list allpolys; /*allpolys[i] will consist list of all monomials in D(1),...,D(n) of degree i-1*/ allpolys[1]=list(1); list minvar; list keepv; minvar[1]=list(1); for (i=1; i<=iupto-1; i++) { allpolys[i+1]=list(); minvar[i+1]=list(); for (k=1; k<=size(allpolys[i]); k++) { for (j=minvar[i][k]; j<=nvars(W) div 2; j++) { allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j); minvar[i+1][size(minvar[i+1])+1]=j; } } } { allpolys[i+1]=list(); minvar[i+1]=list(); for (k=1; k<=size(allpolys[i]); k++) { for (j=minvar[i][k]; j<=nvars(W) div 2; j++) { if (diffforms==0) { allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j); } else { allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*x(j); } minvar[i+1][size(minvar[i+1])+1]=j; } } } list keepformatrix,sizetruncom,fortrun,fst; int count,stc; intvec v,forin; matrix subm; list keepcount; list passendespoly; /*now we compute the truncation*/ for (i=1; i<=size(newcomplex) div 3; i++) { /*truncatedcomplex[2*i-1] will contain all the generators for the truncation of D_n/(x(1),..,x(n))\otimes C[i]*/ truncatedcomplex[2*i-1]=list(); sizetruncom[2*i-1]=list(); sizetruncom[2*i]=list(); /*truncatedcomplex[2*i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i]) ->trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/ truncatedcomplex[2*i]=newcomplex[3*i]; v=0;count=0; sizetruncom[2*i][1]=0; for (j=1; j<=newcomplex[3*i-2]; j++) { if (size(shorten[3*i-1][j])!=0) { fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]); truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1]; count=count+fortrun[2]; fst=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1); sizetruncom[2*i-1][size(sizetruncom[2*i-1])+1]=fst; sizetruncom[2*i][size(sizetruncom[2*i])+1]=count; if (v!=0) { v[size(v)+1]=j; } else { v[1]=j; } } } if (v!=0) { subm=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i])); truncatedcomplex[2*i]=subm; if (i!=1) { i1=1..nrows(truncatedcomplex[2*(i-1)]); subm=submat(truncatedcomplex[2*(i-1)],i1,v); truncatedcomplex[2*(i-1)]=subm; } } else { truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i])); if (i!=1) { nr=nrows(truncatedcomplex[2*(i-1)]); truncatedcomplex[2*(i-1)]=matrix(0,nr,1); } } } { /*truncatedcomplex[2*i-1] will contain all the generators for the truncation of D_n/(x(1),..,x(n))\otimes C[i]*/ truncatedcomplex[2*i-1]=list(); sizetruncom[2*i-1]=list(); sizetruncom[2*i]=list(); passendespoly[i]=list(); /*truncatedcomplex[2*i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i]) ->trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/ truncatedcomplex[2*i]=newcomplex[3*i]; v=0;count=0; sizetruncom[2*i][1]=0; for (j=1; j<=newcomplex[3*i-2]; j++) { if (size(shorten[3*i-1][j])!=0) { fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]); truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1]; for (k=1; k<=size(fortrun[1]); k++) { for (l=1; l<=size(fortrun[1][k]); l++) { passendespoly[i][size(passendespoly[i])+1]=list(fortrun[1][k][l][1],j); } } count=count+fortrun[2]; fst=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1); sizetruncom[2*i-1][size(sizetruncom[2*i-1])+1]=fst; sizetruncom[2*i][size(sizetruncom[2*i])+1]=count; if (v!=0) { v[size(v)+1]=j; } else { v[1]=j; } } } if (v!=0) { keepv[i]=v; subm=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i])); truncatedcomplex[2*i]=subm; if (i!=1) { i1=1..nrows(truncatedcomplex[2*(i-1)]); subm=submat(truncatedcomplex[2*(i-1)],i1,v); truncatedcomplex[2*(i-1)]=subm; } } else { keepv[i]=list(); truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i])); if (i!=1) { nr=nrows(truncatedcomplex[2*(i-1)]); truncatedcomplex[2*(i-1)]=matrix(0,nr,1); } } } list keeptruncatedcomplex=truncatedcomplex; matrix M; int st,pi,pj; poly ptc; int b,d,ideg,kplus,lplus; int z; poly form,lform,nform; /*computation of the maps*/ if (diffforms==1) { def ConvWeyl=makeConverseWeyl(nvars(basering) div 2); setring ConvWeyl; poly form,lform,nform; poly ptc; list truncatedcomplex; matrix M; ideal I=x(1); for (i=2; i<=nvars(basering) div 2; i++) { I=I,var(nvars(basering) div 2 + i); } for (i=1; i<=nvars(basering) div 2; i++) { I=I,var(i); } map transtc=W,I; truncatedcomplex=transtc(truncatedcomplex); } for (i=1; i=0) { nr=ideg+1; st=size(truncatedcomplex[2*(i+1)-1][l][nr]); for (d=1; d<=st;d++) { nr=max(1,sizetruncom[2*i][size(sizetruncom[2*i])]); nc=max(1,sizetruncom[2*i+2][size(sizetruncom[2*i+2])]); M=matrix(0,nr,nc); for (k=1; k<=size(truncatedcomplex[2*i-1]);k++) { for (l=1; l<=size(truncatedcomplex[2*(i+1)-1]); l++) { if (size(sizetruncom[2*i])!=1) { for (j=1; j<=size(truncatedcomplex[2*i-1][k]); j++) { nc=2*(i+1)-1; ptc=truncatedcomplex[nc][l][ideg+1][d][1]; if (leadmonom(lform)==ptc) { nr=2*i-1; pi=truncatedcomplex[nr][k][j][b][2]; pi=pi+sizetruncom[2*i][k]; nc=2*(i+1)-1; nr=ideg+1; pj=truncatedcomplex[nc][l][nr][d][2]; pj=pj+sizetruncom[2*(i+1)][l]; M[pi,pj]=leadcoef(lform); break; } for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++) { form=truncatedcomplex[2*i-1][k][j][b][1]; form=form*truncatedcomplex[2*i][k,l]; for (z=1; z<=nvars(basering) div 2; z++)//neu {// form=subst(form,var(z),0);// }// while (form!=0) { lform=lead(form); v=leadexp(lform); v=v[1..n]; // if (v==(0:n)) //{ ideg=deg(lform)-sizetruncom[2*(i+1)-1][l][1]; if (ideg>=0) { nr=ideg+1; st=size(truncatedcomplex[2*(i+1)-1][l][nr]); for (d=1; d<=st;d++) { nc=2*(i+1)-1; ptc=truncatedcomplex[nc][l][ideg+1][d][1]; if (leadmonom(lform)==ptc) { nr=2*i-1; pi=truncatedcomplex[nr][k][j][b][2]; pi=pi+sizetruncom[2*i][k]; nc=2*(i+1)-1; nr=ideg+1; pj=truncatedcomplex[nc][l][nr][d][2]; pj=pj+sizetruncom[2*(i+1)][l]; M[pi,pj]=leadcoef(lform); break; } } } //        } form=form-lform; } } } } } form=form-lform; } } } } } } truncatedcomplex[2*i]=M; truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; } } } } truncatedcomplex[2*i]=M; truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; } truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; if (truncatedcomplex[2*i-1]!=0) { truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1); { truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1); } if (diffforms==1) { setring W; truncatedcomplex=imap(ConvWeyl,truncatedcomplex); } setring R; list truncatedcomplex=imap(W,truncatedcomplex); /*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex, list truncatedcomplex=imap(W,truncatedcomplex); /*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex, i.e. D^i=C^truncatedcomplex[2*i-1] and d^i=truncatedcomplex[2*i]*/ list derhamhom=findCohomology(truncatedcomplex,le); option(set,saveoptions); return (derhamhom); if (diffforms==0) { list derhamhom=findCohomology(truncatedcomplex,le); option(set,saveoptions); return (derhamhom); } list outall=findCohomologyDiffForms(truncatedcomplex,le); setring W; list dimanddiff=imap(R,outall); list alldiffforms=dimanddiff[2]; while(size(alldiffforms)size(generators)) { generators=insert(generators,list()); } while (size(dimanddiff[1])>size(quasiiso)) { quasiiso=insert(quasiiso,list()); } int keepsign; list derhamdiff; list doublecom=makeDoubleComplex(newcomplexmod,omegacomplex,quasiiso,generators); matrix diffform; int stopping; int p; matrix convert; list interim; list correspondingposition; list allforms=list(); for (i=1; i<=size(newdiffforms); i++) { derhamdiff[i]=list(); allforms[i]=list(); for (j=1; j<=size(newdiffforms[i]); j++) { allforms[i][j]=list(); keepsign=1; derhamdiff[i][j]=0; diffform=newdiffforms[i][j];//Zeilenform correspondingposition=doublecom[i][1];//needed fpr transformation process interim=transferDiffforms(diffform,correspondingposition); if (size(interim)!=0) { allforms[i][j][size(allforms[i][j])+1]=interim; } stopping=0; p=1; for (k=i; k<=size(newdiffforms); k++) { keepsign=(-1)*keepsign; if (stopping==0) { if (size(doublecom[k][p][2])==0) { stopping=1; } else { if (size(doublecom[k+1][p][3])!=0) { diffform=diffform*doublecom[k][p][2];//Spaltenform if (diffform!=matrix(0,nrows(diffform),ncols(diffform))) { diffform=findPreimage(doublecom[k+1][p][3],transpose(diffform));//Zeilenform correspondingposition=doublecom[k+1][p+1];//needed for transformation process interim=transferDiffforms(keepsign*diffform,correspondingposition); if (size(interim)!=0) { allforms[i][j][size(allforms[i][j])+1]=interim; } p=p+1; } else { stopping=1; } } else { stopping=1; } } } } } } setring R; list allforms=fetch(W,allforms); option(set,saveoptions); return (allforms); } -L is a list of non-constant polynomials RETURN: ring W: the nth Weyl algebra @* W contains a list MV, which represents the Mayer-Vietrois complex (C^i,d^i) of the W contains a list MV, which represents the Mayer-Vietrois complex (C^i,d^i) of the polynomials contained in L as follows:@* the C^i are given  by D_n^ncols(C[2*i-1])/im(C[2*i-1]) and the differentials { /* We follow algorithm 3.2.5 in [R],if #!=0 we use also  Remark 3.2.6 in [R] for an additional iterative localization*/ an additional iterative localization*/ def R=basering; int i; int iterative=1; if (size(#)!=0) { iterative=#[1]; } { iterative=#[1]; } for (i=1; i<=size(L); i++) { if (L[i]==0) { print("localization with respect to 0 not possible"); return(); } if (leadcoef(L[i])-L[i]==0) { print("polynomials must be non-constant"); return(); } } { if (L[i]==0) { print("localization with respect to 0 not possible"); return(); } if (leadcoef(L[i])-L[i]==0) { print("polynomials must be non-constant"); return(); } } if (iterative==1) { /*compute the localizations by factorizing the polynomials and iterative localization of the factors */ for (i=1; i<=size(L); i++) { L[i]=factorize(L[i],1); } } { /*compute the localizations by factorizing the polynomials and iterative localization of the factors */ for (i=1; i<=size(L); i++) { L[i]=factorize(L[i],1); } } int r=size(L); int n=nvars(basering); /*construct the ring Ws*/ def W=makeWeyl(n); setring W; setring W; list man=ringlist(W); if (n==1) { man[2][1]="x(1)"; man[2][2]="D(1)"; def Wi=ring(man); setring Wi; kill W; def W=Wi; setring W; list man=ringlist(W); } { man[2][1]="x(1)"; man[2][2]="D(1)"; def Wi=ring(man); setring Wi; kill W; def W=Wi; setring W; list man=ringlist(W); } man[2][size(man[2])+1]="s";; man[3][3]=man[3][2]; matrix M[1][1]; man[6]=transpose(concat(transpose(concat(man[6],M)),M)); def Ws=ring(man); def Ws=ring(man); setring Ws; int j,k,l,c; list L=fetch(R,L); list Cech; int j,k,l,c; list L=fetch(R,L); list Cech; ideal J=var(1+n); for (i=2; i<=n; i++) { J=J,var(i+n); } { J=J,var(i+n); } Cech[1]=list(J); list Theta, remminroots; int rmr; if (iterative==0) {/*computation of the modules of the MV complex*/ for (i=1; i<=r; i++) { findminintroot=list(); Cech[i+1]=list(); Theta[i+1]=list(); k1=1; for (j=1; j<=i; j++) { k1[size(k1)+1]=size(Theta[j+1]); for (k=1; k<=k1[j]; k++) { Theta[j+1][size(Theta[j+1])+1]=list(Theta[j][k][1]+list(i)); Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i]; /*We compute the s-parametric annihilator J(s)  and the b-function of the polynomial L[i] and Cech[i][k] to localize the module D_n/(D(1),...,D(n))[L[i]^(-1)]\otimes D_n^c/im(Cech[i][k]), where c=ncols(Cech[i][k]) and the im(Cech[i][k]) is generated by the rows of the matrix. If we plug the minimal integer root r(or a smaller integer value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic to the above localization*/ rem=SannfsIBM(L[i],Cech[j][k]); Cech[j+1][size(Cech[j+1])+1]=rem[1]; findminintroot[size(findminintroot)+1]=rem[2]; } } /* we compute the minimal root of all b-functions of L[i] computed above, because we want to plug in the same root r in all s-parametric annihilators we computed for L[i]  ->this will ensure  we can compute the maps of the MV complex*/ minroot=minIntRoot0(findminintroot); for (j=1; j<=i; j++) { for (k=1; k<=k1[j]; k++) { sk=size(Cech[j+1])+1-k; Cech[j+1][size(Cech[j+1])+1-k]=subst(Cech[j+1][sk],s,minroot); } } remminroots[i]=minroot; } Cech=delete(Cech,1); Theta=delete(Theta,1); list zw; poly reme; /*computation of the maps of the MV complex*/ for (i=1; ithis will ensure  we can compute the maps of the MV complex*/ minroot=minIntRoot(findminintroot); for (j=1; j<=i; j++) { for (k=1; k<=k1[j]; k++) { sk=size(Cech[j+1])+1-k; Cech[j+1][size(Cech[j+1])+1-k]=subst(Cech[j+1][sk],s,minroot); } } remminroots[i]=minroot; } Cech=delete(Cech,1); Theta=delete(Theta,1); list zw; poly reme; /*computation of the maps of the MV complex*/ for (i=1; ithis will ensure  we can compute the maps of the MV complex*/ minroot=minIntRoot0(findminintroot); for (j=1; j<=i; j++) { for (k=1; k<=k1[j]; k++) { st=size(Cech[j+1])+1-k; Cech[j+1][st]=subst(Cech[j+1][st],s,minroot); } } if (c==1) { remminroots[i]=list(); } remminroots[i][c]=minroot; } } Cech=delete(Cech,1); Theta=delete(Theta,1); list zw; poly reme; /*maps of the MV Complex*/ for (i=1; ithis will ensure  we can compute the maps of the MV complex*/ minroot=minIntRoot(findminintroot); for (j=1; j<=i; j++) { for (k=1; k<=k1[j]; k++) { st=size(Cech[j+1])+1-k; Cech[j+1][st]=subst(Cech[j+1][st],s,minroot); } } if (c==1) { remminroots[i]=list(); } remminroots[i][c]=minroot; } } Cech=delete(Cech,1); Theta=delete(Theta,1); list zw; poly reme; /*maps of the MV Complex*/ for (i=1; iD_n^(r_1)/im(M_1)-> D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i) is generated by the rows of M_i. In particular it hold:@* - The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @* -the image of M_1*f_i is contained in the image of M_(i+1) @* d is an integer between 1 and n. If no value for d is given, it is assumed to be n @* df is an optional int, if df equals 1 a \tilde(V_d)-strict complex will be computed (instead of a V_d-strict one) (for a definition see [W3]) RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @* L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,i_s,g_s,m_s) such that:@* -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs of size i_j@* -D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0  is a V_d-strict complex with differentials m_i that is quasi-isomorphic to the complex given by L@* L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and the n_i are shift vectors such that:@* -coker(H_i) is the ith cohomology group of the complex given by L_1@* -the n_i are the shift vectors of the coker(H_i) THEORY: We follow Proposition 3.2 and Corollary 3.3 in [W3] " { int tilde; if (size(#)!=0) { tilde=#[1]; } def B=basering; int n=nvars(B) div 2 + 1;//+1 müsste stimmen! bitte kontrollieren! int d=nvars(B) div 2; int r=size(L) div 2; int lonc=n+r; int Kiold=0; matrix kerold; // matrix kernew=out[r][2][2]; matrix kernew=diag(1,ncols(L[size(L)-1])); module mL; int i; int k; matrix testm; int Kinew=nrows(kernew); int Jiold=0; int Jinew=0; matrix Niold; matrix Ninew; list newcomplex; int Aiold=Kinew; matrix savediv; newcomplex[3*lonc-2]=Kinew; newcomplex[3*lonc-1]=intvec(0:Kinew); newcomplex[3*lonc]=matrix(0,Kinew,1); list quasiiso; quasiiso[lonc]=diag(1,Kinew); matrix invimage; matrix keralpha; intvec v; int j; matrix sc; matrix fnc; int indk; int indj; int Aiold; list saveres; matrix Liplus; for (i=r-1; i>=0; i--) { indk=0; indj=0; Kiold=Kinew; kerold=kernew; if (i!=0) { // kernew=divdr(L[2*i],L[2*i+1],1); kernew=divdr(L[2*i],L[2*i+1]); mL=slimgb(transpose(L[2*i-1])); for (k=1; k<=nrows(kernew); k++) { testm=reduce(transpose(submat(kernew,k,intvec(1..ncols(kernew)))),mL); if (testm==matrix(0,nrows(testm),ncols(testm))) { kernew=transpose(deletecol(transpose(kernew),k)); k=k-1; } } Kinew=nrows(kernew); if (kernew==matrix(0,nrows(kernew),ncols(kernew))) { Kinew=0; indk=1; } } else { Kinew=0; indk=1; } Jiold=Jinew; Niold=Ninew; keralpha=transpose(syz(transpose(newcomplex[3*(i+n)+3]))); if (i!=0) { invimage=divdr(quasiiso[n+i+1],transpose(concat(transpose(L[2*i]),transpose(L[2*i+1])))); Ninew=vdStrictIntersect(keralpha,invimage,newcomplex[3*(n+i+1)-1],tilde);////////////// } else { invimage=divdr(quasiiso[n+i+1],L[2*i+1]); saveres=vdStrictIntersectPlus(keralpha,invimage,newcomplex[3*(n+i+1)-1],tilde);//////////////////////// ///////////////////BIS HIERHIN VERALLGEMEINERT//////////////////////////////////////////////////////////////////// Ninew=saveres[1]; } Jinew=nrows(Ninew); if (Ninew==matrix(0,nrows(Ninew),ncols(Ninew))) { Jinew=0; indk=1; } newcomplex[3*(n+i)-2]=Kinew+Jinew; v=0; if (indk==0) { v=(0:Kinew); if (indj==0) { fnc=transpose(concat(transpose(matrix(0,Kinew,Kiold+Jiold)),transpose(Ninew))); } else { fnc=matrix(0,Kinew,Kiold+Jiold); } } else { if (indj==0) { fnc=Ninew; } else { fnc=matrix(0,1,Kiold+Jiold); newcomplex[3*(n+i)-2]=1; } } Aiold=Jinew+Kinew; if (Aiold==0) { Aiold=1; } newcomplex[3*(n+i)]=fnc; for (j=1; j<=Jinew; j++) { if (tilde==0) { v[Kinew+j]=VdDeg(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]); } else { v[Kinew+j]=VdDegTilde(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]); } } newcomplex[3*(n+i)-1]=v; if (i==0) { quasiiso[n+i]=matrix(0,Jinew,1); } else { if (indj==0) { sc=submat(fnc,intvec(Kinew+1..nrows(fnc)),intvec(1..ncols(fnc)))*quasiiso[n+i+1]; Liplus=transpose(concat(transpose(L[2*i]),transpose(L[2*i+1]))); sc=matrixLift(Liplus,sc);//stimmt das jetzt sc=submat(sc,intvec(1..nrows(sc)),intvec(1..nrows(L[2*i]))); if (indk==0) { //pi=kernew quasiiso[n+i]=transpose(concat(transpose(kernew),transpose(sc))); } else { quasiiso[n+i]=sc; } } else { if (indk==0) { quasiiso[n+i]=kernew; } else { quasiiso[n+i]=matrix(0,1,ncols(kernew)); } } } } for (i=1; i<=n-1; i++) { quasiiso[n-i]=list(); if (size(saveres[2][i])!=0) { newcomplex[3*(n-i)]=saveres[2][i]; newcomplex[3*(n-i)-2]=nrows(saveres[2][i]); v=0; for (j=1; j<=newcomplex[3*(n-i)-2]; j++) { if (tilde==0) { v[j]=VdDeg(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(n-i)+2]); } else { v[j]=VdDegTilde(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(n-i)+2]); } } newcomplex[3*(n-i)-1]=v; } else { newcomplex[3*(n-i)]=matrix(0,1,1); if (newcomplex[3*(n-i)+1]!=0) { newcomplex[3*(n-i)]=matrix(0,1,newcomplex[3*(n-i)+1]); } newcomplex[3*(n-i)-2]=int(0); newcomplex[3*(n-i)-1]=intvec(0); } } list result; result[1]=newcomplex; result[2]=list(); list forsep; for (i=1; i<=size(L) div 2+1; i++) { forsep[2*i]=newcomplex[3*(n+i-1)]; forsep[2*i-1]=matrix(0,1,nrows(forsep[2*i])); } forsep=shortExactPieces(forsep); list listofHis; matrix forVd; for (i=1; i<=size(L) div 2; i++) { v=0; listofHis[i]=list(forsep[i+1][1][5]); forVd=forsep[i+1][2][2]; for (j=1; j<=nrows(forVd); j++) { if (tilde==0) { v[j]=VdDeg(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)-1]); } else { v[j]=VdDegTilde(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)-1]); } } listofHis[i][2]=v; } result[2]=listofHis; result[3]=quasiiso; return(result); } //////////////////////////////////////////////////////////////////////////////////// static proc vdStrictIntersect(matrix M, matrix N, intvec v, int tilde) { def B=basering; option(returnSB);//                    alternative:erst intersect und dann SB-Berechung mit slimgb if (tilde==0) { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v); } else { def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v); } setring HomWeyl; matrix M=fetch(B,M); matrix N=fetch(B,N); M=nHomogenize(M); N=nHomogenize(N); matrix vdintersection=transpose(intersect(transpose(M),transpose(N))); vdintersection=subst(vdintersection,h,1); setring B; matrix vdintersection=fetch(HomWeyl,vdintersection); option(noreturnSB); return(vdintersection); } //////////////////////////////////////////////////////////////////////////////////// static proc vdStrictIntersectPlus(matrix M, matrix N, intvec v, int tilde) { def B=basering; int n=nvars(B) div 2; matrix vdint=transpose(intersect(transpose(M),transpose(N))); if (tilde==0) { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v); } else { def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v); } setring HomWeyl; matrix vdint=fetch(B,vdint); matrix N=fetch(B,N); vdint=nHomogenize(vdint); intvec i1; intvec i2; int i; int nr; int nc; def ringofSyz=Sres(transpose(vdint),n);//////////////////////////////////////////////////////////////// setring ringofSyz; matrix vdint=transpose(matrix(RES[2])); vdint=subst(vdint,h,1); int logens=ncols(vdint)+1; int omitemptylist; matrix zerom; list rofA; for (i=3; i<=n+3; i++)////////////////////////////////////////////////////////////////////////////n und si müssen noch definiert werden { if (size(RES)>=i) { zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); if (RES[i]!=zerom) { rofA[i-2]=(matrix(RES[i])); if (i==3) { if (nrows(rofA[i-2])-logens+1!=nrows(vdint)) { //build the resolution nr=nrows(vdint)+logens-1; nc=ncols(rofA[i-2]); rofA[i-2]=matrix(rofA[i-2],nr,nc); } } if (i!=3) { if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3])) { nr=nrows(rofA[i-3])+logens-1; nc=ncols(rofA[i-2]); rofA[i-2]=matrix(rofA[i-2],nr,nc); } } i1=intvec(logens..nrows(rofA[i-2])); i2=intvec(1..ncols(rofA[i-2])); rofA[i-2]=transpose(submat(rofA[i-2],i1,i2)); logens=logens+ncols(rofA[i-2]); rofA[i-2]=subst(rofA[i-2],h,1); } else { rofA[i-2]=list(); } } else { rofA[i-2]=list(); } } if(size(rofA[1])==0) { omitemptylist=1; } setring B; vdint=fetch(ringofSyz,vdint); if (omitemptylist!=1) { list rofA=fetch(ringofSyz,rofA); } kill HomWeyl; kill ringofSyz; return(list(vdint,rofA)); } //////////////////////////////////////////////////////////////////////////////////// static proc toVdStrictFreeComplex(list L,string Syzstring,list #) "USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form (M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a "USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form (M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a string, d an optional integer ASSUME: Basering is the Weyl algebra D_n @* (M_1,f_1,...,M_s,f_s) represents a complex 0->D_n^(r_1)/im(M_1)-> D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i) D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i) is generated by the rows of M_i. In particular it hold:@* - The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @* -the image of M_1*f_i is contained in the image of M_(i+1) @* d is an integer between 1 and n. If no value for d is given, it is assumed to be n @* d is an optional integer which indices in the case size(L)=2, whether a V_d-strict or \tilde(V_d)-strict will be computed@* Syzstring is either: @* -'Sres' (computes the resolutions and Groebner bases in the homogenized Weyl algebra using Schreyer's method)@* -'Sres' (computes the resolutions and Groebner bases in the homogenized Weyl algebra using Schreyer's method)@* or @* -'Vdres' (computes the resolutions via V_d-homogenization and without -'Vdres' (computes the resolutions via V_d-homogenization and without Schreyer's method)@* RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @* L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,i_s,g_s,m_s) such that:@* -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs -the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs of size i_j@* -D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0  is a V_d-strict complex with differentials m_i that is quasi-isomorphic to the complex given by L@* L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and the n_i are shift vectors such that:@* -coker(H_i) is the ith cohomology group of the complex given by L_1@* -coker(H_i) is the ith cohomology group of the complex given by L_1@* -the n_i are the shift vectors of the coker(H_i) THEORY: We follow Algorithm 3.8 in [W2] " { def B=basering; def B=basering; int n=nvars(B) div 2+2; int d=nvars(B) div 2; matrix mem; intvec i1,i2; int tilde; if (size(#)!=0) { for (i=1; i<=size(#); i++) { if (typeof(#[i])=="int") { if (#[1]>=1 and #[1]<=n) { d=#[i]; } } } } { for (i=1; i<=size(#); i++) { if (typeof(#[i])=="int") { tilde=#[i]; } } } /* If size(L)=2, our complex consists for only one non-trivial module. Therefore, we just have to compute a V_d-strict resolution of this module.*/ if (size(L)==2) { v=(0:ncols(L[1])); out[3*n-1]=v; out[3*n-2]=ncols(L[1]); out[3*n]=L[2]; if (Syzstring=="Vdres") { /*if Syzstring="Vdres", we compute a V_d-strict Groebner basis of L[1] using F-homogenization (Prop. 3.9 in [OT]); then we compute the syzygies and make them V_d-strict using Prop  3.9[OT] and so on*/ out[3*n-3]=VdStrictGB(L[1],d,v); for (i=n-1; i>=1; i--) { out[3*i-2]=nrows(out[3*i]); v=0; for (j=1; j<=out[3*i-2]; j++) { mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector } out[3*i-1]=v; if (i!=1) { /*next step in the resolution*/ out[3*i-3]=transpose(syz(transpose(out[3*i]))); if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3]))) { /*makes the resolution V_d-strict*/ out[3*i-3]=VdStrictGB(out[3*i-3],d,out[3*i-1]); } else { /*resolution is already computed*/ out[3*i-3]=matrix(0,1,ncols(out[3*i-3])); out[3*i-4]=intvec(0); out[3*i-5]=int(0); for (j=i-2; j>=1; j--) { out[3*j]=matrix(0,1,1); out[3*j-1]=intvec(0); out[3*j-2]=int(0); } break; } } } } else { /*in the case Syzstring!="Vdres" we compute the resolution in the homogenized Weyl algebra using Thm 9.10 in[OT]*/ def HomWeyl=makeHomogenizedWeyl(d); setring HomWeyl; list L=fetch(B,L); L[1]=nHomogenize(L[1]); list out=fetch(B,out); out[3*n-3]=L[1]; /*computes a ring with a list RES; RES is a V_d-strict resolution of L[1]*/ def ringofSyz=Sres(transpose(L[1]),d); setring ringofSyz; int logens=2; matrix mem; list out=fetch(HomWeyl,out); out[3*n-3]=transpose(matrix(RES[2])); out[3*n-3]=subst(out[3*n-3],h,1); for (i=n-1; i>=1; i--) { out[3*i-2]=nrows(out[3*i]); v=0; for (j=1; j<=out[3*i-2]; j++) { mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); v[j]=VdDeg(mem,d, out[3*i+2]); } out[3*i-1]=v;//shift vector such that the resolution RES is V_d-strict if (i!=1) { indi=0; if (size(RES)>=n-i+2) { nr=nrows(matrix(RES[n-i+2])); mem=matrix(0,nr,ncols(matrix(RES[n-i+2]))); if (matrix(RES[n-i+2])!=mem) { indi=1; out[3*i-3]=(matrix(RES[n-i+2])); if (nrows(out[3*i-3])-logens+1!=nrows(out[3*i])) { mem=out[3*i-3]; out[3*i-3]=matrix(mem,nrows(mem)+logens-1,ncols(mem)); } mem=out[3*i-3]; i1=intvec(logens..nrows(mem)); mem=submat(mem,i1,intvec(1..ncols(mem))); out[3*i-3]=transpose(mem); out[3*i-3]=subst(out[3*i-3],h,1); logens=logens+ncols(out[3*i-3]); } } if(indi==0) { out[3*i-3]=matrix(0,1,nrows(out[3*i])); out[3*i-4]=intvec(0); out[3*i-5]=int(0); for (j=i-2; j>=1; j--) { out[3*j]=matrix(0,1,1); out[3*j-1]=intvec(0); out[3*j-2]=int(0); } break; } } } setring B; out=fetch(ringofSyz,out);//contains the V_d-strict resolution kill ringofSyz; } outall[1]=out; outall[2]=list(list(out[3*n-3],out[3*n-1])); return(outall); } /*case size(L)>2: We compute a quasi-isomorphic free complex following Alg 3.8 in [W2]*/ /* We denote the complex given by L as (C^i,d^i). We start by computing in the proc shortExaxtPieces representations for the short exact sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1), where the B^i, Z^i and H^i are coboundaries, cocycles and cohomology groups, respectively.*/ out=shortExactPieces(L); Therefore, we just have to compute a V_d-strict resolution of this module.*/ if (size(L)==2) { v=(0:ncols(L[1])); out[3*n-1]=v; out[3*n-2]=ncols(L[1]); out[3*n]=L[2]; if (Syzstring=="Vdres") { /*if Syzstring="Vdres", we compute a V_d-strict Groebner basis of L[1] using F-homogenization (Prop. 3.9 in [OT]); then we compute the syzygies and make them V_d-strict using Prop  3.9[OT] and so on*/ out[3*n-3]=VdStrictGB(L[1],d,v); for (i=n-1; i>=1; i--) { out[3*i-2]=nrows(out[3*i]); v=0; for (j=1; j<=out[3*i-2]; j++) { mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector } out[3*i-1]=v; if (i!=1) { /*next step in the resolution*/ out[3*i-3]=transpose(syz(transpose(out[3*i]))); if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3]))) { /*makes the resolution V_d-strict*/ out[3*i-3]=VdStrictGB(out[3*i-3],d,out[3*i-1]); } else { /*resolution is already computed*/ out[3*i-3]=matrix(0,1,ncols(out[3*i-3])); out[3*i-4]=intvec(0); out[3*i-5]=int(0); for (j=i-2; j>=1; j--) { out[3*j]=matrix(0,1,1); out[3*j-1]=intvec(0); out[3*j-2]=int(0); } break; } } } } else { /*in the case Syzstring!="Vdres" we compute the resolution in the homogenized Weyl algebra using Thm 9.10 in[OT]*/ if (tilde==0) { def HomWeyl=makeHomogenizedWeyl(d); } else { def HomWeyl=makeHomogenizedWeylTilde(d); } setring HomWeyl; list L=fetch(B,L); L[1]=nHomogenize(L[1]); list out=fetch(B,out); out[3*n-3]=L[1]; /*computes a ring with a list RES; RES is a V_d-strict resolution of L[1]*/ def ringofSyz=Sres(transpose(L[1]),d); setring ringofSyz; int logens=2; matrix mem; list out=fetch(HomWeyl,out); out[3*n-3]=transpose(matrix(RES[2])); out[3*n-3]=subst(out[3*n-3],h,1); for (i=n-1; i>=1; i--) { out[3*i-2]=nrows(out[3*i]); v=0; for (j=1; j<=out[3*i-2]; j++) { mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); if (tilde==0) { v[j]=VdDeg(mem,d, out[3*i+2]); } else { v[j]=VdDegTilde(mem,d, out[3*i+2]); } } out[3*i-1]=v;//shift vector such that the resolution RES is V_d-strict if (i!=1) { indi=0; if (size(RES)>=n-i+2) { nr=nrows(matrix(RES[n-i+2])); mem=matrix(0,nr,ncols(matrix(RES[n-i+2]))); if (matrix(RES[n-i+2])!=mem) { indi=1; out[3*i-3]=(matrix(RES[n-i+2])); if (nrows(out[3*i-3])-logens+1!=nrows(out[3*i])) { mem=out[3*i-3]; out[3*i-3]=matrix(mem,nrows(mem)+logens-1,ncols(mem)); } mem=out[3*i-3]; i1=intvec(logens..nrows(mem)); mem=submat(mem,i1,intvec(1..ncols(mem))); out[3*i-3]=transpose(mem); out[3*i-3]=subst(out[3*i-3],h,1); logens=logens+ncols(out[3*i-3]); } } if(indi==0) { out[3*i-3]=matrix(0,1,nrows(out[3*i])); out[3*i-4]=intvec(0); out[3*i-5]=int(0); for (j=i-2; j>=1; j--) { out[3*j]=matrix(0,1,1); out[3*j-1]=intvec(0); out[3*j-2]=int(0); } break; } } } setring B; out=fetch(ringofSyz,out);//contains the V_d-strict resolution kill ringofSyz; } outall[1]=out; outall[2]=list(list(out[3*n-3],out[3*n-1])); return(outall); } /*case size(L)>2: We compute a quasi-isomorphic free complex following Alg 3.8 in [W2]*/ /* We denote the complex given by L as (C^i,d^i). We start by computing in the proc shortExaxtPieces representations for the short exact sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1), where the B^i, Z^i and H^i are coboundaries, cocycles and cohomology groups, respectively.*/ out=shortExactPieces(L); list rem; /* shortExactpiecesToVdStrict makes the sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1) V_d-strict*/ /* shortExactpiecesToVdStrict makes the sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1) V_d-strict*/ rem=shortExactPiecesToVdStrict(out,d,Syzstring); /*VdStrictDoubleComplexes computes V_d-strict resolutions over the seqeunces from proc shortExactPiecesToVdstrict*/ /*VdStrictDoubleComplexes computes V_d-strict resolutions over the seqeunces from proc shortExactPiecesToVdstrict*/ out=VdStrictDoubleComplexes(rem[1],d,Syzstring); for (i=1;i<=size(out); i++) { rem[2][i][1]=out[i][1][5][1]; rem[2][i][2]=out[i][1][8][1]; } /* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences Z^i->C^i->B^(i+1)) together to obtain a Cartan-Eilenberg resolution of (C^i,d^i)*/ { rem[2][i][1]=out[i][1][5][1]; rem[2][i][2]=out[i][1][8][1]; } /* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences Z^i->C^i->B^(i+1)) together to obtain a Cartan-Eilenberg resolution of (C^i,d^i)*/ out=assemblingDoubleComplexes(out); /*the proc totalComplex takes the total complex of the double complex from the proc assemblingDoubleComplexes*/ proc assemblingDoubleComplexes*/ out=totalComplex(out); outall[1]=out; outall[2]=rem[2];//contains the cohomology groups and their shift vectors outall[2]=rem[2];//contains the cohomology groups and their shift vectors return (outall); } int count; for (i=m; i<=n; i++) { out[size(out)+1]=list(); for (j=1; j<=size(L[i]); j++) { count=count+1; out[size(out)][j]=list(L[i][j],count); } } { out[size(out)+1]=list(); for (j=1; j<=size(L[i]); j++) { count=count+1; out[size(out)][j]=list(L[i][j],count); } } list o=list(out,count); return(o); //////////////////////////////////////////////////////////////////////////////////// static proc LMSubset(list L,list M) static proc LMSubset(list L,list M, list #) { int i; int j=1; list position=(M[size(M)],(-1)^(size(L))); if (size(#)==0) { list position=(M[size(M)],(-1)^(size(L))); } else { list position=(M[size(M)],1); } for (i=1; i<=size(L); i++) { if (L[i]!=M[j]) { if (L[i]!=M[i+1] or j!=i) { return (L[i],0); } else { position=(M[i],(-1)^(i-1)); j=j+1; } } j=j+1; } { if (L[i]!=M[j]) { if (L[i]!=M[i+1] or j!=i) { return (L[i],0); } else { if (size(#)==0) { position=(M[i],(-1)^(i-1)); } else { position=(M[i],(-1)^(size(L)+1-i)); } j=j+1; } } j=j+1; } return (position); } { /*we follow Section 3.3 in [W2]*/ /* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex  C=(C^i,d^i) as in the procedure toVdstrictcomplex*/ /* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex  C=(C^i,d^i) as in the procedure toVdstrictcomplex*/ matrix  Bnew= divdr(L[2],L[3]); matrix Bold=Bnew; list sep; /* the list sep will be of size s such that -sep[i]=(sep[i][1],sep[i][2]) is a list of two lists -sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)->coker(Z^i) ->coker(H^i) represents the short exact seqeuence B^i(C)->Z^i(C)->H^i(C) -sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)->coker(C^i)-> coker(B^(i+1)) represents the short exact seqeuence Z^i(C)->C^i->B^(i+1)(C)*/ -sep[i]=(sep[i][1],sep[i][2]) is a list of two lists -sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)->coker(Z^i) ->coker(H^i) represents the short exact seqeuence B^i(C)->Z^i(C)->H^i(C) -sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)->coker(C^i)-> coker(B^(i+1)) represents the short exact seqeuence Z^i(C)->C^i->B^(i+1)(C)*/ sep[1]=list(bzh,zcb); int i; list out; for (i=3; i<=size(L)-2; i=i+2) { /*the proc bzhzcb computes representations for the short exact seqeunces */ out=bzhzcb(Bold, L[i-1] , L[i], L[i+1], L[i+2]); sep[size(sep)+1]=out[1]; Bold=out[2]; } { /*the proc bzhzcb computes representations for the short exact seqeunces */ out=bzhzcb(Bold, L[i-1] , L[i], L[i+1], L[i+2]); sep[size(sep)+1]=out[1]; Bold=out[2]; } bzh=(divdr(L[size(L)-2], L[size(L)-1]),L[size(L)-2], L[size(L)-1]); bzh[4]=unitmat(ncols(L[size(L)-1])); static proc shortExactPiecesToVdStrict(list C,int d,list #) {/* We transform the short exact pieces from procedure shortExactPieces to V_d- strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in [W2].*/ strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in [W2].*/ /* If we compute our Groebner bases in the homogenized Weyl algebra, we already compute some resolutions it omit additional Groebner basis computations later on.*/ compute some resolutions it omit additional Groebner basis computations later on.*/ int s =size(C);int i; int j; string Syzstring="Sres"; intvec v=0:ncols(C[s][1][5]); intvec v=0:ncols(C[s][1][5]); if (size(#)!=0) { for (i=1; i<=size(#); i++) { if (typeof(#[i])=="string") { Syzstring=#[i]; } if (typeof(#[i])=="intvec") { v=#[i]; } } } { for (i=1; i<=size(#); i++) { if (typeof(#[i])=="string") { Syzstring=#[i]; } if (typeof(#[i])=="intvec") { v=#[i]; } } } list out; list forout; if (Syzstring=="Vdres") { out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s)); } { out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s)); } else { forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s); list resolutionofA=forout[9]; list resolutionofC=forout[10]; forout=delete(forout,10); forout=delete(forout,9); out[s]=list(forout); for (i=1; i<=size(resolutionofC); i++) { out[s][1][5][i+1]=resolutionofC[i];//save the resolutions out[s][1][1][i+1]=resolutionofA[i]; } } { forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s); list resolutionofA=forout[9]; list resolutionofC=forout[10]; forout=delete(forout,10); forout=delete(forout,9); out[s]=list(forout); for (i=1; i<=size(resolutionofC); i++) { out[s][1][5][i+1]=resolutionofC[i];//save the resolutions out[s][1][1][i+1]=resolutionofA[i]; } } out[s][2]=list(list(out[s][1][3][1])); out[s][2][2]=list(unitmat(ncols(out[s][1][3][1]))); list resolutionofF; for (i=s-1; i>=2; i--) { C[i][2][5]=out[i+1][1][1][1]; forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s); if (Syzstring=="Sres") { resolutionofD=forout[3];//save the resolutions resolutionofF=forout[4]; forout=delete(forout,4); forout=delete(forout,3); } out[i]=forout; if(Syzstring=="Sres") { for (j=2; j<=size(out[i+1][1][1]); j++) { out[i][2][5][j]=out[i+1][1][1][j]; } for (j=1; j<=size(resolutionofD);j++) { out[i][1][1][j+1]=resolutionofD[j]; out[i][1][5][j+1]=resolutionofF[j]; } } } { C[i][2][5]=out[i+1][1][1][1]; forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s); if (Syzstring=="Sres") { resolutionofD=forout[3];//save the resolutions resolutionofF=forout[4]; forout=delete(forout,4); forout=delete(forout,3); } out[i]=forout; if(Syzstring=="Sres") { for (j=2; j<=size(out[i+1][1][1]); j++) { out[i][2][5][j]=out[i+1][1][1][j]; } for (j=1; j<=size(resolutionofD);j++) { out[i][1][1][j+1]=resolutionofD[j]; out[i][1][5][j+1]=resolutionofF[j]; } } } out[1]=list(list());//initalize our list C[1][2][5]=out[2][1][1][1]; /*Compute the last V_d-strict seqeunce*/ if (Syzstring=="Vdres") { out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv"); } { out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv"); } else { forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv"); out[1][2]=delete(forout,9); list resolutionofA2=forout[9]; for (i=1; i<=size(out[2][1][1]); i++) { /*put the modules for the resolutions in the right spot*/ out[1][2][5][i]=out[2][1][1][i]; } for (i=1; i<=size(resolutionofA2); i++) { out[1][2][1][i+1]=resolutionofA2[i]; } } { forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv"); out[1][2]=delete(forout,9); list resolutionofA2=forout[9]; for (i=1; i<=size(out[2][1][1]); i++) { /*put the modules for the resolutions in the right spot*/ out[1][2][5][i]=out[2][1][1][i]; } for (i=1; i<=size(resolutionofA2); i++) { out[1][2][1][i+1]=resolutionofA2[i]; } } out[1][1][3]=list(out[1][2][1][1]); out[1][1][5]=list(out[1][2][1][1]); out[1][1][6]=list(list()); if (Syzstring=="Sres") { for (i=1; i<=size(out[1][2][1]); i++) { out[1][1][3][i]=out[1][2][1][i]; out[1][1][5][i]=out[1][2][1][i]; } } { for (i=1; i<=size(out[1][2][1]); i++) { out[1][1][3][i]=out[1][2][1][i]; out[1][1][5][i]=out[1][2][1][i]; } } list Hi; for (i=1; i<=size(out); i++) { Hi[i]=list(out[i][1][5][1],out[i][1][8][1]); } { Hi[i]=list(out[i][1][5][1],out[i][1][8][1]); } list outall; outall[1]=out; //////////////////////////////////////////////////////////////////////////////////// static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #) static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #) { /*this is the Algorithm 3.11 in [W2]*/ def B=basering; matrix bi=slimgb(transpose(C[5])); /* Computation of a V_d-strict Groebner basis of C[5]: -if Syzstring=="Vdres" this is done using the method of weighted homogenization (Prop. 3.9 [OT]) -else we use the homogenized Weyl algebra for Groebner basis computations (Prop 9.9 [OT]), in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit extra Groebner basis computations later on*/ /* Computation of a V_d-strict Groebner basis of C[5]: -if Syzstring=="Vdres" this is done using the method of weighted homogenization (Prop. 3.9 [OT]) -else we use the homogenized Weyl algebra for Groebner basis computations (Prop 9.9 [OT]), in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit extra Groebner basis computations later on*/ int nr,nc; intvec i1,i2; if (Syzstring=="Vdres") { if(size(#)==0) { matrix J_C=VdStrictGB(C[5],n,list(v)); } else { matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis } } { if(size(#)==0) { matrix J_C=VdStrictGB(C[5],n,list(v)); } else { matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis } } else { if (size(#)==0) { matrix MC=C[5]; def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v); setring HomWeyl; matrix J_C=fetch(B,MC); J_C=nHomogenize(J_C); /*computation of V_d-strict resolution of C[5]->needed for proc VdstrictDoubleComplexes*/ def ringofSyz=Sres(groebner(transpose(J_C)), lengthofres); setring ringofSyz; matrix J_C=transpose(matrix(RES[2])); J_C=subst(J_C,h,1); logens=ncols(J_C)+1; matrix zerom; list rofC;//will contain resolution of C for (i=3; i<=n+si+1; i++) { if (size(RES)>=i) { zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); if (RES[i]!=zerom) { rofC[i-2]=(matrix(RES[i])); if (i==3) { if (nrows(rofC[i-2])-logens+1!=nrows(J_C)) { //build the resolution nr=nrows(J_C)+logens-1; nc=ncols(rofC[i-2]); rofC[i-2]=matrix(rofC[i-2],nr,nc); } } if (i!=3) { if (nrows(rofC[i-2])-logens+1!=nrows(rofC[i-3])) { nr=nrows(rofC[i-3])+logens-1; nc=ncols(rofC[i-2]); rofC[i-2]=matrix(rofC[i-2],nr,nc); } } i1=intvec(logens..nrows(rofC[i-2])); i2=intvec(1..ncols(rofC[i-2])); rofC[i-2]=transpose(submat(rofC[i-2],i1,i2)); logens=logens+ncols(rofC[i-2]); rofC[i-2]=subst(rofC[i-2],h,1); } else { rofC[i-2]=list(); } } else { rofC[i-2]=list(); } } if(size(rofC[1])==0) { omitemptylist=1; } setring B; matrix  J_C=fetch(ringofSyz,J_C); if (omitemptylist!=1) { list rofC=fetch(ringofSyz,rofC); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } else { matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis } } { if (size(#)==0) { matrix MC=C[5]; def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v); setring HomWeyl; matrix J_C=fetch(B,MC); J_C=nHomogenize(J_C); /*computation of V_d-strict resolution of C[5]->needed for proc VdstrictDoubleComplexes*/ def ringofSyz=Sres(transpose(J_C),lengthofres); setring ringofSyz; matrix J_C=transpose(matrix(RES[2])); J_C=subst(J_C,h,1); logens=ncols(J_C)+1; matrix zerom; list rofC;//will contain resolution of C for (i=3; i<=n+si+1; i++) { if (size(RES)>=i) { zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); if (RES[i]!=zerom) { rofC[i-2]=(matrix(RES[i])); if (i==3) { if (nrows(rofC[i-2])-logens+1!=nrows(J_C)) { //build the resolution nr=nrows(J_C)+logens-1; nc=ncols(rofC[i-2]); rofC[i-2]=matrix(rofC[i-2],nr,nc); } } if (i!=3) { if (nrows(rofC[i-2])-logens+1!=nrows(rofC[i-3])) { nr=nrows(rofC[i-3])+logens-1; nc=ncols(rofC[i-2]); rofC[i-2]=matrix(rofC[i-2],nr,nc); } } i1=intvec(logens..nrows(rofC[i-2])); i2=intvec(1..ncols(rofC[i-2])); rofC[i-2]=transpose(submat(rofC[i-2],i1,i2)); logens=logens+ncols(rofC[i-2]); rofC[i-2]=subst(rofC[i-2],h,1); } else { rofC[i-2]=list(); } } else { rofC[i-2]=list(); } } if(size(rofC[1])==0) { omitemptylist=1; } setring B; matrix  J_C=fetch(ringofSyz,J_C); if (omitemptylist!=1) { list rofC=fetch(ringofSyz,rofC); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } else { matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis } } /* we compute a V_d-strict Groebner basis for C[3]*/ matrix J_A=C[1]; matrix Pi[1][ncols(J_AC)]; for (i=1; i<=nrows(J_C); i++) { for (j=1; j<=nrows(J_AC);j++) { Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC))); } storePi[i]=Pi; Pi=0; } { for (j=1; j<=nrows(J_AC);j++) { Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC))); } storePi[i]=Pi; Pi=0; } /*we compute the shift vector for C[1]*/ intvec m_a; list findMin; int comMin; int comMin; for (i=1; i<=ncols(J_A); i++) { for (j=1; j<=size(storePi);j++) { if (storePi[j][1,i]!=0) { comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v); comMin=comMin-VdDeg(storePi[j][1,i],n,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { m_a[i]=Min(findMin); findMin=list(); } else { m_a[i]=0; } } { for (j=1; j<=size(storePi);j++) { if (storePi[j][1,i]!=0) { comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v); comMin=comMin-VdDeg(storePi[j][1,i],n,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { m_a[i]=Min(findMin); findMin=list(); } else { m_a[i]=0; } } matrix zero[ncols(J_A)][ncols(J_C)]; matrix g_AB=concat(unitmat(ncols(J_A)),zero); matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C))))); intvec m_b=m_a,v; /* computation of a V_d-strict Groebner basis of C[1] (and resolution if Syzstring=='Vdres') */ /* computation of a V_d-strict Groebner basis of C[1] (and resolution if Syzstring=='Vdres') */ if (Syzstring=="Vdres") { J_A=VdStrictGB(J_A,n,m_a); } { J_A=VdStrictGB(J_A,n,m_a); } else { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a); setring HomWeyl; matrix J_A=fetch(B,J_A); J_A=nHomogenize(J_A); def ringofSyz=Sres(groebner(transpose(J_A)),lengthofres); setring ringofSyz; matrix J_A=transpose(matrix(RES[2])); matrix zerom; J_A=subst(J_A,h,1); logens=ncols(J_A)+1; list rofA; for (i=3; i<=n+si+1; i++) { if (size(RES)>=i) { zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); if (RES[i]!=zerom) { rofA[i-2]=matrix(RES[i]);// resolution for C[1] if (i==3) { if (nrows(rofA[i-2])-logens+1!=nrows(J_A)) { nr=nrows(J_A)+logens-1; nc=ncols(rofA[i-2]); rofA[i-2]=matrix(rofA[i-2],nr,nc); } } if (i!=3) { if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3])) { nr=nrows(rofA[i-3])+logens-1; nc=ncols(rofA[i-2]); rofA[i-2]=matrix(rofA[i-2],nr,nc); } } i1=intvec(logens..nrows(rofA[i-2])); i2=intvec(1..ncols(rofA[i-2])); rofA[i-2]=transpose(submat(rofA[i-2],i1,i2)); logens=logens+ncols(rofA[i-2]); rofA[i-2]=subst(rofA[i-2],h,1); } else { rofA[i-2]=list(); } } else { rofA[i-2]=list(); } } if(size(rofA[1])==0) { omitemptylist=1; } setring B; J_A=fetch(ringofSyz,J_A); if (omitemptylist!=1) { list rofA=fetch(ringofSyz,rofA); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a); setring HomWeyl; matrix J_A=fetch(B,J_A); J_A=nHomogenize(J_A); def ringofSyz=Sres(transpose(J_A),lengthofres); setring ringofSyz; matrix J_A=transpose(matrix(RES[2])); matrix zerom; J_A=subst(J_A,h,1); logens=ncols(J_A)+1; list rofA; for (i=3; i<=n+si+1; i++) { if (size(RES)>=i) { zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); if (RES[i]!=zerom) { rofA[i-2]=matrix(RES[i]);// resolution for C[1] if (i==3) { if (nrows(rofA[i-2])-logens+1!=nrows(J_A)) { nr=nrows(J_A)+logens-1; nc=ncols(rofA[i-2]); rofA[i-2]=matrix(rofA[i-2],nr,nc); } } if (i!=3) { if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3])) { nr=nrows(rofA[i-3])+logens-1; nc=ncols(rofA[i-2]); rofA[i-2]=matrix(rofA[i-2],nr,nc); } } i1=intvec(logens..nrows(rofA[i-2])); i2=intvec(1..ncols(rofA[i-2])); rofA[i-2]=transpose(submat(rofA[i-2],i1,i2)); logens=logens+ncols(rofA[i-2]); rofA[i-2]=subst(rofA[i-2],h,1); } else { rofA[i-2]=list(); } } else { rofA[i-2]=list(); } } if(size(rofA[1])==0) { omitemptylist=1; } setring B; J_A=fetch(ringofSyz,J_A); if (omitemptylist!=1) { list rofA=fetch(ringofSyz,rofA); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } J_AC=transpose(storePi[1]); for (i=2; i<= size(storePi); i++) { J_AC=concat(J_AC, transpose(storePi[i])); } { J_AC=concat(J_AC, transpose(storePi[i])); } J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC)); list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a)); Vdstrict[7]=list(m_b); Vdstrict[8]=list(v); Vdstrict[8]=list(v); if(Syzstring=="Sres") { Vdstrict[9]=rofA; if(size(#)==0) { Vdstrict[10]=rofC; } } { Vdstrict[9]=rofA; if(size(#)==0) { Vdstrict[10]=rofC; } } return (Vdstrict); } //////////////////////////////////////////////////////////////////////////////////// static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL) static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL) { /* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact pieces. We asume that we are given two sequences of the form coker(L[i][1])-> coker(L[i][3])->coker(L[i][5]) with differentials L[i][2] and L[i][4] such that L[1][3]=L[2][1].We are going to transform them to V_d-strict sequences J_D->J_A->J_F and J_A->J_B->J_C*/ /* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact pieces. We asume that we are given two sequences of the form coker(L[i][1])-> coker(L[i][3])->coker(L[i][5]) with differentials L[i][2] and L[i][4] such that L[1][3]=L[2][1].We are going to transform them to V_d-strict sequences J_D->J_A->J_F and J_A->J_B->J_C*/ int omitemptylist; int lengthofres=sizeL+d-1; matrix J_D=L[1][1]; matrix f_FA=L[1][4]; /*We find new presentations coker(J_DF) and coker(J_DFC)  for L[1][4]=L[2][1] and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */ /*We find new presentations coker(J_DF) and coker(J_DFC)  for L[1][4]=L[2][1] and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */ matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA))); matrix J_DF=divdr(f_DFA,L[1][3]);//coker(J_DF) is isomorphic to coker(L[2][1]); matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB))); matrix J_DFC=divdr(f_DFCB,L[2][3]);//coker(J_DFC) are coker(L[2][3)]) isomorphic /* find a shift vector on the range of J_F such that the first sequence is exact*/ /* find a shift vector on the range of J_F such that the first sequence is exact*/ matrix P=matrixLift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C); list storePi; int i; int j; for (i=1; i<=nrows(J_C); i++) { for (j=1; j<=nrows(J_DFC);j++) { Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC))); } storePi[i]=Pi; Pi=0; } { for (j=1; j<=nrows(J_DFC);j++) { Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC))); } storePi[i]=Pi; Pi=0; } intvec m_a; list findMin; intvec i1,i2; for (i=1; i<=ncols(J_DF); i++) { for (j=1; j<=size(storePi);j++) { if (storePi[j][1,i]!=0) { comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v); comMin=comMin-VdDeg(storePi[j][1,i],d,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { m_a[i]=Min(findMin);// shift vector for L[2][1] findMin=list(); noMin[i]=0; } else { noMin[i]=1; } } { for (j=1; j<=size(storePi);j++) { if (storePi[j][1,i]!=0) { comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v); comMin=comMin-VdDeg(storePi[j][1,i],d,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { m_a[i]=Min(findMin);// shift vector for L[2][1] findMin=list(); noMin[i]=0; } else { noMin[i]=1; } } if (size(m_a) < ncols(J_DF)) { m_a[ncols(J_DF)]=0; } { m_a[ncols(J_DF)]=0; } intvec m_f=m_a[ncols(J_D)+1..size(m_a)]; /* Computation of a V_d-strict Groebner basis of J_F=L[1][5]: if Syzstring=="Vdres" this is done using the method of weighted homogenization (Prop. 3.9 [OT]) else we use the homogenized Weyl algerba for Groebner basis computations (Prop 9.9 [OT]), in this case we already compute resolutions (Thm. 9.10 in [OT]) to omit extra Groebner basis  computations later on*/ /* Computation of a V_d-strict Groebner basis of J_F=L[1][5]: if Syzstring=="Vdres" this is done using the method of weighted homogenization (Prop. 3.9 [OT]) else we use the homogenized Weyl algerba for Groebner basis computations (Prop 9.9 [OT]), in this case we already compute resolutions (Thm. 9.10 in [OT]) to omit extra Groebner basis  computations later on*/ if (Syzstring=="Vdres") { J_F=VdStrictGB(J_F,d,m_f); } { J_F=VdStrictGB(J_F,d,m_f); } else { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f); setring HomWeyl; matrix J_F=fetch(B,J_F); J_F=nHomogenize(J_F); def ringofSyz=Sres(groebner(transpose(J_F)),lengthofres); setring ringofSyz; matrix J_F=transpose(matrix(RES[2])); J_F=subst(J_F,h,1); logens=ncols(J_F)+1; list rofF; for (i=3; i<=d+sizeL+1; i++) { if (size(RES)>=i) { if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])))) { rofF[i-2]=(matrix(RES[i]));// resolution for J_F if (i==3) { if (nrows(rofF[i-2])-logens+1!=nrows(J_F)) { nr=nrows(J_F)+logens-1; nc=ncols(rofF[i-2]); rofF[i-2]=matrix(rofF[i-2],nr,nc); } } if (i!=3) { if (nrows(rofF[i-2])-logens+1!=nrows(rofF[i-3])) { nr=nrows(rofF[i-3])+logens-1; rofF[i-2]=matrix(rofF[i-2],nr,ncols(rofF[i-2])); } } i1=intvec(logens..nrows(rofF[i-2])); i2=intvec(1..ncols(rofF[i-2])); rofF[i-2]=transpose(submat(rofF[i-2],i1,i2)); logens=logens+ncols(rofF[i-2]); rofF[i-2]=subst(rofF[i-2],h,1); } else { rofF[i-2]=list(); } } else { rofF[i-2]=list(); } } if(size(rofF[1])==0) { omitemptylist=1; } setring B; J_F=fetch(ringofSyz,J_F); if (omitemptylist!=1) { list rofF=fetch(ringofSyz,rofF); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f); setring HomWeyl; matrix J_F=fetch(B,J_F); J_F=nHomogenize(J_F); def ringofSyz=Sres(transpose(J_F),lengthofres); setring ringofSyz; matrix J_F=transpose(matrix(RES[2])); J_F=subst(J_F,h,1); logens=ncols(J_F)+1; list rofF; for (i=3; i<=d+sizeL+1; i++) { if (size(RES)>=i) { if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])))) { rofF[i-2]=(matrix(RES[i]));// resolution for J_F if (i==3) { if (nrows(rofF[i-2])-logens+1!=nrows(J_F)) { nr=nrows(J_F)+logens-1; nc=ncols(rofF[i-2]); rofF[i-2]=matrix(rofF[i-2],nr,nc); } } if (i!=3) { if (nrows(rofF[i-2])-logens+1!=nrows(rofF[i-3])) { nr=nrows(rofF[i-3])+logens-1; rofF[i-2]=matrix(rofF[i-2],nr,ncols(rofF[i-2])); } } i1=intvec(logens..nrows(rofF[i-2])); i2=intvec(1..ncols(rofF[i-2])); rofF[i-2]=transpose(submat(rofF[i-2],i1,i2)); logens=logens+ncols(rofF[i-2]); rofF[i-2]=subst(rofF[i-2],h,1); } else { rofF[i-2]=list(); } } else { rofF[i-2]=list(); } } if(size(rofF[1])==0) { omitemptylist=1; } setring B; J_F=fetch(ringofSyz,J_F); if (omitemptylist!=1) { list rofF=fetch(ringofSyz,rofF); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } /*find shift vectors on the range of J_D*/ P=matrixLift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F); matrix Pidf[1][ncols(J_DF)]; for (i=1; i<=nrows(J_F); i++) { for (j=1; j<=nrows(J_DF);j++) { Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF))); } storePinew[i]=Pidf; Pidf=0; } { for (j=1; j<=nrows(J_DF);j++) { Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF))); } storePinew[i]=Pidf; Pidf=0; } intvec m_d; for (i=1; i<=ncols(J_D); i++) { for (j=1; j<=size(storePinew);j++) { if (storePinew[j][1,i]!=0) { comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f); comMin=comMin-VdDeg(storePinew[j][1,i],d,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { if (noMin[i]==0) { m_d[i]=Min(insert(findMin,m_a[i])); m_a[i]=m_d[i]; } { for (j=1; j<=size(storePinew);j++) { if (storePinew[j][1,i]!=0) { comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f); comMin=comMin-VdDeg(storePinew[j][1,i],d,intvec(0)); findMin[size(findMin)+1]=comMin; } } if (size(findMin)!=0) { if (noMin[i]==0) { m_d[i]=Min(insert(findMin,m_a[i])); m_a[i]=m_d[i]; } else { m_d[i]=Min(findMin); m_a[i]=m_d[i]; } } else { m_d[i]=Min(findMin); m_a[i]=m_d[i]; } } else { m_d[i]=m_a[i]; } findMin=list(); } /* compute a V_d-strict Groebner basis (and resolution of J_D if Syzstring!='Vdres') for J_D*/ { m_d[i]=m_a[i]; } findMin=list(); } /* compute a V_d-strict Groebner basis (and resolution of J_D if Syzstring!='Vdres') for J_D*/ if (Syzstring=="Vdres") { J_D=VdStrictGB(J_D,d,m_d); } { J_D=VdStrictGB(J_D,d,m_d); } else { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d); setring HomWeyl; matrix J_D=fetch(B,J_D); J_D=nHomogenize(J_D); def ringofSyz=Sres(groebner(transpose(J_D)),lengthofres); setring ringofSyz; matrix J_D=transpose(matrix(RES[2])); J_D=subst(J_D,h,1); logens=ncols(J_D)+1; list rofD; for (i=3; i<=d+sizeL+1; i++) { if (size(RES)>=i) { if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])))) { rofD[i-2]=(matrix(RES[i]));// resolution for J_D if (i==3) { if (nrows(rofD[i-2])-logens+1!=nrows(J_D)) { nr=nrows(J_D)+logens-1; rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2])); } } if (i!=3) { if (nrows(rofD[i-2])-logens+1!=nrows(rofD[i-3])) { nr=nrows(rofD[i-3])+logens-1; rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2])); } } i1=intvec(logens..nrows(rofD[i-2])); i2=intvec(1..ncols(rofD[i-2])); rofD[i-2]=transpose(submat(rofD[i-2],i1,i2)); logens=logens+ncols(rofD[i-2]); rofD[i-2]=subst(rofD[i-2],h,1); } else { rofD[i-2]=list(); } } else { rofD[i-2]=list(); } } if(size(rofD[1])==0) { omitemptylist=1; } setring B; J_D=fetch(ringofSyz,J_D); if (omitemptylist!=1) { list rofD=fetch(ringofSyz,rofD); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } /* compute new matrices for J_A and J_B  such that their rows form a V_d-strict Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and nrows(J_B)=nrows(J_A)+nrows(J_C)*/ { def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d); setring HomWeyl; matrix J_D=fetch(B,J_D); J_D=nHomogenize(J_D); def ringofSyz=Sres(transpose(J_D),lengthofres); setring ringofSyz; matrix J_D=transpose(matrix(RES[2])); J_D=subst(J_D,h,1); logens=ncols(J_D)+1; list rofD; for (i=3; i<=d+sizeL+1; i++) { if (size(RES)>=i) { if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])))) { rofD[i-2]=(matrix(RES[i]));// resolution for J_D if (i==3) { if (nrows(rofD[i-2])-logens+1!=nrows(J_D)) { nr=nrows(J_D)+logens-1; rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2])); } } if (i!=3) { if (nrows(rofD[i-2])-logens+1!=nrows(rofD[i-3])) { nr=nrows(rofD[i-3])+logens-1; rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2])); } } i1=intvec(logens..nrows(rofD[i-2])); i2=intvec(1..ncols(rofD[i-2])); rofD[i-2]=transpose(submat(rofD[i-2],i1,i2)); logens=logens+ncols(rofD[i-2]); rofD[i-2]=subst(rofD[i-2],h,1); } else { rofD[i-2]=list(); } } else { rofD[i-2]=list(); } } if(size(rofD[1])==0) { omitemptylist=1; } setring B; J_D=fetch(ringofSyz,J_D); if (omitemptylist!=1) { list rofD=fetch(ringofSyz,rofD); } omitemptylist=0; kill HomWeyl; kill ringofSyz; } /* compute new matrices for J_A and J_B  such that their rows form a V_d-strict Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and nrows(J_B)=nrows(J_A)+nrows(J_C)*/ J_DF=transpose(storePinew[1]); for (i=2; i<=nrows(J_F); i++) { J_DF=concat(J_DF,transpose(storePinew[i])); } { J_DF=concat(J_DF,transpose(storePinew[i])); } J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF)); J_DFC=transpose(storePi[1]); for (i=2; i<=nrows(J_C); i++) { J_DFC=concat(J_DFC,transpose(storePi[i])); } { J_DFC=concat(J_DFC,transpose(storePi[i])); } J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC)); intvec m_b=m_a,v; matrix g_AB=concat(unitmat(ncols(J_DF)),zero1); matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C)))); list out; list out; out[1]=list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F)); out[1]=out[1]+list(list(m_d),list(m_a),list(m_f)); out[2]=out[2]+list(list(m_a),list(m_b),list(v)); if (Syzstring=="Sres") { out[3]=rofD; out[4]=rofF; } { out[3]=rofD; out[4]=rofF; } return(out); } static proc VdStrictDoubleComplexes(list L,int d,string Syzstring) { /* We compute  V_d-strict resolutions over the V_d-strict short exact pieces from the procedure shortExactPiecesToVdStrict. We use Algorithms 3.14 and 3.15 in [W2]*/ /* We compute  V_d-strict resolutions over the V_d-strict short exact pieces from the procedure shortExactPiecesToVdStrict. We use Algorithms 3.14 and 3.15 in [W2]*/ int i,k,c,j,l,totaldeg,comparedegs,SBcom,verk; intvec fordegs; intvec fordegs; intvec n_b,i1,i2; matrix rem,forML,subm,zerom,unitm,subm2; list forhW; if (Syzstring=="Sres") { /*we already computed some of the resolutions in the procedure shortExactPiecesToVdStrict*/ matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew; for (k=1; k<=(size(L)+d-1); k++) { { /*we already computed some of the resolutions in the procedure shortExactPiecesToVdStrict*/ matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew; for (k=1; k<=(size(L)+d-1); k++) { L[1][1][1][k+1]=list(); L[1][1][2][k+1]=list(); L[1][1][6][k+1]=list(); } L[1][1][6][size(L)+d+1]=list(); matrix mem; for (i=2; i<=d+size(L)+1; i++) {; v=0; if(size(L[1][1][3][i-1])!=0) { if(i!=d+size(L)+1) { /*horizontal differential*/ L[1][1][4][i-1]=unitmat(nrows(L[1][1][3][i-1])); } for (j=1; j<=nrows(L[1][1][3][i-1]); j++) { mem=submat(L[1][1][3][i-1],j,intvec(1..ncols(L[1][1][3][i-1]))); v[j]=VdDeg(mem,d,L[1][1][7][i-1]); } L[1][1][7][i]=v;//new shift vector L[1][1][8][i]=v; L[1][2][6][i]=v; } else { if (i!=d+size(L)+1) { L[1][1][4][i-1]=list(); } L[1][1][7][i]=list(); L[1][1][8][i]=list(); L[1][2][6][i]=list(); } } if (size(L[1][1][3][d+size(L)])!=0) { /*horizontal differential*/ L[1][1][4][d+size(L)]=unitmat(nrows(L[1][1][3][d+size(L)])); } else { L[1][1][4][d+size(L)]=list(); } for (k=1; k coker(L[k][2][3][1])->coker(L[k][2][5][1]) using the resolution obtained for coker(L[k][1][3][1]). L[k][2][i][j] will be the jth module in the resolution of L[k][2][i][1] for i=1,3,5. L[k][2][i+5][j] will be the jth  shift vector in the resolution of L[k][2][i][1](this holds also for the case Syzstring=="Vdres")*/ for (i=2; i<=d+size(L); i++) { v=0; if (size(L[k][2][5][i-1])!=0) { for (j=1; j<=nrows(L[k][2][5][i-1]); j++) { i1=intvec(1..ncols(L[k][2][5][i-1])); mem=submat(L[k][2][5][i-1],j,i1); v[j]=VdDeg(mem,d,L[k][2][8][i-1]); } /*next shift vector in th resolution of coker(L[k][2][5][1])*/ L[k][2][8][i]=v; } else { L[k][2][8][i]=list(); } /* we build step by step a resolution for coker(L[k][2][5][1]) using the resolutions of coker(L[k][2][1][1]) and coker(L[k][2][5][1])*/ if (size(L[k][2][5][i])!=0) { if (size(L[k][2][1][i])!=0 or size(L[k][2][1][i-1])!=0) { L[k][2][3][i]=transpose(syz(transpose(L[k][2][3][i-1]))); nr= nrows(L[k][2][1][i-1]); nc=ncols(L[k][2][5][i]); Pold=matrixLift(L[k][2][3][i]*prodr(nr,nc), L[k][2][5][i]); matrix Pi[1][ncols(L[k][2][3][i])]; for (l=1; l<=nrows(L[k][2][5][i]); l++) { for (j=1; j<=nrows(L[k][2][3][i]); j++) { i2=intvec(1..ncols(L[k][2][3][i])); Pi=Pi+Pold[l,j]*submat(L[k][2][3][i],j,i2); } if (l==1) { Picombined=transpose(Pi); } else { Picombined=concat(Picombined,transpose(Pi)); } Pi=0; } kill Pi; Picombined=transpose(Picombined); if (size(L[k][2][1][i])!=0) { if (i==2) { containsndeg=(0:ncols(L[k][2][1][1])); } containsndeg=nDeg(L[k][2][1][i-1],containsndeg); forhW=list(L[k][2][6][i],containsndeg); def HomWeyl=makeHomogenizedWeyl(n,forhW); setring HomWeyl; list L=fetch(B,L); matrix M=L[k][2][1][i]; module Mmod; list forM=nHomogenize(M,containsndeg,1); M=forM[1]; totaldeg=forM[2]; kill forM; matrix Maorig=fetch(B,Picombined); matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M))); matrix mem,subm,zerom; matrix Pinew; M=transpose(M); SBcom=0; for (l=1; l<=nrows(Ma); l++) { zerom=matrix(0,1,(ncols(Maorig)-ncols(Ma))); i1=(ncols(Ma)+1..ncols(Maorig)); if (submat(Maorig,l,i1)==zerom) { for (cc=1; cc<=ncols(Ma); cc++) { Maorig[l,cc]=0; } } i2=(ncols(Ma)+1..ncols(Maorig)); i1=(1..ncols(Ma)); if (VdDeg(submat(Maorig,l,i1),d,L[k][2][6][i])> VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]) and submat(Maorig,l,i1)!=matrix(0,1,ncols(Ma))) { /*V_d-Grad is to big--> we make it smaller using Vdnormal form computations*/ if (SBcom==0) { Mmod=slimgb(M); M=Mmod; SBcom=1; } //print("Reduzierung des V_d-Grades(Stelle1)"); i2=(ncols(Ma)+1..ncols(Maorig)); vd1=VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]); mem=submat(Ma,l,(1..ncols(Ma))); mem=nHomogenize(mem,containsndeg); mem=h^totaldeg*mem; mem=transpose(mem); mem=reduce(mem,Mod);////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// matrix jt=transpose(subst(mem,h,1)); setring B; matrix jt=fetch(HomWeyl,jt); matrix need=fetch(HomWeyl,Maorig); need=submat(need,l,(1..ncols(need))); i1=L[k][2][6][i]; i2=L[k][2][8][i]; jt=VdNormalForm(need,L[k][2][1][i],d,i1,i2); setring HomWeyl; mem=fetch(B,jt); mem=transpose(mem); if (l==1) { Pinew=mem; } else { Pinew=concat(Pinew,mem); } vd2=VdDeg(transpose(mem),d,L[k][2][6][i]); if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem))) {//should not happen!! //print("Reduzierung fehlgeschlagen!!(Stelle1)"); } } else { if (l==1) { Pinew=transpose(submat(Ma,l,(1..ncols(Ma)))); } else { subm=transpose(submat(Ma,l,(1..ncols(Ma)))); Pinew=concat(Pinew,subm); } } } Pinew=subst(Pinew,h,1); Pinew=transpose(Pinew); setring B; Pinew=fetch(HomWeyl,Pinew); kill HomWeyl; L[k][2][3][i]=concat(Pinew,L[k][2][5][i]); subm=transpose(L[k][2][3][i]); subm=concat(transpose(L[k][2][1][i]),subm); L[k][2][3][i]=transpose(subm); } else { L[k][2][3][i]=Picombined; } L[k+1][1][1][i]=L[k][2][5][i]; nr=nrows(L[k][2][1][i-1]); nc=ncols(L[k][2][5][i]); L[k][2][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nc); v=L[k][2][6][i],L[k][2][8][i]; L[k][2][7][i]=v; L[k+1][1][6][i]=L[k][2][8][i]; } else { L[k][2][3][i]=L[k][2][5][i]; L[k][2][2][i]=list(); L[k][2][7][i]=L[k][2][8][i]; L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1])); L[k+1][1][6][i]=L[k][2][8][i]; L[k+1][1][1][i]=L[k][2][5][i]; } } else { if (size(L[k][2][1][i])!=0) { if (size(L[k][2][5][i-1])!=0) { nr=nrows(L[k][2][5][i-1]); L[k][2][3][i]=concat(L[k][2][1][i],matrix(0,1,nr)); v=L[k][2][6][i],L[k][2][8][i]; L[k][2][7][i]=v; nc=nrows(L[k][2][1][i-1]); L[k][2][2][i]=concat(unitmat(nc),matrix(0,nc,nr)); L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nr); } else { L[k][2][3][i]=L[k][2][1][i]; L[k][2][7][i]=L[k][2][6][i]; L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1])); L[k][2][4][i]=list(); } L[k+1][1][1][i]=L[k][2][5][i]; L[k+1][1][6][i]=L[k][2][8][i]; } else { L[k][2][3][i]=list(); if (size(L[k][2][6][i])!=0) { if (size(L[k][2][8][i])!=0) { v=L[k][2][6][i],L[k][2][8][i]; L[k][2][7][i]=v; nr=nrows(L[k][2][1][i-1]); nc=nrows(L[k][2][5][i-1]); L[k][2][2][i]=concat(unitmat(nc),matrix(0,nr,nc)); L[k][2][4][i]=prodr(nr,nrows(L[k][2][5][i-1])); } else { L[k][2][7][i]=L[k][2][6][i]; L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1])); L[k][2][4][i]=list(); } } else { if (size(L[k][2][8][i])!=0) { L[k][2][7][i]=L[k][2][8][i]; L[k][2][2][i]=list(); L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1])); } else { L[k][2][7][i]=list(); L[k][2][2][i]=list(); L[k][2][4][i]=list(); } } L[k+1][1][1][i]=L[k][2][5][i]; L[k+1][1][6][i]=L[k][2][8][i]; } } } i=d+size(L)+1; v=0; if (size(L[k][2][5][i-1])!=0) { for (j=1; j<=nrows(L[k][2][5][i-1]); j++) { mem=submat(L[k][2][5][i-1],j,intvec(1..ncols(L[k][2][5][i-1]))); v[j]=VdDeg(mem,d,L[k][2][8][i-1]); } L[k][2][8][i]=v; if (size(L[k][2][6][i])!=0) { v=L[k][2][6][i],L[k][2][8][i]; L[k][2][7][i]=v; } else { L[k][2][7][i]=L[k][2][8][i]; } } else { L[k][2][8][i]=list(); L[k][2][7][i]=L[k][2][6][i]; } L[k+1][1][6][i]=L[k][2][8][i]; /* now we build V_d-strict resolutions for the sequences coker(L[k+1][1][1][1])->coker(L[k+1][1][3][1])->coker(L[k+1][1][5][i]) using the resolutions  for coker(L[k][2][5][1]) we just obtained (works exactly the same as above)*/ for (i=2; i<=d+size(L); i++) { v=0; if (size(L[k+1][1][5][i-1])!=0) { for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++) { i1=intvec(1..ncols(L[k+1][1][5][i-1])); mem=submat(L[k+1][1][5][i-1],j,i1); v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]); } L[k+1][1][8][i]=v; } else { L[k+1][1][8][i]=list(); } if (size(L[k+1][1][5][i])!=0) { if (size(L[k+1][1][1][i])!=0 or size(L[k+1][1][1][i-1])!=0) { L[k+1][1][3][i]=transpose(syz(transpose(L[k+1][1][3][i-1]))); nr=nrows(L[k+1][1][1][i-1]); nc=ncols(L[k+1][1][5][i]); Pold=matrixLift(L[k+1][1][3][i]*prodr(nr,nc),L[k+1][1][5][i]); matrix Pi[1][ncols(L[k+1][1][3][i])]; for (l=1; l<=nrows(L[k+1][1][5][i]); l++) { for (j=1; j<=nrows(L[k+1][1][3][i]); j++) { i2=intvec(1..ncols(L[k+1][1][3][i])); Pi=Pi+Pold[l,j]*submat(L[k+1][1][3][i],j,i2); } if (l==1) { Picombined=transpose(Pi); } else { Picombined=concat(Picombined,transpose(Pi)); } Pi=0; } kill Pi; Picombined=transpose(Picombined); if(size(L[k+1][1][1][i])!=0) { if (i==2) { containsndeg=(0:ncols(L[k+1][1][1][i-1])); } containsndeg=nDeg(L[k+1][1][1][i-1],containsndeg); forhW=list(L[k+1][1][6][i], containsndeg); def HomWeyl=makeHomogenizedWeyl(n,forhW); setring HomWeyl; list L=fetch(B,L); matrix M=L[k+1][1][1][i]; module Mmod; list forM=nHomogenize(M,containsndeg,1); M=forM[1]; totaldeg=forM[2]; kill forM; matrix Maorig=fetch(B,Picombined); matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M))); Ma=nHomogenize(Ma,containsndeg); matrix mem,subm,zerom,subm2; matrix Pinew; M=transpose(M); SBcom=0; for (l=1; l<=nrows(Ma); l++) { i2=(ncols(Ma)+1..ncols(Maorig)); nc=ncols(Maorig)-ncols(Ma); if (submat(Maorig,l,i2)==matrix(0,1,nc)) { for (cc=1; cc<=ncols(Ma); cc++) { Maorig[l,cc]=0; } } i1=(1..ncols(Ma)); i2=L[k+1][1][8][i]; subm=submat(Maorig,l,i1); subm2=submat(Maorig,l,(ncols(Ma)+1..ncols(Maorig))); if (VdDeg(subm,d,L[k+1][1][6][i])>VdDeg(subm2,d,i2) and subm!=matrix(0,1,ncols(Ma))) { //print("Reduzierung des Vd-Grades (Stelle2)"); if (SBcom==0) { Mmod=slimgb(M); M=Mmod; SBcom=1; } vd1=VdDeg(subm2,d,L[k+1][1][8][i]); mem=submat(Ma,l,(1..ncols(Ma))); mem=nHomogenize(mem,containsndeg); mem=h^totaldeg*mem; mem=transpose(mem); mem=reduce(mem,Mmod); if (l==1) { Pinew=mem; } else { Pinew=concat(Pinew,mem); } vd2=VdDeg(transpose(mem),d,L[k+1][1][6][i]); if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem))) {//should not happen //print("Reduzierung fehlgeschlagen!!!!(Stelle2)"); } } else { if (l==1) { Pinew=transpose(submat(Ma,l,(1..ncols(Ma)))); } else { subm=transpose(submat(Ma,l,(1..ncols(Ma)))); Pinew=concat(Pinew,subm); } } } Pinew=subst(Pinew,h,1); Pinew=transpose(Pinew); setring B; Pinew=fetch(HomWeyl,Pinew); kill HomWeyl; L[k+1][1][3][i]=concat(Pinew,L[k+1][1][5][i]); subm=transpose(L[k+1][1][1][i]); subm2=transpose(L[k+1][1][3][i]); L[k+1][1][3][i]=transpose(concat(subm,subm2)); } else { L[k+1][1][3][i]=Picombined; } L[k+1][2][1][i]=L[k+1][1][3][i]; nr=nrows(L[k+1][1][1][i-1]); nc=ncols(L[k+1][1][5][i]); L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); L[k+1][1][4][i]=prodr(nr,nc); v=L[k+1][1][6][i],L[k+1][1][8][i]; L[k+1][1][7][i]=v; L[k+1][2][6][i]=L[k+1][1][7][i]; } else { L[k+1][1][3][i]=L[k+1][1][5][i]; L[k+1][1][2][i]=list(); L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1])); L[k+1][1][7][i]=L[k+1][1][8][i]; L[k+1][2][6][i]=L[k+1][1][7][i]; L[k+1][2][1][i]=L[k+1][1][3][i]; } } else { if (size(L[k+1][1][1][i])!=0) { if (size(L[k+1][1][5][i-1])!=0) { zerom=matrix(0,1,nrows(L[k+1][1][5][i-1])); L[k+1][1][3][i]=concat(L[k+1][1][1][i],zerom); v=L[k+1][1][6][i],L[k+1][1][8][i]; L[k+1][1][7][i]=v; nr=nrows(L[k+1][1][1][i-1]); nc=nrows(L[k+1][1][5][i-1]); L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); L[k+1][1][4][i]=prodr(nr,nc); } else { L[k+1][1][3][i]=L[k+1][1][1][i]; &