From 87975bc73befb8a8429794fc9b6c7a03c77c7054 Mon Sep 17 00:00:00 2001 From: Mike Stillman Date: Wed, 21 Aug 2024 21:01:15 -0400 Subject: [PATCH] change `complex` to a method with Options => true. Add some Complex functions: status, rank, transpose, nullhomotopy (synonym of nullHomotopy). Allow complex to take a single matrix. Do not compute remainder in one version of nullHomotopy(serious performance degradation). Moved Ext(Module, Module) functionality from the Core to Complexes. Added minimalBetti fixes. --- M2/Macaulay2/packages/Complexes.m2 | 8 +- .../packages/Complexes/ChainComplex.m2 | 51 +++++++-- .../packages/Complexes/ChainComplexDoc.m2 | 77 ++++++++++++- .../packages/Complexes/ChainComplexMap.m2 | 9 +- .../packages/Complexes/ChainComplexMapDoc.m2 | 7 ++ .../packages/Complexes/ChainComplexTests.m2 | 4 +- M2/Macaulay2/packages/Complexes/Ext.m2 | 108 ++++++++++++++++++ .../packages/Complexes/FreeResolutionTests.m2 | 7 ++ .../packages/Complexes/FreeResolutions.m2 | 64 +++++++++++ 9 files changed, 315 insertions(+), 20 deletions(-) create mode 100644 M2/Macaulay2/packages/Complexes/Ext.m2 diff --git a/M2/Macaulay2/packages/Complexes.m2 b/M2/Macaulay2/packages/Complexes.m2 index b2e8a047ee6..1bdbd033473 100644 --- a/M2/Macaulay2/packages/Complexes.m2 +++ b/M2/Macaulay2/packages/Complexes.m2 @@ -46,9 +46,10 @@ export { "isNullHomotopyOf", "isShortExactSequence", "liftMapAlongQuasiIsomorphism", +-- "minimalBetti", "minimizingMap", "nullHomotopy", - "nullhomotopy" => "nullHomotopy", + --"nullhomotopy" => "nullHomotopy", "naiveTruncation", "randomComplexMap", -- "res" => "resolution", @@ -116,6 +117,7 @@ load "Complexes/ChainComplex.m2" load "Complexes/FreeResolutions.m2" load "Complexes/ChainComplexMap.m2" load "Complexes/Tor.m2" +load "Complexes/Ext.m2" -------------------------------------------------------------------- -- interface code to legacy types ---------------------------------- @@ -129,7 +131,7 @@ chainComplex Complex := ChainComplex => (cacheValue symbol ChainComplex) (C -> ( D )) -complex ChainComplex := Complex => opts -> (cacheValue symbol Complex)(D -> ( +complex ChainComplex := Complex => {} >> opts -> (cacheValue symbol Complex)(D -> ( (lo,hi) := (min D, max D); while lo < hi and (D_lo).numgens == 0 do lo = lo+1; while lo < hi and (D_hi).numgens == 0 do hi = hi-1; @@ -150,7 +152,7 @@ chainComplex ComplexMap := ChainComplexMap => f -> ( g ) -complex ChainComplexMap := ComplexMap => opts -> g -> ( +complex ChainComplexMap := ComplexMap => {} >> opts -> g -> ( map(complex target g, complex source g, i -> g_i, Degree => degree g) ) -------------------------------------------------------------------- diff --git a/M2/Macaulay2/packages/Complexes/ChainComplex.m2 b/M2/Macaulay2/packages/Complexes/ChainComplex.m2 index 3ae5d7ebef5..c3a0ce551c4 100644 --- a/M2/Macaulay2/packages/Complexes/ChainComplex.m2 +++ b/M2/Macaulay2/packages/Complexes/ChainComplex.m2 @@ -53,8 +53,10 @@ concentration ComplexMap := Sequence => f -> ( max Complex := ZZ => C -> max concentration C min Complex := ZZ => C -> min concentration C -complex = method(Options => {Base=>0}) -complex HashTable := Complex => opts -> maps -> ( +complexOptions = {Base => 0} +--complex = method(Options => {Base=>0}) +complex = method(Options => true) +complex HashTable := Complex => complexOptions >> opts -> maps -> ( spots := sort keys maps; if #spots === 0 then error "expected at least one matrix"; @@ -80,7 +82,7 @@ complex HashTable := Complex => opts -> maps -> ( C.dd = map(C,C,maps,Degree=>-1); C ) -complex List := Complex => opts -> L -> ( +complex List := Complex => complexOptions >> opts -> L -> ( -- L is a list of matrices or a list of modules if not instance(opts.Base, ZZ) then error "expected Base to be an integer"; @@ -104,7 +106,10 @@ complex List := Complex => opts -> L -> ( ); error "expected a list of matrices or a list of modules"; ) -complex Module := Complex => opts -> (M) -> ( +complex Matrix := Complex => complexOptions >> opts -> M -> ( + complex({M}, opts) + ) +complex Module := Complex => complexOptions >> opts -> (M) -> ( if not instance(opts.Base, ZZ) then error "complex: expected base to be an integer"; if M.cache.?Complex and opts.Base === 0 then return M.cache.Complex; @@ -118,9 +123,9 @@ complex Module := Complex => opts -> (M) -> ( C.dd = map(C,C,0,Degree=>-1); C ) -complex Ring := Complex => opts -> R -> complex(R^1, opts) -complex Ideal := Complex => opts -> I -> complex(module I, opts) -complex Complex := Complex => opts -> C -> ( +complex Ring := Complex => complexOptions >> opts -> R -> complex(R^1, opts) +complex Ideal := Complex => complexOptions >> opts -> I -> complex(module I, opts) +complex Complex := Complex => complexOptions >> opts -> C -> ( -- all this does is change the homological degrees -- so the concentration begins at opts.Base (lo,hi) := concentration C; @@ -133,7 +138,7 @@ complex Complex := Complex => opts -> C -> ( complex(L, Base=>opts.Base) ) ) -complex ComplexMap := Complex => opts -> f -> ( +complex ComplexMap := Complex => complexOptions >> opts -> f -> ( if degree f === -1 then ( if source f =!= target f then error "expected a differential"; (lo,hi) := concentration source f; @@ -483,6 +488,8 @@ defaultLengthLimit = (R, baselen, len) -> ( len ) +-- MES: note, this list of options is all of the ones from resolution, in the Core, +-- except FastNonminimal is not present (use instead: Strategy => Nonminimal). freeResolution = method(Options => { StopBeforeComputation => false, LengthLimit => infinity, -- (infinity means numgens R) @@ -533,6 +540,29 @@ freeResolution Matrix := ComplexMap => opts -> f -> extend( matrix f ) +-- TODO: reinstate these once we remove all uses of ChainComplex... +-- resolution Module := Complex => opts -> M -> ( +-- o := pairs opts; +-- o2 := new OptionTable from select(pairs opts, x -> x#0 =!= FastNonminimal); +-- if opts.FastNonminimal then ( +-- o2 = o2 ++ {Strategy => Nonminimal}; +-- << "warning: `FastNonminimal => true` is deprecated. Use: res(..., Strategy => Nonminimal) instead" << endl; +-- ); +-- freeResolution(M, o2) +-- ) +-- resolution Ideal := Complex => opts -> I -> resolution(comodule I, opts) +-- resolution MonomialIdeal := Complex => opts -> I -> resolution(comodule ideal I, opts) +-- resolution Matrix := ComplexMap => opts -> f -> extend( +-- resolution(target f, opts), +-- resolution(source f, opts), +-- matrix f +-- ) + +complete Complex := C -> C +complete ComplexMap := F -> F +nullhomotopy ComplexMap := F -> nullHomotopy F +status Complex := C -> << "resolution status of a Complex needs to be implemented" << endl; + isHomogeneous Complex := (C) -> isHomogeneous dd^C -- These next two local functions are lifted from previous code in m2/chaincomplexes.m2 @@ -592,6 +622,11 @@ poincareN Complex := C -> ( f ) +rank Complex := ZZ => C -> ( + (lo, hi) := concentration C; + sum for i from lo to hi list (-1)^i * rank C_i + ) + minimalPresentation Complex := prune Complex := Complex => opts -> (cacheValue symbol minimalPresentation)(C -> ( -- opts is ignored here diff --git a/M2/Macaulay2/packages/Complexes/ChainComplexDoc.m2 b/M2/Macaulay2/packages/Complexes/ChainComplexDoc.m2 index 6d303e275d1..a9a78e4b447 100644 --- a/M2/Macaulay2/packages/Complexes/ChainComplexDoc.m2 +++ b/M2/Macaulay2/packages/Complexes/ChainComplexDoc.m2 @@ -408,9 +408,10 @@ doc /// doc /// Key - complex (complex, List) - [complex, Base] + complex + (complex, Matrix) + [(complex, List), Base] Base Headline make a chain complex @@ -418,7 +419,7 @@ doc /// complex L Inputs L:List - of maps + of maps, or a single matrix. Base => ZZ the index of the target of the first map in the differential. @@ -483,6 +484,14 @@ doc /// HH C prune HH C prune HH_1 C + Text + Having the input be a matrix is equivalent to having it be + a singleton list containing this matrix. + Example + C' = complex F1 + assert isWellDefined C' + C'' = complex(F1, Base => 3) + assert isWellDefined C'' Caveat This constructor minimizes computation and does very little error checking. To verify that a complex @@ -4363,6 +4372,66 @@ doc /// freeResolution /// +doc /// + Key + (Ext,Module,Module) + (Ext,Ideal,Ideal) + (Ext,Ideal,Module) + (Ext,Ideal,Ring) + (Ext,Module,Ideal) + (Ext,Module,Ring) + Headline + total Ext module + Usage + Ext(M, N) + Inputs + M:Module + or ofClass{Ideal,Ring}, that is homogeneous + N:Module + or ofClass{Ideal,Ring} over the same ring as $M$, that is also homogeneous + Outputs + :Module + the $\operatorname{Ext}$ module of $M$ and $N$, as a + multigraded module, with the modules + $\operatorname{Ext}^i(M,N)$ for all values of $i$ + appearing simultaneously. + Description + Text + The computation of the total Ext module is possible for modules over the + ring $R$ of a complete intersection, according to the algorithm + of Shamash-Eisenbud-Avramov-Buchweitz. The result is provided as a finitely + presented module over a new ring with one additional variable of degree + $\{-2,-d\}$ for each equation of degree $d$ defining $R$. The + variables in this new ring have degree length $1$ more than the degree length of + the original ring. In other words, it is multigraded with the + degree $d$ part of $\operatorname{Ext}^n(M,N)$ appearing as the degree + $\{-n,d\}$ part of $\operatorname{Ext}(M,N)$. + Text + We illustrate this in the following example. + Example + R = QQ[x,y]/(x^3,y^2); + N = cokernel matrix {{x^2, x*y}} + H = Ext(N,N); + ring H + S = ring H; + H + isHomogeneous H + rank source basis( { -2,-3 }, H) + rank source basis( { -3 }, Ext^2(N,N) ) + rank source basis( { -4,-5 }, H) + rank source basis( { -5 }, Ext^4(N,N) ) + hilbertSeries H + hilbertSeries(H,Order=>11) + Text + For more information, see the chapter {\it Resolutions and cohomology over complete intersections} + by Luchezar L. Avramov and Daniel R. Grayson, in the book + @HREF("https://macaulay2.com/Book/ComputationsBook/book/book.pdf", "Computatations in Algebraic Geometry with Macaulay2")@. + Text + The result of the computation is cached for future reference. + SeeAlso + "computing with Ext" +/// + /// Key Headline @@ -4375,3 +4444,5 @@ doc /// Caveat SeeAlso /// + + diff --git a/M2/Macaulay2/packages/Complexes/ChainComplexMap.m2 b/M2/Macaulay2/packages/Complexes/ChainComplexMap.m2 index 86e8b358e1e..d0aa0d9747b 100644 --- a/M2/Macaulay2/packages/Complexes/ChainComplexMap.m2 +++ b/M2/Macaulay2/packages/Complexes/ChainComplexMap.m2 @@ -634,6 +634,7 @@ Hom(Matrix, ComplexMap) := ComplexMap => opts -> (f,g) -> Hom(map(complex target f, complex source f, i -> if i === 0 then f), g, opts) dual ComplexMap := ComplexMap => {} >> o -> f -> Hom(f, (ring f)^1) +transpose ComplexMap := ComplexMap => f -> dual f homomorphism ComplexMap := ComplexMap => (h) -> ( -- h should be a homomorphism of complexes from R^1[-i] --> E = Hom(C,D) @@ -1116,13 +1117,13 @@ nullHomotopyFreeSource = f -> ( (lo,hi) := concentration f; for i from lo to hi do ( if hs#?(i-1) then ( - rem := (f_i - hs#(i-1) * dd^C_i) % (dd^D_(i+deg)); - if rem != 0 then return null; -- error "can't construct homotopy"; + --rem := (f_i - hs#(i-1) * dd^C_i) % (dd^D_(i+deg)); + --if rem != 0 then return null; -- error "can't construct homotopy"; hs#i = (f_i - hs#(i-1) * dd^C_i) // (dd^D_(i+deg)) ) else ( - rem = f_i % dd^D_(i+deg); - if rem != 0 then return null; -- error "can't construct homotopy"; + --rem = f_i % dd^D_(i+deg); + --if rem != 0 then return null; -- error "can't construct homotopy"; hs#i = f_i // dd^D_(i+deg) ) ); diff --git a/M2/Macaulay2/packages/Complexes/ChainComplexMapDoc.m2 b/M2/Macaulay2/packages/Complexes/ChainComplexMapDoc.m2 index 6795b53ce06..ed084385121 100644 --- a/M2/Macaulay2/packages/Complexes/ChainComplexMapDoc.m2 +++ b/M2/Macaulay2/packages/Complexes/ChainComplexMapDoc.m2 @@ -1271,6 +1271,7 @@ doc /// doc /// Key (dual, ComplexMap) + (transpose, ComplexMap) Headline the dual of a map of complexes Usage @@ -1299,10 +1300,16 @@ doc /// D' = (freeResolution coker matrix{{a^2,a*b,c^3}})[-1] f' = randomComplexMap(D', D) dual(f' * f) == dual f * dual f' + Text + For backwards compatibility, one can also use {\tt transpose}. + Example + h' = transpose f + assert(h == h') SeeAlso (Hom, ComplexMap, ComplexMap) (randomComplexMap, Complex, Complex) (dual, Matrix) + (transpose, Matrix) /// doc /// diff --git a/M2/Macaulay2/packages/Complexes/ChainComplexTests.m2 b/M2/Macaulay2/packages/Complexes/ChainComplexTests.m2 index ea4f0509c28..a506951c497 100644 --- a/M2/Macaulay2/packages/Complexes/ChainComplexTests.m2 +++ b/M2/Macaulay2/packages/Complexes/ChainComplexTests.m2 @@ -2280,8 +2280,8 @@ TEST /// beta = map(E, D, {d,e,f}) isWellDefined alpha isWellDefined beta - phi = complex chainComplex alpha - psi = complex chainComplex beta + phi = complex alpha + psi = complex beta prune HH phi prune HH psi isNullHomotopic phi diff --git a/M2/Macaulay2/packages/Complexes/Ext.m2 b/M2/Macaulay2/packages/Complexes/Ext.m2 new file mode 100644 index 00000000000..a5d27844a9c --- /dev/null +++ b/M2/Macaulay2/packages/Complexes/Ext.m2 @@ -0,0 +1,108 @@ +-- total ext over complete intersections + +-- Code originally written by Dan Grayson +-- change 2009/1/12: +-- some modifications contributed 23 Mar 2007 by "Paul S. Aspinwall" , +-- but we also get rid of the fudge factor entirely, depending instead on automatic +-- computation of the heft vector + +-- TODO: Ext(R, S) should work as well +Ext(Ring, Ring) := +Ext(Ring, Ideal) := +Ext(Ring, Module) := +Ext(Ideal, Ring) := +Ext(Ideal, Ideal) := +Ext(Ideal, Module) := +Ext(Module, Ring) := +Ext(Module, Ideal) := Module => opts -> (M, N) -> Ext(module M, module N, opts) +Ext(Module, Module) := Module => opts -> (M, N) -> ( + -- c.f. caching in Hom(Module,Module) + cacheModule := youngest(M.cache.cache, N.cache.cache); + cacheKey := (Ext,M,N); + if cacheModule#?cacheKey then return cacheModule#cacheKey; + B := ring M; + if B =!= ring N + then error "expected modules over the same ring"; + if not isCommutative B + then error "'Ext' not implemented yet for noncommutative rings."; + if not isHomogeneous B + then error "'Ext' received modules over an inhomogeneous ring"; + if not isHomogeneous N or not isHomogeneous M + then error "'Ext' received an inhomogeneous module"; + if N == 0 or M == 0 then return cacheModule#cacheKey = B^0; + p := presentation B; + A := ring p; + I := ideal mingens ideal p; + n := numgens A; + c := numgens I; + if c =!= codim B + then error "total Ext available only for complete intersections"; + f := I_*; -- apply(c, i -> I_i); + pM := lift(presentation M,A); + pN := lift(presentation N,A); + M' := cokernel ( pM | p ** id_(target pM) ); + N' := cokernel ( pN | p ** id_(target pN) ); + assert isHomogeneous M'; + assert isHomogeneous N'; + C := freeResolution M'; + X := getSymbol "X"; + K := coefficientRing A; + S := K(monoid [X_1 .. X_c, toSequence A.generatorSymbols, + Degrees => { + apply(0 .. c-1, i -> prepend(-2, - degree f_i)), + apply(0 .. n-1, j -> prepend( 0, degree A_j)) + }]); + -- make a monoid whose monomials can be used as indices + Rmon := monoid [X_1 .. X_c,Degrees=>{c:{2}}]; + -- make group ring, so 'basis' can enumerate the monomials + R := K Rmon; + -- make a hash table to store the blocks of the matrix + blks := new MutableHashTable; + blks#(exponents 1_Rmon) = C.dd; + scan(0 .. c-1, i -> + blks#(exponents Rmon_i) = nullHomotopy (- f_i*id_C)); + -- a helper function to list the factorizations of a monomial + factorizations := (gamma) -> ( + -- Input: gamma is the list of exponents for a monomial + -- Return a list of pairs of lists of exponents showing the + -- possible factorizations of gamma. + if gamma === {} then { ({}, {}) } + else ( + i := gamma#-1; + splice apply(factorizations drop(gamma,-1), + (alpha,beta) -> apply (0..i, + j -> (append(alpha,j), append(beta,i-j)))))); + scan(4 .. length C + 1, d -> if even d then ( + scan( flatten \ exponents \ leadMonomial \ first entries basis(d,R), + gamma -> ( + s := - sum(factorizations gamma, + (alpha,beta) -> ( + if blks#?alpha and blks#?beta + then blks#alpha * blks#beta + else 0)); + -- compute and save the nonzero nullhomotopies + if s != 0 then blks#gamma = nullhomotopy s)))); + -- make a free module whose basis elements have the right degrees + (loC, hiC) := concentration C; + Cstar := S^(apply(toList(loC..hiC), + i -> toSequence apply(degrees C_i, d -> prepend(i,d)))); + -- assemble the matrix from its blocks. + -- We omit the sign (-1)^(n+1) which would ordinarily be used, + -- which does not affect the homology. + toS := map(S,A,apply(toList(c .. c+n-1), i -> S_i), + DegreeMap => prepend_0); + Delta := map(Cstar, Cstar, + transpose sum(keys blks, m -> S_m * toS sum blks#m), + Degree => { -1, degreeLength A:0 }); + DeltaBar := Delta ** (toS ** N'); + if debugLevel > 10 then ( + assert isHomogeneous DeltaBar; + assert(DeltaBar * DeltaBar == 0); + stderr << describe ring DeltaBar < 2) /// + +TEST /// + R = QQ[x, Degrees => {{}}] + M = image x + F = freeResolution M -- gives: "error: map with index 0 has inconsistent source" + assert(isWellDefined F) +/// diff --git a/M2/Macaulay2/packages/Complexes/FreeResolutions.m2 b/M2/Macaulay2/packages/Complexes/FreeResolutions.m2 index 6ffef2fcb0d..7e528fd28fd 100644 --- a/M2/Macaulay2/packages/Complexes/FreeResolutions.m2 +++ b/M2/Macaulay2/packages/Complexes/FreeResolutions.m2 @@ -601,7 +601,67 @@ truncate(BettiTally, InfiniteNumber, InfiniteNumber) := BettiTally => {} >> opts else continue ) +----------------------------------------------------------------------------- +-- minimalBetti +----------------------------------------------------------------------------- + +-- TODO: place this here, not in Core +-- minimalBetti = method( +-- TypicalValue => BettiTally, +-- Options => { +-- DegreeLimit => null, +-- LengthLimit => infinity, +-- Weights => null, +-- ParallelizeByDegree => false -- currently: only used over primes fields of positive characteristic +-- }) + +--- version 1.24.05 version of minimalBetti. -* +minimalBetti Module := BettiTally => opts -> M -> ( + R := ring M; + degreelimit := resolutionDegreeLimit(R, opts.DegreeLimit); + lengthlimit := resolutionLengthLimit(R, opts.LengthLimit); + -- check to see if a cached resolution is sufficient + cacheKey := ResolutionContext{}; + if M.cache#?cacheKey and isComputationDone(C := M.cache#cacheKey, + DegreeLimit => degreelimit, LengthLimit => lengthlimit) + then return betti(C.Result.Resolution, Weights => opts.Weights); + -- if not, compute a fast non-minimal resolution + -- the following line is because we need to make sure we have the resolution + -- either complete, or one more than the desired minimal betti numbers. + + -- We see if we can now compute a non-minimal resolution. + -- If not, we compute a usual resolution. + -- TODO: this isn't quite correct. + useFastNonminimal := not isQuotientRing R and + char R > 0 and char R < (1<<15); + + if not useFastNonminimal then + return betti resolution(M, DegreeLimit => degreelimit, LengthLimit => lengthlimit); + -- At this point, we think we are good to use the faster algorithm. + -- First, we need to comppute the non-minimal resolution to one further step. + if instance(opts.LengthLimit, ZZ) then lengthlimit = lengthlimit + 1; + C = resolution(M, + StopBeforeComputation => true, FastNonminimal => true, ParallelizeByDegree => opts.ParallelizeByDegree, + DegreeLimit => degreelimit, LengthLimit => lengthlimit); + rC := if C.?Resolution and C.Resolution.?RawComputation then C.Resolution.RawComputation + -- TODO: when can this error happen? + else error "cannot use 'minimalBetti' with this input. Input must be an ideal or module in a + polynomial ring or skew commutative polynomial ring over a finite field, which is singly graded. + These restrictions might be removed in the future."; + -- + B := unpackEngineBetti rawMinimalBetti(rC, + if opts.DegreeLimit =!= null then {opts.DegreeLimit} else {}, + if opts.LengthLimit =!= infinity then {opts.LengthLimit} else {} + ); + betti(B, Weights => heftvec(opts.Weights, heft R)) + ) +minimalBetti Ideal := BettiTally => opts -> I -> minimalBetti( + if I.cache.?quotient then I.cache.quotient + else I.cache.quotient = cokernel generators I, opts + ) +*- + minimalBetti Module := BettiTally => opts -> M -> ( R := ring M; degreelimit := opts.DegreeLimit; @@ -640,6 +700,10 @@ minimalBetti Module := BettiTally => opts -> M -> ( betti(B, Weights => heftvec(opts.Weights, heft R)) ) +minimalBetti Ideal := BettiTally => opts -> I -> minimalBetti(comodule I, opts) + +-* +-- older version of Core version of minimalBetti. Can't use here? minimalBetti(Module, Thing) := BettiTally => opts -> (M, junk) -> ( R := ring M; degreelimit := resolutionDegreeLimit(R, opts.DegreeLimit);