"fraction" fixes

M2/Macaulay2/e/interface/factory.cpp

@@ -48,20 +48,22 @@ enum factoryCoeffMode {
   modeGF,
   modeUnknown
 };
-static enum factoryCoeffMode coeffMode(const PolynomialRing *P)
+static enum factoryCoeffMode coeffMode(const PolynomialRing *P,
+                                       bool reportErrors = true)
 {
   const Ring *F = P->getCoefficientRing();
   // if (F->cast_to_QQ()) return modeQQ;
   if (F->is_QQ()) return modeQQ;
   if (F->cast_to_RingZZ()) return modeZZ;
   // factory will abort if the characteristic is too large
   if (F->characteristic() > 536870909) {
-    ERROR("characteristic is too large (max is 2^29)");
+    if (reportErrors) ERROR("characteristic is too large (max is 2^29)");
     return modeError;
   }
   if (F->isFinitePrimeField()) return modeZn;
   if (F->isGaloisField()) return modeGF;
-  ERROR("expected coefficient ring of the form ZZ/n, ZZ, QQ, or GF");
+  if (reportErrors)
+    ERROR("expected coefficient ring of the form ZZ/n, ZZ, QQ, or GF");
   return modeError;
 }
 
@@ -478,11 +480,10 @@ void displayCF(const PolynomialRing *R, const CanonicalForm &h)  // for debuggin
   emit(o.str());
 }
 
-// TODO: figure out where this should be used
-bool factoryGoodRing(const PolynomialRing *P)
+// Check coefficient support without setting the engine error state.
+M2_bool factoryGoodRing(const PolynomialRing *P)
 {
-  struct enter_factory foo(P);
-  return foo.mode != modeError;
+  return coeffMode(P, false) != modeError;
 }
 
 const RingElement /* or null */ *rawGCDRingElement(const RingElement *f,

M2/Macaulay2/e/interface/factory.h

@@ -6,9 +6,11 @@
 // TODO: fix this
 #  if defined(__cplusplus)
 class Matrix;
+class PolynomialRing;
 class RingElement;
 #  else
 typedef struct Matrix Matrix;
+typedef struct PolynomialRing PolynomialRing;
 typedef struct RingElement RingElement;
 #  endif
 
@@ -30,6 +32,8 @@ const RingElement *rawExtendedGCDRingElement(const RingElement *f,
                                              const RingElement **A,
                                              const RingElement **B);
 
+M2_bool factoryGoodRing(const PolynomialRing *P);
+
 const RingElement *rawPseudoRemainder(const RingElement *f,
                                       const RingElement *g);
 

M2/Macaulay2/e/rings/frac.cpp

@@ -99,6 +99,31 @@ ring_elem FractionField::set_non_unit_frac(ring_elem top) const
   return zero();
 }
 
