added new algorithm for Hom(Module,Module)

Implemented with Devlin Mallory and David Eisenbud,
this algorithm takes an optional argument DegreeLimit,
which is passed on to syz.

M2/Macaulay2/m2/modules.m2

@@ -385,7 +385,7 @@ super(Module) := Module => (M) -> (
 -- Homomorphisms and Endomorphisms
 -----------------------------------------------------------------------------
 
-Hom = method(Options => { MinimalGenerators => false })
+Hom = method(Options => { DegreeLimit => null, MinimalGenerators => false })
 
 End = method(Options => options Hom)
 End(Thing) := o -> M -> Hom(M, M, o)

M2/Macaulay2/m2/modules2.m2

@@ -271,14 +271,25 @@ Hom(Module, Ring)   :=
 Hom(Module, Ideal)  := Module => opts -> (M, N) -> Hom(M, module N, opts)
 
 Hom(Module, Module) := Module => opts -> (M, N) -> (
-     Y := youngest(M.cache.cache,N.cache.cache);
-     if Y#?(Hom,M,N) then return Y#(Hom,M,N);
-     trim' := if opts.MinimalGenerators then trim else identity;
-     H := trim' kernel (transpose presentation M ** N);
-     H.cache.homomorphism = (f) -> map(N,M,adjoint'(f,M,N), Degree => first degrees source f + degree f);
-     Y#(Hom,M,N) = H; -- a hack: we really want to type "Hom(M,N) = ..."
-     H.cache.formation = FunctionApplication { Hom, (M,N) };
-     H)
+    -- TODO: take advantage of cached results with higher e
+    e := opts.DegreeLimit;
+    Y := youngest(M.cache.cache, N.cache.cache);
+    if Y#?(Hom, M, N, e) then return Y#(Hom, M, N, e);
+    -- Previously Hom(M,N) was computed simply as:
+    -- H := trim kernel(transpose presentation M ** N);
+    -- This algorithm is more flexible, and hopefully not slower
+    A := presentation M; (G, F) := (target A, source A); -- M <-- G <-- F
+    B := presentation N; (L, K) := (target B, source B); -- N <-- L <-- K
+    piN := inducedMap(N, L, generators N);
+    psi := (Hom(A, N) * Hom(G, piN)) // Hom(F, piN);
+    p := id_(Hom(G, L)) | map(Hom(G, L), Hom(F, K), 0);
+    r := syz(psi | -Hom(F, B), DegreeLimit => e);
+    trim' := if opts.MinimalGenerators then trim else identity;
+    Y#(Hom, M, N, e) =
+    H := trim' image(Hom(G, piN) * p * r);
+    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, e) };
+    H)
 
 adjoint' = method()
 adjoint'(Matrix,Module,Module) := Matrix => (m,G,H) -> (