[flang][runtime] Better real MOD/MODULO results (llvm#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.

Author: klausler

flang/docs/Extensions.md

@@ -102,6 +102,11 @@ end
   in F'2023 subclause 19.4 paragraphs 6 & 8 should apply.  Since this
   compiler properly scopes these names, violations of these restrictions
   elicit only portability warnings by default.
+* 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
 

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 &,

flang/lib/Evaluate/real.cpp

@@ -401,48 +401,59 @@ 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) {
+        result.value =
+            result.value.Subtract(adj, rounding).AccumulateFlags(result.flags);
+        if (result.value.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;
 }

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,43 @@ 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) {
+        tmp -= adj;
+        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 +270,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)) {

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