-- Routines for computing Orlik-Solomon algebras and related ... -- Version of 12/15/00 -- David Eisenbud & Sorin Popescu -- Available routines are : -- isMonomial -- monomials -- randomNonZero -- orlikSolomon -- orlikSolomonDiff -- symExt -- FA -- -- Most of the speedup comes from the idea that to -- compute circuits it is equivalent to finding the ideal -- generated by the monomials in a given ideal is the -- exterior algebra. isMonomial = method() isMonomial Ideal := (I) -> ( if not isPolynomialRing(ring I) then error "expected an ideal in a polynomial ring or an exterior algebra"; M:=mingens I; rank source compress(M-leadTerm(M)) == 0) isMonomial Module := (I) -> ( if not isPolynomialRing(ring I) then error "expected an ideal in a polynomial ring or an exterior algebra"; M:=mingens I; rank source compress(M-leadTerm(M)) == 0) document { isMonomial, TT "isMonomial", " -- checks wether something is a monomial ideal or module"} TEST /// R=ZZ/101[x_0 .. x_4] I=ideal(x_0*x_2, x_0^7*x_3-x_3^(10), x_1^5*x_4, x_1*x_3, x_2^2*x_4) isMonomial(I) -- J=ideal(x_0*x_2, x_0^7*x_3, x_1^5*x_4, x_1*x_3, x_2^2*x_4) isMonomial(J) /// randomNonZero = method() randomNonZero Ring := (R) -> ( x:=0_R; while (x==0_R) do ( x = random(R)); x) monomials1 = (I) -> ( R:=ring I; K := I; J := ideal(1_R); while (not isMonomial(K)) do ( J = ideal leadTerm gens gb K; K = intersect(I,J) ); ideal mingens K) monomials2 = (I) -> ( R:=ring I; K := I; while (not isMonomial(K)) do ( fran := map(R,R,apply(flatten entries vars R, i-> randomNonZero(coefficientRing R)*i)); K = intersect(K,fran(I)); ); ideal mingens K) monomials3= (I) -> ( R:=ring I; K := I; while (not isMonomial(K)) do ( fran1 := map(R,R,apply(flatten entries vars R, i-> randomNonZero(coefficientRing R)*i)); fran2 := map(R,R,apply(flatten entries vars R, i-> randomNonZero(coefficientRing R)*i)); K = intersect(K,fran1(I), fran2(I)) ); ideal mingens K) monomials4 = (I) -> ( K := I; R := ring I; J := ideal 1_R; while not isSubset(J,K) do ( J = ideal leadTerm gens gb K; K = intersect(I,J) ); K) monomials = method(Options => { Strategy => null} ) monomials(Ideal) := Ideal => options -> (I) -> ( if isCommutative(ring I) then error "expected an ideal in an exterior algebra"; strategy := options.Strategy; if strategy === null then strategy = GBLeadTerms; local f; if strategy === GBLeadTerms then ( f = monomials1; ) else if strategy === TorusDeformation then ( f = monomials2; ) else if strategy === TorusDeformation2 then ( f = monomials3; ); f(I) ) document { monomials, TT "monomials(I)", " -- computes the ideal generated by all monomials in the given ideal in the exterior algebra"} -- document { monomials => Strategy, -- HEADLINE "which algorithm to use", -- TT "An option for ", TO "monomials", " which specifies -- which strategy to use. Strategys are specified by the -- following tokens", BR, -- MENU{ -- (TT "Strategy => GBLeadTerms", "-- Iterates the intersection between -- the ideal and its initial ideal with respect to the current -- monomial order in the exterior algebra. This is the default -- since it appears to be on the average the fastest."), -- (TT "Strategy => TorusDeformation", "--Iterates the intersection -- between the ideal and a random torus deformation of it. This -- avoids computing GBs and the improvement can be very dramatic -- for very dense ideals."), -- (TT "Strategy => TorusDeformation2", "--Iterates the intersection -- between the ideal and two random torus deformation of it. This -- avoids also computing GBs and can sometimes be marginally -- faster then TorusDeformation.")} -- } TEST /// R=ZZ/101[a,b,c,d, SkewCommutative=> true] I=ideal(a*b,b*c,a*c+b*d) time monomials(I) time monomials(I, Strategy=>GBLeadTerms) time monomials(I, Strategy=>TorusDeformation) time monomials(I, Strategy=>TorusDeformation2) -- R=ZZ/101[a,b,c,d,e,f,g,h,k,l, SkewCommutative=> true] I=ideal(a*c+b*d+e*g, a*d+b*f+f*h, a*e+b*g+g*k, a*f+b*h+h*l, a*g+b*k, a*h+b*l) time monomials(I) time monomials(I, Strategy=>GBLeadTerms) -- used 1.35 seconds time monomials(I, Strategy=>TorusDeformation) -- used 1.55 seconds time monomials(I, Strategy=>TorusDeformation2) -- used 1.18 seconds --time emonomials(I,2) -- used 0.91 seconds --time emonomials(I,1) -- used 0.77 seconds -- R=ZZ/101[a,b,c,d,e,f,g, SkewCommutative=> true] I=ideal(a*b,b*c,a*c+b*d, a*d+b*f, a*e+b*g) time monomials(I) time monomials(I, Strategy=>GBLeadTerms) time monomials(I, Strategy=>TorusDeformation) time monomials(I, Strategy=>TorusDeformation2) -- R=ZZ/101[a,b,c,d,e,f,g,h,k,l, SkewCommutative=> true] I=ideal(a*c+b*d+e*g, a*d+b*f+f*h)+ideal random(R^1, R^{1:-1}) time monomials(I) -- used 2.36 seconds time monomials(I, Strategy=>GBLeadTerms) -- used 2.36 seconds time monomials(I, Strategy=>TorusDeformation) -- used 29.91 seconds time monomials(I, Strategy=>TorusDeformation2) -- n=8 m=3 R=ZZ/101[x_0..x_(n),SkewCommutative=> true] I= ideal random(R^1, R^{m:-1}) time monomials(I) -- used 8.57 seconds time monomials(I, Strategy=>GBLeadTerms) -- used 8.57 seconds time monomials(I, Strategy=>TorusDeformation) -- used 90.47 seconds time monomials(I, Strategy=>TorusDeformation2) -- used 54.99 seconds --time monomials4(I) -- used 7.77 seconds --time emonomials(I,2) -- used 62.18 seconds -- so only for generic linear forms "TorusDeformation" wins ! /// orlikSolomonDiff = (m) -> ( R:=ring m; d:=(degree m)_0; (compress diff(vars R, m))* (transpose map(R^1,R^d,{{d:1}})) ) document { orlikSolomonDiff, TT "orlikSolomonDiff(m)", " -- given a monomial m in the exterior algebra returns its Orlik-Solomon Differential"} TEST /// R=ZZ/101[a,b,c,d, SkewCommutative=> true] orlikSolomonDiff(a*b*c) /// orlikSolomon = method(Options => { Strategy => null} ) orlikSolomon(Matrix, Ring) := Ideal => options -> (M, E) -> ( if (not(isPolynomialRing(E)) or isCommutative(E)) then error "expected an exterior algebra"; if ( numgens E != rank source M) then error "the given exterior algebra doesn't have the expected number of variables"; N:=substitute(syz jacobian(M), E); mI:=monomials(ideal((vars E)*N), options); trim ideal matrix apply(flatten entries gens mI, i-> flatten entries(orlikSolomonDiff(i))) ) orlikSolomon(Matrix) := Ideal => options -> (M) -> ( ncols:= rank source M; x:=symbol x; E:= coefficientRing(ring M)[x_0..x_(ncols-1), SkewCommutative=> true]; orlikSolomon(M, E, options) ) document { orlikSolomon, TT "orlikSolomon", " -- computes the Orlik-Solomon ideal of the complement of a complex arrangement of hyperplanes", PARA, EXAMPLE { "n=4", "R=ZZ/101[x_1..x_n]", "M = gens ideal toList flatten apply(1..(n-1), i-> flatten toList apply(i+1..n, j-> {x_i-x_j}))", "--these are the hyps of the braid arrangement", "E=ZZ/101[y_1..y_(binomial(n,2)), SkewCommutative=>true]", "time (j=orlikSolomon(M,E))", "time (j=orlikSolomon(M,E, Strategy=>GBLeadTerms))", "time (j=orlikSolomon(M,E, Strategy=>TorusDeformation))", "time (j=orlikSolomon(M,E, Strategy=>TorusDeformation2))", } } -- document { orlikSolomon => Strategy, -- HEADLINE "specify which strategy to use when computing the -- ideal generated by all monomials in a given ideal in -- the exterior algebra", -- TT "orlikSolomon(M,S, Strategy=>v)", " -- The value of the option -- is simply passed to ", TO "monomials", ", so see the -- documentation there for details."} TEST/// n=4 m=6 R=ZZ/101[x_0..x_(n)] M= random(R^1, R^{m:-1}) time (j=orlikSolomon(M);) time (j=orlikSolomon(M, Strategy=>GBLeadTerms)) time (j=orlikSolomon(M, Strategy=>TorusDeformation)) time (j=orlikSolomon(M, Strategy=>TorusDeformation2)) -- n=5 R=ZZ/101[x_1..x_n] M=ideal(0_R) scan(1..(n-1), i-> scan(i+1..n, j-> M=M + ideal(x_i-x_j))) M=compress gens M --these are the hyps of the braid arrangement time (j=orlikSolomon(M)) time (j=orlikSolomon(M, Strategy=>GBLeadTerms)) time (j=orlikSolomon(M, Strategy=>TorusDeformation)) time (j=orlikSolomon(M, Strategy=>TorusDeformation2)) /// ---- the following needs to be tuned up symExt= (m,R) ->( --this function converts between the a linear presentation map --of a module over a symmetric algebra and the first map --in the corresponding complex of modules over the exterior --algebra, and vice-versa. The matrix m is the given --** linear ** -- presentation map, and R is the target ring. --the script returns a linear map over R. If --coker m has linear resolution (regularity 0) --then the resolution of coker p will be linear, --and will be the complex corresponding to coker m. if (not(isPolynomialRing(R))) then error "expected a polynomial ring or an exterior algebra"; if (numgens R != numgens ring m) then error "the given ring has a wrong number of variables"; ev := map(R,ring m,vars R); mt := transpose jacobian m; jn := gens kernel mt; q := vars(ring m)**id_(target m); n := ev(q*jn)) FA = (j, R) -> ( -- here j is assumed to be an Orlik-Solomon ideal; -- R is a polynomial ring with as many vars as the -- exterior algebra where j is contained. -- This is computing the module F(A) modT := (ring j)^1/(j*(ring j^1)); F:=res(prune modT, LengthLimit=>3); g:=transpose F.dd_2; G:=res(coker g,LengthLimit=>4); symExt(G.dd_4, R)) TEST/// --- Braid arrangement n = 4 R = ZZ/101[x_1..x_n] M = gens ideal toList flatten apply(1..(n-1), i-> flatten toList apply(i+1..n, j-> {x_i-x_j})) E=ZZ/101[y_1..y_(binomial(n,2)), SkewCommutative=>true] j=orlikSolomon(M,E) modT = (ring j)^1/(j*(ring j^1)) betti (F=res(modT, LengthLimit=>4)) g=transpose F.dd_2; betti (G=res (coker g,LengthLimit=>10)) S=ZZ/101[y_0..y_(binomial(n,2)-1)] Fa=FA(j,S) betti (FF=res coker Fa) betti FF.dd_2 fit = ideal mingens minors(5,FF.dd_2); fit =saturate(fit, ideal vars S) degree fit -- (for n=4 this is 5) codim fit llist=decompose fit --union of linear subspaces scan(llist, i-> print i) --The following is the output for n=4 --Comp1 (3,4,5) -- --ideal (y , y + y + y , y , y ) -- 0 3 4 5 2 1 --Comp2 (non-local component) --ideal (y - y , y + y + y , y + y + y , y - y ) -- 0 5 3 4 5 2 4 5 1 4 --Comp3 (0,2,4) -- --ideal y + y + y , y , y , y ) -- 0 2 4 3 1 5 --comp4 (1,2,5) -- --ideal (y , y , y + y + y , y ) -- 0 3 1 2 5 4 --Comp5 (0,1,3) -- --ideal (y + y + y , y , y , y ) -- 0 1 3 2 4 5 --Comp 1, 3, 4, 5 are the local comp of the char variety. Comp2 is --the only nonlocal one. The local ones correspond to the triple --points, while the non-local one corresponds to the partition --(05,14,23) (notation a la Falk) /// TEST/// -- deleted D3 arangement R=ZZ/101[x_0..x_2] M=gens ideal(x_0-x_2, x_1-x_2, x_0, x_1, x_0-x_1+x_2, x_2, x_0-x_1-x_2, x_0-x_1) j=orlikSolomon(M) modT = (ring j)^1/(j*(ring j^1)) betti (F=res(modT, LengthLimit=>3)) g=transpose F.dd_2; betti (G=res (coker g,LengthLimit=>6)) S=ZZ/101[y_0..y_7] Fa=FA(j,S) betti (FF=res coker Fa) betti FF.dd_2 fit = ideal mingens minors(7,FF.dd_2); --this takes for ever! a bug in minors ? degree fit codim fit decompose fit --union of linear subspaces scan(oo, i-> print i) /// TEST/// --- Reflection arrangement of type B3 R = ZZ/101[x,y,z] M = gens ideal(x,y,z,x-y,x-z,y-z,x-y-z,x-y+z) j=orlikSolomon(M) modT = (ring j)^1/(j*(ring j^1)) betti (F=res(modT, LengthLimit=>2)) g=transpose F.dd_2; betti (G=res (coker g,LengthLimit=>5)) S=ZZ/101[Y_0..Y_7] Fa=FA(j,S); betti (FF=res coker Fa) betti FF.dd_2 fit = ideal mingens minors(7,FF.dd_2); --this takes for ever! a bug in minors ? fit =saturate(fit, ideal vars S) degree fit -- codim fit llist=decompose fit --union of linear subspaces scan(llist, i-> print i) ///