[builtins] Divide shouldn't underflow if rounded result would be normal.

[efriedma]

Summary

We were treating certain edge cases that are actually normal as denormal results, and flushing them to zero; we shouldn't do that. Not sure this is the cleanest way to implement this edge case, but I wanted to avoid adding any code on the common path.

(This doesn't touch the behavior for results that are actually denormal; they're still flushed to zero.)

[compnerd]

Would be nice if @scanon would also take a look, but, this seems like a good thing to fix.

[scanon]

In the parlance of IEEE 754, there are two ways to "detect tininess": "before rounding" and "after rounding". The standard doesn't define how to flush subnormal results, but in practice most HW flushes results that are "tiny". The existing code flushes as though tininess is detected before rounding. This proposed update flushes as though tininess were detected after rounding.

Of mainstream platforms supported by LLVM, only x86 detects tininess after rounding, and these functions should approximately never be used on x86, because HW floating-point is always available there. So I would lean towards keeping the existing behavior, so that e.g. soft-float and hard-float ARM match their behavior.

[efriedma]

Okay, that makes sense, in general.

Do you have any suggestions for fixing the algorithm so it doesn't consider the result "tiny" when it actually isn't?

[scanon]

These results *are* tiny in the before rounding sense.

[efriedma]

Rereading my question, it isn't really clear; I'll try to explain in more detail.

The current algorithm for computing the "unrounded" result does not attempt to ensure that it is actually between the two nearest floating-point values. Instead, it only tries to produce a value such that a round-to-even step will result in the correct answer. This means that the "writtenExponent < 1" check does not reliably compute whether a value is tiny: there are cases where the check fails, but the mathematical result of the division is actually in the normal range. The most obvious example of this is my testcase: DBL_MIN is normal, so (DBL_MIN * 2) / 2 should also be normal.

Given the intermediate results the code currently, is there some reasonable way to detect this case?

[scanon]

Ah, now I see what you're talking about. And in fact, because of the way divide works out, there's a little gap of results that are even possible to achieve just below each binade boundary, so the code you have here will work out fine. We *should* add a comment to clarify this somewhat, but I'm happy to do that in a separate commit. LGTM.

lib/builtins/divdf3.c

