[flang-commits] [flang] [flang][runtime] Use std::fmod for most MOD/MODULO (PR #78745)

[Peter Klausler via flang-commits]

================
@@ -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);
----------------
klausler wrote:

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

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