[Draft] Compute Hom / Ext via pushforward

M2/Macaulay2/m2/Hom.m2

@@ -7,12 +7,13 @@ needs "modules.m2"
 needs "matrix1.m2"
 
 Hom = method(Options => {
-	DegreeLimit       => null,
-	MinimalGenerators => true,
-	Strategy          => null,
-	})
+    DegreeLimit         => null,
+    MinimalGenerators   => true,
+    Strategy            => null
+})
 
-End = options Hom >> o -> M -> Hom(M, M, o)
+End = method(Options => options Hom)
+End(Thing) := Thing => o -> X -> Hom(X, X, o)
 
 -- TODO: Hom(R, S) should work as well
 Hom(Ring, Ring)   :=
@@ -32,16 +33,19 @@ Hom(Module, Module) := Module => opts -> (M, N) -> (
     -- M.cache is a hashless (hence ageless) CacheTable, but
     -- M.cache.cache is a MutableHashTable, hence has an age.
     Y := youngest(M.cache.cache, N.cache.cache);
-    if Y#?(Hom, M, N, e) then return Y#(Hom, M, N, e);
-    H := runHooks((Hom, Module, Module), (opts, M, N), Strategy => opts.Strategy);
-    if H === null then error "Hom: no strategy found for the given input";
-    trim' := if opts.MinimalGenerators then trim else identity;
-    -- a hack: we really want to type "Hom(M, N) = ..."
-    H = Y#(Hom, M, N, e) = if opts.MinimalGenerators then trim H else H;
-    H.cache.homomorphism = f -> map(N, M, adjoint'(f, M, N), Degree => first degrees source f + degree f);
-    H.cache.formation = FunctionApplication { Hom, (M, N, DegreeLimit => e) };
-    H)
-
+    Y#(Hom, M, N, e) ?? (
+        H := runHooks((Hom, Module, Module), (opts, M, N), Strategy => opts.Strategy);
+        if H === null then error "Hom: no strategy found for the given input";
+        trim' := if opts.MinimalGenerators then trim else identity;
+        -- a hack: we really want to type "Hom(M, N) = ..."
+        H = Y#(Hom, M, N, e) = if opts.MinimalGenerators then trim H else H;
+        H.cache.homomorphism = f -> map(N, M, adjoint'(f, M, N), Degree => first degrees source f + degree f);
+        H.cache.formation = FunctionApplication { Hom, (M, N, DegreeLimit => e) };
+        H
+    )
+)
+
+rawTranspose = m -> transpose matrix for row in entries m list for x in row list antipode(x);
 basicHom = (M, N) -> kernel(transpose presentation M ** N)
 addHook((Hom, Module, Module), Strategy => Default,  (opts, M, N) -> basicHom(M, N))
 addHook((Hom, Module, Module), Strategy => Syzygies, (opts, M, N) -> (

M2/Macaulay2/m2/centerRing.m2

@@ -0,0 +1,52 @@
+-*
+computes the center of a ring R and returns the map from center(R) -> R
+
+maybe this should just be called "center" but the term is ubiquitous across
+mathematics and appears in several different. trying to merge them all into a
+single method runs into the awkwardness of method options leaking out to all
+the different overrides.
+*-
+
+centerRing = method(Options => {Strategy => null})
+centerRing(Ring) := RingMap => opts -> (R) -> R.cache.centerRing ??= (
+    z := runHooks((centerRing, Ring), (opts, R), Strategy => opts.Strategy);
+    if z =!= null then z else error "centerRing: no strategy for computing center of the given ring"
+)
+
+addHook((centerRing, Ring), Strategy => "isCommutative", (opts, R) -> if isCommutative R then map(R, R) else null)
+addHook((centerRing, Ring), Strategy => "isSkewAffineRing", (opts, R) -> (
+    -- David Eisenbud observed this strategy for computing the center of a skew polynomial ring
+    if not isPolynomialRing R then return null;
+    if not isSkewAffineRing R then return null;
+
+    (R', phi) := flattenRing R;
+    -- get the natural prime coefficient ring
+    kk := coefficientRing R';
+
+    x := symbol x; -- for pairs of skew variables
+    y := symbol y; -- for commuting variables
+    skews := R'.SkewCommutative;
+    n := #skews;
+    xvars := flatten for i from 0 to n - 1 list
+        for j from i + 1 to n - 1 list
+            x_{i, j};
+    ximages := flatten for i from 0 to n - 1 list
+        for j from i + 1 to n - 1 list
+            R'_(skews_i) * R'_(skews_j);
+
+    commuters := toList (0..numgens R' - 1) - set skews;
+    yvars := toList(y_0..y_(#commuters - 1));
+    yimages := apply(commuters, i -> R'_i);
+    S := kk[xvars, yvars, Degrees => apply(ximages| yimages, degree)];
+
+    -- hack to lazy-load PushForward package only when this hook needs it
+    PushForward := needsPackage "PushForward";
+    pushFwd := value PushForward.Dictionary#"pushFwd";
+
+    -- the ideal defining the center is the annihilator of R^1 as a module over the polynomial ring S
+    K := annihilator pushFwd(map(R', S, ximages | yimages), R'^1);
+    Z := S/K;
+    Z.formation = FunctionApplication {centerRing, R};
+
+    map(R, Z, phi^-1 matrix {ximages | yimages})
+))

M2/Macaulay2/m2/exports.m2

@@ -539,6 +539,7 @@ export {
 	"capture",
 	"catch",
 	"ceiling",
+        "centerRing",
 	"centerString",
 	"changeBase",
 	"changeDirectory",
@@ -704,6 +705,7 @@ export {
 	"flagLookup",
 	"flatten",
 	"flattenRing",
+	"flattenDegreeMap",
 	"flip",
 	"floor",
 	"flush",

M2/Macaulay2/m2/loadsequence

@@ -80,6 +80,7 @@ monomcurve.m2
 fano.m2
 schubert.m2
 local.m2
+centerRing.m2
 
 packages.m2
 document.m2

M2/Macaulay2/m2/pushforward.m2

@@ -111,7 +111,36 @@ ordertab := new HashTable from {
     Eliminate    => (nR, nS) -> Eliminate nR,
     ProductOrder => (nR, nS) -> ProductOrder{nR, nS},
     Lex          => (nR, nS) -> Lex,
-    }
+}
+
+monomialOrderFromOpts = (opts) -> (
+    monorder := opts.MonomialOrder;
+    monorder = if ordertab#?monorder then ordertab#monorder
+    else error("pushForward: MonomialOrder option expected one of ", demark_", " \\ toString \ keys ordertab)
+)
+
+-- a splitting of f.cache.DegreeMap to use as a DegreeMap when constructing the
+-- projection from graphRing f to source f.
+makeSectionForDegreeMap = (f) -> (
+    if degreeLength source f === 0 then return null;
+
+    G := degreeGroup source f;
+    -- the linear map that f induces on DegreeGroups as a matrix
+    M := transpose matrix apply(entries G_{0..numgens G - 1}, f.cache.DegreeMap);
+
+    -- if M is not injective then we cannot produce a section
+    if (kernel M != 0) then return null;
+
+    D := image M;
+    (d) -> (
+        -- d' is the element of degreeGroup target f corresponding to d
+        d' := matrix transpose {d};
+        -- get coefficients for d' against generators for G ~ D.
+        -- this produces mild nonsense if d' is not in D but that doesn't matter
+        -- since we only apply the resulting projecting to image f.
+        flatten entries(d' // inducedMap(ambient D, D) // coverMap D)
+    )
+)
 
 pushNonLinear := (opts, f0, M) -> (
     -- given f: R --> S, and M an S-module, finite over R,
@@ -129,19 +158,13 @@ pushNonLinear := (opts, f0, M) -> (
     if isHomogeneous M  then assert isHomogeneous m;
     if isHomogeneous f0 then assert isHomogeneous f;
 
-    -- set up some variables that are used throughout
-    deglenS := degreeLength S;
-    numgensS := numgens S;
-    numgensR := numgens R;
-
-    monorder := opts.MonomialOrder;
-    monorder  = if ordertab#?monorder then (ordertab#monorder)(numgensR, numgensS)
-    else error("pushForward: MonomialOrder option expected one of ",
-	demark_", " \\ toString \ keys ordertab);
+    s := numgens S;
+    r := numgens R;
+    monorder := (monomialOrderFromOpts opts)(r, s);
 
     J := graphIdeal(f, MonomialOrder => monorder, VariableBaseName => local X);
     G := ring J;
-    xvars := map(G, R, submatrix(vars G, toList(0..numgensR - 1)));
+    xvars := map(G, R, submatrix(vars G, toList(0..r - 1)));
     m1 := presentation (cokernel xvars m  **  cokernel generators J);
 
     if opts.UseHilbertFunction and all({f, m}, isHomogeneous) then (
@@ -150,26 +173,25 @@ pushNonLinear := (opts, f0, M) -> (
 	T := degreesRing G;
 	hf = hf * product(degrees source generators J, d -> 1 - T_d);
 	-- cache poincare
-	poincare cokernel m1 = hf);
-
-    mapbackdeg := d -> take(d, -deglenS);
-    -- that choice of degree map was chosen to make the symmetricPower functor homogeneous, but it doesn't have much
-    -- else to recommend it.
-    -- we should really be *lifting* the result to S along the natural map S ---> G
-    mapback := map(S, G, map(S^1, S^numgensR, 0) | vars S, DegreeMap => mapbackdeg );
-
-    -- let's at least check it splits f's degree map:
-    for i from 0 to deglenS-1 do (
-	e := for j from 0 to deglenS-1 list if i === j then 1 else 0;
-	if mapbackdeg f.cache.DegreeMap e =!= e
-	then error "not implemented yet: unexpected degree map of ring map");
+	poincare cokernel m1 = hf
+    );
 
     g := gb(m1,
 	StopBeforeComputation => opts.StopBeforeComputation,
 	DegreeLimit           => opts.DegreeLimit,
-	PairLimit             => opts.PairLimit);
+	PairLimit             => opts.PairLimit
+    );
+
+    -- todo - what should happen here if the map on degree groups induced by f is not injective?
+    mapback := if (degmapback := makeSectionForDegreeMap(f)) =!= null then
+        map(S, G, map(S^1, S^r, 0) | vars S, DegreeMap => degmapback)
+    else
+        map(S, G, map(S^1, S^r, 0) | vars S);
+
+    result := phiS^-1 mapback selectInSubring(if r > 0 then 1 else 0, generators g);
     -- MES: check if the monomial order restricts to R.  If so, then do `` forceGB result ''
-    phiS^-1 mapback selectInSubring(if numgensR > 0 then 1 else 0, generators g))
+    result
+)
 
 -*
 pushLinear := opts -> (f,M) -> (