fix(libc): preserve nextafter floating-point environment behavior

core/libs/prx/libc/src/Math.cpp

@@ -1273,6 +1273,8 @@ float APS5_VABI cbrtf_nid_postfix(float x) {
     return static_cast<float>(T);
 }
 float APS5_VABI remainderf_nid_postfix(float x, float y) { return std::remainder(x, y); }
+void APS5_VABI _Feraise_nid_postfix(int exceptions);
+
 float APS5_VABI nextafterf_nid_postfix(float x, float y) {
     const std::uint32_t ux = std::bit_cast<std::uint32_t>(x);
     const std::uint32_t uy = std::bit_cast<std::uint32_t>(y);
@@ -1285,17 +1287,20 @@ float APS5_VABI nextafterf_nid_postfix(float x, float y) {
         return y;
     }
     if ((ux & 0x7fffffffu) == 0) {
-        return std::bit_cast<float>((uy & 0x80000000u) | 1u);
+        const volatile float smallest = std::bit_cast<float>((uy & 0x80000000u) | 1u);
+        volatile float product = smallest * smallest;
+        return product == smallest ? product : smallest;
     }
-    const bool towardZero = (x > 0.0f) == (x > y);
+    const bool towardZero = ((ux ^ uy) & 0x80000000u) != 0 || ux > uy;
     const std::uint32_t next = towardZero ? ux - 1 : ux + 1;
     const std::uint32_t exponent = next & 0x7f800000u;
     if (exponent == 0x7f800000u) {
-        errno = ERANGE;
+        _Feraise_nid_postfix(0x08);
         return x + x;
     }
     if (exponent == 0) {
-        errno = ERANGE;
+        volatile float product = x * x;
+        if (product != x) _Feraise_nid_postfix(0x10);
     }
     return std::bit_cast<float>(next);
 }
@@ -1311,17 +1316,20 @@ double APS5_VABI nextafter_nid_postfix(double x, double y) {
         return y;
     }
     if ((ux & 0x7fffffffffffffffull) == 0) {
-        return std::bit_cast<double>((uy & 0x8000000000000000ull) | 1u);
+        const volatile double smallest = std::bit_cast<double>((uy & 0x8000000000000000ull) | 1u);
+        volatile double product = smallest * smallest;
+        return product == smallest ? product : smallest;
     }
-    const bool towardZero = (x > 0.0) == (x > y);
+    const bool towardZero = ((ux ^ uy) & 0x8000000000000000ull) != 0 || ux > uy;
     const std::uint64_t next = towardZero ? ux - 1 : ux + 1;
     const std::uint64_t exponent = next & 0x7ff0000000000000ull;
     if (exponent == 0x7ff0000000000000ull) {
-        errno = ERANGE;
+        _Feraise_nid_postfix(0x08);
         return x + x;
     }
     if (exponent == 0) {
-        errno = ERANGE;
+        volatile double product = x * x;
+        if (product != x) _Feraise_nid_postfix(0x10);
     }
     return std::bit_cast<double>(next);
 }

core/libs/tests/GuestMath.cpp

@@ -671,6 +671,64 @@ static void CheckRoundingNaNConsoleBits() {
     }
 }
 