__divdf3(fp_t a, fp_t b) {
     const int writtenExponent = quotientExponent + exponentBias;
 
     if (writtenExponent >= maxExponent) {
         // If we have overflowed the exponent, return infinity.
         return fromRep(infRep | quotientSign);
     }
 
     else if (writtenExponent < 1) {
+        if (writtenExponent == 0) {
+            // Check whether the rounded result is normal.
+            const bool round = (residual << 1) > bSignificand;
+            // Clear the implicit bit.
+            rep_t absResult = quotient & significandMask;
+            // Round.
+            absResult += round;
+            if (absResult & ~significandMask) {
+                // The rounded result is normal; return it.
+                return fromRep(absResult | quotientSign);
+            }
+        }
         // Flush denormals to zero.  In the future, it would be nice to add
         // code to round them correctly.
         return fromRep(quotientSign);
     }
 
     else {
         const bool round = (residual << 1) > bSignificand;
         // Clear the implicit bit

lib/builtins/divsf3.c

__divsf3(fp_t a, fp_t b) {
     const int writtenExponent = quotientExponent + exponentBias;
 
     if (writtenExponent >= maxExponent) {
         // If we have overflowed the exponent, return infinity.
         return fromRep(infRep | quotientSign);
     }
 
     else if (writtenExponent < 1) {
+        if (writtenExponent == 0) {
+            // Check whether the rounded result is normal.
+            const bool round = (residual << 1) > bSignificand;
+            // Clear the implicit bit.
+            rep_t absResult = quotient & significandMask;
+            // Round.
+            absResult += round;
+            if (absResult & ~significandMask) {
+                // The rounded result is normal; return it.
+                return fromRep(absResult | quotientSign);
+            }
+        }
         // Flush denormals to zero.  In the future, it would be nice to add
         // code to round them correctly.
         return fromRep(quotientSign);
     }
 
     else {
         const bool round = (residual << 1) > bSignificand;
         // Clear the implicit bit

lib/builtins/divtf3.c

COMPILER_RT_ABI fp_t __divtf3(fp_t a, fp_t b) {
 
     const int writtenExponent = quotientExponent + exponentBias;
 
     if (writtenExponent >= maxExponent) {
         // If we have overflowed the exponent, return infinity.
         return fromRep(infRep | quotientSign);
     }
     else if (writtenExponent < 1) {
+        if (writtenExponent == 0) {
+            // Check whether the rounded result is normal.
+            const bool round = (residual << 1) > bSignificand;
+            // Clear the implicit bit.
+            rep_t absResult = quotient & significandMask;
+            // Round.
+            absResult += round;
+            if (absResult & ~significandMask) {
+                // The rounded result is normal; return it.
+                return fromRep(absResult | quotientSign);
+            }
+        }
         // Flush denormals to zero.  In the future, it would be nice to add
         // code to round them correctly.
         return fromRep(quotientSign);
     }
     else {
         const bool round = (residual << 1) >= bSignificand;
         // Clear the implicit bit
         rep_t absResult = quotient & significandMask;

test/builtins/Unit/divdf3_test.c

+// RUN: %clang_builtins %s %librt -o %t && %run %t
+//===--------------- divdf3_test.c - Test __divdf3 ------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+//
+// This file tests __divdf3 for the compiler_rt library.
+//
+//===----------------------------------------------------------------------===//
+
+#include "int_lib.h"
+#include <stdio.h>
+
+#include "fp_test.h"
+
+// Returns: a / b
+COMPILER_RT_ABI double __divdf3(double a, double b);
+
+int test__divdf3(double a, double b, uint64_t expected)
+{
+    double x = __divdf3(a, b);
+    int ret = compareResultD(x, expected);
+
+    if (ret){
+        printf("error in test__divdf3(%.20e, %.20e) = %.20e, "
+               "expected %.20e\n", a, b, x,
+               fromRep64(expected));
+    }
+    return ret;
+}
+
+int main()
+{
+    // 1/3
+    if (test__divdf3(1., 3., 0x3fd5555555555555ULL))
+      return 1;
+    // smallest normal result
+    if (test__divdf3(4.450147717014403e-308, 2., 0x10000000000000ULL))
+      return 1;
+
+    return 0;
+}

test/builtins/Unit/divsf3_test.c

+// RUN: %clang_builtins %s %librt -o %t && %run %t
+//===--------------- divsf3_test.c - Test __divsf3 ------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+//
+// This file tests __divsf3 for the compiler_rt library.
+//
+//===----------------------------------------------------------------------===//
+
+#include "int_lib.h"
+#include <stdio.h>
+
+#include "fp_test.h"
+
+// Returns: a / b
+COMPILER_RT_ABI float __divsf3(float a, float b);
+
+int test__divsf3(float a, float b, uint32_t expected)
+{
+    float x = __divsf3(a, b);
+    int ret = compareResultF(x, expected);
+
+    if (ret){
+        printf("error in test__divsf3(%.20e, %.20e) = %.20e, "
+               "expected %.20e\n", a, b, x,
+               fromRep32(expected));
+    }
+    return ret;
+}
+
+int main()
+{
+    // 1/3
+    if (test__divsf3(1.f, 3.f, 0x3EAAAAABU))
+      return 1;
+    // smallest normal result
+    if (test__divsf3(2.3509887e-38, 2., 0x00800000U))
+      return 1;
+
+    return 0;
+}

test/builtins/Unit/divtf3_test.c

 
 int test__divtf3(long double a, long double b,
                  uint64_t expectedHi, uint64_t expectedLo)
 {
     long double x = __divtf3(a, b);
     int ret = compareResultLD(x, expectedHi, expectedLo);
 
     if (ret){
-        printf("error in test__divtf3(%.20Lf, %.20Lf) = %.20Lf, "
-               "expected %.20Lf\n", a, b, x,
+        printf("error in test__divtf3(%.20Le, %.20Le) = %.20Le, "
+               "expected %.20Le\n", a, b, x,
                fromRep128(expectedHi, expectedLo));
     }
     return ret;
 }
 
 char assumption_1[sizeof(long double) * CHAR_BIT == 128] = {0};
 
 #endif
if (test__divtf3(0x1.2d3456f789ba6322bc665544edefp-234L,
                      UINT64_C(0x507b38442b539266),
                      UINT64_C(0x22ce0f1d024e1252)))
         return 1;
     if (test__divtf3(0x1.2345f6b77b7a8953365433abcdefp+234L,
                      0x1.edcba987d6bb3aa467754354321fp-4055L,
                      UINT64_C(0x50bf2e02f0798d36),
                      UINT64_C(0x5e6fcb6b60044078)))
         return 1;
+    if (test__divtf3(6.72420628622418701252535563464350521E-4932L,
+                     2.L,
+                     UINT64_C(0x0001000000000000),
+                     UINT64_C(0)))
+        return 1;
 
 #else
     printf("skipped\n");
 
 #endif
     return 0;
 }