[builtins] Make divXf3 handle denormal results

[atrosinenko]

Summary

This is the last patch from the patchset introducing proper support for subnormal results to div[sdt]f3.

[atrosinenko]

Re-upload after amending previous commit

[atrosinenko]

Re-upload after D85731: [NFC][builtins] Make softfloat-related errors less noisy to get rid of "error: unknown type name 'fp_t'" and similar clang-tidy diagnostics for fp_div_impl.inc.

[atrosinenko]

Ping.

[atrosinenko]

Rebase the entire patch stack against the up-to-date master and re-upload.

[atrosinenko]

Re-upload after parent commits were changed.

[atrosinenko]

No-change re-upload: rebase onto current master branch.

[atrosinenko]

Re-upload after changing parent diff.

[atrosinenko]

Re-upload after amending parent diff.

[sepavloff]

LGTM.

[atrosinenko]

Re-upload after amending parent diff + add minor clarification.

compiler-rt/lib/builtins/fp_div_impl.inc

 
 #define HW (typeWidth / 2)
 #define loMask (REP_C(-1) >> HW)
 
 #define NUMBER_OF_ITERATIONS \
   (NUMBER_OF_HALF_ITERATIONS + NUMBER_OF_FULL_ITERATIONS)
 
 #if NUMBER_OF_FULL_ITERATIONS < 1
 #   error At least one full iteration is required
 #endif
 
 static __inline fp_t __divXf3(fp_t a, fp_t b) {
 
   const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
   const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
   const rep_t quotientSign = (toRep(a) ^ toRep(b)) & signBit;
 
   rep_t aSignificand = toRep(a) & significandMask;
   rep_t bSignificand = toRep(b) & significandMask;
   int scale = 0;
 
   // Detect if a or b is zero, denormal, infinity, or NaN.
   if (aExponent - 1U >= maxExponent - 1U ||
       bExponent - 1U >= maxExponent - 1U) {
 
     const rep_t aAbs = toRep(a) & absMask;
     const rep_t bAbs = toRep(b) & absMask;
 
     // NaN / anything = qNaN
     if (aAbs > infRep)
       return fromRep(toRep(a) | quietBit);
     // anything / NaN = qNaN
     if (bAbs > infRep)
       return fromRep(toRep(b) | quietBit);
 
     if (aAbs == infRep) {
       // infinity / infinity = NaN
       if (bAbs == infRep)
         return fromRep(qnanRep);
static __inline fp_t __divXf3(fp_t a, fp_t b) {
   }
 
   // Set the implicit significand bit.  If we fell through from the
   // denormal path it was already set by normalize( ), but setting it twice
   // won't hurt anything.
   aSignificand |= implicitBit;
   bSignificand |= implicitBit;
 
   int writtenExponent = (aExponent - bExponent + scale) + exponentBias;
 
   const rep_t b_UQ1 = bSignificand << (typeWidth - significandBits - 1);
 
   // Align the significand of b as a UQ1.(n-1) fixed-point number in the range
   // [1.0, 2.0) and get a UQ0.n approximate reciprocal using a small minimax
   // polynomial approximation: x0 = 3/4 + 1/sqrt(2) - b/2.
   // In infinite precision, for b in [1, 2):
   //   0 < x0(b) < 1
   //   abs(x0(b) - 1/b) <= 3/4 - 1/sqrt(2)
 
#endif
   wideMultiply(x_UQ0, aSignificand << 1, &quotient_UQ1, &dummy);
   // Now, a/b - 4*P * 2^-W < q < a/b
 
   // quotient_UQ1 is in [0.5, 2.0) as UQ1.(SB+1), adjust it to be in [1.0, 2.0) as UQ1.SB
   rep_t residualLo;
   if (quotient_UQ1 < (implicitBit << 1)) {
     residualLo = (aSignificand << (significandBits + 1)) - quotient_UQ1 * bSignificand;
     writtenExponent -= 1;
+    aSignificand <<= 1;
 
     // the error is doubled
   } else {
     quotient_UQ1 >>= 1;
     residualLo = (aSignificand << significandBits) - quotient_UQ1 * bSignificand;
   }
   // Now, q cannot be greater than a/b and can differ by at most 8*P * 2^-W + 2^-SB
   // Each NextAfter() increments the floating point value by at least 2^-SB
#endif
 
   // If we have overflowed the exponent, return infinity
   if (writtenExponent >= maxExponent)
     return fromRep(infRep | quotientSign);
 
   // Now, quotient_UQ1_SB <= the correctly-rounded result
   // and may need taking NextAfter() up to 3 times (see error estimates above)
   // r = a - b * q
-
-  if (writtenExponent < 0) {
-    // Result is definitely subnormal, flushing to zero
-    return fromRep(quotientSign);
-  }
-
+  rep_t absResult;
+  if (writtenExponent > 0) {
     // Clear the implicit bit
-  rep_t absResult = quotient_UQ1 & significandMask;
+    absResult = quotient_UQ1 & significandMask;
     // Insert the exponent
     absResult |= (rep_t)writtenExponent << significandBits;
+    residualLo <<= 1;
+  } else {
+    // Prevent shift amount from being negative
+    if (significandBits + writtenExponent < 0)
+      return fromRep(quotientSign);
+
+    absResult = quotient_UQ1 >> (-writtenExponent + 1);
+
+    // multiplied by two to prevent shift amount to be negative
+    residualLo = (aSignificand << (significandBits + writtenExponent)) - (absResult * bSignificand << 1);
+  }
 
   // Round
-  residualLo <<= 1;
   residualLo += absResult & 1; // tie to even
   absResult += residualLo > bSignificand;
 #if defined(QUAD_PRECISION) || (defined(SINGLE_PRECISION) && NUMBER_OF_HALF_ITERATIONS > 0)
   // Do not round Infinity to NaN
   absResult += absResult < infRep && residualLo > (2 + 1) * bSignificand;
 #endif
 #if defined(QUAD_PRECISION)
   absResult += absResult < infRep && residualLo > (4 + 1) * bSignificand;
 #endif
-
-  if ((absResult & ~significandMask) == 0) {
-    // Result is subnormal, flushing to zero
-    return fromRep(quotientSign);
-  }
-  // Result is normal, insert the sign and return
   return fromRep(absResult | quotientSign);
 }

compiler-rt/test/builtins/Unit/divdf3_test.c

int main()
 
     // divisor is exactly 1.0
     if (test__divdf3(0x1.0p+0, 0x1.0p+0, UINT64_C(0x3ff0000000000000)))
       return 1;
     // divisor is 1.0 as UQ1.31
     if (test__divdf3(0x1.0p+0, 0x1.00000001p+0, UINT64_C(0x3fefffffffe00000)))
       return 1;
 
+    // smallest normal value divided by 2.0
+    if (test__divdf3(0x1.0p-1022, 2., UINT64_C(0x0008000000000000)))
+      return 1;
+    // smallest subnormal result
+    if (test__divdf3(0x1.0p-1022, 0x1.0p+52, UINT64_C(0x0000000000000001)))
+      return 1;
+
     // some misc test cases obtained by fuzzing against h/w implementation
     if (test__divdf3(0x1.fdc239dd64735p-658, -0x1.fff9364c0843fp-948, UINT64_C(0xd20fdc8fc0ceffb1)))
       return 1;
     if (test__divdf3(-0x1.78abb261d47c8p+794, 0x1.fb01d537cc5aep+266, UINT64_C(0xe0e7c6148ffc23e3)))
       return 1;
     if (test__divdf3(-0x1.da7dfe6048b8bp-875, 0x1.ffc7ea3ff60a4p-610, UINT64_C(0xaf5dab1fe0269e2a)))
       return 1;
+    if (test__divdf3(0x1.0p-1022, 0x1.9p+5, UINT64_C(0x000051eb851eb852)))
+      return 1;
+    if (test__divdf3(0x1.0p-1022, 0x1.0028p+41, UINT64_C(0x00000000000007ff)))
+      return 1;
+    if (test__divdf3(0x1.0p-1022, 0x1.0028p+52, UINT64_C(0x1)))
+      return 1;
 
     return 0;
 }

compiler-rt/test/builtins/Unit/divsf3_test.c

int main()
 
     // divisor is exactly 1.0
     if (test__divsf3(0x1.0p+0F, 0x1.0p+0F, UINT32_C(0x3f800000)))
       return 1;
     // divisor is 1.0 as UQ1.15
     if (test__divsf3(0x1.0p+0F, 0x1.0001p+0F, UINT32_C(0x3f7fff00)))
       return 1;
 
+    // smallest normal value divided by 2.0
+    if (test__divsf3(0x1.0p-126F, 2.0F, UINT32_C(0x00400000)))
+      return 1;
+    // smallest subnormal result
+    if (test__divsf3(0x1.0p-126F, 0x1p+23F, UINT32_C(0x00000001)))
+      return 1;
+
+    // some misc test cases obtained by fuzzing against h/w implementation
+    if (test__divsf3(-0x1.3e75e6p-108F, -0x1.cf372p+38F, UINT32_C(0x00000006)))
+      return 1;
+    if (test__divsf3(0x1.e77c54p+81F, -0x1.e77c52p-47F, UINT32_C(0xff800000)))
+      return 1;
+    if (test__divsf3(0x1.fffffep-126F, 2.F, UINT32_C(0x00800000)))
+      return 1;
+
     return 0;
 }

compiler-rt/test/builtins/Unit/divtf3_test.c

if (test__divtf3(0x1.0p+0L,
         return 1;
     // divisor is 1.0 as UQ1.63
     if (test__divtf3(0x1.0p+0L,
                      0x1.0000000000000001p+0L,
                      UINT64_C(0x3ffeffffffffffff),
                      UINT64_C(0xfffe000000000000)))
         return 1;
 
+    // smallest normal value divided by 2.0
+    if (test__divtf3(0x1.0p-16382L, 2.L, UINT64_C(0x0000800000000000), UINT64_C(0x0)))
+      return 1;
+    // smallest subnormal result
+    if (test__divtf3(0x1.0p-1022L, 0x1p+52L, UINT64_C(0x0), UINT64_C(0x1)))
+      return 1;
+
     // any / any
     if (test__divtf3(0x1.a23b45362464523375893ab4cdefp+5L,
                      0x1.eedcbaba3a94546558237654321fp-1L,
                      UINT64_C(0x4004b0b72924d407),
                      UINT64_C(0x0717e84356c6eba2)))
         return 1;
     if (test__divtf3(0x1.a2b34c56d745382f9abf2c3dfeffp-50L,
                      0x1.ed2c3ba15935332532287654321fp-9L,