+bool FractionField::simplify_unit_denominator(frac_elem *f) const
+{
+  if (R_->is_equal(f->denom, R_->one())) return true;
+
+  ring_elem denom_inverse;
+  // TODO uniformise behaviour of invert for noninvertible elements
+  if (R_->is_quotient_ring())
+    {
+      denom_inverse = R_->invert(f->denom); // for quotient rings, don't call is_unit since it calls invert internally
+      if (R_->is_zero(denom_inverse)) return false; // for non invertible elements, returns zero
+    }
+  else
+    {
+      if (!R_->is_unit(f->denom)) return false; // for polynomial rings, test with is_unit first
+      denom_inverse = R_->invert(f->denom); // because invert throws an error for noninvertible elements
+    }
+
+  ring_elem numer = R_->mult(f->numer, denom_inverse);
+  R_->remove(f->numer);
+  R_->remove(f->denom);
+  f->numer = numer;
+  f->denom = R_->one();
+  return true;
+}
+
 ring_elem FractionField::fraction(const ring_elem top,
                                   const ring_elem bottom) const
 {
@@ -108,10 +133,10 @@ ring_elem FractionField::fraction(const ring_elem top,
 void FractionField::simplify(frac_elem *f) const
 {
   ring_elem x, y;
+  if (simplify_unit_denominator(f)) return;
   if (use_gcd_simplify)
     {
       y = f->denom;
-      if (R_->is_equal(y, R_->one())) return;
       x = f->numer;
       const RingElement *a = RingElement::make_raw(R_, x);
       const RingElement *b = RingElement::make_raw(R_, y);
@@ -159,6 +184,7 @@ void FractionField::simplify(frac_elem *f) const
           f->numer = R_->divide_by_given_content(f->numer, ct);
           f->denom = R_->divide_by_given_content(f->denom, ct);
         }
+      simplify_unit_denominator(f);
     }
   else
     {
@@ -176,6 +202,7 @@ void FractionField::simplify(frac_elem *f) const
       R_->remove(f->denom);
       f->numer = y;
       f->denom = x;
+      simplify_unit_denominator(f);
     }
 }
 
@@ -267,12 +294,11 @@ ring_elem FractionField::from_int(mpz_srcptr n) const
 
 bool FractionField::from_rational(mpq_srcptr n, ring_elem &result) const
 {
-  frac_elem *f = new_frac_elem();
-  f->numer = R_->from_int(mpq_numref(n));
-  f->denom = R_->from_int(mpq_denref(n));
-  bool ok = not R_->is_zero(f->denom);
-  if (ok) result = FRAC_RINGELEM(f);
-  return ok;
+  ring_elem numer = R_->from_int(mpq_numref(n));
+  ring_elem denom = R_->from_int(mpq_denref(n));
+  if (R_->is_zero(denom)) return false;
+  result = FRAC_RINGELEM(make_elem(numer, denom));
+  return true;
 }
 
 ring_elem FractionField::var(int v) const

M2/Macaulay2/e/rings/frac.hpp

@@ -22,6 +22,7 @@ class FractionField : public Ring
 
   ring_elem set_non_unit_frac(ring_elem top) const;
   frac_elem *new_frac_elem() const;
+  bool simplify_unit_denominator(frac_elem *f) const;
   void simplify(frac_elem *f) const;
   frac_elem *make_elem(ring_elem a, ring_elem b) const;
 

M2/Macaulay2/e/rings/localring.cpp

@@ -60,6 +60,20 @@ local_elem *LocalRing::make_elem(ring_elem a, ring_elem b) const
 
 local_elem *LocalRing::new_local_elem() const { return newitem(local_elem); }
 
+bool LocalRing::simplify_unit_denominator(local_elem *f) const
+{
+  if (mRing->is_equal(f->denom, mRing->one())) return true;
+  if (!mRing->is_unit(f->denom)) return false;
+
+  ring_elem denom_inverse = mRing->invert(f->denom);
+  ring_elem numer = mRing->mult(f->numer, denom_inverse);
+  mRing->remove(f->numer);
+  mRing->remove(f->denom);
+  f->numer = numer;
+  f->denom = mRing->one();
+  return true;
+}
+
 bool LocalRing::is_in_prime(const ring_elem f) const
 {
   MatrixConstructor mat(mRing->make_FreeModule(1), 1);
@@ -73,10 +87,10 @@ bool LocalRing::is_in_prime(const ring_elem f) const
 void LocalRing::simplify(local_elem *f) const
 {
   ring_elem x, y;
+  if (simplify_unit_denominator(f)) return;
   if (use_gcd_simplify)
     {
       y = f->denom;
-      if (mRing->is_equal(y, mRing->one())) return;
       x = f->numer;
       const RingElement *a = RingElement::make_raw(mRing, x);
       const RingElement *b = RingElement::make_raw(mRing, y);
@@ -124,6 +138,7 @@ void LocalRing::simplify(local_elem *f) const
           f->numer = mRing->divide_by_given_content(f->numer, ct);
           f->denom = mRing->divide_by_given_content(f->denom, ct);
         }
+      simplify_unit_denominator(f);
     }
   else
     {
@@ -141,6 +156,7 @@ void LocalRing::simplify(local_elem *f) const
       mRing->remove(f->denom);
       f->numer = y;
       f->denom = x;
+      simplify_unit_denominator(f);
     }
 }
 
@@ -280,6 +296,12 @@ bool LocalRing::lift(const Ring *Rg, const ring_elem f, ring_elem &result) const
           result = mRing->copy(h->numer);
           return true;
         }
