-
Notifications
You must be signed in to change notification settings - Fork 10.8k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[flang][runtime] Better real MOD/MODULO results #77167
Conversation
@llvm/pr-subscribers-flang-semantics @llvm/pr-subscribers-flang-runtime Author: Peter Klausler (klausler) ChangesThe Fortran standard defines real MOD and MODULO with expressions like MOD(a,p) = a - AINT(a/p)*p. Unfortunately, these definitions have poor accuracy when a is much larger in magnitude than p, and every Fortran compiler uses better algorithms instead. Fixes llvm-test-suite/Fortran/gfortran/regression/mod_large_1.f90. Full diff: https://github.com/llvm/llvm-project/pull/77167.diff 5 Files Affected:
diff --git a/flang/docs/Extensions.md b/flang/docs/Extensions.md
index 16eb67f2e27c81..e102ffcd389b81 100644
--- a/flang/docs/Extensions.md
+++ b/flang/docs/Extensions.md
@@ -96,6 +96,11 @@ end
* `NULL()` without `MOLD=` is not allowed to be associated as an
actual argument corresponding to an assumed-rank dummy argument;
its rank in the called procedure would not be well-defined.
+* The standard defines the intrinsic functions `MOD` and `MODULO`
+ for real arguments using expressions in terms of `AINT` and `FLOOR`.
+ These definitions yield fairly poor results due to floating-point
+ cancellation, and every Fortran compiler (including this one)
+ uses better algorithms.
## Extensions, deletions, and legacy features supported by default
diff --git a/flang/include/flang/Evaluate/real.h b/flang/include/flang/Evaluate/real.h
index a30d0dbfa8e1eb..62c99cebc31684 100644
--- a/flang/include/flang/Evaluate/real.h
+++ b/flang/include/flang/Evaluate/real.h
@@ -170,8 +170,8 @@ class Real : public common::RealDetails<PREC> {
// DIM(X,Y) = MAX(X-Y, 0)
ValueWithRealFlags<Real> DIM(const Real &,
Rounding rounding = TargetCharacteristics::defaultRounding) const;
- // MOD(x,y) = x - AINT(x/y)*y
- // MODULO(x,y) = x - FLOOR(x/y)*y
+ // MOD(x,y) = x - AINT(x/y)*y (in the standard)
+ // MODULO(x,y) = x - FLOOR(x/y)*y (in the standard)
ValueWithRealFlags<Real> MOD(const Real &,
Rounding rounding = TargetCharacteristics::defaultRounding) const;
ValueWithRealFlags<Real> MODULO(const Real &,
diff --git a/flang/lib/Evaluate/real.cpp b/flang/lib/Evaluate/real.cpp
index cb8be98bb6f40f..eb0848e50a58d1 100644
--- a/flang/lib/Evaluate/real.cpp
+++ b/flang/lib/Evaluate/real.cpp
@@ -401,48 +401,63 @@ ValueWithRealFlags<Real<W, P>> Real<W, P>::HYPOT(
return result;
}
-// MOD(x,y) = x - AINT(x/y)*y
+// MOD(x,y) = x - AINT(x/y)*y in the standard; unfortunately, this definition
+// can be pretty inaccurate when x is much larger than y in magnitude due to
+// cancellation. Implement instead with (essentially) arbitrary precision
+// long division, discarding the quotient and returning the remainder.
+// See runtime/numeric.cpp for more details.
template <typename W, int P>
ValueWithRealFlags<Real<W, P>> Real<W, P>::MOD(
- const Real &y, Rounding rounding) const {
+ const Real &p, Rounding rounding) const {
ValueWithRealFlags<Real> result;
- auto quotient{Divide(y, rounding)};
- if (quotient.value.IsInfinite() && IsFinite() && y.IsFinite() &&
- !y.IsZero()) {
- // x/y overflowed -- so it must be an integer in this representation and
- // the result must be a zero.
+ if (IsNotANumber() || p.IsNotANumber() || IsInfinite()) {
+ result.flags.set(RealFlag::InvalidArgument);
+ result.value = NotANumber();
+ } else if (p.IsZero()) {
+ result.flags.set(RealFlag::DivideByZero);
+ result.value = NotANumber();
+ } else if (p.IsInfinite()) {
+ result.value = *this;
+ } else {
+ result.value = ABS();
+ auto pAbs{p.ABS()};
+ Real half, adj;
+ half.Normalize(false, exponentBias - 1, Fraction::MASKL(1)); // 0.5
+ for (adj.Normalize(false, Exponent(), pAbs.GetFraction());
+ result.value.Compare(pAbs) != Relation::Less;
+ adj = adj.Multiply(half).value) {
+ if (result.value.Compare(adj) != Relation::Less) {
+ Real adjusted{
+ result.value.Subtract(adj, rounding).AccumulateFlags(result.flags)};
+ if (adjusted.Compare(result.value) == Relation::Equal) {
+ break;
+ }
+ result.value = adjusted;
+ if (adjusted.IsZero()) {
+ break;
+ }
+ }
+ }
if (IsNegative()) {
- result.value = Real{}.Negate(); // -0.
+ result.value = result.value.Negate();
}
- } else {
- Real toInt{quotient.AccumulateFlags(result.flags)
- .ToWholeNumber(common::RoundingMode::ToZero)
- .AccumulateFlags(result.flags)};
- Real product{toInt.Multiply(y, rounding).AccumulateFlags(result.flags)};
- result.value = Subtract(product, rounding).AccumulateFlags(result.flags);
}
return result;
}
-// MODULO(x,y) = x - FLOOR(x/y)*y
+// MODULO(x,y) = x - FLOOR(x/y)*y in the standard; here, it is defined
+// in terms of MOD() with adjustment of the result.
template <typename W, int P>
ValueWithRealFlags<Real<W, P>> Real<W, P>::MODULO(
- const Real &y, Rounding rounding) const {
- ValueWithRealFlags<Real> result;
- auto quotient{Divide(y, rounding)};
- if (quotient.value.IsInfinite() && IsFinite() && y.IsFinite() &&
- !y.IsZero()) {
- // x/y overflowed -- so it must be an integer in this representation and
- // the result must be a zero.
- if (y.IsNegative()) {
- result.value = Real{}.Negate(); // -0.
+ const Real &p, Rounding rounding) const {
+ ValueWithRealFlags<Real> result{MOD(p, rounding)};
+ if (IsNegative() != p.IsNegative()) {
+ if (result.value.IsZero()) {
+ result.value = result.value.Negate();
+ } else {
+ result.value =
+ result.value.Add(p, rounding).AccumulateFlags(result.flags);
}
- } else {
- Real toInt{quotient.AccumulateFlags(result.flags)
- .ToWholeNumber(common::RoundingMode::Down)
- .AccumulateFlags(result.flags)};
- Real product{toInt.Multiply(y, rounding).AccumulateFlags(result.flags)};
- result.value = Subtract(product, rounding).AccumulateFlags(result.flags);
}
return result;
}
diff --git a/flang/runtime/numeric.cpp b/flang/runtime/numeric.cpp
index 6cbf00e0c36c78..b7fb04ace82f0b 100644
--- a/flang/runtime/numeric.cpp
+++ b/flang/runtime/numeric.cpp
@@ -101,6 +101,25 @@ template <typename T> inline RT_API_ATTRS T Fraction(T x) {
RT_DIAG_POP
+// SET_EXPONENT (16.9.171)
+template <typename T> inline RT_API_ATTRS T SetExponent(T x, std::int64_t p) {
+ if (std::isnan(x)) {
+ return x; // NaN -> same NaN
+ } else if (std::isinf(x)) {
+ return std::numeric_limits<T>::quiet_NaN(); // +/-Inf -> NaN
+ } else if (x == 0) {
+ return x; // return negative zero if x is negative zero
+ } else {
+ int expo{std::ilogb(x) + 1};
+ auto ip{static_cast<int>(p - expo)};
+ if (ip != p - expo) {
+ ip = p < 0 ? std::numeric_limits<int>::min()
+ : std::numeric_limits<int>::max();
+ }
+ return std::ldexp(x, ip); // x*2**(p-e)
+ }
+}
+
// MOD & MODULO (16.9.135, .136)
template <bool IS_MODULO, typename T>
inline RT_API_ATTRS T IntMod(T x, T p, const char *sourceFile, int sourceLine) {
@@ -121,14 +140,47 @@ inline RT_API_ATTRS T RealMod(
Terminator{sourceFile, sourceLine}.Crash(
IS_MODULO ? "MODULO with P==0" : "MOD with P==0");
}
- T quotient{a / p};
- if (std::isinf(quotient) && std::isfinite(a) && std::isfinite(p)) {
- // a/p overflowed -- so it must be an integer, and the result
- // must be a zero of the same sign as one of the operands.
- return std::copysign(T{}, IS_MODULO ? p : a);
+ if (std::isnan(a) || std::isnan(p) || std::isinf(a)) {
+ return std::numeric_limits<T>::quiet_NaN();
+ } else if (std::isinf(p)) {
+ return a;
+ } else {
+ // The standard defines MOD(a,p)=a-AINT(a/p)*p and
+ // MODULO(a,p)=a-FLOOR(a/p)*p, but those definitions lose
+ // precision badly due to cancellation when ABS(a) is
+ // much larger than ABS(p).
+ // Insights:
+ // - MOD(a,p)=MOD(a-n*p,p) when a>0, p>0, integer n>0, and a>=n*p
+ // - when n is a power of two, n*p is exact
+ // - as a>=n*p, a-n*p does not round.
+ // So repeatedly reduce a by all n*p in decreasing order of n;
+ // what's left is the desired remainder. This is basically
+ // the same algorithm as arbitrary precision binary long division,
+ // discarding the quotient.
+ T tmp{std::abs(a)};
+ T pAbs{std::abs(p)};
+ for (T adj{SetExponent(pAbs, Exponent<int>(tmp))}; tmp >= pAbs; adj /= 2) {
+ if (tmp >= adj) {
+ T adjusted{tmp - adj};
+ if (adjusted == tmp) {
+ break;
+ }
+ tmp = adjusted;
+ if (tmp == 0) {
+ break;
+ }
+ }
+ }
+ if (a < 0) {
+ tmp = -tmp;
+ }
+ if constexpr (IS_MODULO) {
+ if ((a < 0) != (p < 0)) {
+ tmp += p;
+ }
+ }
+ return tmp;
}
- T toInt{IS_MODULO ? std::floor(quotient) : std::trunc(quotient)};
- return a - toInt * p;
}
// RRSPACING (16.9.164)
@@ -222,25 +274,6 @@ inline RT_API_ATTRS CppTypeFor<TypeCategory::Integer, 4> SelectedRealKind(
return error ? error : kind;
}
-// SET_EXPONENT (16.9.171)
-template <typename T> inline RT_API_ATTRS T SetExponent(T x, std::int64_t p) {
- if (std::isnan(x)) {
- return x; // NaN -> same NaN
- } else if (std::isinf(x)) {
- return std::numeric_limits<T>::quiet_NaN(); // +/-Inf -> NaN
- } else if (x == 0) {
- return x; // return negative zero if x is negative zero
- } else {
- int expo{std::ilogb(x) + 1};
- auto ip{static_cast<int>(p - expo)};
- if (ip != p - expo) {
- ip = p < 0 ? std::numeric_limits<int>::min()
- : std::numeric_limits<int>::max();
- }
- return std::ldexp(x, ip); // x*2**(p-e)
- }
-}
-
// SPACING (16.9.180)
template <int PREC, typename T> inline RT_API_ATTRS T Spacing(T x) {
if (std::isnan(x)) {
diff --git a/flang/test/Evaluate/fold-mod.f90 b/flang/test/Evaluate/fold-mod.f90
index 96c35cb6e4708a..fa9132300e251a 100644
--- a/flang/test/Evaluate/fold-mod.f90
+++ b/flang/test/Evaluate/fold-mod.f90
@@ -39,4 +39,8 @@ module m1
logical, parameter :: test_modulo_r12b = sign(1., modulo(huge(0.), -tiny(0.))) == -1.
logical, parameter :: test_modulo_r13a = modulo(huge(0.), tiny(0.)) == 0.
logical, parameter :: test_modulo_r13b = sign(1., modulo(-huge(0.), -tiny(0.))) == -1.
+
+ logical, parameter :: test_mod_r14 = mod(1e22, 1.7) == .99592876
+ logical, parameter :: test_modulo_r14 = modulo(1e22, -1.7) == -.7040713
+
end module
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
LGTM
The Fortran standard defines real MOD and MODULO with expressions like MOD(a,p) = a - AINT(a/p)*p. Unfortunately, these definitions have poor accuracy when a is much larger in magnitude than p, and every Fortran compiler uses better algorithms instead. Fixes llvm-test-suite/Fortran/gfortran/regression/mod_large_1.f90.
The Fortran standard defines real MOD and MODULO with expressions like MOD(a,p) = a - AINT(a/p)*p. Unfortunately, these definitions have poor accuracy when a is much larger in magnitude than p, and every Fortran compiler uses better algorithms instead. Fixes llvm-test-suite/Fortran/gfortran/regression/mod_large_1.f90.
The Fortran standard defines real MOD and MODULO with expressions like MOD(a,p) = a - AINT(a/p)*p. Unfortunately, these definitions have poor accuracy when a is much larger in magnitude than p, and every Fortran compiler uses better algorithms instead.
Fixes llvm-test-suite/Fortran/gfortran/regression/mod_large_1.f90.