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

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

Author: lntue

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;

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;
 }

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.fffffep-1F, UINT32_C(0x3f800001)))
+      return 1;
+
     return 0;
 }

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");