[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))`.
This commit is contained in:
lntue 2024-12-17 09:24:57 -08:00 committed by GitHub
parent e8a6563768
commit a57f4c7009
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
4 changed files with 19 additions and 1 deletions

View File

@ -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 // Suppose 1/b - P * 2^-W < x < 1/b + P * 2^-W
x_UQ0 -= RECIPROCAL_PRECISION; x_UQ0 -= RECIPROCAL_PRECISION;
// Now 1/b - (2*P) * 2^-W < x < 1/b // Now 1/b - (2*P) * 2^-W < x < 1/b
// FIXME Is x_UQ0 still >= 0.5?
rep_t quotient_UQ1, dummy; rep_t quotient_UQ1, dummy;
wideMultiply(x_UQ0, aSignificand << 1, &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. // adjust it to be in [1.0, 2.0) as UQ1.SB.
rep_t residualLo; rep_t residualLo;
if (quotient_UQ1 < (implicitBit << 1)) { 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, // Highest bit is 0, so just reinterpret quotient_UQ1 as UQ1.SB,
// effectively doubling its value as well as its error estimation. // effectively doubling its value as well as its error estimation.
residualLo = (aSignificand << (significandBits + 1)) - quotient_UQ1 * bSignificand; residualLo = (aSignificand << (significandBits + 1)) - quotient_UQ1 * bSignificand;

View File

@ -113,5 +113,9 @@ int main()
if (test__divdf3(0x1.0p-1022, 0x1.0028p+52, UINT64_C(0x1))) if (test__divdf3(0x1.0p-1022, 0x1.0028p+52, UINT64_C(0x1)))
return 1; 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; return 0;
} }

View File

@ -107,5 +107,9 @@ int main()
if (test__divsf3(0x1.fffffep-126F, 2.F, UINT32_C(0x00800000))) if (test__divsf3(0x1.fffffep-126F, 2.F, UINT32_C(0x00800000)))
return 1; 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; return 0;
} }

View File

@ -156,6 +156,11 @@ int main() {
UINT64_C(0x0001000000000000), UINT64_C(0))) UINT64_C(0x0001000000000000), UINT64_C(0)))
return 1; 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 #else
printf("skipped\n"); printf("skipped\n");