[flang] Cope with overflow in real MOD/MODULO

[klausler]

Summary

In folding and in the runtime library for real MOD/MODULO(A,P),
detect overflow from the division A/P and return a properly signed
zero result. (When A/P overflows and both A and P are finite numbers
with nonzero P, the quotient would be a large integer when rounded to
the precision of the floating-point representation.)

[vzakhari]

Sorry for the delay. LGTM.

flang/lib/Evaluate/real.cpp

ValueWithRealFlags<Real<W, P>> Real<W, P>::HYPOT(
   return result;
 }
 
 // MOD(x,y) = x - AINT(x/y)*y
 template <typename W, int P>
 ValueWithRealFlags<Real<W, P>> Real<W, P>::MOD(
     const Real &y, Rounding rounding) const {
   ValueWithRealFlags<Real> result;
-  Real quotient{Divide(y, rounding).AccumulateFlags(result.flags)};
-  Real toInt{quotient.ToWholeNumber(common::RoundingMode::ToZero)
+  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 (IsNegative()) {
+      result.value = Real{}.Negate(); // -0.
+    }
+  } 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
 template <typename W, int P>
 ValueWithRealFlags<Real<W, P>> Real<W, P>::MODULO(
     const Real &y, Rounding rounding) const {
   ValueWithRealFlags<Real> result;
-  Real quotient{Divide(y, rounding).AccumulateFlags(result.flags)};
-  Real toInt{quotient.ToWholeNumber(common::RoundingMode::Down)
+  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.
+    }
+  } 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;
 }
 
 template <typename W, int P>
 ValueWithRealFlags<Real<W, P>> Real<W, P>::DIM(
     const Real &y, Rounding rounding) const {
   ValueWithRealFlags<Real> result;
   if (IsNotANumber() || y.IsNotANumber()) {

flang/runtime/numeric.cpp

inline T IntMod(T x, T p, const char *sourceFile, int sourceLine) {
   }
   auto mod{x - (x / p) * p};
   if (IS_MODULO && (x > 0) != (p > 0)) {
     mod += p;
   }
   return mod;
 }
 template <bool IS_MODULO, typename T>
-inline T RealMod(T x, T p, const char *sourceFile, int sourceLine) {
+inline T RealMod(T a, T p, const char *sourceFile, int sourceLine) {
   if (p == 0) {
     Terminator{sourceFile, sourceLine}.Crash(
         IS_MODULO ? "MODULO with P==0" : "MOD with P==0");
   }
-  if constexpr (IS_MODULO) {
-    return x - std::floor(x / p) * p;
-  } else {
-    return x - std::trunc(x / p) * p;
+  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);
   }
+  T toInt{IS_MODULO ? std::floor(quotient) : std::trunc(quotient)};
+  return a - toInt * p;
 }
 
 // RRSPACING (16.9.164)
 template <int PREC, typename T> inline T RRSpacing(T x) {
   if (std::isnan(x)) {
     return x; // NaN -> same NaN
   } else if (std::isinf(x)) {
     return std::numeric_limits<T>::quiet_NaN(); // +/-Inf -> NaN

flang/test/Evaluate/fold-mod.f90

 ! RUN: %python %S/test_folding.py %s %flang_fc1
 ! Tests folding of integer and real MOD and MODULO
 module m1
   logical, parameter :: test_mod_i1 = mod(8, 5) == 3
   logical, parameter :: test_mod_i2 = mod(-8, 5) == -3
   logical, parameter :: test_mod_i3 = mod(8, -5) == 3
   logical, parameter :: test_mod_i4 = mod(-8, -5) == -3
 
   logical, parameter :: test_mod_r1 = mod(3., 2.) == 1.
   logical, parameter :: test_mod_r2 = mod(8., 5.) == 3.
   logical, parameter :: test_mod_r3 = mod(-8., 5.) == -3.
   logical, parameter :: test_mod_r4 = mod(8., -5.) == 3.
   logical, parameter :: test_mod_r5 = mod(-8., -5.) == -3.
 
+  logical, parameter :: test_mod_r10a = mod(huge(0.), tiny(0.)) == 0.
+  logical, parameter :: test_mod_r10b = sign(1., mod(huge(0.), tiny(0.))) == 1.
+  logical, parameter :: test_mod_r11a = mod(-huge(0.), tiny(0.)) == 0.
+  logical, parameter :: test_mod_r11b = sign(1., mod(-huge(0.), tiny(0.))) == -1.
+  logical, parameter :: test_mod_r12a = mod(huge(0.), -tiny(0.)) == 0.
+  logical, parameter :: test_mod_r12b = sign(1., mod(huge(0.), -tiny(0.))) == 1.
+  logical, parameter :: test_mod_r13a = mod(huge(0.), tiny(0.)) == 0.
+  logical, parameter :: test_mod_r13b = sign(1., mod(-huge(0.), -tiny(0.))) == -1.
+
   logical, parameter :: test_modulo_i1 = modulo(8, 5) == 3
   logical, parameter :: test_modulo_i2 = modulo(-8, 5) == 2
   logical, parameter :: test_modulo_i3 = modulo(8, -5) == -2
   logical, parameter :: test_modulo_i4 = modulo(-8, -5) == -3
 
   logical, parameter :: test_modulo_r1 = modulo(8., 5.) == 3.
   logical, parameter :: test_modulo_r2 = modulo(-8., 5.) == 2.
   logical, parameter :: test_modulo_r3 = modulo(8., -5.) == -2.
   logical, parameter :: test_modulo_r4 = modulo(-8., -5.) == -3.
+
+  logical, parameter :: test_modulo_r10a = modulo(huge(0.), tiny(0.)) == 0.
+  logical, parameter :: test_modulo_r10b = sign(1., modulo(huge(0.), tiny(0.))) == 1.
+  logical, parameter :: test_modulo_r11a = modulo(-huge(0.), tiny(0.)) == 0.
+  logical, parameter :: test_modulo_r11b = sign(1., modulo(-huge(0.), tiny(0.))) == 1.
+  logical, parameter :: test_modulo_r12a = modulo(huge(0.), -tiny(0.)) == 0.
+  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.
 end module