Merge pull request #9668 from rooklz/fix/libc-acosf-msun

fix(libc): compute acosf with FreeBSD msun's algorithm, as the console's libc.prx does

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

@@ -69,6 +69,13 @@ float SinhfExpm1(float x) {
     return (oneMinusTwoPowMinusK - (e - x)) * twoPowK;
 }
 
+volatile float AcosfPio2Low = 7.5497894159e-08f;
+
+float AcosfRational(float z) {
+    constexpr float P0 = 1.6666586697e-01f, P1 = -4.2743422091e-02f, P2 = -8.6563630030e-03f, Q1 = -7.0662963390e-01f;
+    return z * (P0 + z * (P1 + z * P2)) / (1.0f + z * Q1);
+}
+
 float SinhfExp(float x) {
     constexpr float Ln2High = 6.9314575195e-01f, Ln2Low = 1.4286067653e-06f, InverseLn2 = 1.4426950216e+00f;
     constexpr float P1 = 1.6666625440e-1f, P2 = -2.7667332906e-3f;
@@ -310,10 +317,23 @@ float APS5_VABI acosf_nid_postfix(float x) {
     if (std::fabs(x) > 1.0f) {
         return std::bit_cast<float>(0xffc00000u);
     }
-    const int saved_errno = errno;
-    const float result = std::acos(x);
-    errno = saved_errno;
-    return result;
+    constexpr float Pi = 3.1415925026e+00f, Pio2High = 1.5707962513e+00f;
+    const std::uint32_t magnitude = std::bit_cast<std::uint32_t>(x) & 0x7fffffffu;
+    if (magnitude == 0x3f800000u) return std::signbit(x) ? Pi + 2.0f * AcosfPio2Low : 0.0f;
+    if (magnitude < 0x3f000000u) {
+        if (magnitude <= 0x32800000u) return Pio2High + AcosfPio2Low;
+        return Pio2High - (x - (AcosfPio2Low - x * AcosfRational(x * x)));
+    }
+    if (std::signbit(x)) {
+        const float z = (1.0f + x) * 0.5f;
+        const float s = std::sqrt(z);
+        return Pi - 2.0f * (s + (AcosfRational(z) * s - AcosfPio2Low));
+    }
+    const float z = (1.0f - x) * 0.5f;
+    const float s = std::sqrt(z);
+    const float high = std::bit_cast<float>(std::bit_cast<std::uint32_t>(s) & 0xfffff000u);
+    const float correction = (z - high * high) / (s + high);
+    return 2.0f * (high + (AcosfRational(z) * s + correction));
 }
 float APS5_VABI atan2f_nid_postfix(float y, float x) { return std::atan2(y, x); }
 float APS5_VABI hypotf_nid_postfix(float x, float y) { return static_cast<float>(std::hypot(static_cast<double>(x), static_cast<double>(y))); }

core/libs/tests/GuestMath.cpp

@@ -1315,6 +1315,15 @@ static bool Near(long double actual, long double expected) {
     return std::fabs(actual - expected) <= std::fabs(expected) * 1e-17L;
 }
 
+static void CheckAcosfMsun() {
+    Require(std::bit_cast<std::uint32_t>(acosf_nid_postfix(std::bit_cast<float>(0x3d58d031u))) == 0x3fc2488au);
+    Require(std::bit_cast<std::uint32_t>(acosf_nid_postfix(std::bit_cast<float>(0x3f00b99au))) == 0x3f859f4fu);
+    Require(std::bit_cast<std::uint32_t>(acosf_nid_postfix(std::bit_cast<float>(0xbf00b99au))) == 0x40064033u);
+    Require(std::bit_cast<std::uint32_t>(acosf_nid_postfix(1.0f)) == 0u);
+    Require(std::bit_cast<std::uint32_t>(acosf_nid_postfix(-1.0f)) == 0x40490fdbu);
+    Require(std::bit_cast<std::uint32_t>(acosf_nid_postfix(std::bit_cast<float>(0x32800000u))) == 0x3fc90fdbu);
+}
+
 static void CheckLongDoubleFunctions() {
     const long double infinity = std::numeric_limits<long double>::infinity();
     const long double nan = std::numeric_limits<long double>::quiet_NaN();
@@ -1740,6 +1749,7 @@ int main() {
     CheckPowfSpecialValue();
     Require(abs_nid_postfix(0) == 0);
     CheckAcosfErrno();
+    CheckAcosfMsun();
     Require(abs_nid_postfix(42) == 42);
     Require(abs_nid_postfix(-42) == 42);
     Require(abs_nid_postfix(std::numeric_limits<int>::max()) == std::numeric_limits<int>::max());

docs/dev/debt/libc.md

@@ -116,6 +116,7 @@
 - [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
+- [acosf](../../../core/libs/prx/libc/src/Math.cpp) (libc) - for |x| > 1 and NaN the console stores `EDOM` and raises `FE_INVALID` before returning `(x - x) / (x - x)`; here those return the same NaNs and leave `errno` as it was (#3688). For |x| <= 1 it follows the `acosf` of the `libc.prx` bundled with PPSA28319 (disassembly): FreeBSD msun's [`e_acosf.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/e_acosf.c), with the same constants, `pio2_lo` read at run time and the `0xfffff000` mask on `sqrtf(z)`, which the binary holds. Titles ship different builds of that `libc.prx`: PPSA13934's is FMA-contracted and differs from this one by 1 ulp on 129,170 inputs (#9668), so no single port matches every title's last bit.
 - [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
 - [j0, j1, jn, y0, y1 and yn](../../../core/libs/prx/libc/src/Bessel.cpp) (libc) call the host's `sin`, `cos`, `sincos` and `log`, so the last bits can differ from the console. The rest follows the system `libc.prx` built as `J03093718` (disassembly), which given the console's own `sin`, `cos`, `sincos` and `log` gave the same bits and exception flags for 90 million arguments: msun's [`e_j0.c`, `e_j1.c` and `e_jn.c`](https://github.com/freebsd/freebsd-src/tree/release/10.0.0/lib/msun/src) with each multiply-add the console fuses done as one FMA. As in that release, before FreeBSD commit `186f620727`, `y0`, `y1` and `yn` return `-inf` for ±0 and `0x7ff8000000000000` for a negative argument without raising an exception. The order wraps at 32 bits as on the console: `jn(1 << 30, x)` never returns for 2^-29 <= |x| < 2^30
 - [getenv, setenv, unsetenv, putenv](../../../core/libs/prx/libc/src/Environment.cpp) (libc) keep one guest environment that starts empty, like the `environ` libkernel exports, and never reads or changes the host's. A title sees only the variables it sets itself; host variables such as `HOME` and `PATH` name host paths a guest cannot use (PS5MC, PPSA99420, exits with `XBMC_Run returned -1` when it gets the host's `HOME`). Whether the console starts a title with any variables set is not measured