-- Available routines are : -- isMonomial -- emonomials -- randomNonZero -- circuitsSupports -- circuitsIdeal -- orlikSolomon -- SP/2002/12/21 -- save this as "circuits.m2". Use in M2 with load "circuits.m2". 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) emonomials = method(Options => { Strategy => null} ) emonomials(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 { emonomials, TT "emonomials(I)", " -- computes the ideal generated by all monomials in the given ideal in the exterior algebra"} document { emonomials => Strategy, TT "which algorithm to use", TT "An option for ", TO "emonomials", " 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 emonomials(I) time emonomials(I, Strategy=>GBLeadTerms) time emonomials(I, Strategy=>TorusDeformation) time emonomials(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 emonomials(I) -- used 0.62 seconds time emonomials(I, Strategy=>GBLeadTerms) -- used 0.61 seconds time emonomials(I, Strategy=>TorusDeformation) -- used 0.71 seconds time emonomials(I, Strategy=>TorusDeformation2) -- used 0.64 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 emonomials(I) time emonomials(I, Strategy=>GBLeadTerms) time emonomials(I, Strategy=>TorusDeformation) time emonomials(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 emonomials(I) -- used 1.21 seconds time emonomials(I, Strategy=>GBLeadTerms) -- used 1.2 seconds time emonomials(I, Strategy=>TorusDeformation) -- used 15.56 seconds time emonomials(I, Strategy=>TorusDeformation2) -- takes 21.14 seconds -- n=8 m=3 R=ZZ/101[x_0..x_(n),SkewCommutative=> true] I= ideal random(R^1, R^{m:-1}) time emonomials(I) -- used 4.48 seconds time emonomials(I, Strategy=>GBLeadTerms) -- used 4.68 seconds time emonomials(I, Strategy=>TorusDeformation) -- used 47.42 seconds time emonomials(I, Strategy=>TorusDeformation2) -- used 31.78 seconds -- so only for generic linear forms "TorusDeformation" wins ! /// support = (m) -> ( sort unique join toSequence (keys \ keys standardForm (m)) ) circuitsSupports = method(Options => {Strategy => null}) circuitsSupports (Matrix) := List => options -> (M) -> ( Gale := syz M; X := symbol X; E := ZZ/31991[X_1..X_(rank source M), SkewCommutative=>true]; monsI := emonomials(ideal((vars E)*Gale), options); sups:= apply(flatten entries gens monsI, i-> support(i)); sups ) circuitsIdeal = method(Options => {Strategy => null}) circuitsIdeal (Matrix, Ring) := Ideal => options -> (M,R) -> ( Gale := syz M; X := symbol X; E := ZZ/31991[X_1..X_(rank source M), SkewCommutative=>true]; monsI := emonomials(ideal((vars E)*Gale), options); move := map (R, E, vars R); move(monsI)) TEST /// A=d->matrix{{1,1,1,1,1},{0,1,1,0,0},{0,0,1,1,d}} time C=circuitsSupports(A(30)) R=ZZ/101[x_0..x_4] time I=circuitsIdeal(A(2),R) -- M= transpose matrix{{1, -1, 0, 0}, {-1, 0, 1, 0}, {1, 0, 0, -1}, {0, -1, 0, 1}, {0, 0, 1, -1}} R=ZZ/101[x_0..x_4] time I=circuitsIdeal(M,R) time I=circuitsIdeal(M,R, Strategy=>GBLeadTerms) -- -- a 9 x 12 example M= matrix {{0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1}, {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1}, {0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 1}, {0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 1, 1}, {0, 0, 1, 0, 0, 0, 0, -2, 0, -1, 1, 2}, {0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0}} R=ZZ/101[x_0..x_11] time I=circuitsIdeal(M,R) -- used 50.15 seconds -- -- a 16 x 20 example M= matrix {{0, -1, 0, 0, 0, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0}, {1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, {-1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0}, {0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0}, {1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 1, 1, -1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0}, {0, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 1, 0, 0, -1, 0, 0, -1, 0, 0, -1, -1, 0, 0, 0, -1, 0, 0, 0}, {0, 0, 0, 1, -1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0}, {1, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1}, {1, 1, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, -1}, {-1, -1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, {0, 0, 0, -2, 0, -1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0}} R=ZZ/101[x_0..x_19] time I=circuitsIdeal(M,R) /// 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:=emonomials(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, TT "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=4 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)) /// 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) /// 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) /// 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) ///