--------------------------------------------------------------------------- -- PURPOSE : Computation of quotient, saturation, and annihilator -- -- UPDATE HISTORY : created 14 April 2018 at M2@UW; -- updated November 2020 -- -- ISSUE TRACKER : https://github.com/orgs/Macaulay2/projects/3 -- -- TODO : 1. why are there shadowed symbols? -- 2. any strategies to incorporate from Shimoyama-Yokoyama.m2 -- or greedySat in GTZ.m2 or getSaturation in newGTZ.m2 -- 3. which algorithms for Ideals can be adapted to submodules? -- 4. move radical here, perhaps intersect as well -- 5. employ storefun! Maybe call it declare --------------------------------------------------------------------------- newPackage( "Saturation", Version => "0.3", Date => "November 29, 2021", Headline => "quotient, saturation, and annihilator routines for ideals and modules", Authors => { {Name => "Justin Chen", Email => "justin.chen@math.gatech.edu"}, {Name => "Mahrud Sayrafi", Email => "mahrud@umn.edu", HomePage => "https://math.umn.edu/~mahrud"}, {Name => "Mike Stillman", Email => "mike@math.cornell.edu", HomePage => "http://www.math.cornell.edu/~mike"}}, Keywords => {"Commutative Algebra"}, PackageExports => { "Elimination" }, AuxiliaryFiles => true, DebuggingMode => false ) export { "isSupportedInZeroLocus" } exportFrom_Core { "saturate", "annihilator" } importFrom_Core { "nonnull", "raw", "rawColon", "rawSaturate", "newMonomialIdeal", "eliminationInfo", "resolutionDegreeLimit" } importFrom_Core { "isComputationDone", "cacheComputation", "fetchComputation", "updateComputation", "cacheHit", "Context", "Computation" } -- TODO: where should these be placed? trim MonomialIdeal := MonomialIdeal => opts -> (cacheValue (symbol trim => opts)) ((I) -> monomialIdeal trim(module I, opts)) -- TODO: pick a better name and move to interpreter for speed? -- this is similar to override, except: -- keys in opts with value null are ignored (TODO: is this an okay way to remove options?) -- keys in opts that don't exist in def don't produce an error override' := (def, opts) -> nonnull apply(keys def, key -> if opts#?key and opts#key =!= null then key => opts#key) -- TODO: is this the right name for this function? ambient' = method() ambient' Module := Module => ambient ambient' Ideal := Ideal => I -> (if instance(I, MonomialIdeal) then monomialIdeal else ideal) 1_(ring I) -- TODO: remove this once the Ideal vs. MonomialIdeal dichotomy is resolved uniform' := L -> all(L, l -> instance(l, Ideal)) or uniform L -- this is a light version of isSubset for modules and ideals, which doesn't compute a gb -- TODO: move these improvements to isSubset, for instance isSubset(0, A) should not compute gb of A -- TODO: add heuristics to decide when is it okay to compute a gb of B, e.g. 2 vars? isSubset' := (A, B) -> A == 0 or isSubset(entries transpose gens A, entries transpose gens B) --Ideal % Matrix := --remainder(Ideal, Matrix) := Matrix => (I, m) -> remainder(gens I, m) --Module % Matrix := --remainder(Module, Matrix) := Matrix => (M, m) -> remainder(gens M, m) -- This is a map from method keys to strategy hash tables algorithms := new MutableHashTable from {} -------------------------------------------------------------------- -- Helpers -------------------------------------------------------------------- -- given a ring R, determines if R is a poly ring over ZZ or a field isFlatPolynomialRing = R -> isPolynomialRing R and (isField(kk := coefficientRing R) or kk === ZZ) -- TODO: can this work with multigraded ideals? isGRevLexRing = (R) -> ( -- returns true if the monomial order in the polynomial ring R -- is graded reverse lexicographic order, w.r.t. the first degree -- vector in the ring. mo := (options monoid R).MonomialOrder; if #mo === 0 then return false; mo = select(mo, x -> x#0 =!= MonomialSize and x#0 =!= Position); isgrevlex := mo#0#0 === GRevLex and mo#0#1 === apply(degrees R, first); #mo === 1 and isgrevlex and all(mo, x -> x#0 =!= Weights and x#0 =!= Lex)) -- Helper for Linear strategies -- currently this is false for x+1, but in all use cases homogeneity is required anyway isLinearForm := f -> degreeLength ring f > 0 and sum degree f === 1 and same(degree \ terms f) -- Return (R1, R1<-R, R<-R1), where generators i and n are switched -- TODO: can this be simplified using newRing? grevLexRing = method() grevLexRing(ZZ, Ring) := (i, R) -> ( X := local X; n := numgens R; degs := degrees R; if i === n - 1 then return (R, identity, identity); perm := toList splice(0..i-1, n-1, i+1..n-2, i); R1 := (coefficientRing R)[X_1..X_n, Degrees => degs_perm, MonomialSize => 16]; fto := map(R1, R, (generators R1)_perm); fback := map(R, R1, (generators R)_perm); (R1, fto, fback)) -------------------------------------------------------------------- -- Quotients -------------------------------------------------------------------- -- quotient methods: -- 1. syzygies -- 2. use elimination methods? I forget how? -- 3. case: x is a variable, I is homogeneous -- case: x is a polynomial -- case: x is an ideal -- 4. DegreeLimit: stop at a certain degree -- 5. BasisElementLimit: stop after one element -- 6. PairLimit: stop after 100 S-pairs -- keys: the second object in the quotient QuotientContext = new SelfInitializingType of Context QuotientContext.synonym = "quotient context" QuotientComputation = new Type of Computation QuotientComputation.synonym = "quotient computation" -- TODO: try to find a compatible context that can be used in the computation new QuotientContext from Sequence := (C, S) -> ( (A, B) := S; QuotientContext{ mingens B } ) new QuotientComputation from Sequence := (C, S) -> new QuotientComputation from { BasisElementLimit => (S#2).BasisElementLimit, DegreeLimit => (S#2).DegreeLimit, PairLimit => (S#2).PairLimit, Result => null} isComputationDone QuotientComputation := Boolean => options quotient >> opts -> container -> ( -- this function determines whether we can use the cached result, or further computation is necessary not instance(container.Result, Nothing) and opts.BasisElementLimit <= container.BasisElementLimit and(opts.DegreeLimit <= container.DegreeLimit and opts.DegreeLimit =!= {} or container.DegreeLimit === {}) and opts.PairLimit <= container.PairLimit) updateComputation(QuotientComputation, Thing) := Thing => options quotient >> opts -> (container, result) -> ( container.BasisElementLimit = opts.BasisElementLimit; container.DegreeLimit = opts.DegreeLimit; container.PairLimit = opts.PairLimit; container.Result = result) --quotient = method(...) -- defined in m2/quotient.m2 quotient(Ideal, Ideal) := Ideal => opts -> (I, J) -> quotientHelper(I, J, (quotient, Ideal, Ideal), opts) quotient(Ideal, Number) := quotient(Ideal, RingElement) := Ideal => opts -> (I, f) -> quotient(I, ideal f, opts) Ideal : Ideal := Ideal => (I, J) -> quotient(I, J) Ideal : Number := Ideal : RingElement := Ideal => (I, f) -> quotient(I, f) -- TODO: why is this the right thing to do? quotient(MonomialIdeal, RingElement) := MonomialIdeal => opts -> (I, f) -> ( quotient(I, if size f === 1 and leadCoefficient f == 1 then monomialIdeal f else ideal f, opts)) MonomialIdeal : RingElement := MonomialIdeal => (I, f) -> quotient(I, f) quotient(Module, Ideal) := Module => opts -> (M, I) -> quotientHelper(M, I, (quotient, Module, Ideal), opts) quotient(Module, Number) := quotient(Module, RingElement) := Module => opts -> (M, f) -> quotient(M, ideal f, opts) Module : Ideal := Module => (M, I) -> quotient(M, I) Module : Number := Module : RingElement := Module => (M, f) -> quotient(M, f) -- annihilator of (B+A)/A, where A and B have a common ambient module -- note: if A is an ideal and B=f, then this is isomorphic to R/(A:f) quotient(Module, Module) := Ideal => opts -> (M, N) -> quotientHelper(M, N, (quotient, Module, Module), opts) Module : Module := Ideal => (M, N) -> quotient(M, N) -- TODO: the methods with "Number" should be made unnecessary via https://github.com/Macaulay2/M2/issues/1519 -- Helper for quotient methods quotientHelper = (A, B, key, opts) -> ( R := ring A; if ring B === ZZ then B = sub(B, R); strategy := opts.Strategy; doTrim := if opts.MinimalGenerators then trim else identity; cast := if instance(A, MonomialIdeal) then monomialIdeal else ideal; opts = opts ++ { DegreeLimit => resolutionDegreeLimit(R, opts.DegreeLimit) }; -- this logic runs the strategies in order, or the specified strategy computation := (opts, container) -> ( if R =!= ring B then error "expected objects in the same ring"; if instance(B, RingElement) then B = ideal B; if uniform' {A, B} and ambient' A != ambient' B then error "expected objects to be contained in the same ambient object"; -- note: if B \subset A then A:B should be "everything", but computing -- a gb for A can be slow, so isSubset' doesn't compute a gb if isSubset'(B, A) then return if uniform' {A, B} then cast 1_R else ambient' A; -- note: ideal(..A..) : f = A <==> f is nzd / A -- note: ideal(..A..) : ideal(..B..) = A <==> -- note: module(.A.) : ideal(..B..) = A <==> B is not contained in any associated primes of A -- TODO: can either of the above be efficiently checked? -- TODO: module(.A.) : module(.B.) = ? <==> A \subset B if instance(B, Ideal) then -- see the above TODO item if isSubset'(ambient' A, B) then return A; -- TODO: add speedup for when isSubset'(A, B), being cautious of (a3):(a3,bc) in kk[abc]/ac -- TODO: what would adding {SyzygyLimit => opts.BasisElementLimit, BasisElementLimit => null} do? opts = new OptionTable from override'(options gb, opts ++ {Strategy => null}); runHooks(key, (opts, A, B), Strategy => strategy)); -- this is the logic for caching partial quotient computations. A.cache contains an option: -- QuotientContext{ mingens B } => QuotientComputation{ Result } container := fetchComputation(QuotientComputation, A, (A, B, opts), new QuotientContext from (A, B)); -- the actual computation of quotient occurs here C := (cacheComputation(opts, container)) computation; if C =!= null then doTrim C else if strategy === null then error("no applicable strategy for ", toString key) else error("assumptions for quotient strategy ", toString strategy, " are not met")) -------------------------------------------------------------------- -- Algorithms for Ideal : Ideal algorithms#(quotient, Ideal, Ideal) = new MutableHashTable from { Iterate => (opts, I, J) -> ( R := ring I; -- TODO: extend this into a strategy for any PID if R === ZZ then return ideal sub(gcd I_* / gcd flatten(I_*, J_*), ZZ); fold(first entries mingens J, ideal 1_R, (f, M1) -> if generators(f * M1) % generators I == 0 then M1 else intersect(M1, quotient(I, ideal f, opts, Strategy => Quotient)))), -- FIXME: potential bug with nonstandard fields, see the last test in quotient-test.m2 Quotient => (opts, I, J) -> ( R := ring I; -- FIXME: this line computes a gb for I!!! mR := transpose generators J ** (R / I); -- if J is a single element, this is the same as -- computing syz gb(matrix{{f}} | generators I, ...) g := syz gb(mR, opts, Strategy => LongPolynomial, Syzygies => true, SyzygyRows => 1); -- The degrees of g are not correct, so we fix that here: -- g = map(R^1, null, g); lift(ideal g, R)), -- TODO Linear => (opts, I, J) -> ( -- assumptions: J is a single linear element, and everything is homogeneous if not isHomogeneous I or not isHomogeneous J or not isLinearForm J_0 then return null; stderr << "warning: quotient strategy Linear is not yet implemented" << endl; null), Monomial => (opts, I, J) -> ( R := ring I; if not isMonomialIdeal I or not isMonomialIdeal J or not isPolynomialRing R or not isCommutative R then return null; cast := if instance(I, MonomialIdeal) then identity else ideal; -- TODO: make sure (monomialIdeal, MonomialIdeal) isn't forgetful cast newMonomialIdeal(R, rawColon(raw monomialIdeal I, raw monomialIdeal J))), } -- Installing hooks for Ideal : Ideal scan({Quotient, Iterate-*, Linear*-, Monomial}, strategy -> addHook(key := (quotient, Ideal, Ideal), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- -- Algorithms for Module : Ideal algorithms#(quotient, Module, Ideal) = new MutableHashTable from { Iterate => (opts, M, J) -> ( -- This is the iterative version, where M is a -- submodule of F/K, or ideal, and J is an ideal. M1 := super M; m := generators M | relations M; scan(numgens J, i -> ( f := J_i; if generators(f*M1) % m != 0 then ( M2 := quotient(M, f, opts, Strategy => Quotient); M1 = intersect(M1, M2); ))); M1), Quotient => (opts, M, J) -> ( m := generators M; F := target m; if M.?relations then m = m | M.relations; j := transpose generators J; g := (j ** F) | (target j ** m); -- We would like to be able to inform the engine that -- it is not necessary to compute various of the pairs -- of the columns of the matrix g. h := syz gb(g, opts, Strategy => LongPolynomial, Syzygies => true, SyzygyRows => numgens F); if M.?relations then subquotient(h % M.relations, M.relations) else image h ), -- TODO: Linear => (opts, M, J) -> ( -- assumptions: J is a single linear element, and everything is homogeneous if not isHomogeneous M or not isHomogeneous J or not isLinearForm J_0 then return null; stderr << "warning: quotient strategy Linear is not yet implemented" << endl; null), } -- Installing hooks for Module : Ideal scan({Quotient, Iterate-*, Linear*-}, strategy -> addHook(key := (quotient, Module, Ideal), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- -- Algorithms for Module : Module algorithms#(quotient, Module, Module) = new MutableHashTable from { Iterate => (opts, I, J) -> ( R := ring I; M1 := ideal 1_R; m := generators I | relations I; scan(numgens J, i -> ( f := image (J_{i}); -- it used to say f ** M1, but that can't have been right. -- I'm just guessing that M1 * f is better. (drg) if generators(M1 * f) % m != 0 then ( M2 := quotient(I, f, opts, Strategy => Quotient); M1 = intersect(M1, M2); ))); M1), Quotient => (opts, M, J) -> ( m := generators M; if M.?relations then m = m | M.relations; j := adjoint(generators J, (ring J)^1, source generators J); F := target m; g := j | (dual source generators J ** m); -- << g << endl; -- We would like to be able to inform the engine that -- it is not necessary to compute various of the pairs -- of the columns of the matrix g. h := syz gb(g, opts, Strategy => LongPolynomial, Syzygies => true, SyzygyRows => 1); ideal h), } -- Installing hooks for Module : Module scan({Quotient, Iterate}, strategy -> addHook(key := (quotient, Module, Module), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- -- Saturations -------------------------------------------------------------------- -- TODO: -- - should saturate(I) use the irrelevant ideal when multigraded? -- keys: the second ideal of saturation SaturateContext = new SelfInitializingType of QuotientContext SaturateContext.synonym = "saturate context" SaturateComputation = new Type of QuotientComputation SaturateComputation.synonym = "saturate computation" -- TODO: try to find other compatible cacheKeys, perhaps of type QuotientContext, that can be used in the computation new SaturateContext from Sequence := (C, S) -> ( (A, B) := S; SaturateContext{ mingens B } ) new SaturateComputation from Sequence := (C, S) -> new SaturateComputation from { BasisElementLimit => (S#2).BasisElementLimit, DegreeLimit => (S#2).DegreeLimit, PairLimit => (S#2).PairLimit, Result => null} -- Note: isComputationDone and updateComputation inherit from QuotientComputation -- TODO: is there something smarter that can be done for saturation in particular? -- saturate = method(Options => options saturate) -- defined in m2/quotient.m2 saturate(Ideal, Ideal) := Ideal => opts -> (I, J) -> saturateHelper(I, J, (saturate, Ideal, Ideal), opts) saturate(Ideal, Number) := saturate(Ideal, RingElement) := Ideal => opts -> (I, f) -> saturateHelper(I, f, (saturate, Ideal, RingElement), opts) saturate Ideal := Ideal => opts -> I -> saturate(I, ideal vars ring I, opts) saturate(MonomialIdeal, RingElement) := MonomialIdeal => opts -> (I, f) -> ( saturate(I, if size f === 1 and leadCoefficient f == 1 then monomialIdeal f else ideal f, opts)) saturate(Module, Ideal) := Module => opts -> (M, J) -> saturateHelper(M, J, (saturate, Module, Ideal), opts) saturate(Module, Number) := saturate(Module, RingElement) := Module => opts -> (M, f) -> saturateHelper(M, f, (saturate, Module, RingElement), opts) saturate Module := Module => opts -> M -> saturate(M, ideal vars ring M, opts) -- TODO: is M / saturate 0_M a correct computation of saturation of M? saturate(Vector, Ideal) := Module => opts -> (v, J) -> saturate(image matrix {v}, J, opts) saturate(Vector, Number) := saturate(Vector, RingElement) := Module => opts -> (v, f) -> saturate(image matrix {v}, f, opts) saturate Vector := Module => opts -> v -> saturate(image matrix {v}, opts) -- used when P = decompose irr -- TODO: do a one-step caching here saturate(Ideal, List) := Ideal => opts -> (I, L) -> fold(L, I, (J, I) -> saturate(I, J, opts)) saturate(Module, List) := Module => opts -> (M, L) -> fold(L, M, (J, M) -> saturate(M, J, opts)) saturate(Vector, List) := Module => opts -> (v, L) -> saturate(image matrix {v}, L, opts) -- TODO: the methods with "Number" should be unnecessary via https://github.com/Macaulay2/M2/issues/1519 -- Helper for saturation methods saturateHelper = (A, B, key, opts) -> ( R := ring A; if ring B === ZZ then B = sub(B, R); -- this is only here because some methods are not implemented for RingElement B' := if instance(B, RingElement) or instance(B, Number) then ideal B else B; strategy := opts.Strategy; doTrim := if opts.MinimalGenerators then trim else identity; opts = opts ++ { DegreeLimit => resolutionDegreeLimit(R, opts.DegreeLimit) }; -- this logic runs the strategies in order, or the specified strategy computation := (opts, container) -> ( if R =!= ring B then error "expected objects in the same ring"; -- note: if B \subset A then A:B^infty should be "everything", but computing -- a gb for A can be slow, so isSubset' doesn't compute a gb -- TODO: if radical A is cached and B \subset radical A then A : B^infty = ambient A -- alternatively, can radical containment be efficiently checked? if isSubset'(B', A) then return ambient' A; -- note: ideal(..A..) : f^infty = A <==> f is nzd /A -- note: ideal(..A..) : ideal(..B..)^infty = A <==> B is not contained in any associated primes of A -- TODO: can either of the above be efficiently checked? if isSubset'(ambient' A, B') then return A; -- TODO: add speedup for when isSubset'(A, B'), being cautious of (a3):(a3,bc) in kk[abc]/ac opts = new OptionTable from override'(options gb, opts ++ {Strategy => null}); runHooks(key, (opts, A, B), Strategy => strategy)); -- this is the logic for caching partial saturation computations. A.cache contains an option: -- SaturateContext{ mingens B } => SaturateComputation{ Result } container := fetchComputation(SaturateComputation, A, (A, B, opts), new SaturateContext from (A, B')); -- the actual computation of saturation occurs here C := (cacheComputation(opts, container)) computation; if C =!= null then doTrim C else if strategy === null then error("no applicable strategy for ", toString key) else error("assumptions for saturation strategy ", toString strategy, " are not met")) -- Helper for GRevLex strategy saturationByGRevLexHelper := (I, v, opts) -> ( R := ring I; (R1, fto, fback) := grevLexRing(index v, R); g1 := groebnerBasis(fto I, Strategy => "F4"); (g1', maxpower) := divideByVariable(g1, R1_(numgens R1 - 1)); if maxpower == 0 then (I, 0) else (ideal fback g1', maxpower)) -------------------------------------------------------------------- -- Algorithms for Module : Ideal^infinity algorithms#(saturate, Module, Ideal) = new MutableHashTable from { Iterate => (opts, M, I) -> ( M' := quotient(M, I, opts); while M' != M do ( M = M'; M' = quotient(M, I, opts) ); M ), -- TODO: if decompose is already cached, do this strategy instead Decompose => (opts, M, I) -> saturate(M, decompose I, Strategy => Iterate), } -- Installing hooks for Module : Ideal^infinity scan({Decompose, Iterate}, strategy -> addHook(key := (saturate, Module, Ideal), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- -- Algorithms for Module : RingElement^infinity algorithms#(saturate, Module, RingElement) = new MutableHashTable from { Iterate => (opts, M, f) -> ( M' := quotient(M, f, opts); while M' != M do ( M = M'; M' = quotient(M, f, opts) ); M ), } -- Installing hooks for Module : RingElement^infinity scan({Iterate}, strategy -> addHook(key := (saturate, Module, RingElement), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- -- Algorithms for Ideal : Ideal^infinity algorithms#(saturate, Ideal, Ideal) = new MutableHashTable from { -- TODO: this is sometimes faster than Eliminate Iterate => (opts, I, J) -> ( R := ring I; m := transpose generators J; while ( S := (ring I)/I; m = m ** S; I = ideal syz gb(m, Syzygies => true); I != 0) do (); ideal (presentation ring I ** R)), Eliminate => (opts, I, J) -> intersect apply(J_*, g -> saturate(I, g, opts)), GRevLex => (opts, I, J) -> ( -- FIXME: this might not be necessary, but the code isn't designed for this case. if not isFlatPolynomialRing ring I or not isHomogeneous I or not isHomogeneous J or not isGRevLexRing ring I then return null; -- First check that all generators are variables of the ring -- TODO: can this strategy work with generators of the irrelevant ideal? if any(index \ J_*, v -> v === null) then return null; -- Saturate with respect to each variable separately L := for g in J_* list saturationByGRevLexHelper(I, g, opts); -- Intersect them all -- TODO: when exactly is I returned? if any(last \ L, x -> x == 0) then I else intersect(first \ L)), Monomial => (opts, I, J) -> ( R := ring I; if not isMonomialIdeal I or not isMonomialIdeal J or not isPolynomialRing R or not isCommutative R then return null; cast := if instance(I, MonomialIdeal) then identity else ideal; -- TODO: make sure (monomialIdeal, MonomialIdeal) isn't forgetful cast newMonomialIdeal(R, rawSaturate(raw monomialIdeal I, raw monomialIdeal J))), } -- Installing hooks for Ideal : Ideal^infinity scan({Iterate, Eliminate, GRevLex, Monomial}, strategy -> addHook(key := (saturate, Ideal, Ideal), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- -- Algorithms for Ideal : RingElement^infinity algorithms#(saturate, Ideal, RingElement) = new MutableHashTable from { Iterate => (opts, I, f) -> saturate(I, ideal f, opts ++ {Strategy => Iterate}), -- backwards compatibility Linear => (opts, I, f) -> ( -- assumptions for this case: -- (1) the ring is of the form k[x1..xn]. No quotients, k a field or ZZ, grevlex order -- (2) all variables have degree 1. -- (3) I is homogeneous -- (4) f = homog linear form R := ring I; if not isFlatPolynomialRing R or not isGRevLexRing R or not all(gens R, g -> sum degree g == 1) or not isHomogeneous I or not isHomogeneous f or not isLinearForm f then return null; -- TODO: what does this do? res := newCoordinateSystem(R, matrix{{f}}); fto := res#1; fback := res#0; v := R_(numgens R - 1); g := gens gb(fto I, opts); ideal fback first divideByVariable(g, v)), Bayer => (opts, I, f) -> ( -- Bayer method. This may be used if I, f are homogeneous. -- Basic idea: in a ring R[z]/(f - z), with the RevLex order, compute GB of I. -- assumptions for this case: -- (1) the ring is of the form k[x1..xn]. No quotients, k a field or ZZ -- (2) I is homogeneous -- (3) f = homog, generated by one element R := ring I; if not isFlatPolynomialRing R or not isHomogeneous I or not isHomogeneous f or sum degree f === 0 then return null; n := numgens R; degs := append(degrees R, degree f); X := local X; R1 := (coefficientRing R)[X_0 .. X_n, Degrees => degs, MonomialSize => 16]; i := map(R1, R, (vars R1)_{0..n-1}); f1 := i f; I1 := ideal (i generators I); A := R1/(f1 - R1_n); -- TODO: add to ideal instead of quotient? iback := map(R, A, vars R | f); IA := generators I1 ** A; g := groebnerBasis(IA, Strategy => "F4"); -- TODO: compare with MGB (g1, notused) := divideByVariable(g, A_n); ideal iback g1), Eliminate => (opts, I, f) -> ( -- Eliminate(t, (I, t * f - 1)) -- assumptions for this case: -- I is an ideal in a flat polynomial ring (ring of the form k[x1..xn], no quotients, k a field or ZZ) -- f is an ideal, generated by one elem R := ring I; if not isFlatPolynomialRing R then return null; (R1, fto, fback) := eliminationInfo R; J := ideal(R1_0 * fto f - 1) + fto I; g := groebnerBasis(J, Strategy => "F4"); -- TODO: compare with MGB p1 := selectInSubring(1, g); ideal fback p1), GRevLex => (opts, I, v) -> ( -- FIXME: this might not be necessary, but the code isn't designed for this case. if not isHomogeneous I or not isFlatPolynomialRing ring I or not isGRevLexRing ring I then return null; -- First check that v is a variable of the ring -- TODO: can this strategy work with generators of the irrelevant ideal? if index v === null then return null; -- Saturate with respect to each variable separately first saturationByGRevLexHelper(I, v, opts)), "Unused" => (opts, I, f) -> ( -- NOT USED; TODO: assumptions? R := ring I; I1 := ideal 1_R; while I1 != I do ( I1 = I; I = ideal syz gb(matrix{{f}} | generators I, Syzygies => true, SyzygyRows => 1) ); I) } -- Installing hooks for Ideal : RingElement^infinity scan({"Unused", Iterate, Eliminate, GRevLex, Bayer, Linear}, strategy -> addHook(key := (saturate, Ideal, RingElement), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- -- isSupportedInZeroLocus -------------------------------------------------------------------- -- Note: this function isn't cached because the usecase in VirtualResolutions doesn't -- require it but it does take advantage of precomputed saturation of the annihilator isSupportedInZeroLocus = method() isSupportedInZeroLocus(Ideal, Ideal) := (B, I) -> isSupportedInZeroLocus(B, comodule I) isSupportedInZeroLocus(Ideal, Module) := (B, M) -> ( -- Returns true if M is supported only on the zero locus of B S := ring B; -- There are two ways to check this: -- 1. saturate the annihilator of M with respect to B and compare it to ideal(1) -- In general, computing the saturation is unnecessary, but when the saturation of the -- annihilator is known, it's faster to check whether saturate(annihilator M, B) == ideal 1 if M.cache.?annihilator then ( N := annihilator M; container := fetchComputation(SaturateComputation, N, (N, B, options saturate), new SaturateContext from (N, B)); if isComputationDone container then ( cacheHit container; return saturate(N, B) == ideal 1_S ) else true); -- 2. check that a high enough power of elements of B annihilates M n := numgens S; all(first entries mingens B, g -> ( supp := index \ support g; perm := toList(set(0..n-1) - set supp) | supp; R := (coefficientRing S)[(S_*)_perm, MonomialSize => 16, MonomialOrder => {Position => Up, n - #supp, #supp}]; -- TODO: can we compute gb once at the top instead? G := groebnerBasis(sub(presentation M, R), Strategy => "F4"); -- TODO: try "MGB" -- TODO: is "ann coker" necessary? 0 != ann coker selectInSubring(1, leadTerm G)))) -- check if a graded module sheafifies to zero isSupportedInZeroLocus(Ideal, GradedModule) := (B, G) -> ( all(min G .. max G, deg -> isSupportedInZeroLocus(B, G_deg))) -------------------------------------------------------------------- -- Annihilators -------------------------------------------------------------------- -- annihilator = method(Options => {Strategy => null}) -- defined in m2/quotient.m2 annihilator RingElement := Ideal => opts -> f -> annihilator(ideal f, opts) annihilator Ideal := Ideal => opts -> I -> annihilator(module I, opts) annihilator Module := Ideal => opts -> (cacheValue symbol annihilator) ( M -> annihilatorHelper(M, (annihilator, Module), opts)) -- Helper for annihilator methods annihilatorHelper = (A, key, opts) -> ( R := ring A; if isWeylAlgebra R then error "annihilator has no meaning for objects over a Weyl algebra"; -- TODO: add more instant checks f := presentation A; -- TODO: is this any different from A == 0? F := target f; if numgens F === 0 then return ideal 1_R; strategy := opts.Strategy; C := runHooks(key, (opts, A), Strategy => strategy); if C =!= null then C else if strategy === null then error("no applicable strategy for ", toString key) else error("assumptions for annihilator strategy ", toString strategy, " are not met")) -- Algorithms for annihilator Module algorithms#(annihilator, Module) = new MutableHashTable from { -- TODO: annihilator routines for MonomialIdeals? Quotient => (opts, M) -> ( f := presentation M; image f : target f), Intersection => (opts, M) -> ( f := presentation M; F := target f; intersect apply(numgens F, i -> ideal modulo(F_{i}, f))), } -- Installing hooks for annihilator Module scan({Quotient, Intersection}, strategy -> addHook(key := (annihilator, Module), algorithms#key#strategy, Strategy => strategy)) -------------------------------------------------------------------- ----- Development section -------------------------------------------------------------------- saturationByGRevLex = (I,J) -> saturate(I, J, Strategy => GRevLex) saturationByElimination = (I,J) -> saturate(I, J, Strategy => Eliminate) intersectionByElimination = L -> intersect(L, Strategy => Eliminate) -------------------------------------------------------------------- ----- Tests section -------------------------------------------------------------------- load "./Saturation/quotient-test.m2" load "./Saturation/saturate-test.m2" load "./Saturation/annihilator-test.m2" -------------------------------------------------------------------- ----- Documentation section -------------------------------------------------------------------- beginDocumentation() -- TODO: review load "./Saturation/doc.m2" load "./Saturation/quotient-doc.m2" load "./Saturation/saturate-doc.m2" load "./Saturation/annihilator-doc.m2" end-- restart debugLevel = 1 debug needsPackage "Saturation" kk = ZZ/32003 R = kk(monoid[x_0, x_1, x_2, x_3, x_4, Degrees => {2:{1, 0}, 3:{0, 1}}, Heft => {1,1}]) B0 = ideal(x_0,x_1) B1 = ideal(x_2,x_3,x_4) I = ideal(x_0^2*x_2^2*x_3^2+44*x_0*x_1*x_2^2*x_3^2+2005*x_1^2*x_2^2*x_3^2+12870 *x_0^2*x_2*x_3^3-725*x_0*x_1*x_2*x_3^3-15972*x_1^2*x_2*x_3^3-7768*x_0^2*x_2 ^2*x_3*x_4-13037*x_0*x_1*x_2^2*x_3*x_4-14864*x_1^2*x_2^2*x_3*x_4+194*x_0^2* x_2*x_3^2*x_4-2631*x_0*x_1*x_2*x_3^2*x_4-2013*x_1^2*x_2*x_3^2*x_4-15080*x_0 ^2*x_3^3*x_4-9498*x_0*x_1*x_3^3*x_4+5151*x_1^2*x_3^3*x_4-12401*x_0^2*x_2^2* x_4^2+4297*x_0*x_1*x_2^2*x_4^2-13818*x_1^2*x_2^2*x_4^2+7330*x_0^2*x_2*x_3*x _4^2-13947*x_0*x_1*x_2*x_3*x_4^2-12602*x_1^2*x_2*x_3*x_4^2-14401*x_0^2*x_3^ 2*x_4^2+8101*x_0*x_1*x_3^2*x_4^2-1534*x_1^2*x_3^2*x_4^2+8981*x_0^2*x_2*x_4^ 3-11590*x_0*x_1*x_2*x_4^3+1584*x_1^2*x_2*x_4^3-13638*x_0^2*x_3*x_4^3-5075*x _0*x_1*x_3*x_4^3-14991*x_1^2*x_3*x_4^3,x_0^7*x_2-6571*x_0^6*x_1*x_2+13908*x _0^5*x_1^2*x_2+11851*x_0^4*x_1^3*x_2+14671*x_0^3*x_1^4*x_2-14158*x_0^2*x_1^ 5*x_2-15190*x_0*x_1^6*x_2+6020*x_1^7*x_2+5432*x_0^7*x_3-8660*x_0^6*x_1*x_3- 3681*x_0^5*x_1^2*x_3+11630*x_0^4*x_1^3*x_3-4218*x_0^3*x_1^4*x_3+6881*x_0^2* x_1^5*x_3-6685*x_0*x_1^6*x_3+12813*x_1^7*x_3-11966*x_0^7*x_4+7648*x_0^6*x_1 *x_4-10513*x_0^5*x_1^2*x_4+3537*x_0^4*x_1^3*x_4+2286*x_0^3*x_1^4*x_4+733*x_ 0^2*x_1^5*x_4+11541*x_0*x_1^6*x_4+660*x_1^7*x_4); -- B0 B1 -- GRevLex 25.95s 0.18s -- Eliminate 28.35s 0.29s -- Iterate 60.02s 0.05s for B in {B0, B1} do ( for strategy in {GRevLex, Eliminate, Iterate} do print(strategy, (try elapsedTime saturate(I, B, Strategy => strategy);))) ans1 = elapsedTime saturationByGRevLex(saturationByGRevLex(I, B0), B1); -- 25.53s ans2 = elapsedTime saturationByGRevLex(saturationByGRevLex(I, B1), B0); -- 22.93s elapsedTime saturationByGRevLex(I, x_0); -- 9.01s elapsedTime saturationByGRevLex(I, x_1); -- 8.77s -- TODO: what a discrepancy ans3 = elapsedTime saturationByElimination(saturationByElimination(I, B0), B1); -- 49.22s ans4 = elapsedTime saturationByElimination(saturationByElimination(I, B1), B0); -- 28.63 elapsedTime J1 = saturationByElimination(I, x_0); elapsedTime J2 = saturationByElimination(I, x_1); elapsedTime J = intersectionByElimination(J1,J2); elapsedTime J' = intersectionByElimination(J2,J1); elapsedTime J'' = intersect(J1,J2); elapsedTime J''' = intersect(J2,J1); J == J' J == J'' time gens gb I; J2 = elapsedTime saturationByElimination(I, x_0); assert isHomogeneous J2 J2' = elapsedTime saturationByElimination(I, x_1); J2 = elapsedTime saturationByElimination(I, ideal(x_0,x_1)); J2' = elapsedTime saturationByElimination(J2, ideal(x_2,x_3,x_4)); J1 = elapsedTime saturate(I, x_0); J1' = elapsedTime saturate(I, x_1); J1 == J2 J1' == J2' betti J2 betti J1 restart load "./Saturation/badsaturations.m2" -- TODO: how was this so fast before?? J = paramRatCurve({2,2},{3,3},{4,2}); elapsedTime genSat(J,2) -- 200 sec elapsedTime genSat2(J,2) -- 50 sec elapsedTime genSat3(J,2) -- 35 sec J = paramRatCurve({2,2},{3,3},{5,2}); elapsedTime genSat(J,2) -- 691 sec elapsedTime genSat2(J,2) -- 104 sec elapsedTime genSat3(J,2) -- 71 sec J = paramRatCurve({2,2},{3,4},{4,3}); elapsedTime genSat(J,2) -- sec elapsedTime genSat2(J,2) -- sec elapsedTime genSat3(J,2) -- 75 sec I = ideal J_*_{5,13} use ring I elapsedTime I1 = saturate(I, x_0); elapsedTime (I2,pow) = saturationByGRevLex(I,x_0); I1 == I2 elapsedTime I1 = saturate(I, x_1); elapsedTime (I2,pow) = saturationByGRevLex(I,x_1); I1 == I2 elapsedTime J1 = intersectionByElimination(I1,I2); elapsedTime I1 = saturationByGRevLex(I, B0); elapsedTime I2 = saturationByGRevLex(I1, B1); elapsedTime saturationByGRevLex(saturationByGRevLex(I, B0), B1); elapsedTime saturationByGRevLex(saturationByGRevLex(I, B1), B0); elapsedTime saturationByElimination(saturationByElimination(I, B0), B1); elapsedTime saturationByElimination(saturationByElimination(I, B1), B0); elapsedTime J0a = saturationByGRevLex(I,x_0); elapsedTime J0b = saturationByGRevLex(I,x_1); --J1 = elapsedTime intersectionByElimination(first J0a,first J0b); La = elapsedTime trim first J0a; Lb = elapsedTime trim first J0b; J1 = elapsedTime intersectionByElimination(La, Lb); J1a = elapsedTime saturationByGRevLex(J1,x_2); J1b = elapsedTime saturationByGRevLex(J1,x_3); J1c = elapsedTime saturationByGRevLex(J1,x_4); J1a#1, J1b#1, J1c#1 J1ab = elapsedTime intersectionByElimination(J1a,J1b); elapsedTime J2 = intersectionByElimination{first J1a, first J1b, first J1c}; elapsedTime saturationByGRevLex(I,B0); saturationByElimination(I,x_0); (R1,fto,fback) = grevLexRing(0,S) L = fto I; satL = ideal first divideByVariable(gens gb L, R1_4); fback satL oo == I1 leadTerm oo ideal oo (R1,fto,fpack) = grevLexRing(1,S) use S R = ZZ/101[a..d] I = ideal"ab-ac,b2-cd" I1 = saturate(I,a) elapsedTime (I2,pow) = saturationByGRevLex(I,a); I1 == I2 pow (R1,fto,fback) = grevLexRing(0,R) fto I fto ---------------------------- -- Benchmarking example: restart needsPackage "Saturation" R = ZZ/101[vars(0..14)] M = genericMatrix(R, a, 3, 5) I = minors(3, M); codim I d = 4 J = ideal((gens I) * random(R^10, R^d)); -- algorithm; d = 2 3 4 5 -- null 0.45 40 -- Linear N/A N/A N/A N/A -- Iterate 0.41 40 -- Quotient 22 271 elapsedTime J' = quotient(J, I); for strategy in {Linear, Iterate, Quotient} do print(strategy, (try (elapsedTime J' === quotient(J, I, Strategy => strategy)) else "not applicable")) -- algorithm; d = 2 3 4 5 -- null 0.45 430 -- GRevLex N/A N/A N/A N/A -- Eliminate 2.87 378 -- Iterate 20 575 elapsedTime J'' = saturate(J, I); for strategy in {GRevLex, Eliminate, Iterate} do print(strategy, (try (elapsedTime J'' === saturate(J, I, Strategy => strategy)) else "not applicable")) elapsedTime quotient(J, I, Strategy => Iterate); elapsedTime saturate(J, I, Strategy => Eliminate); degree I elapsedTime(J : I_0);