+template<typename T, typename UInt>
+static void CheckNextafterEnvironment(T (APS5_VABI *next)(T, T), UInt sign, UInt one, UInt minimum, UInt maximum, UInt infinity) {
+    const RestoreFloatingPointState saved;
+    const auto savedControl = _mm_getcsr();
+    const struct {
+        UInt x, y, normalResult, dazResult;
+        unsigned normalFlags, dazFlags;
+        int normalError, dazError;
+    } cases[] = {
+        {0, one, 1, 0, 0x32, 0, 97, 97},
+        {sign, sign | one, sign | 1, 0, 0x32, 0, 97, 97},
+        {1, 0, 0, 0, 0x32, 0, ERANGE, 97},
+        {1, one, 2, 2, 0x32, 0, ERANGE, 97},
+        {sign | 1, sign | one, sign | 2, sign | 2, 0x32, 0, ERANGE, 97},
+        {minimum, 0, minimum - 1, minimum - 1, 0x30, 0x30, ERANGE, ERANGE},
+        {maximum, infinity, infinity, infinity, 0x28, 0x28, ERANGE, ERANGE},
+        {1, 1, 1, 1, 0x02, 0, 97, 97},
+        {one, one + 1, one + 1, one + 1, 0, 0, 97, 97},
+    };
+    for (unsigned denorm = 0; denorm != 4; ++denorm) {
+        const unsigned control = 0x1f80u | ((denorm & 1u) ? 0x40u : 0u) | ((denorm & 2u) ? 0x8000u : 0u);
+        const bool daz = (denorm & 1u) != 0;
+        for (const auto& item : cases) {
+            Require(std::fesetround(FE_TONEAREST) == 0);
+            Require(std::feclearexcept(FE_ALL_EXCEPT) == 0);
+            errno = 97;
+            _mm_setcsr(control);
+            const auto result = std::bit_cast<UInt>(next(std::bit_cast<T>(item.x), std::bit_cast<T>(item.y)));
+            const auto flags = _mm_getcsr() & 0x3fu;
+            const int exceptions = std::fetestexcept(FE_ALL_EXCEPT);
+            const int error = errno;
+            const auto expectedFlags = daz ? item.dazFlags : item.normalFlags;
+            Require(result == (daz ? item.dazResult : item.normalResult));
+            Require(flags == expectedFlags);
+            Require(exceptions == static_cast<int>(expectedFlags & FE_ALL_EXCEPT));
+            Require(error == (daz ? item.dazError : item.normalError));
+        }
+    }
+    constexpr int rounding[] = {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO};
+    for (const int mode : rounding) {
+        Require(std::fesetround(mode) == 0);
+        const unsigned control = (_mm_getcsr() & 0x6000u) | 0x1f80u;
+        errno = 97;
+        Require(std::feclearexcept(FE_ALL_EXCEPT) == 0);
+        _mm_setcsr(control);
+        Require(std::bit_cast<UInt>(next(std::bit_cast<T>(maximum), std::bit_cast<T>(infinity))) == (mode == FE_UPWARD ? infinity : maximum));
+        Require((_mm_getcsr() & 0x3fu) == 0x28u);
+        Require(errno == ERANGE);
+        errno = 97;
+        Require(std::feclearexcept(FE_ALL_EXCEPT) == 0);
+        _mm_setcsr(control);
+        Require(std::bit_cast<UInt>(next(std::bit_cast<T>(UInt{1}), std::bit_cast<T>(one))) == 2);
+        Require((_mm_getcsr() & 0x3fu) == 0x32u);
+        Require(errno == (mode == FE_UPWARD ? 97 : ERANGE));
+    }
+    _mm_setcsr(savedControl);
+}
+
 static void CheckAsinl() {
     const struct { std::uint64_t inputMantissa; std::uint16_t inputSignExponent; std::uint64_t outputMantissa; std::uint16_t outputSignExponent; } cases[] = {
         {0, 0x0000, 0, 0x0000}, {0, 0x8000, 0, 0x8000},
@@ -1816,6 +1874,8 @@ int main() {
     CheckNextafterConsoleBits();
     CheckTanhConsoleBits();
     CheckRoundingNaNConsoleBits();
+    CheckNextafterEnvironment(nextafterf_nid_postfix, 0x80000000u, 0x3f800000u, 0x00800000u, 0x7f7fffffu, 0x7f800000u);
+    CheckNextafterEnvironment(nextafter_nid_postfix, std::uint64_t{0x8000000000000000ull}, std::uint64_t{0x3ff0000000000000ull}, std::uint64_t{0x0010000000000000ull}, std::uint64_t{0x7fefffffffffffffull}, std::uint64_t{0x7ff0000000000000ull});
     CheckAsinl();
     Require(scalbn_nid_postfix(0.75, -2) == 0.1875);
     Require(scalbnf_nid_postfix(0.75f, 4) == 12.f);

docs/dev/debt/libc.md

@@ -117,7 +117,7 @@
 - [coshf](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows FreeBSD msun's [`e_coshf.c`](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_coshf.c), using the `expm1f` and `expf` paths already inlined for `sinhf` below; the Windows UCRT differed from the `libc.prx` copies in PPSA07923 and PPSA13934 on 500,554 and 500,614 of 20 million inputs in the #9707 survey. The exact system `libSceLibcInternal` behavior and the FMA-contracted `libc.prx` results remain unverified
 - [expm1](../../../core/libs/prx/libc/src/Math.cpp) (libc) does not compare floating-point exceptions, and its double inputs are checked by sampling, not exhaustively. It is FreeBSD msun's [`s_expm1.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_expm1.c) with `huge` read at run time, and sets `ERANGE` only when a finite argument overflows to infinity, as the titles' `libc.prx` does. It gives the same bits and `errno` as the `expm1` bundled with PPSA07923 and PPSA05916 on 320 million random and 47,222 structured arguments on a Windows host; the FMA-contracted build bundled with PPSA13934 differs on 43,433 of the random ones.
 - [log1p](../../../core/libs/prx/libc/src/Math.cpp) (libc) does not compare floating-point exceptions, and its double inputs are checked by sampling, not exhaustively. It is FreeBSD msun's [`s_log1p.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_log1p.c). As the titles' `libc.prx` does, it sets `EDOM` at msun's NaN return (arguments below -1, `-inf` and negative NaNs) and leaves `errno` alone for `log1p(-1)`, where it returns `-inf`. It gives the same bits, NaN payloads included, and `errno` as the `log1p` bundled with PPSA07923 and PPSA05916 on 320 million random and 68,216 structured arguments on a Windows host. The FMA-contracted build bundled with PPSA13934 differs on 73,682 of the random ones.
-- [nextafter and nextafterf](../../../core/libs/prx/libc/src/Math.cpp) (libc) do not compare floating-point exceptions, and their inputs are checked by sampling, not exhaustively. They are FreeBSD msun's [`s_nextafter.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_nextafter.c) and [`s_nextafterf.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_nextafterf.c), with `ERANGE` at msun's overflow and underflow returns as the titles' `libc.prx` sets it; a step away from zero returns early and sets none. The copies bundled with eleven titles give the same bits and `errno` on 4,000,784 sampled pairs per function on a Windows host.
+- [nextafter and nextafterf](../../../core/libs/prx/libc/src/Math.cpp) (libc) have no physical-console or unmasked-exception validation. A Linux CPU probe compares result bits, `errno` and masked floating-point flags with the unchanged math code in FAR: Changing Tides PPSA03766 bundled `libc.prx` (`nextafter` 0x51b30, `nextafterf` 0x51c80, `_Feraise` 0x11da0): 16,864 calls across both precisions, four rounding modes and all DAZ/FTZ combinations, not exhaustive inputs. Earlier Windows sampling of eleven bundled libraries compares only bits and `errno`; exception behavior outside the FAR module remains unverified.
 - [tanh](../../../core/libs/prx/libc/src/Math.cpp) (libc) does not compare floating-point exceptions, and its double inputs are checked by sampling, not exhaustively. It is FreeBSD msun's [`s_tanh.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_tanh.c) over this file's msun `Expm1`, with `tiny` read at run time; it sets no `errno`, as the titles' `libc.prx` does not. It gives the same bits as the `tanh` bundled with PPSA07923 and PPSA05916 on 16 million random and 68,211 structured arguments on a Windows host; the FMA-contracted build bundled with PPSA13934 differs on 1,054 of them.
 - [sinhf](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows the `sinhf` of the `libc.prx` bundled with PPSA28319 (disassembly): FreeBSD msun's [`e_sinhf.c`](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_sinhf.c) from before `__ldexp_expf`, which squares `expf(|x| / 2)` between `logf(FLT_MAX)` and the overflow threshold, calling msun's [`expm1f`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/s_expm1f.c) and [`expf`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/e_expf.c), whose paths for the arguments it passes are inlined here. Past the threshold Sony's error helper sets `errno` to `ERANGE` before returning `x * 1e37`; msun itself sets no `errno`. That code, run on a Windows host, gives the same bits as this implementation for every finite argument up to the threshold, of both signs; it was not run on Linux. That the system `libSceLibcInternal`, which Ruffle (PPSA68091) and RetroArch (PPSA99169) import while bundling a zero-return `libc.prx`, has the same code is assumed
 - [acoshf](../../../core/libs/prx/libc/src/Math.cpp) (libc) exact system `libSceLibcInternal` behavior and FMA-contracted copies remain unverified. Positive finite inputs use [msun](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_acoshf.c) range handling and this file's `logf`; the middle range divides -1 before addition under directed rounding. Bits, errno and SSE/x87 status and controls match the supplied PPSA03766/PPSA07923 copies on 236,592 positive finite calls (4,929 inputs, four rounding modes, DAZ/FTZ combinations and sticky flags); these samples are not exhaustive. Invalid-domain errno and SSE/x87 invalid flags are checked against two supplied bundled copies; those samples are not exhaustive either.