Merge pull request #8023 from Ondra97/codex/fix-scaled-hyperbolic

fix(libc): avoid premature scaled hyperbolic overflow
This commit is contained in:
Dean authored and GitHub committed 2026-10-21 06:26:06 +13:00
commit 542c65d74c
3 files changed
+131 -2

No files matched your search

+95
View File
@@ -116,6 +116,8 @@ double APS5_VABI sinh_nid_postfix(double);
double APS5_VABI cosh_nid_postfix(double);
float APS5_VABI coshf_nid_postfix(float);
float APS5_VABI sinhf_nid_postfix(float);
float APS5_VABI _FSinh_nid_postfix(float, float);
float APS5_VABI _FCosh_nid_postfix(float, float);
double APS5_VABI asinh_nid_postfix(double);
float APS5_VABI asinhf_nid_postfix(float);
double APS5_VABI acosh_nid_postfix(double);
@@ -791,6 +793,98 @@ static void CheckHyperbolic() {
Require(Near(acosh_nid_postfix(1e300), 691.4686750787736, 1e-15));
}
static void CheckScaledHyperbolic() {
const int rounding = std::fegetround();
Require(std::fesetround(FE_TONEAREST) == 0);
unsigned checks = 0;
for (bool sine : {true, false}) {
const auto function = sine ? _FSinh_nid_postfix : _FCosh_nid_postfix;
const char* name = sine ? "_FSinh" : "_FCosh";
const auto checkValue = [&](float x, float y) {
const long double value = sine ? std::sinh(static_cast<long double>(x)) : std::cosh(static_cast<long double>(x));
const float expected = static_cast<float>(value * static_cast<long double>(y));
std::feclearexcept(FE_ALL_EXCEPT);
errno = 12345;
const float result = function(x, y);
const auto actualBits = std::bit_cast<std::uint32_t>(result);
const auto expectedBits = std::bit_cast<std::uint32_t>(expected);
const auto distance = actualBits > expectedBits ? actualBits - expectedBits : expectedBits - actualBits;
const bool exact = expected == 0.0f || std::isinf(expected);
if ((exact ? actualBits != expectedBits : !std::isfinite(result) || distance > 2) ||
(std::isfinite(expected) && (std::fetestexcept(FE_OVERFLOW | FE_INVALID) != 0 || errno != 12345))) {
std::fprintf(stderr, "%s(%a, %a): got %08x, expected %08x, errno %d, flags %x\n",
name, static_cast<double>(x), static_cast<double>(y), actualBits, expectedBits, errno, std::fetestexcept(FE_ALL_EXCEPT));
std::abort();
}
++checks;
};
checkValue(90.0f, 1e-20f);
checkValue(100.0f, 1e-40f);
const float magnitudes[] = {0.0f, 0x1p-126f, 0.001f, 0.5f, 1.0f, 4.0f, 7.0f, 8.0f, 10.0f, 20.0f, 50.0f, 88.0f, 90.0f, 100.0f, 150.0f, 180.0f, 190.0f, 193.0f};
const float scales[] = {0.0f, 0x1p-149f, 0x1p-140f, 0x1p-126f, 1e-20f, 0.5f, 1.0f, 1e20f, std::numeric_limits<float>::max()};
for (float magnitude : magnitudes) {
for (float scale : scales) {
for (float xSign : {1.0f, -1.0f}) {
for (float ySign : {1.0f, -1.0f}) checkValue(magnitude * xSign, scale * ySign);
}
}
}
for (float scale : scales) {
if (scale == 0.0f) continue;
for (long double target : {static_cast<long double>(std::numeric_limits<float>::max()), 0x1p-150L}) {
if (!sine && target < scale) continue;
const long double ratio = target / static_cast<long double>(scale);
const float boundary = static_cast<float>(sine ? std::asinh(ratio) : std::acosh(ratio));
float magnitude = std::nextafter(std::nextafter(boundary, -INFINITY), -INFINITY);
for (unsigned neighbor = 0; neighbor != 5; ++neighbor) {
for (float xSign : {1.0f, -1.0f}) {
for (float ySign : {1.0f, -1.0f}) checkValue(magnitude * xSign, scale * ySign);
}
magnitude = std::nextafter(magnitude, INFINITY);
}
}
}
const auto checkSpecial = [&](float x, float y, std::uint32_t expectedBits, int error, int flags) {
std::feclearexcept(FE_ALL_EXCEPT);
std::feraiseexcept(FE_DIVBYZERO);
errno = 12345;
const float result = function(x, y);
const int raised = std::fetestexcept(FE_ALL_EXCEPT);
if (std::bit_cast<std::uint32_t>(result) != expectedBits || errno != error || raised != (FE_DIVBYZERO | flags)) {
std::fprintf(stderr, "%s special (%a, %a): got %08x, expected %08x, errno %d, flags %x\n",
name, static_cast<double>(x), static_cast<double>(y), std::bit_cast<std::uint32_t>(result), expectedBits, errno, raised);
std::abort();
}
++checks;
};
const float infinity = std::numeric_limits<float>::infinity();
for (float sign : {1.0f, -1.0f}) {
checkSpecial(sign * infinity, 0.0f, 0x7fc00000u, EDOM, FE_INVALID);
checkSpecial(sign * 100.0f, 0.0f, sine && sign < 0 ? 0x80000000u : 0u, 12345, 0);
checkSpecial(sign * 100.0f, -0.0f, sine && sign < 0 ? 0u : 0x80000000u, 12345, 0);
checkSpecial(sign * 1.0f, infinity, sine && sign < 0 ? 0xff800000u : 0x7f800000u, 12345, 0);
}
for (std::uint32_t bits : {0x7fc12345u, 0xffc54321u, 0x7fa12345u, 0xffa54321u}) {
const float nan = std::bit_cast<float>(bits);
checkSpecial(nan, 0.0f, bits, 12345, 0);
checkSpecial(infinity, nan, bits, 12345, 0);
checkSpecial(nan, std::bit_cast<float>(0x7fc99999u), bits, 12345, 0);
}
checkSpecial(0.0f, infinity, sine ? 0x7fc00000u : 0x7f800000u, sine ? EDOM : 12345, sine ? FE_INVALID : 0);
checkSpecial(1.0f, std::numeric_limits<float>::max(), 0x7f800000u, 12345, FE_OVERFLOW | FE_INEXACT);
checkSpecial(100.0f, 1.0f, 0x7f800000u, ERANGE, FE_OVERFLOW | FE_INEXACT);
for (int mode : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO}) {
Require(std::fesetround(mode) == 0);
checkSpecial(std::numeric_limits<float>::max(), 0x1p-149f, 0x7f800000u, ERANGE, FE_OVERFLOW | FE_INEXACT);
}
Require(std::fesetround(FE_TONEAREST) == 0);
if (sine) checkSpecial(0x1p-126f, 0x1p-149f, 0u, 12345, FE_UNDERFLOW | FE_INEXACT);
}
std::feclearexcept(FE_ALL_EXCEPT);
Require(std::fesetround(rounding) == 0);
std::printf("Scaled hyperbolic checks: %u\n", checks);
}
static void CheckSinhf() {
const float infinity = std::numeric_limits<float>::infinity();
const struct { std::uint32_t input; std::uint32_t output; } cases[] = {
@@ -1355,6 +1449,7 @@ int main() {
Require(std::isinf(__powidf2_nid_postfix(-0.0, -3)) && std::signbit(__powidf2_nid_postfix(-0.0, -3)));
Require(__powidf2_nid_postfix(2.0, -127) == 0x1p-127);
CheckHyperbolic();
CheckScaledHyperbolic();
CheckSinhf();
CheckAtanhf();
CheckHypotSpecialValue();