From ac3b199860766a5b781e24b18dce32208235328e Mon Sep 17 00:00:00 2001 From: JJenkx Date: Fri, 23 Oct 2026 13:09:21 -0400 Subject: [PATCH 1/2] fix(libc): compute atanh with FreeBSD msun's algorithm, as titles' libc.prx does atanh returned the host's atanh, whose last bit differs from the titles' copy on Windows. The titles' copy is FreeBSD msun's e_atanh.c over msun's log1p, and it sets EDOM for |x| > 1, the infinities and every NaN, and ERANGE for +-1. atanh now computes that algorithm over this file's log1p and sets errno there. Co-Authored-By: Claude Opus 5.5 --- core/libs/prx/libc/src/Math.cpp | 29 ++++++++++++++++++++++++++++- core/libs/tests/GuestMath.cpp | 26 ++++++++++++++++++++++++++ docs/dev/debt/libc.md | 1 + 3 files changed, 55 insertions(+), 1 deletion(-) diff --git a/core/libs/prx/libc/src/Math.cpp b/core/libs/prx/libc/src/Math.cpp index 9c34d18ed7..c8ea682e14 100644 --- a/core/libs/prx/libc/src/Math.cpp +++ b/core/libs/prx/libc/src/Math.cpp @@ -1234,7 +1234,34 @@ float APS5_VABI acoshf_nid_postfix(float x) { const float t = x - 1.0f; return AtanhfLog1p(t + std::sqrt(2.0f * t + t * t)); } -double APS5_VABI atanh_nid_postfix(double x) { return std::atanh(x); } +double APS5_VABI atanh_nid_postfix(double x) { + const double huge = 1e300; + static volatile double zero = 0.0; + const std::uint64_t bits = std::bit_cast(x); + const std::int32_t hx = static_cast(bits >> 32); + const std::uint32_t lx = static_cast(bits); + const std::int32_t ix = hx & 0x7fffffff; + if ((static_cast(ix) | ((lx | (0u - lx)) >> 31)) > 0x3ff00000u) { + errno = EDOM; + return (x - x) / (x - x); + } + if (ix == 0x3ff00000) { + errno = ERANGE; + return x / zero; + } + if (ix < 0x3e300000 && huge + x > 0.0) { + return x; + } + x = std::bit_cast((static_cast(static_cast(ix)) << 32) | lx); + double t; + if (ix < 0x3fe00000) { + t = x + x; + t = 0.5 * log1p_nid_postfix(t + t * x / (1.0 - x)); + } else { + t = 0.5 * log1p_nid_postfix((x + x) / (1.0 - x)); + } + return hx >= 0 ? t : -t; +} float APS5_VABI atanhf_nid_postfix(float x) { const std::uint32_t magnitude = std::bit_cast(x) & 0x7fffffffu; diff --git a/core/libs/tests/GuestMath.cpp b/core/libs/tests/GuestMath.cpp index 80f1217e86..34a0a628d9 100644 --- a/core/libs/tests/GuestMath.cpp +++ b/core/libs/tests/GuestMath.cpp @@ -1038,6 +1038,31 @@ static void CheckExpm1ScalingRounding() { } _mm_setcsr(savedCsr); } +static void CheckAtanhConsoleBits() { + constexpr struct { + std::uint64_t input; + std::uint64_t expected; + int error; + } cases[] = { + {0x0000000000000000ull, 0x0000000000000000ull, 0}, {0x8000000000000000ull, 0x8000000000000000ull, 0}, + {0x3e2fffffffffffffull, 0x3e2fffffffffffffull, 0}, {0x3e30000000000000ull, 0x3e30000000000000ull, 0}, + {0xbe30000000000000ull, 0xbe30000000000000ull, 0}, {0x3fdfffffffffffffull, 0x3fe193ea7aad030aull, 0}, + {0x3fe0000000000000ull, 0x3fe193ea7aad030aull, 0}, {0xbfe0000000000000ull, 0xbfe193ea7aad030aull, 0}, + {0x3fefffffffffffffull, 0x4032b708872320e2ull, 0}, {0xbfefffffffffffffull, 0xc032b708872320e2ull, 0}, + {0x3ff0000000000000ull, 0x7ff0000000000000ull, ERANGE}, {0xbff0000000000000ull, 0xfff0000000000000ull, ERANGE}, + {0x3ff0000000000001ull, 0xfff8000000000000ull, EDOM}, {0xc000000000000000ull, 0xfff8000000000000ull, EDOM}, + {0x7ff0000000000000ull, 0xfff8000000000000ull, EDOM}, {0xfff0000000000000ull, 0xfff8000000000000ull, EDOM}, + {0x7ff8000000000001ull, 0x7ff8000000000001ull, EDOM}, {0xfff8000000000002ull, 0xfff8000000000002ull, EDOM}, + {0x3f03a67d5567227bull, 0x3f03a67d558ea76eull, 0}, {0x3faa28c4621dd8d7ull, 0x3faa2e9a77f5a83full, 0}, + {0x3f4742fce55398bbull, 0x3f4742fd26e249ceull, 0}, {0x3f2cf9c68a297d4aull, 0x3f2cf9c69214cc60ull, 0}, + {0x3fc796d8c6ddd135ull, 0x3fc7dca409c1ba11ull, 0}, {0x3f73b85b718bdc38ull, 0x3f73b8656ddda145ull, 0}, + }; + for (const auto& item : cases) { + errno = 1234; + Require(std::bit_cast(atanh_nid_postfix(std::bit_cast(item.input))) == item.expected); + Require(errno == (item.error != 0 ? item.error : 1234)); + } +} static void CheckAsinl() { const struct { std::uint64_t inputMantissa; std::uint16_t inputSignExponent; std::uint64_t outputMantissa; std::uint16_t outputSignExponent; } cases[] = { @@ -2369,6 +2394,7 @@ int main() { 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}); CheckExpm1ScalingRounding(); + CheckAtanhConsoleBits(); CheckAsinl(); Require(scalbn_nid_postfix(0.75, -2) == 0.1875); Require(scalbnf_nid_postfix(0.75f, 4) == 12.f); diff --git a/docs/dev/debt/libc.md b/docs/dev/debt/libc.md index 568e1931c1..3ef9660417 100644 --- a/docs/dev/debt/libc.md +++ b/docs/dev/debt/libc.md @@ -142,6 +142,7 @@ - [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) 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. +- [atanh](../../../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 [`e_atanh.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/e_atanh.c) over this file's msun `log1p`. As the titles' `libc.prx` does, it sets `EDOM` for |x| > 1, the infinities and every NaN, and `ERANGE` for ±1, msun's two error returns. It gives the same bits, NaN payloads included, and `errno` as the `atanh` 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 2,778 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. - [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`). `tgammaf` delegates to the host's double `tgamma` as FreeBSD 11 does ([`s_tgammaf.c`](https://github.com/freebsd/freebsd-src/blob/stable/11/lib/msun/src/s_tgammaf.c)), then restores the sign of a zero result for negative non-integer inputs from the parity of `floor(-x)`; nine title-bundled `libc.prx` builds return that sign for underflowing inputs, where Windows UCRT can return the opposite sign ([#9707](https://github.com/boykopovar/AnyPS5/issues/9707)). 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 From 72a6683e64531e4fc948e6115af96037b367ab5e Mon Sep 17 00:00:00 2001 From: JJenkx Date: Wed, 28 Oct 2026 06:01:54 -0400 Subject: [PATCH 2/2] fix(libc): raise the titles' flags from atanh at its special inputs Titles' libc.prx raises FE_INVALID wherever atanh stores EDOM, a NaN argument included, where (x - x) / zero raised nothing for a NaN operand. The domain path now raises through _Feraise, which stores EDOM and raises the flag in both units. Results and errno are unchanged. CheckAtanhEdgeFlags covers it. Co-Authored-By: Claude Opus 5.5 --- core/libs/prx/libc/src/Math.cpp | 2 +- core/libs/tests/GuestMath.cpp | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/core/libs/prx/libc/src/Math.cpp b/core/libs/prx/libc/src/Math.cpp index c8ea682e14..1fb4fca6ee 100644 --- a/core/libs/prx/libc/src/Math.cpp +++ b/core/libs/prx/libc/src/Math.cpp @@ -1242,7 +1242,7 @@ double APS5_VABI atanh_nid_postfix(double x) { const std::uint32_t lx = static_cast(bits); const std::int32_t ix = hx & 0x7fffffff; if ((static_cast(ix) | ((lx | (0u - lx)) >> 31)) > 0x3ff00000u) { - errno = EDOM; + _Feraise_nid_postfix(0x01); return (x - x) / (x - x); } if (ix == 0x3ff00000) { diff --git a/core/libs/tests/GuestMath.cpp b/core/libs/tests/GuestMath.cpp index 34a0a628d9..220f103fdd 100644 --- a/core/libs/tests/GuestMath.cpp +++ b/core/libs/tests/GuestMath.cpp @@ -1038,6 +1038,13 @@ static void CheckExpm1ScalingRounding() { } _mm_setcsr(savedCsr); } +static void CheckAtanhEdgeFlags() { + for (const double input : {2.0, -2.0, std::numeric_limits::infinity(), std::numeric_limits::quiet_NaN(), -std::numeric_limits::quiet_NaN()}) { + errno = 0; + std::feclearexcept(FE_ALL_EXCEPT); + Require(std::isnan(atanh_nid_postfix(input)) && errno == EDOM && std::fetestexcept(FE_INVALID) != 0); + } +} static void CheckAtanhConsoleBits() { constexpr struct { std::uint64_t input; @@ -2395,6 +2402,7 @@ int main() { CheckNextafterEnvironment(nextafter_nid_postfix, std::uint64_t{0x8000000000000000ull}, std::uint64_t{0x3ff0000000000000ull}, std::uint64_t{0x0010000000000000ull}, std::uint64_t{0x7fefffffffffffffull}, std::uint64_t{0x7ff0000000000000ull}); CheckExpm1ScalingRounding(); CheckAtanhConsoleBits(); + CheckAtanhEdgeFlags(); CheckAsinl(); Require(scalbn_nid_postfix(0.75, -2) == 0.1875); Require(scalbnf_nid_postfix(0.75f, 4) == 12.f);