Merge pull request #10665 from JJenkx/fix/libc-acoshf-positive-rounding-20261024

fix(libc): match acoshf positive-range evaluation

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

@@ -1045,10 +1045,10 @@ float APS5_VABI acoshf_nid_postfix(float x) {
     if ((bits & 0x80000000u) != 0 || magnitude < 0x3f800000u) return (x - x) / (x - x);
     if (magnitude >= 0x7f800000u) return x + x;
     if (magnitude == 0x3f800000u) return 0.0f;
-    if (magnitude >= 0x4d800000u) return std::log(x) + ln2;
+    if (magnitude >= 0x4d800000u) return logf_nid_postfix(x) + ln2;
     if (magnitude > 0x40000000u) {
         const float square = x * x;
-        return std::log(2.0f * x - 1.0f / (x + std::sqrt(square - 1.0f)));
+        return logf_nid_postfix(2.0f * x + (-1.0f / (x + std::sqrt(square - 1.0f))));
     }
     const float t = x - 1.0f;
     return AtanhfLog1p(t + std::sqrt(2.0f * t + t * t));

core/libs/tests/GuestMath.cpp

@@ -937,6 +937,53 @@ static void CheckAcoshfMsunBoundaries() {
     Require(std::isinf(acoshf_nid_postfix(INFINITY)) && acoshf_nid_postfix(INFINITY) > 0.0f);
     Require(std::isnan(acoshf_nid_postfix(-INFINITY)));
     Require(std::isnan(acoshf_nid_postfix(std::numeric_limits<float>::quiet_NaN())));
+    const RestoreFloatingPointState saved;
+    const unsigned savedCsr = _mm_getcsr();
+    const struct {
+        std::uint32_t input;
+        std::uint32_t result[4];
+        unsigned raised;
+    } cases[] = {
+        {0x3f800000u, {0x00000000u, 0x00000000u, 0x00000000u, 0x00000000u}, 0u},
+        {0x3f800001u, {0x39ffffffu, 0x39ffffffu, 0x3a000001u, 0x39ffffffu}, 32u},
+        {0x3fffffffu, {0x3fa89213u, 0x3fa89212u, 0x3fa89214u, 0x3fa89212u}, 32u},
+        {0x40000000u, {0x3fa89214u, 0x3fa89213u, 0x3fa89214u, 0x3fa89213u}, 32u},
+        {0x40000001u, {0x3fa89215u, 0x3fa89214u, 0x3fa89215u, 0x3fa89214u}, 32u},
+        {0x40000002u, {0x3fa89216u, 0x3fa89215u, 0x3fa89217u, 0x3fa89215u}, 32u},
+        {0x40000003u, {0x3fa89217u, 0x3fa89217u, 0x3fa89218u, 0x3fa89217u}, 32u},
+        {0x4d7fffffu, {0x41a0cf66u, 0x41a0cf65u, 0x41a0cf66u, 0x41a0cf65u}, 32u},
+        {0x4d800000u, {0x41a0cf66u, 0x41a0cf64u, 0x41a0cf66u, 0x41a0cf64u}, 32u},
+        {0x4d800001u, {0x41a0cf66u, 0x41a0cf65u, 0x41a0cf67u, 0x41a0cf65u}, 32u},
+        {0x4d8050f2u, {0x41a0d473u, 0x41a0d472u, 0x41a0d474u, 0x41a0d472u}, 32u},
+        {0x7f7fffffu, {0x42b2d4fcu, 0x42b2d4fbu, 0x42b2d4fdu, 0x42b2d4fbu}, 32u},
+    };
+    const int rounding[] = {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO};
+    float (APS5_VABI * volatile call)(float) = acoshf_nid_postfix;
+    for (unsigned mode = 0; mode != 4; ++mode) {
+        for (unsigned denorm = 0; denorm != 4; ++denorm) {
+            const unsigned control = 0x1f80u | (mode << 13u) | ((denorm & 1u) ? 0x40u : 0u) | ((denorm & 2u) ? 0x8000u : 0u);
+            for (std::uint16_t sticky : {0u, 5u, 32u}) {
+                for (const auto& item : cases) {
+                    Require(std::fesetenv(FE_DFL_ENV) == 0 && std::fesetround(rounding[mode]) == 0);
+                    Require(fesetexceptflag_nid_postfix(&sticky, 0x3f) == 0);
+                    _mm_setcsr(control | sticky);
+                    GuestFloatingPointEnvironment before;
+                    Require(fegetenv_nid_postfix(&before) == 0 && (before.status & 0x3fu) == sticky);
+                    errno = 1234;
+                    const auto result = std::bit_cast<std::uint32_t>(call(std::bit_cast<float>(item.input)));
+                    const int error = errno;
+                    GuestFloatingPointEnvironment after;
+                    Require(fegetenv_nid_postfix(&after) == 0);
+                    if (result != item.result[mode] || after.mxcsr != (control | sticky | item.raised) || (after.status & 0x3fu) != sticky || after.control != before.control || error != 1234) {
+                        std::fprintf(stderr, "acoshf positive: input=%08llx mode=%u denorm=%u sticky=%u bits=%08llx expected=%08llx csr=%x x87=%x errno=%d\n", static_cast<unsigned long long>(item.input), mode, denorm, sticky, static_cast<unsigned long long>(result), static_cast<unsigned long long>(item.result[mode]), after.mxcsr, after.status, error);
+                        Require(false);
+                    }
+                }
+            }
+        }
+    }
+    _mm_setcsr(savedCsr);
+
 }
 
 static void CheckCoshfMsunBoundaries() {

docs/dev/debt/libc.md

@@ -113,7 +113,7 @@
 - [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.
 - [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) uses FreeBSD msun's [`e_acoshf.c`](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_acoshf.c) range handling; the Windows UCRT differed from the `libc.prx` copies in PPSA07923 and PPSA13934 on 216,017 and 216,024 of 20 million inputs in the #9707 survey. The invalid-domain path raises `FE_INVALID` without setting `errno` as the source does. The exact system `libSceLibcInternal` behavior and FMA-contracted results remain unverified
+- [acoshf](../../../core/libs/prx/libc/src/Math.cpp) (libc) does not reproduce bundled invalid-domain errno or x87 exception flags (#10610); 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.
 - [acoshl, asinhl, atanhl, atan2l, cbrtl, erfl, erfcl, expm1l, hypotl, log1pl, log2l, lgammal, tgammal, j0f, j1f, jnf, y0f, y1f and ynf](../../../core/libs/prx/libc/src/Math.cpp) (libc) return the host C library's results (the float Bessel functions through its double ones), not those of the console's own code, so the last bits can differ (`y0f`, `y1f`, `ynf` follow the console at `+-0`). As in the console's `libSceLibcInternal.sprx` (SDK 4.5x build of `J03950001`, 0x56410), `lgammal` does not store the sign in `signgam`; `hypotl` returns infinity for an infinite argument even when the other is a NaN, as C11 F.10.4.3 requires. [hypot3, hypot3f and hypot3l](../../../core/libs/prx/libc/src/Math.cpp) follow that module (0x51D30, 0x520D0, 0x52420): the first NaN of x, y and z is returned unchanged, any infinity gives positive infinity, a zero argument gives the two-argument hypot of the other two, and otherwise the largest magnitude is combined with the hypot of the other two
 - [asinl](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows the `asinl` of the `libc.prx` bundled with PPSA28319 (disassembly), FreeBSD msun's [`e_asinl.c`](https://github.com/freebsd/freebsd-src/blob/9f0432041b8fb69832f10c32da50a67f34feeb61/lib/msun/src/e_asinl.c) with the ld80 [`invtrig.c`](https://github.com/freebsd/freebsd-src/blob/9f0432041b8fb69832f10c32da50a67f34feeb61/lib/msun/ld80/invtrig.c) coefficients, which the binary holds as the same 80-bit constants. For |x| > 1, the infinities and NaN, Sony's error helper raises `FE_INVALID` and sets `errno` to `EDOM` before returning `(x - x) / (x - x)`, so a NaN argument also sets `EDOM`; msun itself sets no `errno`. 20 million random arguments in (-1, 1) give the same bits as that code run on a Windows host. The result needs the x87 precision control at 64 bits, as FreeBSD amd64 starts threads (control word `0x37f`), but guest threads do not set it: on Windows a thread started by `_beginthreadex` runs at 53 bits (`0x27f`), where `asinl(0.75)` returns mantissa `0xd91a98ae3406e000` instead of `0xd91a98ae3406e041`. The system `libSceLibcInternal` that RetroArch (PPSA99169) imports is assumed to have the same code
 - [atanhf](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows the `atanhf` of the `libc.prx` bundled with PPSA28319 (disassembly): FreeBSD msun's [`e_atanhf.c`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/e_atanhf.c) calling msun's [`log1pf`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/s_log1pf.c), whose paths for the positive arguments it passes are inlined here; the exported `log1pf` is still the host's. For |x| > 1, the infinities and NaN, Sony's error helper raises `FE_INVALID` and sets `errno` to `EDOM` before returning `(x - x) / (x - x)`, and for |x| = 1 it raises `FE_DIVBYZERO` and sets `errno` to `ERANGE` before returning `x / 0`; msun itself sets no `errno`. That code, run on a Windows host, gives the same bits, `errno` and floating-point exceptions as this implementation for all 2^32 arguments; it was not run on Linux. That the system `libSceLibcInternal`, which Dusklight (PPSA99305) imports while bundling a zero-return `libc.prx`, has the same code is assumed