[compiler-rt] [compiler-rt] Fix a bug in fp_div_impl when an intermediate result is out of expected range. (PR #119449)

[via llvm-commits]

https://github.com/lntue created https://github.com/llvm/llvm-project/pull/119449

Before this fix, `1.0L / (1.0L - 0x1.0p-113L)1 will return `2 * (1 + eps(1))`.

>From da6c19585c3041d3b656323157fcf20c72b89600 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue at google.com>
Date: Tue, 10 Dec 2024 15:16:36 -0500
Subject: [PATCH] [compiler-rt] Fix a bug in fp_div_impl when an intermediate
 result is out of expected range.

Before this fix, 1.0L / (1.0L - 0x1.0p-113L) will return 2 * (1 + eps(1)).
---
 compiler-rt/lib/builtins/fp_div_impl.inc     | 7 ++++++-
 compiler-rt/test/builtins/Unit/divdf3_test.c | 4 ++++
 compiler-rt/test/builtins/Unit/divsf3_test.c | 4 ++++
 compiler-rt/test/builtins/Unit/divtf3_test.c | 5 +++++
 4 files changed, 19 insertions(+), 1 deletion(-)

diff --git a/compiler-rt/lib/builtins/fp_div_impl.inc b/compiler-rt/lib/builtins/fp_div_impl.inc
index 29bcd1920edfb4..de61e55cd083be 100644
--- a/compiler-rt/lib/builtins/fp_div_impl.inc
+++ b/compiler-rt/lib/builtins/fp_div_impl.inc
@@ -334,7 +334,6 @@ static __inline fp_t __divXf3__(fp_t a, fp_t b) {
   // Suppose 1/b - P * 2^-W < x < 1/b + P * 2^-W
   x_UQ0 -= RECIPROCAL_PRECISION;
   // Now 1/b - (2*P) * 2^-W < x < 1/b
-  // FIXME Is x_UQ0 still >= 0.5?
 
   rep_t quotient_UQ1, dummy;
   wideMultiply(x_UQ0, aSignificand << 1, &quotient_UQ1, &dummy);
@@ -344,6 +343,12 @@ static __inline fp_t __divXf3__(fp_t a, fp_t b) {
   // adjust it to be in [1.0, 2.0) as UQ1.SB.
   rep_t residualLo;
   if (quotient_UQ1 < (implicitBit << 1)) {
+    if (quotient_UQ1 < implicitBit) {
+      // In a rare case where quotient is < 0.5, we can adjust the quotient and
+      // the written exponent, and then treat them the same way as in [0.5, 1.0)
+      quotient_UQ1 <<= 1;
+      writtenExponent -= 1;
+    }
     // Highest bit is 0, so just reinterpret quotient_UQ1 as UQ1.SB,
     // effectively doubling its value as well as its error estimation.
     residualLo = (aSignificand << (significandBits + 1)) - quotient_UQ1 * bSignificand;
diff --git a/compiler-rt/test/builtins/Unit/divdf3_test.c b/compiler-rt/test/builtins/Unit/divdf3_test.c
index 1b8f2b398a91ff..f1f3f3b4f7e096 100644
--- a/compiler-rt/test/builtins/Unit/divdf3_test.c
+++ b/compiler-rt/test/builtins/Unit/divdf3_test.c
@@ -113,5 +113,9 @@ int main()
     if (test__divdf3(0x1.0p-1022, 0x1.0028p+52, UINT64_C(0x1)))
       return 1;
 
+    // test 1 / (1 - eps(0.5)) = 1 + eps(1)
+    if (test__divdf3(1.0, 0x1.fffffffffffffp-1, UINT64_C(0x3ff0000000000001)))
+      return 1;
+
     return 0;
 }
diff --git a/compiler-rt/test/builtins/Unit/divsf3_test.c b/compiler-rt/test/builtins/Unit/divsf3_test.c
index 7a783cdfb863cb..12ed4396aa891a 100644
--- a/compiler-rt/test/builtins/Unit/divsf3_test.c
+++ b/compiler-rt/test/builtins/Unit/divsf3_test.c
@@ -107,5 +107,9 @@ int main()
     if (test__divsf3(0x1.fffffep-126F, 2.F, UINT32_C(0x00800000)))
       return 1;
 
+    // test 1 / (1 - eps(0.5)) = 1 + eps(1)
+    if (test__divsf3(1.0F, 0x1.ffffep-1F, UINT32_C(0x3f800001)))
+      return 1;
+
     return 0;
 }
diff --git a/compiler-rt/test/builtins/Unit/divtf3_test.c b/compiler-rt/test/builtins/Unit/divtf3_test.c
index d46fcc04776938..3c70e43741afb2 100644
--- a/compiler-rt/test/builtins/Unit/divtf3_test.c
+++ b/compiler-rt/test/builtins/Unit/divtf3_test.c
@@ -156,6 +156,11 @@ int main() {
                    UINT64_C(0x0001000000000000), UINT64_C(0)))
     return 1;
 
+  // test 1 / (1 - eps(0.5)) = 1 + eps(1).
+  if (test__divtf3(1.0L, TF_C(0x1.ffffffffffffffffffffffffffffp-1),
+                   UINT64_C(0x3FFF000000000000), UINT64_C(1)))
+    return 1;
+
 #else
   printf("skipped\n");