[flang][runtime] Use std::fmod for most MOD/MODULO

[klausler]

The new accurate algorithm for real MOD and MODULO in the runtime is not as fast as std::fmod(), which is also accurate. So use std::fmod() for those floating-point types that it supports.

Fixes #78641.

[llvmbot]

@llvm/pr-subscribers-flang-runtime

Author: Peter Klausler (klausler)

Changes

The new accurate algorithm for real MOD and MODULO in the runtime is not as fast as std::fmod(), which is also accurate. So use std::fmod() for those floating-point types that it supports.

Fixes #78641.

---

Full diff: https://github.com/llvm/llvm-project/pull/78745.diff

1 Files Affected:

• (modified) flang/runtime/numeric.cpp (+25-17)

diff --git a/flang/runtime/numeric.cpp b/flang/runtime/numeric.cpp
index 3f6f553e7bb554d..ad6b0e854522496 100644
--- a/flang/runtime/numeric.cpp
+++ b/flang/runtime/numeric.cpp
@@ -145,25 +145,33 @@ inline RT_API_ATTRS T RealMod(
   } 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 (tmp < pAbs) {
+    } else if constexpr (std::is_same_v<T, float> ||
+        std::is_same_v<T, double> || std::is_same_v<T, long double>) {
+      tmp = std::fmod(tmp, pAbs);
+    } 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) and the values are not
+      // integers
+      // 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.
+      for (T adj{SetExponent(pAbs, Exponent<int>(tmp))}; tmp >= pAbs;
+           adj /= 2) {
+        if (tmp >= adj) {
+          tmp -= adj;
+          if (tmp == 0) {
+            break;
+          }
         }
       }
     }

[vzakhari]

Maybe we should just use fmod(a, p) for these types, and avoid the abs computations, tmp < pAbs evaluation, and the tmp = -tmp fixup below. I think this may remove some overhead for the likely cases (e.g. when a >= p). What do you think?

[klausler]

Getting the result right for both MOD and MODULO when either a or p or both are negative is tricky. Will try.

[klausler]

Updated.

[vzakhari]

Thank you, Peter! It looks good to me, though, there might be some unneeded overhead (see my comment).
I suggest merging this, and I will check out performance after. It looks like using fmod is the standard way in other compilers as well.

[vzakhari]

Thank you!