Improvements to Simplicial Toric Varieties

M2/Macaulay2/m2/basis.m2

@@ -65,6 +65,7 @@ findHeftandVars := (R, varlist, ndegs) -> (
 
 -- TODO: can this be done via a tensor or push forward?
 -- c.f. https://github.com/Macaulay2/M2/issues/1522
+-- TODO: think about how this compares with part
 liftBasis = (M, phi, B, offset) -> (
     -- lifts a basis B of M via a ring map phi
     (R, S) := (phi.target, phi.source);
@@ -221,6 +222,25 @@ basis(Sequence, RingMap) := Matrix => opts -> (degs, f) -> (
 
 -----------------------------------------------------------------------------
 
+-- Note: f needs to be homogeneous, otherwise returns nonsense
+basis(ZZ,       RingMap) :=
+basis(List,     RingMap) := Matrix => opts -> (degs, f) -> basis((degs, degs), f, opts)
+basis(Sequence, RingMap) := Matrix => opts -> (degs, f) -> (
+    -- if not isHomogeneous f then error "expected a graded ring map (try providing a DegreeMap)";
+    if #degs != 2          then error "expected a sequence (tardeg, srcdeg) of degrees for target and source rings";
+    (T, S) := (target f, source f);
+    (tardeg, srcdeg) := degs;
+    if tardeg === null then tardeg = f.cache.DegreeMap  srcdeg;
+    if srcdeg === null then srcdeg = f.cache.DegreeLift tardeg;
+    targens := basis(tardeg, T);
+    srcgens := basis(srcdeg, S);
+    -- TODO: should matrix RingMap return this map instead?
+    -- mon := map(module T, image matrix(S, { S.FlatMonoid_* }), f, matrix f)
+    mat := last coefficients(f cover srcgens, Monomials => targens);
+    map(image targens, image srcgens, f, mat))
+
+-----------------------------------------------------------------------------
+
 basisHelper = (opts, lo, hi, M) -> (
     R := ring M;
     n := degreeLength R;
@@ -300,7 +320,7 @@ algorithms#(basis, List, List, Module) = new MutableHashTable from {
 	G := degreeGroup(R := ring M);
 	if G == 0 or isFreeModule G then return null;
 	(S, phi, psi, p, f) := permuteDegreeGroup R;
-	B := psi map_S basisDefaultStrategy(opts, f lo, f hi, phi M);
+	B := psi map_S basisDefaultStrategy(opts, f lo, f hi, phi ** M);
 	L := positions(degrees source B, deg -> lo <= deg and deg <= hi);
 	raw submatrix(B, , L)),
     -- TODO: add separate strategies for skew commutative rings, vector spaces, and ZZ-modules
@@ -309,3 +329,26 @@ algorithms#(basis, List, List, Module) = new MutableHashTable from {
 -- Installing hooks for resolution
 scan({Default, Torsion}, strategy ->
     addHook(key := (basis, List, List, Module), algorithms#key#strategy, Strategy => strategy))
+
+-----------------------------------------------------------------------------
+-- part
+-----------------------------------------------------------------------------
+-- returns the graded component of an object in a specific degree
+
+-- gives a map back to the coefficient ring
+residueMap = (cacheValue symbol residueMap) (R -> map(A := coefficientRing R, R, DegreeMap => d -> take(d, - degreeLength A)))
+
+-- TODO: document these
+part(ZZ,   Ring) :=
+part(List, Ring) := Module => (deg, R) -> (residueMap R) source basis(deg, deg, module R)
+
+part(ZZ, Ideal)  :=
+part(ZZ, Module) :=
+part(List, Ideal)  :=
+part(List, Module) := Module => (deg, M) -> (residueMap ring M) source basis(deg, deg, M)
+
+part(ZZ,   Matrix) :=
+part(List, Matrix) := Matrix => (deg, f) -> (residueMap ring f) cover  basis(deg, f)
+
+-- TODO: currently doesn't work unless source f === target f
+part(Sequence, RingMap) := Matrix => (degs, f) -> (residueMap target f) cover basis(degs, f)

M2/Macaulay2/m2/chaincomplexes.m2

@@ -14,6 +14,17 @@ needs "Hom.m2"
 union := (x,y) -> keys(set x + set y)
 intersection := (x,y) -> keys(set x * set y)
 
+-----------------------------------------------------------------------------
+-- Categories for which complexes are currently implemented in Macaulay2
+-----------------------------------------------------------------------------
+-- SheafMap and CoherentSheaf are added later in Varieties
+
+isMorphism Thing    := x -> false
+isMorphism Matrix   := f -> true
+
+isAbelianCategory Thing         := x -> false
+isAbelianCategory Module        := M -> true
+
 -----------------------------------------------------------------------------
 -- ChainComplex, ChainComplexMap, and Resolution type declarations and basic constructors
 -----------------------------------------------------------------------------

M2/Macaulay2/m2/classes.m2

@@ -10,6 +10,19 @@ ancestors' = T -> unique join({T}, while (T = class  T) =!= Type  list T, {Type}
 
 -- TODO: make this TT toString X later?
 synonym = X -> if X.?synonym then X.synonym else "object of class " | toString X
+-- TODO: find a more permanent solution
+plurals = new MutableHashTable from {
+    "body"       => "bodies",
+    "dictionary" => "dictionaries",
+    "matrix"     => "matrices",
+    "sheaf"      => "sheaves",
+    "variety"    => "varieties",
+    }
+pluralize = s -> demark_" " append(
+    drop(ws := separate_" " s, -1),
+    if  plurals#?(last ws)
+    then plurals#(last ws) else last ws | "s")
+pluralsynonym = T -> if T.?synonym then pluralize T.synonym else "objects of class " | toString T
 
 Time.synonym = "timing result"
 Boolean.synonym = "Boolean value"

M2/Macaulay2/m2/matrix.m2

@@ -245,6 +245,10 @@ isHomogeneous Matrix := (cacheValue symbol isHomogeneous) ( m -> ( isHomogeneous
 
 isWellDefined Matrix := f -> matrix f * presentation source f % presentation target f == 0
 
+-----------------------------------------------------------------------------
+-- directSum and friends
+-----------------------------------------------------------------------------
+
 ggConcatCols := (tar,src,mats) -> (
      map(tar,src,if mats#0 .?RingMap then mats#0 .RingMap,rawConcatColumns (raw\mats),Degree => if same(degree \ mats) then degree mats#0)
      )
@@ -318,28 +322,29 @@ indices HashTable := X -> (
      else error "expected an object with components"
      )
 
+-- This is used for methods like directSum or pullback/pushout
+-- which accept an arbitrary list of objects of the same type.
+applyUniformMethod = (symb, name) -> args -> (
+    if #args === 0 then error("expected at least one argument for ", name);
+    type := if uniform args then class args#0 else error("expected uniform objects for ", name);
+    meth := lookup(symb, type) ?? error("no method for ", name, " of ", pluralsynonym type);
+    if (Y := youngest args) =!= null and Y.?cache
+    then Y.cache#(symb, args) ??= meth args else meth args)
+
 directSum List := args -> directSum toSequence args
-directSum Sequence := args -> (
-    if #args === 0 then error "expected more than 0 arguments";
-    type := if uniform args then class args#0 else error "incompatible objects in direct sum";
-    meth := lookup(symbol directSum, type);
-    if meth === null then error "no method for direct sum";
-    Y := youngest args;
-    key := (directSum, args);
-    if Y === null or not Y.?cache then meth args else
-    if Y.cache#?key then Y.cache#key else Y.cache#key = meth args)
+directSum Sequence := applyUniformMethod(symbol directSum, "direct sum")
 
 -- Number.directSum = v -> directSum apply(v, a -> matrix{{a}})
 
 Option ++ Option := directSum
 directSum Option := o -> directSum(1 : o)
 Option.directSum = args -> (
-     if #args === 0 then error "expected more than 0 arguments";
+     if #args === 0 then error "expected at least one argument";
      objects := apply(args,last);
      labels  := toList args/first;
      type := if uniform objects then class objects#0 else error "incompatible objects in direct sum";
      meth := lookup(symbol directSum, type);
-     if meth === null then error "no method for direct sum";
+     if meth === null then error("no method for direct sum of ", pluralsynonym type);
      M := meth objects;
      M.cache.indices = labels;
      ic := M.cache.indexComponents = new HashTable from apply(#labels, i -> labels#i => i);

M2/Macaulay2/m2/matrix1.m2

@@ -4,9 +4,6 @@ needs "matrix.m2"
 needs "modules.m2"
 needs "quotient.m2"
 
-plurals := hashTable { "matrix" => "matrices" }
-pluralize := s -> if plurals#?s then plurals#s else s|"s"
-pluralsynonym := T -> try pluralize T.synonym else "objects of class "|toString T;
 notsamering := (X,Y) -> (
      if X === Y then error("expected ",pluralsynonym X, " for the same ring")
      else error("expected ",X.synonym," and ",Y.synonym," for the same ring"))
@@ -178,25 +175,22 @@ concatBlocks = mats -> (
 Matrix.matrix = options -> (f) -> concatBlocks f
 
 commonRing = method()
-trythese := rings -> ring try sum apply(rings, R -> 0_R) else error "common ring not found"
-commonRing List := f -> (
-     f = flatten f;
-     if #f == 0 then return ZZ;
-     types := unique apply(f, class);
-     if # types === 1 and instance(types#0,Ring) then return types#0;
-     rings := if all(types, T -> instance(T,Ring)) then types else unique apply(f,ring);
-     ring try sum apply(rings, R -> 0_R) else error "common ring not found")
-
-matrixTable := opts -> (f) -> (
+commonRing List := L -> (
+    rings := unique apply(flatten L, ring);
+    if #rings === 0 then ZZ      else
+    if #rings === 1 then rings#0 else
+    ring try sum apply(rings, R -> 0_R) else error "common ring not found")
+
+matrixTable = opts -> (f) -> (
      R := commonRing f;
      havemat := false;
-     f = applyTable(f, x -> (
-	       if instance(x,Matrix) then havemat = true
-	       else if not ( instance(x,RingElement) or instance(x,Number) )
-	       then error "expected numbers, ring elements, and matrices";
-	       promote(x,R)));
-     if not havemat then return map(R^#f,, f, opts);
+     f = applyTable(f,
+	 x -> if isMorphism x then (havemat = true; promote(x, R))
+	 else if instance(x, RingElement) or instance(x, Number) then promote(x, R)
+	 else error "expected numbers, ring elements, and matrices");
+     if not havemat then return map(R^#f, , f, opts);
      types := unique apply(flatten f, class);
+     -- Note: we use Matrix.matrix here, which is different from Matrix#matrix
      if # types === 1 and types#0 .?matrix then return ( types#0 .matrix opts)(f);
      f = apply(f, row -> new MutableList from row);
      m := #f;
@@ -274,20 +268,16 @@ matrix(Matrix) := Matrix => opts -> (m) -> (
 
 matrix RingElement := matrix Number := opts -> r -> matrix({{r}}, opts)
 
-matrix(List) := Matrix => opts -> (m) -> (
-     if #m === 0 then return matrix(ZZ, {});
-     mm := apply(splice m,splice);
-     if #mm === 0 then error "expected nonempty list";
-     types := unique apply(mm,class);
-     if #types === 1 then (
-	  type := types#0;
-	  if instance(type,Module) then matrix { apply(mm, v -> new Matrix from v) }
-	  else if ancestor(List, type) then (
-	       if isTable mm then (matrixTable opts)(mm)
-	       else error "expected rows all to be the same length"
-	       )
-	  else error "expected a table of ring elements or matrices, or a list of elements of the same module")
-     else error "expected a table of ring elements or matrices, or a list of elements of the same module")
+matrix List := Matrix => opts -> L -> (
+    if #L === 0 then return matrix(ZZ, {});
+    L = apply(splice L, splice); -- TODO: is this deepSplice?
+    if not uniform L and not isTable L
+    then error "expected a list of vectors or a table of entries or maps";
+    -- construct a matrix by concatenating column vectors
+    if instance(L#0, Vector) then matrix { apply(L, matrix) } else
+    -- construct a matrix from a table of entries or maps
+    if isTable L then (matrixTable opts)(L)
+    else error "expected a rows all to be the same length")
 
 align := g -> (
      -- generator and relation maps can just as well have a nonzero degree

M2/Macaulay2/m2/matrix2.m2

@@ -325,32 +325,40 @@ quotientRemainder(Matrix,Matrix) := Matrix => (f,g) -> (
 	  map(M, L, rem)
      ))
 
-leftQuotientWarn = true
-
-Matrix // Matrix := Matrix => (f,g) -> quotient(f,g)
-Matrix \\ Matrix := (g,f) -> (if leftQuotientWarn then (leftQuotientWarn = false; printerr "Warning: 'm \\\\ n' is deprecated; use 'n // m' instead."); f // g)
-quotient'(Matrix,Matrix) := Matrix => (f,g) -> (
-     if not isFreeModule source f or not isFreeModule target f
-     or not isFreeModule source g or not isFreeModule target g then error "expected maps between free modules";
-     dual quotient(dual f, dual g))
-quotient(Matrix,Matrix) := Matrix => opts -> (f,g) -> (
-     if target f != target g then error "quotient: expected maps with the same target";
-     c := runHooks((quotient, Matrix, Matrix), (opts, f, g));
-     if c =!= null then c else error "quotient: no method implemented for this type of input")
-
-addHook((quotient, Matrix, Matrix), Strategy => Default, (opts, f, g) -> (
+--------------------------------------------------------------------------------
+-- factoring of a matrix through another
+
+-- Note: when f and g are endomorphisms, the sources and targets all agree,
+-- so we need both functions quotient and quotient' to distinguish them.
+
+-- factor matrices with same targets
+Matrix // Matrix := Matrix => (f, g) -> quotient(f, g)
+Number // Matrix := RingElement // Matrix := Matrix => (r, g) -> map(source g, target g, map(target g, cover target g, r) // g)
+Matrix // Number := Matrix // RingElement := Matrix => (f, r) -> map(target f, source f, f // map(target f, cover target f, r))
+
+-- factor matrices with same sources
+Matrix \\ Matrix := Matrix => (g, f) -> quotient'(f, g)
+-- commented because they don't seem very meaningful
+--Matrix \\ Number := Matrix \\ RingElement := Matrix => (g, r) -> map(source g, target g, g \\ map(cover source g, source g, r))
+--Number \\ Matrix := RingElement \\ Matrix := Matrix => (r, f) -> map(target g, source g, map(cover source g, source g, r) \\ f)
+
+quotient(Matrix, Matrix) := Matrix => opts -> (f, g) -> (
+    -- given f: A-->C and g: B-->C, then find (f//g): A-->B such that g o (f//g) + r = f
+    if target f != target g then error "quotient: expected maps with the same target";
+    c := runHooks((quotient, Matrix, Matrix), (opts, f, g), Strategy => opts.Strategy);
+    if c =!= null then c else error "quotient: no method implemented for this type of input")
+
+addHook((quotient, Matrix, Matrix), Strategy => Default,
+    -- Note: this strategy only works if the remainder is zero, i.e.:
+    -- homomorphism' f % image Hom(source f, g) == 0
+    -- TODO: should we pass MinimalGenerators => false to Hom and homomorphism'?
+    (opts, f, g) -> map(source g, source f, homomorphism(homomorphism' f // Hom(source f, g))))
+
+-- FIXME: this is still causing unreasonable slow downs, e.g. for (large m) // (scalar)
+addHook((quotient, Matrix, Matrix), Strategy => "Reflexive", (opts, f, g) -> if f == 0 or isFreeModule source f then (
      L := source f;	     -- result may not be well-defined if L is not free
      M := target f;
      N := source g;
-     if instance(ring M, InexactField)
-       and numRows g === numColumns g
-       and isFreeModule source g and isFreeModule source f
-       then return solve(g,f);	   	     	       	    
-     local retVal;
-     if isQuotientOf(ZZ,ring target f)
-       and isFreeModule source g and isFreeModule source f
-       then retVal = solve(g,f);
-     if retVal =!= null then return retVal;
      if M.?generators then (
 	  M = cokernel presentation M;	    -- this doesn't change the cover
 	  );
@@ -368,16 +376,36 @@ addHook((quotient, Matrix, Matrix), Strategy => Default, (opts, f, g) -> (
 	  Degree => degree f' - degree g'  -- set the degree in the engine instead
 	  )))
 
-Number // Matrix :=
-RingElement // Matrix := (r,f) -> (r * id_(target f)) // f
-Matrix \\ Number :=
-Matrix \\ RingElement := (f,r) -> r // f
-
-Matrix // Number :=
-Matrix // RingElement := (f,r) -> f // (r * id_(target f))
-Number \\ Matrix :=
-RingElement \\ Matrix      := (r,f) -> f // r
+addHook((quotient, Matrix, Matrix), Strategy => InexactField, (opts, f, g) ->
+    if instance(ring target f, InexactField)
+    -- TODO: c.f. https://github.com/Macaulay2/M2/issues/3252
+    and numRows g === numColumns g
+    and isFreeModule source g
+    and isFreeModule source f
+    then solve(g, f))
+
+addHook((quotient, Matrix, Matrix), Strategy => ZZ, (opts, f, g) ->
+    if isQuotientOf(ZZ, ring target f)
+    and isFreeModule source g
+    and isFreeModule source f
+    then solve(g, f))
+
+quotient'(Matrix, Matrix) := Matrix => opts -> (f, g) -> (
+    -- given f: A-->C and g: A-->B, then finds (g\\f): B-->C such that (g\\f) o g + r = f
+    if source f != source g then error "quotient': expected maps with the same source";
+    c := runHooks((quotient', Matrix, Matrix), (opts, f, g), Strategy => opts.Strategy);
+    if c =!= null then c else error "quotient': no method implemented for this type of input")
+
+addHook((quotient', Matrix, Matrix), Strategy => Default,
+    -- Note: this strategy only works if the remainder is zero, i.e.:
+    -- homomorphism' f % image Hom(g, target f) == 0
+    -- TODO: should we pass MinimalGenerators => false to Hom and homomorphism'?
+    (opts, f, g) -> map(target f, target g, homomorphism(homomorphism' f // Hom(g, target f))))
+
+addHook((quotient', Matrix, Matrix), Strategy => "Reflexive",
+    (opts, f, g) -> if all({source f, target f, source g, target g}, isFreeModule) then dual quotient(dual f, dual g, opts))
 
+--------------------
 
 remainder'(Matrix,Matrix) := Matrix => (f,g) -> (
      if not isFreeModule source f or not isFreeModule source g

M2/Macaulay2/m2/methods.m2

@@ -271,7 +271,7 @@ setupMethods((), {
 	  koszul, target, source,
 	  getChangeMatrix, cover, coverMap, super, terms,
 	  cokernel, coimage, comodule, image, someTerms, scanKeys, scanValues,
-	  substitute, complete, ambient, remainder, quotientRemainder, remainder', quotientRemainder', quotient',
+	  substitute, complete, ambient, remainder, quotientRemainder, remainder', quotientRemainder',
 	  coefficients, monomials, size, sum, product, nullhomotopy, module, raw,
 	  content, leadTerm, leadCoefficient, leadMonomial, components,
 	  assign, realPart, imaginaryPart, conjugate,