Skip to content
Open
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 51 additions & 16 deletions M2/Macaulay2/e/rings/frac.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include "groebner-computations/gbring.hpp"
#include "ring-elements/ring-element.hpp"
#include "rings/polyring.hpp"
#include "rings/polyquotient.hpp"
#include "exceptions.hpp"

#define FRAC_VAL(f) (reinterpret_cast<frac_elem *>((f).poly_val))
Expand Down Expand Up @@ -99,6 +100,30 @@ 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;
if (dynamic_cast<const PolyRingQuotient *>(R_) != nullptr)
Comment thread
pzinn marked this conversation as resolved.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think R_->is_quotient_ring() would be clearer

{
denom_inverse = R_->invert(f->denom);
if (R_->is_zero(denom_inverse)) return false;
}
else
{
if (!R_->is_unit(f->denom)) return false;
denom_inverse = R_->invert(f->denom);
}

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
{
Expand All @@ -108,15 +133,17 @@ 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;
Comment thread
pzinn marked this conversation as resolved.
x = f->numer;
const RingElement *a = RingElement::make_raw(R_, x);
const RingElement *b = RingElement::make_raw(R_, y);
const RingElement *c = rawGCDRingElement(a, b, nullptr, false);
if (!c) return;
if (!R_->is_equal(y, R_->one()))

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If $y = 1$, then we would have already returned a few lines ago (since simplify_unit_denominator(f) would have returned true, so this if statement isn't necessary. (And the same comment for the local ring version below.)

{
x = f->numer;
const RingElement *a = RingElement::make_raw(R_, x);
const RingElement *b = RingElement::make_raw(R_, y);
const RingElement *c = rawGCDRingElement(a, b, nullptr, false);
if (!c) return;

#if 0
// Debugging code
Expand All @@ -131,10 +158,11 @@ void FractionField::simplify(frac_elem *f) const
o << newline;
emit(o.str());
#endif
if (!R_->is_equal(c->get_value(), R_->one()))
{
f->numer = R_->divide(f->numer, c->get_value());
f->denom = R_->divide(f->denom, c->get_value());
if (!R_->is_equal(c->get_value(), R_->one()))
{
f->numer = R_->divide(f->numer, c->get_value());
f->denom = R_->divide(f->denom, c->get_value());
}
}
// Now, let's take the content of the denominator, and divide the
// numerator
Expand All @@ -159,6 +187,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
{
Expand All @@ -176,6 +205,7 @@ void FractionField::simplify(frac_elem *f) const
R_->remove(f->denom);
f->numer = y;
f->denom = x;
simplify_unit_denominator(f);
}
}

Expand Down Expand Up @@ -267,12 +297,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
Expand Down Expand Up @@ -340,6 +369,12 @@ bool FractionField::lift(const Ring *Rg,
result = R_->copy(h->numer);
return true;
}
else if (R_->is_unit(h->denom))
{
ring_elem hinv = R_->invert(h->denom);
result = R_->mult(hinv, h->numer);
return true;
}
else
{
if (R_->is_field())
Expand Down
1 change: 1 addition & 0 deletions M2/Macaulay2/e/rings/frac.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
34 changes: 28 additions & 6 deletions M2/Macaulay2/e/rings/localring.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -73,6 +87,7 @@ 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;
Expand Down Expand Up @@ -124,6 +139,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
{
Expand All @@ -141,6 +157,7 @@ void LocalRing::simplify(local_elem *f) const
mRing->remove(f->denom);
f->numer = y;
f->denom = x;
simplify_unit_denominator(f);
}
}

Expand Down Expand Up @@ -280,6 +297,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())
Expand Down Expand Up @@ -317,12 +340,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
Expand Down
1 change: 1 addition & 0 deletions M2/Macaulay2/e/rings/localring.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
7 changes: 7 additions & 0 deletions M2/Macaulay2/tests/normal/frac.m2
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,13 @@ 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)
-----------------------------------------------------------------------------
Comment thread
pzinn marked this conversation as resolved.
A = ZZ/101[a,b]/(a*b)
L = toField A
assert try (1/(a+b); false) else true
Expand Down
7 changes: 7 additions & 0 deletions M2/Macaulay2/tests/normal/matrix-equality.m2
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
Loading