+      else if (mRing->is_unit(h->denom))
+        {
+          ring_elem hinv = mRing->invert(h->denom);
+          result = mRing->mult(hinv, h->numer);
+          return true;
+        }
       else
         {
           if (mRing->is_field())
@@ -317,12 +339,11 @@ bool LocalRing::promote(const Ring *Rf,
 
 bool LocalRing::from_rational(mpq_srcptr n, ring_elem &result) const
 {
-  local_elem *f = new_local_elem();
-  f->numer = mRing->from_int(mpq_numref(n));
-  f->denom = mRing->from_int(mpq_denref(n));
-  bool ok = not mRing->is_zero(f->denom);
-  if (ok) result = ring_elem(f);
-  return ok;
+  ring_elem numer = mRing->from_int(mpq_numref(n));
+  ring_elem denom = mRing->from_int(mpq_denref(n));
+  if (mRing->is_zero(denom)) return false;
+  result = ring_elem(make_elem(numer, denom));
+  return true;
 }
 
 ring_elem LocalRing::from_long(long n) const

M2/Macaulay2/e/rings/localring.hpp

@@ -29,6 +29,7 @@ class LocalRing : public Ring
   local_elem *new_local_elem() const;
 
   bool is_in_prime(const ring_elem f) const;
+  bool simplify_unit_denominator(local_elem *f) const;
   void simplify(local_elem *f) const;
 
   // FIXME remove:

M2/Macaulay2/e/rings/polyquotient.cpp

@@ -162,21 +162,17 @@ ring_elem PolyRingQuotient::invert(const ring_elem f) const
 {
   if (nvars_ == 1 && n_quotients() == 1 && K_->is_field() && ! K_->is_fraction_field())
     {
+      // Avoid setting the engine error state for unsupported coefficient rings.
+      if (!factoryGoodRing(getAmbientRing())) return from_long(0);
+
       ring_elem g = quotient_element(0);
 
       RingElement *f1 = RingElement::make_raw(getAmbientRing(), f);
       RingElement *g1 = RingElement::make_raw(getAmbientRing(), g);
       const RingElement *u1;
       const RingElement *v1;
       const RingElement *ret = rawExtendedGCDRingElement(f1, g1, &u1, &v1);
-      if (ret == nullptr)
-        {
-          // one reason this might return nullptr is if the coefficient ring is not
-          // ZZ/n, ZZ, or QQ
-          // now what do we do?
-          // we can't return nullptr
-          INTERNAL_ERROR("ring element gcd computation failed");
-        }
+      if (ret == nullptr) return from_long(0);
       if (!getAmbientRing()->is_unit(ret->get_value())) return from_long(0);
       return u1->get_value();
     }

M2/Macaulay2/packages/LocalRings/tests.m2

@@ -872,6 +872,25 @@ TEST /// -- dim and char of a LocalRing reflect the underlying prime
   assert(degreeLength Rmax == 1)
 ///
 
+TEST /// -- normalize ambient-unit denominators
+  R = QQ[x]
+  RP = localRing(R, ideal x)
+  f = sub(1/2, RP)
+  assert(numerator f == 1/2_R)
+  assert(denominator f == 1_R)
+  assert(lift(f, R) == 1/2_R)
+///
+
+TEST /// -- promoting/lifting to/from fraction field
+  S = QQ[x]
+  p = ideal x
+  R = S_p
+  F = frac R
+  assert(promote(x_R, F) === x_F)
+  assert(lift(x_F, R) === x_R)
+  assert not liftable(1/x, R)
+///
+
 end--
 
 --============================ Tests Under Development ===================================--
@@ -957,17 +976,6 @@ end--
   gens gb Iloc
 ///
 
-TEST ///
--- promoting/lifting to/from fraction field
-S = QQ[x]
-p = ideal x
-R = S_p
-F = frac R
-assert(promote(x_R, F) === x_F)
-assert(lift(x_F, R) === x_R)
-assert not liftable(1/x, R)
-///
-
 end--
 
 --  Development stuff

M2/Macaulay2/tests/normal/frac.m2

@@ -80,6 +80,20 @@ assert try (a/b; false) else true
 getNonUnit B
 isField B
 -----------------------------------------------------------------------------
+R = QQ[a,b]/(a*b-1)
+F = frac R
+assert(numerator(1/b) == a)
+assert(denominator(1/b) == 1_R)
+assert(lift(1/b,R) == a)
+assert(lift((b+1)/b,R) == a+1)
+-----------------------------------------------------------------------------
+K = toField(QQ[a]/(a^2-2))
+B = K[t]/(t^2-a)
+F = frac B
+f = 1/(t+1)
+assert(numerator f == 1_B)
+assert(denominator f == t+1)
+-----------------------------------------------------------------------------
 A = ZZ/101[a,b]/(a*b)
 L = toField A
 assert try (1/(a+b); false) else true

M2/Macaulay2/tests/normal/matrix-equality.m2

@@ -34,6 +34,13 @@ debug Core
   assert isWellDefined f
   assert isWellDefined g
   assert(f == g)  -- fails in version 1.13
+
+  -- issue #4461 fix
+  F = frac(QQ[y])
+  a = matrix{{1/2}} * matrix{{y}}
+  b = matrix{{y/2}}
+  assert(a_(0,0) === b_(0,0))
+  assert(a == b)
 end--
 
 -- Local Variables:

M2/Macaulay2/tests/normal/toField.m2

@@ -8,6 +8,11 @@ assert( dim I == 2 )					    -- fails in 1.2
 assert( dim Proj (R/I) == 1 )				    -- fails in 1.2
 assert( hilbertPolynomial I == - 15 * hilbertPolynomial (QQ[x]) + 6 * hilbertPolynomial (QQ[x,y]) ) -- fails in 1.2
 
+-- Extended GCD is unavailable over some nested algebraic coefficient fields.
+F = toField(QQ[i]/(i^2+1))
+K = toField(F[r]/(r^4-2))
+assert try (r^-1; false) else true
+
 
 A = ZZ[a]/(a^2+3);
 L = toField A