Git issues 3937 1935

M2/Macaulay2/e/hilb.cpp

@@ -709,33 +709,26 @@ int hilb_comp::coeff_of(const RingElement *h, int deg)
   // exp[0]=deg.
   const PolynomialRing *P = h->get_ring()->cast_to_PolynomialRing();
 
-  exponents_t exp = newarray_atomic(int, P->n_vars());
+  exponents_t exp = new int[P->n_vars()];
   int result = 0;
   for (Nterm& f : h->get_value())
     {
       P->getMonoid()->to_expvector(f.monom, exp);
       if (exp[0] < deg)
         {
-          ERROR("incorrect Hilbert function given");
-          fprintf(
-              stderr,
-              "internal error: incorrect Hilbert function given, aborting\n");
-          fprintf(
-              stderr,
-              "exp[0]: %d   deg: %d\n", exp[0], deg);    
-          abort();
+          throw exc::engine_error("incorrect Hilbert function given");
         }
       else if (exp[0] == deg)
         {
           std::pair<bool, long> res =
               P->getCoefficientRing()->coerceToLongInteger(f.coeff);
-          assert(res.first &&
-                 std::abs(res.second) < std::numeric_limits<int>::max());
+          if (not res.first or std::abs(res.second) > std::numeric_limits<int>::max())
+            throw exc::engine_error("Hilbert function value too large to use with Groebner basis computation");
           int n = static_cast<int>(res.second);
           result += n;
         }
     }
-  freemem(exp);
+  delete [] exp;
   return result;
 }
 

M2/Macaulay2/m2/exports.m2

@@ -529,6 +529,7 @@ export {
 	"cache",
 	"cacheValue",
 	"cancelTask",
+    "canUseHilbertHint",
 	"capture",
 	"catch",
 	"ceiling",

M2/Macaulay2/m2/gb.m2

@@ -257,32 +257,43 @@ gbGetPartialComputation := (m, type) -> (
 	null))
 
 -- handle the Hilbert numerator later, which might be here:
-checkHilbertHint = m -> (
-    R := ring m;
-    -- Needed for using Hilbert functions to aid in Groebner basis computation:
-    --    Ring is poly ring over a field (or skew commutative, or quotient ring of such, or both)
-    --    Ring is singly graded, every variable is positive
-    --    Ring is homogeneous in this grading
-    --    Matrix is homogeneous in this grading
-    isHomogeneous m
-    and degreeLength R === 1
+canUseHilbertHint = method()
+
+canUseHilbertHint Ring := Boolean => R -> (
+    degreeLength R === 1
     and (instance(R, PolynomialRing) or isQuotientOf(PolynomialRing, R))
     and isField coefficientRing R
     and (isCommutative R or isSkewCommutative R)
-    and all(degree \ generators(R, CoefficientRing => ZZ), deg -> deg#0 > 0)
+    and all(degree \ generators R, deg -> deg#0 > 0)
     )
+canUseHilbertHint Ideal :=
+canUseHilbertHint Module :=
+canUseHilbertHint Matrix := Boolean => m -> canUseHilbertHint ring m and isHomogeneous m
+
+-- checkHilbertHint = m -> (
+--     -- Needed for using Hilbert functions to aid in Groebner basis computation:
+--     --    Ring is poly ring over a field (or skew commutative, or quotient ring of such, or both)
+--     --    Ring is singly graded, every variable is positive
+--     --    Ring is homogeneous in this grading
+--     --    Matrix is homogeneous in this grading
+--     isHomogeneous m and canUseHilbertHint ring m
+--     )
 
 gbGetHilbertHint := (m, opts) -> (
     if opts.Hilbert =!= null then (
+        if not canUseHilbertHint m then (
+            --<< "warning: Hilbert hint given, but is being ignored" << endl;
+            return null;
+            );
 	if ring opts.Hilbert === degreesRing ring m then opts.Hilbert
 	else error "expected Hilbert option to be in the degrees ring of the ring of the matrix")
     else if m.cache.?cokernel and m.cache.cokernel.cache.?poincare
-    and   checkHilbertHint m then m.cache.cokernel.cache.poincare
+    and canUseHilbertHint m then m.cache.cokernel.cache.poincare
     else if m.?generators then (
 	g := m.generators;
 	if g.cache.?image then (
 	    M := g.cache.image;
-	    if M.cache.?poincare and checkHilbertHint m
+	    if M.cache.?poincare and canUseHilbertHint m
 	    then poincare target g - poincare M)))
 
 checkArgGB := m -> (

M2/Macaulay2/m2/ringmap.m2

@@ -304,7 +304,7 @@ algorithms#(kernel, RingMap) = new MutableHashTable from {
 	       graph := generators graphIdeal f;
 	       assert( not isHomogeneous f or isHomogeneous graph );
 	       SS := ring graph;
-	       chh := checkHilbertHint graph;
+	       chh := canUseHilbertHint graph;
 	       if chh then (
 		   -- compare with pushNonLinear
 		   hf := poincare module target f;

M2/Macaulay2/packages/Macaulay2Doc/functions.m2

@@ -21,6 +21,7 @@ load "./functions/Beta-doc.m2"
 load "./functions/betti-doc.m2"
 load "./functions/between-doc.m2"
 load "./functions/binomial-doc.m2"
+load "./functions/canUseHilbertHint-doc.m2"
 load "./functions/ceiling-doc.m2"
 load "./functions/changeBase-doc.m2"
 load "./functions/char-doc.m2"

M2/Macaulay2/packages/Macaulay2Doc/functions/canUseHilbertHint-doc.m2

@@ -0,0 +1,61 @@
+doc ///
+  Key
+    canUseHilbertHint
+    (canUseHilbertHint, Ring)
+    (canUseHilbertHint, Ideal)
+    (canUseHilbertHint, Module)
+    (canUseHilbertHint, Matrix)
+  Headline
+    whether certain Groebner computations can make use of the Hilbert function
+  Usage
+    canUseHilbertHint m
+  Inputs
+    m:{Ring,Ideal,Module,Matrix}
+  Outputs
+    :Boolean
+  Description
+    Text
+      Certain Groebner basis computations, when given homogeneous input, can run faster if
+      the Hilbert function is know: the algorithm can use this information to bypass
+      computing zero reductions in many cases.
+
+      In order to be able to use a Hilbert function hint, currently the object must be
+      in a singly graded polynomial ring (or quotient of such), over a field,
+      with all generators of positive degree, and it must be either commutative, or skew-commutative.
+
+      Given a ring, this function returns whether this holds for the ring.  For
+      a module, matrix or ideal, this function additionally checks whether the input is
+      homogeneous.
+    Example
+      R1 = ZZ/101[x,y,z, Degrees => {1,2,3}];
+      assert canUseHilbertHint R1
+      R2 = ZZ[x,y,z];
+      assert not canUseHilbertHint R2
+      R3 = (GF 8)[x,y,z]
+      assert canUseHilbertHint R3
+      use R1
+      R4 = R1[u,v, Join => false]/(x*u+y*v)
+      R5 = first flattenRing R4
+      canUseHilbertHint R4
+      canUseHilbertHint R5
+    Text
+      Here is one way to use Hilbert hints. If one knows the Hilbert function, one can use it.
+      If one computes a Hilbert function (via poincare), one can use that in a ring
+      which is identical, except the monomial order is more computationally intensive (e.g. Lex).
+    Example
+      R = ZZ/101[x,y,z,w]
+      I = ideal random(R^1, R^{-5,-5,-5})
+      hf = poincare I
+      codim I == 3 -- a complete intersection
+      Rlex = newRing(R, MonomialOrder => Lex)
+      Ilex = sub(I, Rlex)
+      elapsedTime g1 = gens gb Ilex;
+      Ilex = ideal(Ilex_*) -- clear out the previous Groebner basis
+      elapsedTime g2 = gens gb(Ilex, Hilbert => hf);
+      g1 == g2
+  Caveat
+    One cannot currently use Hilbert hints for multigraded input
+  SeeAlso
+    poincare
+    newRing
+///

M2/Macaulay2/packages/Macaulay2Doc/functions/gb-doc.m2

@@ -111,6 +111,7 @@ document {
 	 TO gbRemove,
 	 TO "gbTrace",
 	 TO LongPolynomial,
+         TO canUseHilbertHint
          }
      }
 

M2/Macaulay2/packages/MinimalPrimes/splitIdeals.m2

@@ -521,7 +521,7 @@ splitFunction#IndependentSet = (I,opts) -> (
     J := I.Ideal;
     if J == 1 then error "Internal error: Input should not be unit ideal.";
     R := ring J;
-    hf := if isHomogeneous J then poincare J else null;
+    hf := if canUseHilbertHint J then poincare J else null;
     indeps := independentSets(J, Limit=>1);
     basevars := support first indeps;
     if opts.Verbosity >= 3 then
@@ -540,7 +540,7 @@ splitFunction#IndependentSet = (I,opts) -> (
     -- otherwise compute over the fraction field.
     -- in the next test, I am pretty sure that whenever J is homogeneous, so is J (and hf is then defined).
     -- But I'm keeping the test here anyway
-    if isHomogeneous JS and hf =!= null then gb(JS, Hilbert=>hf) else gb JS;
+    if hf =!= null and canUseHilbertHint JS then gb(JS, Hilbert=>hf) else gb JS;
     (JSF, coeffs) := minimalizeOverFrac(JS, SF);
     if coeffs == {} then (
         I.IndependentSet = (basevars,S,SF);

M2/Macaulay2/tests/normal/gb-hilbert.m2

@@ -0,0 +1,27 @@
+-- test of canUseHilbertHint, Hilbert hint code with M2.
+restart
+R1 = ZZ/101[x,y,z, Degrees => {1,2,3}];
+assert canUseHilbertHint R1
+R2 = ZZ[x,y,z];
+assert not canUseHilbertHint R2
+R3 = (GF 8)[x,y,z]
+assert canUseHilbertHint R3
+use R1
+R4 = R1[u,v, Join => false]/(x*u+y*v)
+R5 = first flattenRing R4
+assert not canUseHilbertHint R4 -- we might want to change this behavior
+assert canUseHilbertHint R5
+
+-- test in presence of multi-degrees
+R = ZZ/101[a,b,c,d, Degrees => {2:{1,0}, 2:{0,1}}]
+I = ideal for i from 1 to 3 list random({2,3}, R)
+hf = poincare I
+assert not canUseHilbertHint R
+assert not canUseHilbertHint I
+
+Rlex = ZZ/101[a,b,c,d, Degrees => {2:{1,0}, 2:{0,1}}, MonomialOrder => Lex]
+Ilex = sub(I, Rlex)
+gblex = gens gb(Ilex, Hilbert => hf) -- gives warning, but no error
+assert(numgens ideal gblex == 37) -- happens with given random seed, at least....
+
+