[flang-commits] [flang] [flang][runtime] Better real MOD/MODULO results (PR #77167)
Peter Klausler via flang-commits
flang-commits at lists.llvm.org
Fri Jan 5 17:19:54 PST 2024
https://github.com/klausler created https://github.com/llvm/llvm-project/pull/77167
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.
>From 0c45846308234de6ef264088b49eadc654823a6b Mon Sep 17 00:00:00 2001
From: Peter Klausler <pklausler at nvidia.com>
Date: Fri, 5 Jan 2024 17:16:12 -0800
Subject: [PATCH] [flang][runtime] Better real MOD/MODULO results
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.
---
flang/docs/Extensions.md | 5 ++
flang/include/flang/Evaluate/real.h | 4 +-
flang/lib/Evaluate/real.cpp | 75 +++++++++++++++----------
flang/runtime/numeric.cpp | 85 ++++++++++++++++++++---------
flang/test/Evaluate/fold-mod.f90 | 4 ++
5 files changed, 115 insertions(+), 58 deletions(-)
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
More information about the flang-commits
mailing list