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

[Peter Klausler via flang-commits]

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

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 https://github.com/llvm/llvm-project/issues/78641.

>From 27b58346a29a6e689daaaf121b04caca56d84245 Mon Sep 17 00:00:00 2001
From: Peter Klausler <pklausler at nvidia.com>
Date: Fri, 19 Jan 2024 09:01:20 -0800
Subject: [PATCH] [flang][runtime] Use std::fmod for most MOD/MODULO

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 https://github.com/llvm/llvm-project/issues/78641.
---
 flang/runtime/numeric.cpp | 42 +++++++++++++++++++++++----------------
 1 file changed, 25 insertions(+), 17 deletions(-)

diff --git a/flang/runtime/numeric.cpp b/flang/runtime/numeric.cpp
index 3f6f553e7bb554..ad6b0e85452249 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;
+          }
         }
       }
     }