//////////////////////////////////////////////////////////////////////////////////// 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 (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, 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) PROCEDURES: deRhamCohomology(list[,opt]); computes the de Rham cohomology MVComplex(list); computes the Mayer-Vietoris complex "; LIB "nctools.lib"; LIB "matrix.lib"; LIB "qhmoduli.lib"; LIB "general.lib"; LIB "dmod.lib"; LIB "bfun.lib"; LIB "dmodapp.lib"; LIB "poly.lib"; LIB "schreyer.lib"; //////////////////////////////////////////////////////////////////////////////////// proc deRhamCohomology(list L,list #) "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 and without using Schreyer's method @* -'Sres': compute the Cartan-Eilenberg resolutions in the homogenized Weyl algebra using Schreyer's method@* one of the strings@* -'iterativeloc': compute localizations by factorizing the polynomials and sucessive localization of the factors @* -'no iterativeloc': compute localizations by directly localizing the product@* and one of the strings -'onlybounds': computes bounds for the minimal and maximal interger roots of the global b-function -'exactroots' computes the minimal and maximal integer root of the global b-function The default is 'Sres', '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 EXAMPLE:example deRhamCohomology; shows an example " { intvec saveoptions=option(get); intvec i1,i2; option(none); int recursiveloc=1; int i,j,nr,nc; def R=basering; poly islcm, forlcm; int n=nvars(R); int le=size(L)+n; string Syzstring="Sres"; int onlybounds=1; for (i=1; i<=size(#); i++) { if (#[i]=="Vdres") { Syzstring="Vdres"; } if (#[i]=="noiterativeloc") { recursiveloc=0; } if (#[i]=="exactroots") { onlybounds=0; } } for (i=1; i<=size(L); i++) { if (L[i]==0) { L=delete(L,i); i=i-1; } } if (size(L)==0) { return (list(0)); } for (i=1; i<= size(L); i++) { if (leadcoef(L[i])-L[i]==0) { return(list(1)); } } if (size(L)==0) { /*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(); } } /* 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); 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); } 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]*/ 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; } } } } int iupto=Max(upto);//maximal degree +1 of the polynomials we have to consider if (iupto<=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; 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; } } } list keepformatrix,sizetruncom,fortrun,fst; int count,stc; intvec v,forin; matrix subm; /*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); } } } matrix M; int st,pi,pj; poly ptc; int b,d,ideg,kplus,lplus; poly form,lform,nform; /*computation of the maps*/ for (i=1; i=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; } } } } } } 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); /*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); } example { "EXAMPLE:"; ring r = 0,(x,y,z),dp; list L=(xy,xz); deRhamCohomology(L); } //////////////////////////////////////////////////////////////////////////////////// //COMPUTATION OF THE MAYER-VIETORIS COMPLEX //////////////////////////////////////////////////////////////////////////////////// proc MVComplex(list L,list #) "USAGE:MVComplex(L); L a list of polynomials ASSUME: -Basering is a polynomial ring with n vwariables and rational coefficients -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 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 d^i are given by C[2*i] EXAMPLE:example MVComplex; shows an example " { /* We follow algorithm 3.2.5 in [R],if #!=0 we use also Remark 3.2.6 in [R] for an additional iterative localization*/ def R=basering; int i; int iterative=1; if (size(#)!=0) { 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 (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); } } int r=size(L); int n=nvars(basering); int le=size(L)+n; /*construct the ring Ws*/ 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 Ws; 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); } Cech[1]=list(J); list Theta, remminroots; Theta[1]=list(list(list(),1,1)); list rem,findminintroot,diffmaps; int minroot,st,sk; intvec k1; poly fred,forfetch; matrix subm; 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=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; 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 @* Syzstring is either: @* -'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 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 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 Algorithm 3.8 in [W2] " { def B=basering; int n=nvars(B) div 2+2; int d=nvars(B) div 2; intvec v; list out, outall; int i,j,k,indi,nc,nr; matrix mem; intvec i1,i2; if (size(#)!=0) { for (i=1; i<=size(#); i++) { if (typeof(#[i])=="int") { if (#[1]>=1 and #[1]<=n) { d=#[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); list rem; /* 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*/ 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)*/ out=assemblingDoubleComplexes(out); /*the proc totalComplex takes the total complex of the double complex from the proc assemblingDoubleComplexes*/ out=totalComplex(out); outall[1]=out; outall[2]=rem[2];//contains the cohomology groups and their shift vectors return (outall); } //////////////////////////////////////////////////////////////////////////////////// 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 LMSubset(list L,list 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[i+1] or j!=i) { return (L[i],0); } else { position=(M[i],(-1)^(i-1)); j=j+1; } } j=j+1; } return (position); } //////////////////////////////////////////////////////////////////////////////////// static proc shortExactPieces(list L) { /*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*/ matrix Bnew= divdr(L[2],L[3]); matrix Bold=Bnew; matrix Z=divdr(Bnew,L[1]); list bzh,zcb; bzh=list(list(),list(),Z,unitmat(ncols(Z)),Z); zcb=(Z, Bnew, L[1], unitmat(ncols(L[1])), 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[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]; } 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])); bzh[5]=transpose(concat(transpose(L[size(L)-2]),transpose(L[size(L)-1]))); zcb=(L[size(L)-1], unitmat(ncols(L[size(L)-1])), L[size(L)-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); matrix H=transpose(concat(transpose(lift1),transpose(Z))); list bzh=(Bold, lift1, Z, unitmat(ncols(Z)),H); list zcb=(Z, Bnew, C1, unitmat(ncols(C1)),Bnew); list out=(list(bzh, zcb), Bnew); return(out); } //////////////////////////////////////////////////////////////////////////////////// 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].*/ /* If we compute our Groebner bases in the homogenized Weyl algebra, we already 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]); if (size(#)!=0) { 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)); } 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]; } } out[s][2]=list(list(out[s][1][3][1])); out[s][2][2]=list(unitmat(ncols(out[s][1][3][1]))); out[s][2][3]=list(out[s][1][3][1]); out[s][2][4]=list(list()); out[s][2][5]=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()); list resolutionofD; 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]; } } } 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"); } 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]; } } 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()); 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]; } } 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; outall[2]=Hi; return(outall); } //////////////////////////////////////////////////////////////////////////////////// static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #) { /*this is the Algorithm 3.11 in [W2]*/ int omitemptylist; int lengthofres=si+n-1; int i,j,logens; 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*/ 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 } } 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 } } /* we compute a V_d-strict Groebner basis for C[3]*/ 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)]; 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; } /*we compute the shift vector for C[1]*/ 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); 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') */ if (Syzstring=="Vdres") { 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; } 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)); Vdstrict[7]=list(m_b); Vdstrict[8]=list(v); if(Syzstring=="Sres") { Vdstrict[9]=rofA; if(size(#)==0) { Vdstrict[10]=rofC; } } return (Vdstrict); } //////////////////////////////////////////////////////////////////////////////////// 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*/ int omitemptylist; int lengthofres=sizeL+d-1; int logens; def B=basering; matrix J_F=L[1][5]; 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]) */ 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 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]);//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*/ 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; int nr,nc; 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; } } 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)]; /* 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); } 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; } /*find shift vectors on the range of J_D*/ P=matrixLift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F); 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); 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]=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); } 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)*/ 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; 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]=list(list(J_DF),list(g_AB),list(J_DFC),list(g_BC),list(J_C)); out[2]=out[2]+list(list(m_a),list(m_b),list(v)); if (Syzstring=="Sres") { 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]*/ int i,k,c,j,l,totaldeg,comparedegs,SBcom,verk; intvec fordegs; intvec n_b,i1,i2; matrix rem,forML,subm,zerom,unitm,subm2; matrix J_B; list store; int t=size(L)+d; int vd1,vd2,nr,nc; def B=basering; int n=nvars(B) div 2; intvec v; 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++) { 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]; 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) { M=slimgb(M); 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,groebner(M)); 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]; 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) { M=slimgb(M); 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, groebner(M)); 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]; L[k+1][1][7][i]=L[k+1][1][6][i]; L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1])); L[k+1][1][4][i]=list(); } L[k+1][2][1][i]=L[k+1][1][3][i]; L[k+1][2][6][i]=L[k+1][1][7][i]; } else { L[k+1][1][3][i]=list(); if (size(L[k+1][1][6][i])!=0) { if (size(L[k+1][1][8][i])!=0) { 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,nrows(L[k+1][1][5][i-1])); } else { L[k+1][1][7][i]=L[k+1][1][6][i]; L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1])); L[k+1][1][4][i]=list(); } } else { if (size(L[k+1][1][8][i])!=0) { L[k+1][1][7][i]=L[k+1][1][8][i]; L[k+1][1][2][i]=list(); L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1])); } else { L[k+1][1][7][i]=list(); L[k+1][1][2][i]=list(); L[k+1][1][4][i]=list(); } } L[k+1][2][1][i]=L[k+1][1][3][i]; L[k+1][2][6][i]=L[k+1][1][7][i]; } } } i=size(L)+d+1; 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; if (size(L[k+1][1][6][i])!=0) { v=L[k+1][1][6][i],L[k+1][1][8][i]; L[k+1][1][7][i]=v; } else { L[k+1][1][7][i]=L[k+1][1][8][i]; } } else { L[k+1][1][8][i]=list(); L[k+1][1][7][i]=L[k+1][1][8][i]; } L[k+1][2][6][i]=L[k+1][1][7][i]; } for (k=1; k<=(size(L)+d); k++) { L[size(L)][2][5][k]=list(); L[size(L)][2][4][k]=list(); L[size(L)][2][8][k]=list(); L[size(L)][2][3][k]=L[size(L)][2][1][k]; L[size(L)][2][7][k]=L[size(L)][2][6][k]; } L[size(L)][2][7][size(L)+d+1]=L[size(L)][2][6][size(L)+d+1]; L[size(L)][2][8][size(L)+d+1]=list(); /* building the resolution of the last short exact piece*/ for (i=2; i<=d+size(L); i++) { v=0; if(size(L[size(L)][2][1][i-1])!=0) { L[size(L)][2][2][i]=unitmat(nrows(L[size(L)][2][1][i-1])); } else { L[size(L)][2][2][i-1]=list(); } } return(L); } /*case Syzstring=="Vdres"*/ list forVd; for (k=1; k<=(size(L)+d); k++)//????? { /* we compute a V_d-strict resolution for the first short exact piece*/ 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(); } /* we compute step by step V_d-strict resolutions over coker(L[i][2][1][1])->coker(L[i][2][3][1])->coker(L[i][2][1][5]) and coker(L[i+1][1][1][1])->coker(L[i+1][1][3][1])->coker(L[i+1][1][1][5]) using the already computed resolutions for coker(L[i][2][1][1])= coker(L[i][1][3][1]) and coker(L[i+1][1][1][1])=coker(L[i][2][5][1])*/ for (i=1; i=0) for the seqeunce coker(L[1][2][3][1])->...->coker(L[size(L)][2][3][1])*/ list out; int i,j,k,l,oldj,newj,nr,nc; for (i=1; i<=size(L); i++) { out[i]=list(list()); out[i][1][1]=ncols(L[i][2][3][1]);//rank of module P^i_0 if (size(L[i][2][5][1])!=0) { /*horizontal differential P^i_0->P^(i+1)_0*/ nc=ncols(L[i][2][5][1]); out[i][1][4]=prodr(ncols(L[i][2][3][1])-ncols(L[i][2][5][1]),nc); } else { /*horizontal differential P^i_0->0*/ out[i][1][4]=matrix(0,ncols(L[i][2][3][1]),1); } oldj=newj; for (j=1; j<=size(L[i][2][3]);j++) { out[i][j][2]=L[i][2][7][j];//shift vector of P^i_{j-1} if (size(L[i][2][3][j])==0) { newj =j; break; } out[i][j+1]=list(); out[i][j+1][1]=nrows(L[i][2][3][j]);//rank of the module P^i_j out[i][j+1][3]=L[i][2][3][j];//vertical differential P^i_j->P^(i+1)_j if (size(L[i][2][5][j])!=0) { //horizonal differential P^i_j->P^(i-1)_j nr=nrows(L[i][2][3][j])-nrows(L[i][2][5][j]); out[i][j+1][4]=(-1)^j*prodr(nr,nrows(L[i][2][5][j])); } else { /*horizontal differential P^i_j->P^(i-1)_j*/ out[i][j+1][4]=matrix(0,nrows(L[i][2][3][j]),1); } if(j==size(L[i][2][3])) { out[i][j+1][2]=L[i][2][7][j+1];//shift vector of P^i_j newj=j+1; } } if (i>1) { for (k=1; k<=Min(list(oldj,newj)); k++) { /*horizonal differential P^(i-1)_(k-1)->P^i_(k-1)*/ nr=nrows(out[i-1][k][4]); out[i-1][k][4]=matrix(out[i-1][k][4],nr,out[i][k][1]); } for (k=newj+1; k<=oldj; k++) { /*no differential needed*/ out[i-1][k]=delete(out[i-1][k],4); } } } return (out); } //////////////////////////////////////////////////////////////////////////////////// static proc totalComplex(list L); { /* Input is the output of assemblingDoubleComplexes. We obtain a complex C^1[m^1]->...->C^(r)[m^r] with differentials d^i and shift vectors m^i (where C^r is placed in degree size(L)-1). This complex is dercribed in the list out as follows: rank(C^i)=out[3*i-2]; m_i=out[3*i-1] and d^i=out[3*i]*/ 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+1; 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]; 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); } //////////////////////////////////////////////////////////////////////////////////// //COMPUTATION OF THE BLOBAL B-FUNCTION //////////////////////////////////////////////////////////////////////////////////// static proc globalBFun(list L,list #) { /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]), where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an intvec of size n_i. We compute bounds for the minimal and maximal integer roots of the b-functions of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def. 6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal intersection (cf. Remark 6.1.7 in [R] 2012). This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this proc is only called from the procedure deRhamCohomology and the input comes originally from the procedure toVdstrictFreeComplex*/ if (size(#)==0)//# may contain the Syzstring { string Syzstring="Sres"; } else { string Syzstring=#[1]; } int i,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; ideal I; poly f; intvec i1; for (j=1; j<=size(L); j++) { /*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1] is already a V_d-strict Groebner basis, because it was obtained by the procedure toVdStrictFreeComplex*/ if (L[j][2]!=intvec(0:size(L[j][2]))) { if (Syzstring=="Vdres") { L[j][1]=VdStrictGB(L[j][1],n); } else { def HomWeyl=makeHomogenizedWeyl(n); setring HomWeyl; list L=fetch(W,L); L[j][1]=nHomogenize(L[j][1]); L[j][1]=transpose(matrix(slimgb(transpose(L[j][1])))); L[j][1]=subst(L[j][1],h,1); setring W; L=fetch(HomWeyl,L); kill HomWeyl; } } for (i=1; i<=ncols(L[j][1]); i++) { G0[j][i]=I; } for (i=1; i<=nrows(L[j][1]); i++) { /*computes the terms of maximal V_d-degree of the biggest non-zero component of submat(L[j][1],i,(1..ncols(L[j][1])))*/ i1=(1..ncols(L[j][1])); out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1); f=L[j][1][i,out[2]]; G0[j][out[2]]=G0[j][out[2]],f; G0[j][out[2]]=compress(G0[j][out[2]]); } } list save; int l; list weights; /*bFctIdealModified computes the intersection of G0[j][i] and x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/ for (j=1; j<=size(G0); j++) { for (i=1; i<=size(G0[j]); i++) { G0[j][i]=bFctIdealModified(G0[j][i]); } for (i=1; i<=size(G0[j]); i++) { weights=list(); if (size(G0[j][i])!=0) { for (l=i; l<=size(G0[j]); l++) { if (size(G0[j][l])!=0) { weights[size(weights)+1]=L[j][2][l]; } } G0[j][i]=list(G0[j][i][1]+Min(weights),G0[j][i][2]+Max(weights)); } } } list allmin; list allmax; for (j=1; j<=size(G0); j++) { for (i=1; i<=size(G0[j]); i++) { if (size(G0[j][i])!=0) { allmin[size(allmin)+1]=G0[j][i][1]; allmax[size(allmax)+1]=G0[j][i][2]; } } } list minmax=list(Min(allmin),Max(allmax)); return(minmax); } //////////////////////////////////////////////////////////////////////////////////// static proc exactGlobalBFun(list L,list #) { /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]), where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an intvec of size n_i. We compute bounds for the minimal and maximal integer roots of the b-functions of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def. 6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal intersection (cf. Remark 6.1.7 in [R] 2012). This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this proc is only called from the procedure deRhamCohomology and the input comes originally from the procedure toVdstrictFreeComplex*/ if (size(#)==0)//# may contain the Syzstring { string Syzstring="Sres"; } else { string Syzstring=#[1]; } int i,j,k; 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; matrix M; ideal I; poly f; intvec i1; for (j=1; j<=size(L); j++) { M=L[j][1]; /*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1] is already a V_d-strict Groebner basis, because it was obtained by the procedure toVdStrictFreeComplex*/ for (k=1; k<=ncols(L[j][1]); k++) { L[j][1]=permcol(M,1,k); if (Syzstring=="Vdres") { L[j][1]=VdStrictGB(L[j][1],n); } else { def HomWeyl=makeHomogenizedWeyl(n); setring HomWeyl; list L=fetch(W,L); L[j][1]=nHomogenize(L[j][1]); L[j][1]=transpose(matrix(slimgb(transpose(L[j][1])))); L[j][1]=subst(L[j][1],h,1); setring W; L=fetch(HomWeyl,L); kill HomWeyl; } for (i=1; i<=nrows(L[j][1]); i++) { /*computes the terms of maximal V_d-degree of the biggest non-zero component of submat(L[j][1],i,(1..ncols(L[j][1])))*/ i1=(1..ncols(L[j][1])); out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1); f=L[j][1][i,out[2]]; if (out[2]==1) { G0[j][k]=G0[j][k],f; G0[j][k]=compress(G0[j][k]); } } } } list save; int l; list weights; /*bFctIdealModified computes the intersection of G0[j][i] and x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/ for (j=1; j<=size(G0); j++) { for (i=1; i<=size(G0[j]); i++) { G0[j][i]=bFctIdealModified(G0[j][i]); } for (i=1; i<=size(G0[j]); i++) { if (size(G0[j][i])!=0) { G0[j][i]=list(G0[j][i][1]+L[j][2][i],G0[j][i][2]+L[j][2][i]); } } } list allmin; list allmax; for (j=1; j<=size(G0); j++) { for (i=1; i<=size(G0[j]); i++) { if (size(G0[j][i])!=0) { allmin[size(allmin)+1]=G0[j][i][1]; allmax[size(allmax)+1]=G0[j][i][2]; } } } list minmax=list(Min(allmin),Max(allmax)); return(minmax); } //////////////////////////////////////////////////////////////////////////////////// static proc bFctIdealModified (ideal I) {/*modified version of the procedure bfunIdeal from bfun.lib*/ def B= basering; int n = nvars(B) div 2; intvec w=(1:n); I= initialIdealW(I,-w,w); poly s; int i; for (i=1; i<=n; i++) { s=s+x(i)*D(i); } /*pIntersect computes the intersection on s and I*/ vector b = pIntersect(s,I); list RL = ringlist(B); RL = RL[1..4]; RL[2] = list(safeVarName("s")); RL[3] = list(list("dp",intvec(1)),list("C",intvec(0))); def @S = ring(RL); setring @S; vector b = imap(B,b); poly bs = vec2poly(b); ring r=0,s,dp; poly bs=imap(@S,bs); /*find minimal and maximal integer root*/ ideal allfac=factorize(bs,1); list allfacs; for (i=1; i<=ncols(allfac); i++) { allfacs[i]=allfac[i]; } number testzero; list zeros; for (i=1; i<=size(allfacs); i++) { if (deg(allfacs[i])==1) { testzero=number(subst(allfacs[i],s,0))/leadcoef(allfacs[i]); if (testzero-int(testzero)==0) { zeros[size(zeros)+1]=int(-1)*int(testzero); } } } if (size(zeros)!=0) { list minmax=(Min(zeros),Max(zeros)); } else { list minmax=list(); } setring B; return(minmax); } //////////////////////////////////////////////////////////////////////////////////// static proc safeVarName (string s) {/* from the library "bfun.lib"*/ string S = "," + charstr(basering) + "," + varstr(basering) + ","; s = "," + s + ","; while (find(S,s) <> 0) { s[1] = "@"; s = "," + s; } s = s[2..size(s)-1]; return(s) } //////////////////////////////////////////////////////////////////////////////////// static proc globalBFunOT(list L,list #) { /*this proc is currently not used since globalBFun computes the same output and is faster, however globalBFun works only for non-zero b-functions!*/ /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]), where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an intvec of size n_i. We compute bounds for the minimal and maximal integer roots of the b-functions of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def. 6.1.1 in [R]) using Algorithm 6.1.6 in [R].*/ if (size(#)==0) { string Syzstring="Sres"; } else { string Syzstring=#[1]; } int i; int j; def W=basering; int n=nvars(W) div 2; list G0; ideal I; intvec i1; 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++) { if (L[j][2]!=intvec(0:size(L[j][2]))) { if (Syzstring=="Vdres") { L[j][1]=VdStrictGB(L[j][1],n); } else { def HomWeyl=makeHomogenizedWeyl(n); setring HomWeyl; list L=fetch(W,L); L[j][1]=nHomogenize(L[j][1]); L[j][1]=transpose(matrix(slimgb(transpose(L[j][1])))); L[j][1]=subst(L[j][1],h,1); setring W; L=fetch(HomWeyl,L); kill HomWeyl; } } for (i=1; i<=nrows(L[j][1]); i++) { i1=(1..ncols(L[j][1])); out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(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; 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,G2,LL; intvec e,f; int kapp,k,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]=slimgb(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]; LL=list(lexp,leadcoef(lterm)); G3[l][i][size(G3[l][i])][size(G3[l][i][size(G3[l][i])])+1]=LL; } } } } } ring r=0,(s(1..n)),dp; ideal I; map G3forr=Wuv,I; list G3=G3forr(G3); poly fs,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); J=intersect(J,G4[l][i]); } G4[l]=poly(std(J)[1]); } list minmax; list mini=list(); list maxi=list(); list L=fetch(W,L); for (i=1; i<=size(G4); i++) { minmax[i]=minIntRoot0(G4[i],1); if (size(minmax[i])!=0) { mini=insert(mini,minmax[i][1]+Min(L[i][2])); maxi=insert(maxi,minmax[i][2]+Max(L[i][2])); } } mini=Min(mini); maxi=Max(maxi); minmax=list(mini[1],maxi[1]); option(none); return(minmax); } //////////////////////////////////////////////////////////////////////////////////// //COMPUTATION OF THE COHOMOLOGY //////////////////////////////////////////////////////////////////////////////////// static proc findCohomology(list L,int le) { /*computes the cohomology of the complex (D^i,d^i) given by D^i=C^L[2*i-1] and d^i=L[2*i]*/ def R=basering; ring r=0,(x),dp; list L=imap(R,L); list out; int i, ker, im; matrix S; option(returnSB); option(redSB); for (i=2; i<=size(L); i=i+2) { if (L[i-1]==0) { out[i div 2]=0; im=0; } else { S=matrix(syz(transpose(L[i]))); if (S!=matrix(0,nrows(S),ncols(S))) { ker=ncols(S); out[i div 2]=ker-im; im=L[i-1]-ker; } else { out[i-1]=0; im=L[i-1]; } } } option(none); while (size(out)>le) { out=delete(out,1); } setring R; return(out); } //////////////////////////////////////////////////////////////////////////////////// //AUXILIARY PROCEDURES //////////////////////////////////////////////////////////////////////////////////// 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) { intvec v=option(get); option(none); matrix l=transpose(lift(transpose(M),transpose(N))); option(set,v); return(l); } //////////////////////////////////////////////////////////////////////////////////// static proc VdStrictGB (matrix M,int d,list #) "USAGE:VdStrictGB(M,d[,v]); M a matrix, d an integer, v an optional intvec ASSUME:-basering is the nth Weyl algebra D_n @* -1<=d<=n @* -v (if given) is the shift vector on the range of M (in particular, size(v)=ncols(M)); otherwise v is assumed to be the zero shift vector RETURN:matrix N; the rows of N form a V_d-strict Groebner basis with respect to v for the module generated by the rows of M " { if (M==matrix(0,nrows(M),ncols(M))) { return (matrix(0,1,ncols(M))); } intvec op=option(get); 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]=list("dp",intvec(1)); matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); Data[5]=re; int k,l; Data[6]=transpose(concat(matrix(0,1,1),transpose(concat(matrix(0,1,1),Data[6])))); def Whom=ring(Data);// D_n[nhv] with the new commuative variable nhv 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);//homogenization of M with respect to the new variable Mnew=transpose(Mnew); Mnew=slimgb(Mnew);// computes a Groebner basis of the homogenzition of M Mnew=subst(Mnew,nhv,1);// substitution of 1 gives V_d-strict Groebner basis of M Mnew=compress(Mnew); Mnew=transpose(Mnew); setring W; M=imap(Whom,Mnew); option(set,op); return(M); } //////////////////////////////////////////////////////////////////////////////////// static proc VdNormalForm(matrix F,matrix M,int d,intvec v,list #) "USAGE:VdNormalForm(F,M,d,v[,w]); F and M matrices, d int, v intvec, w an optional intvec ASSUME:-basering is the nth Weyl algebra D_n @* -F a n_1 x n_2-matrix and M a m_1 x m_2-matrix with m_2<=n_2 @* -d is an integer between 1 and n @* -v is a shift vector for D_n^(m_2) and hence size(v)=m_2 @* -w is a shift vector for D_n^(m_1-m_2) and hence size(v)=m_1-m_2 @* RETURN:a n_1 x n_2-matrix N such that:@* -If no optional intvec w is given:(N[i,1],..,N[i,m_2]) is a V_d-strict normal form of (F[i,1],...,F[i,m_2]) with respect to a V_d-strict Groebner basis of the rows of M and the shift vector v -If w is given:(N[i,1],..,N[i,m_2]) is chosen such that Vddeg((N[i,1],...,N[i,m_2])[v])<=Vddeg((F[i,m_2+1],...,F[i,m_1])[v]); -N[i,j]=F[i,j] for j>m_2 " { int SBcom; def W =basering; int c=ncols(M); matrix keepF=F; if (size(#)!=0) { intvec w=#[1]; } 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]=list("dp",intvec(1)); matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); Data[5]=re; int k,l,nr,nc; 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);//new ring D_n[nvh] this new commuative variable nhv setring Whom; matrix Mnew=imap(W,M); list forMnew=homogenize(Mnew,d,v,1);//commputes homogenization of M; Mnew=forMnew[1]; int rightexp=forMnew[2]; matrix Fnew=imap(W,F); matrix keepF=imap(W,keepF); matrix Fb; int cc; intvec i1,i2; matrix zeromat,subm1,subm2,zeromat2; for (l=1; l<=nrows(Fnew); l++) { if (size(#)!=0) { subm2=submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF))); zeromat2=matrix(0,1,ncols(subm2)); if (submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF)))==zeromat2) { for (cc=1; cc<=ncols(Fnew); c++) { Fnew[l,cc]=0; } } i1=intvec(1..ncols(Fnew)); subm1=submat(Fnew,l,i1); subm2=submat(keepF,l,(ncols(Fnew)+1)..ncols(keepF)); zeromat=matrix(0,1,ncols(Fnew)); if (VdDegnhv(subm1,d,v)>VdDegnhv(subm2,d,w) and submat(Fnew,l,intvec(1..ncols(Fnew)))!=zeromat) { //print("Reduzierung des V_d-Grades nötig"); /*We need to reduce the V_d-degree. First we homogenize the lth row of Fnew*/ Fb=homogenize(subm1,d,v)*(nhv^rightexp); if (SBcom==0) { /*computes a V_d-strict standard basis*/ Mnew=slimgb(transpose(Mnew));// SBcom=1; } /*computes a V_d-strict normal form for FB*/ Fb=transpose(reduce(transpose(Fb),groebner(Mnew))); if (VdDegnhv(Fb,d,v)> VdDegnhv(subm2,d,w) and Fb!=matrix(0,nrows(Fb),ncols(Fb)))//should not happen { //print("Reduzierung fehlgeschlagen!!!!!!!!!!!!!!!!"); } } else { /*condition on V_ddeg already satisfied -> no normal form computation is needed*/ Fb=submat(Fnew,l,intvec(1..ncols(Fnew))); } } else { Fb=homogenize(submat(Fnew,l,intvec(1..ncols(Fnew))),d,v); if (SBcom==0) { Mnew=slimgb(transpose(Mnew));// computes a V_d-strict Groebner basis SBcom=1; } Fb=transpose(reduce(transpose(Fb),groebner(Mnew)));//normal form } for (k=1; k<=ncols(Fnew);k++) { Fnew[l,k]=Fb[1,k]; } } Fnew=subst(Fnew,nhv,1);//obtain normal form in D_n setring W; F=imap(Whom,Fnew); return(F); } //////////////////////////////////////////////////////////////////////////////////// static proc homogenize (matrix M,int d,intvec v,list #) { /* we compute the F[v]-homogenization of each row of M (cf. Def. 3.4 in [OT])*/ if (M==matrix(0,nrows(M),ncols(M))) { return(M); } int i,l,s, kmin, nhvexp; poly f; intvec vnm; list findmin,maxnhv,rempoly,remk,rem1,rem2; int n=(nvars(basering)-1) div 2; for (int k=1; k<=nrows(M); k++) { 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++) { nhvexp=remk[l][i]-kmin; M[k,l]=M[k,l]+nhv^(nhvexp)*rempoly[l][i]; maxnhv[size(maxnhv)+1]=nhvexp; } } } rempoly=list(); remk=list(); findmin=list(); } maxnhv=Max(maxnhv); nhvexp=maxnhv[1]; if (size(#)!=0) { return(list(M,nhvexp));//only needed for normal form computations } return(M); } //////////////////////////////////////////////////////////////////////////////////// static proc soldr (matrix M,matrix N) { /* We compute a ncols(M) x nrows(M)-matrix C such that C[i,1]M_1+...+C[i,nrows(M)]M_(nrows(M))= e_i mod im(N), where e_i is the ith basis element on the range of M, M_j denotes the jth row of M and im(N) is generated by the rows of 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,Smodnew; option(returnSB); int i,j; for (i=1;i<=q;i++) { E[i,1]=1; Smod2=concat(E,Smod); 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); option(none); return (Snew); } //////////////////////////////////////////////////////////////////////////////////// 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); } //////////////////////////////////////////////////////////////////////////////////// static proc VdDeg(matrix M,int d,intvec v,list #) { /* We assume that the basering it the nth Weyl algebra and that M is a 1 x r- matrix. We compute the V_d-deg of M with respect to the shift vector v, i.e V_ddeg(M)=max (V_ddeg(M_i)+v[i]), where k=V_ddeg(M_i) if k is the minimal integer, such that M_i can be expressed as a sum of operators x(1)^(a_1)*...*x(n)^(a_n)*D(1)^(b_1)*...*D(n)^(b_n) with a_1+..+a_d+k>=b_1+..+b_d*/ int i, j, etoint; int n=nvars(basering) div 2; intvec e; list findmax; int c=ncols(M); poly l; list positionpoly,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) { /*returns the largest component that has maximal V_d-degree and its terms of maximal V_d-deg (needed for globalBFun)*/ return (list(l,i)); } } } //////////////////////////////////////////////////////////////////////////////////// static proc VdDegnhv(matrix M,int d,intvec v,list #) { /* As the procedure VdDeg, but the basering is the nth Weyl algebra with a commutative variable nhv*/ int i,j,etoint; 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[2..d+1]+e[n+2..n+d+1]; 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)); } } //////////////////////////////////////////////////////////////////////////////////// static proc deletecol(matrix M,int l) { int s=ncols(M); if (l==1) { M=submat(M,(1..nrows(M)),(2..ncols(M))); return(M); } if (l==s) { M=submat(M,(1..nrows(M)),(1..(ncols(M)-1))); return(M); } intvec v=(1..(l-1)),((l+1)..s); M=submat(M,(1..nrows(M)),v); return(M); } //////////////////////////////////////////////////////////////////////////////////// static proc mHom(poly f) {/*for globalBFunOT*/ 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)=2) { if (typeof(#[2])=="intvec") { intvec n=#[2]; for (i=1; i<=size(n); i++) { w[size(w)+1]=n[i]; } RL[3]=insert(RL[3],list("am",w)); } else { RL[3]=insert(RL[3],list("a",w)); } } else { RL[3]=insert(RL[3],list("a",w)); } /*this ordering is needed for globalBFun and globalBFunOT*/ list saveord=RL[3][3]; RL[3][3]=RL[3][4]; RL[3][4]=saveord; intvec notforh=(1:(size(RL[3][4][2])-1)); RL[3][4][2]=notforh; RL[3][5]=list("dp",1); def @@R=ring(RL); return(@@R); } //////////////////////////////////////////////////////////////////////////////////// static proc nHomogenize (matrix M,list #) { /* # may contain an intvec v, if no intvec is given, we assume that v=(0:ncols(M)) We compute the h[v]-homogenization of the rows of M as in Definition 9.2 [OT]*/ int l; poly f; int s; int i; intvec vnm;int kmin; list findmax; int n=(nvars(basering)-1) div 2; list rempoly; list remk; list rem1; list rem2; list maxhexp; int hexp; intvec v=(0:ncols(M)); if (size(#)!=0) { if (typeof(#[1])=="intvec") { v=#[1]; } } if (size(v)1) { list forreturn=M,hexp; return(forreturn); } return(M); } //////////////////////////////////////////////////////////////////////////////////// static proc max(int i,int j) { if(i>j){return(i);} return(j); } //////////////////////////////////////////////////////////////////////////////////// static proc nDeg (matrix M,intvec m) {/*we compute an intvec n such that n[i]=max(deg(M[i,j])+m[j]|M[i,j]!=0) (where deg stands for the total degree) if (M[i,j]!=0 for some j) and n[i]=0 else*/ int i; int j; intvec n; list L; for (i=1; i<=nrows(M); i++) { L=list(); for (j=1; j<=ncols(M); j++) { if (M[i,j]!=0) { L=insert(L,deg(M[i,j])+m[j]); } } if (size(L)==0) { n[i]=0; } else { n[i]=Max(L); } } return(n); } //////////////////////////////////////////////////////////////////////////////////// static proc minIntRoot0(list L,list #) "USAGE:minIntRoot0(L [,M]); L list, M optinonal list ASSUME:L a list of univariate polynomials with rational coefficients @* the variable of the polynomial is s if size(#)==0 (needed for proc MVComplex) and t else (needed for globalBFun) RETURN:-if size(#)==0: int i, where i is an integer root of one of the polynomials and it is minimal with respect to that property@* -if size(#)!=0: list L=(i,j), where i is as above and j is an integer root of one of the polynomials and is maximal with respect to that property (if an integer root exists) or L=list() else " { def B=basering; if (size(#)==0) { ring rnew=0,s,dp; } else { ring rnew=0,t,dp; } list L=imap(B,L); int i; int j; number isint; list possmin; ideal allfac; list allfacs; for (i=1; i<=size(L); i++) { allfac=factorize(L[i],1); for (j=1; j<=ncols(allfac); j++) { allfacs[j]=allfac[j]; } for (j=1; j<=size(allfacs); j++) { if (deg(allfacs[j])==1) { isint=number(subst(allfacs[j],var(1),0)/leadcoef(allfacs[j])); if (isint-int(isint)==0) { possmin[size(possmin)+1]=int(isint); } } } allfacs=list(); } int zerolist; if (size(possmin)!=0) { int miniroot=(-1)*Max(possmin); int maxiroot=(-1)*Min(possmin); } else { zerolist=1; } setring B; if (size(#)==0) { return(miniroot); } else { if (zerolist==0) { return(list(miniroot,maxiroot)); } else { return(list()); } } } //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////// /* //////////////////////////////////////////////////////////////////////////////////// FURTHER EXAMPLES FOR TESTING THE PROCEDURES //////////////////////////////////////////////////////////////////////////////////// LIB "derham.lib"; //---------------------------------------- //EXAMPLE 1 //---------------------------------------- ring r=0,(x,y),dp; poly f=y2-x3-2x+3; list L=deRhamCohomology(f); L; kill r; //---------------------------------------- //EXAMPLE 2 //---------------------------------------- ring r=0,(x,y),dp; poly f=y2-x3-x; list L=deRhamCohomology(f); L; kill r; //---------------------------------------- //EXAMPLE 3 //---------------------------------------- ring r=0,(x,y),dp; list C=list(x2-1,(x+1)*y,y*(y2+2x+1)); list L=deRhamCohomology(C); L; kill r; //---------------------------------------- //EXAMPLE 4 //---------------------------------------- ring r=0,(x,y,z),dp; list C=list(x*(x-1),y,z*(z-1),z*(x-1)); list L=deRhamCohomology(C); L; kill r; //---------------------------------------- //EXAMPLE 5 //---------------------------------------- ring r=0,(x,y,z),dp; list C=list(x*y,y*z); list L=deRhamCohomology(C,"Vdres"); L; kill r; //---------------------------------------- //EXAMPLE 6 //---------------------------------------- ring r=0,(x,y,z,u),dp; list C=list(x,y,z,u); list L=deRhamCohomology(C); L; kill r; //---------------------------------------- //EXAMPLE 7 //---------------------------------------- ring r=0,(x,y,z),dp; poly f=x3+y3+z3; list L=deRhamCohomology(f); L; kill r; //---------------------------------------- //EXAMPLE 8 //---------------------------------------- ring r=0,(x,y,z),dp; poly f=x2+y2+z2; list L=deRhamCohomology(f,"Vdres"); L; kill r; //---------------------------------------- //EXAMPLE 9 //---------------------------------------- ring r=0,(x,y,z,u),dp; list C=list(x2+y2+z2,u); list L=deRhamCohomology(C); L; kill r; //---------------------------------------- //EXAMPLE 10 //---------------------------------------- ring r=0,(x,y,z),dp; list C=list((x*(y-1),y2-1)); list L=deRhamCohomology(C); L; kill r; */