fix(libc): compute acosl with FreeBSD msun's ld80 code as the bundled libc.prx does

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

@@ -178,6 +178,11 @@ long double ClearLowMantissaWord(long double value) {
     return value;
 }
 
+long double SumAtCallerPrecision(long double high, long double low) {
+    volatile long double sum = high;
+    return sum + low;
+}
+
 constexpr long double AsinNumeratorCoefficients[] = {
     1.66666666666666666631e-01L, -4.16313987993683104320e-01L, 3.69068046323246813704e-01L, -1.36213932016738603108e-01L,
     1.78324189708471965733e-02L, -2.19216428382605211588e-04L, -7.10526623669075243183e-06L,
@@ -189,6 +194,7 @@ constexpr long double AsinDenominatorCoefficients[] = {
 constexpr long double HalfPiHigh = 1.57079632679489661926e+00L;
 constexpr long double HalfPiLow = -2.50827880633416601173e-20L;
 constexpr long double QuarterPiHigh = 7.85398163397448309628e-01L;
+const volatile long double Pi = 3.14159265358979323846264338327950280e+00L;
 
 long double AsinNumerator(long double t) {
     const auto& p = AsinNumeratorCoefficients;
@@ -549,7 +555,46 @@ float APS5_VABI fmaxf_nid_postfix(float x, float y) { return GuestFmax(x, y); }
 long double APS5_VABI fmaxl_nid_postfix(long double x, long double y) { return GuestFmax(x, y); }
 double APS5_VABI fdim_nid_postfix(double x, double y) { return std::fdim(x, y); }
 float APS5_VABI fdimf_nid_postfix(float x, float y) { return std::fdim(x, y); }
-long double APS5_VABI acosl_nid_postfix(long double x) { return std::acos(x); }
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+/* Adapted from FreeBSD release/11.4.0 lib/msun/src/e_acosl.c: https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/e_acosl.c */
+long double APS5_VABI acosl_nid_postfix(long double x) {
+    const ExtendedParts parts = SplitExtended(x);
+    const int exponent = parts.signExponent & 0x7fff;
+    const bool negative = (parts.signExponent & 0x8000) != 0;
+    if (exponent >= 0x3fff) {
+        if (exponent == 0x3fff && (parts.mantissa & 0x7fffffffffffffffull) == 0) return negative ? SumAtCallerPrecision(Pi, 2.0L * HalfPiLow) : 0.0L;
+        std::feraiseexcept(FE_INVALID);
+        errno = EDOM;
+        return (x - x) / (x - x);
+    }
+    if (exponent < 0x3ffe) {
+        if (exponent < 0x3fff - 65) return SumAtCallerPrecision(HalfPiHigh, HalfPiLow);
+        const long double z = x * x;
+        const long double r = AsinNumerator(z) / AsinDenominator(z);
+        return HalfPiHigh - (x - (HalfPiLow - x * r));
+    }
+    if (negative) {
+        const long double z = (1.0L + x) * 0.5L;
+        const long double s = std::sqrt(z);
+        const long double r = AsinNumerator(z) / AsinDenominator(z);
+        return Pi - 2.0L * (s + (r * s - HalfPiLow));
+    }
+    const long double z = (1.0L - x) * 0.5L;
+    const long double s = std::sqrt(z);
+    const long double high = ClearLowMantissaWord(s);
+    const long double correction = (z - high * high) / (s + high);
+    const long double r = AsinNumerator(z) / AsinDenominator(z);
+    return 2.0L * (high + (r * s + correction));
+}
 long double APS5_VABI acoshl_nid_postfix(long double x) { return std::acosh(x); }
 long double APS5_VABI asinhl_nid_postfix(long double x) { return std::asinh(x); }
 long double APS5_VABI atanhl_nid_postfix(long double x) { return std::atanh(x); }

core/libs/tests/GuestMath.cpp

@@ -1073,6 +1073,70 @@ static void CheckAsinl() {
     *__error_nid_postfix() = 0;
 }
 
+static bool ExtendedBitsAt53Bits(long double (APS5_VABI* function)(long double), std::uint64_t inputMantissa, std::uint16_t inputSignExponent,
+                                 std::uint64_t outputMantissa, std::uint16_t outputSignExponent) {
+    const long double input = MakeExtended(inputMantissa, inputSignExponent);
+    const std::uint16_t saved = X87ControlWord();
+    const auto control = static_cast<std::uint16_t>((saved & ~0x0300u) | 0x0200u);
+    asm volatile("fldcw %0" : : "m"(control) : "memory");
+    const long double result = function(input);
+    const std::uint16_t after = X87ControlWord();
+    asm volatile("fldcw %0" : : "m"(saved) : "memory");
+    return after == control && ExtendedBits(result, outputMantissa, outputSignExponent);
+}
+
+static void CheckAcosl() {
+    const struct { std::uint64_t inputMantissa; std::uint16_t inputSignExponent; std::uint64_t outputMantissa; std::uint16_t outputSignExponent; } cases[] = {
+        {0, 0x0000, 0xc90fdaa22168c235ull, 0x3fff}, {0, 0x8000, 0xc90fdaa22168c235ull, 0x3fff},
+        {0x8000000000000000ull, 0x3fbd, 0xc90fdaa22168c235ull, 0x3fff}, {0x8000000000000000ull, 0x3fbe, 0xc90fdaa22168c235ull, 0x3fff},
+        {0x8000000000000000ull, 0xbfbe, 0xc90fdaa22168c235ull, 0x3fff},
+        {0x8000000000000000ull, 0x3ffd, 0xa8b807060a5927f8ull, 0x3fff}, {0x8000000000000000ull, 0xbffd, 0xe967ae3e38785c72ull, 0x3fff},
+        {0xcccccccccccccccdull, 0x3ffb, 0xbc3d91671fac8055ull, 0x3fff},
+        {0x8000000000000000ull, 0x3ffe, 0x860a91c16b9b2c23ull, 0x3fff}, {0x8000000000000000ull, 0xbffe, 0x860a91c16b9b2c23ull, 0x4000},
+        {0xc000000000000000ull, 0x3ffe, 0xb9051c960ecaa429ull, 0x3ffe}, {0xc000000000000000ull, 0xbffe, 0x9ace937c9db6192bull, 0x4000},
+        {0xe666666666666666ull, 0x3ffe, 0xe6ecfc7c9174ca46ull, 0x3ffd}, {0xfd70a3d70a3d70a4ull, 0xbffe, 0xc000df0de008d77cull, 0x4000},
+        {0xffffffffffffffffull, 0x3ffe, 0xb504f333f9de6484ull, 0x3fdf}, {0xffffffffffffffffull, 0xbffe, 0xc90fdaa1c6e6489bull, 0x4000},
+        {0x8000000000000000ull, 0x3fff, 0, 0x0000}, {0x8000000000000000ull, 0xbfff, 0xc90fdaa22168c235ull, 0x4000},
+        {0xc1d51c773e6f53e1ull, 0x3ffd, 0x975d54d5f258ea8aull, 0x3fff}, {0xcdd3484f06240dc9ull, 0xbffe, 0xa04e976af0cf72f5ull, 0x4000},
+        {0xaabe336043197207ull, 0x3ffd, 0x9d8ae5f7320d25aaull, 0x3fff}, {0x848f5ad972ed8af4ull, 0x3ffe, 0x83647f790cc185e8ull, 0x3fff},
+    };
+    for (const auto& test : cases) {
+        *__error_nid_postfix() = 12345;
+        const long double result = acosl_nid_postfix(MakeExtended(test.inputMantissa, test.inputSignExponent));
+        if (!ExtendedBits(result, test.outputMantissa, test.outputSignExponent) || *__error_nid_postfix() != 12345) {
+            std::fprintf(stderr, "Guest acosl failed for %016llx:%04x\n", static_cast<unsigned long long>(test.inputMantissa), test.inputSignExponent);
+            std::abort();
+        }
+    }
+    const long double domainErrors[] = {
+        1.5L, -2.0L, MakeExtended(0x8000000000000001ull, 0x3fff), MakeExtended(0x8000000000000001ull, 0xbfff),
+        std::numeric_limits<long double>::max(), std::numeric_limits<long double>::infinity(), -std::numeric_limits<long double>::infinity(),
+        std::numeric_limits<long double>::quiet_NaN(),
+    };
+    for (const long double input : domainErrors) {
+        std::feclearexcept(FE_ALL_EXCEPT);
+        *__error_nid_postfix() = 0;
+        Require(std::isnan(acosl_nid_postfix(input)) && *__error_nid_postfix() == EDOM && std::fetestexcept(FE_INVALID) != 0);
+    }
+    const int rounding = std::fegetround();
+    std::fesetround(FE_DOWNWARD);
+    const long double downward[] = {acosl_nid_postfix(MakeExtended(0xa6f97bbf48a499f5ull, 0xbffe)), acosl_nid_postfix(-0.75L), acosl_nid_postfix(-1.0L)};
+    std::fesetround(rounding);
+    Require(ExtendedBits(downward[0], 0x9201704b764da194ull, 0x4000) && ExtendedBits(downward[1], 0x9ace937c9db6192aull, 0x4000) && ExtendedBits(downward[2], 0xc90fdaa22168c234ull, 0x4000));
+    *__error_nid_postfix() = 0;
+    const struct { std::uint64_t inputMantissa; std::uint16_t inputSignExponent; std::uint64_t outputMantissa; std::uint16_t outputSignExponent; } at53Bits[] = {
+        {0x8000000000000000ull, 0x3fbd, 0xc90fdaa22168c000ull, 0x3fff}, {0x8000000000000000ull, 0x3ffd, 0xa8b807060a592800ull, 0x3fff},
+        {0x8000000000000000ull, 0xbffd, 0xe967ae3e38786000ull, 0x3fff}, {0xc000000000000000ull, 0x3ffe, 0xb9051c960ecaa800ull, 0x3ffe},
+        {0xc000000000000000ull, 0xbffe, 0x9ace937c9db61800ull, 0x4000}, {0x8000000000000000ull, 0xbfff, 0xc90fdaa22168c000ull, 0x4000},
+    };
+    for (const auto& test : at53Bits) {
+        if (!ExtendedBitsAt53Bits(acosl_nid_postfix, test.inputMantissa, test.inputSignExponent, test.outputMantissa, test.outputSignExponent)) {
+            std::fprintf(stderr, "Guest acosl at 53 bits failed for %016llx:%04x\n", static_cast<unsigned long long>(test.inputMantissa), test.inputSignExponent);
+            std::abort();
+        }
+    }
+}
+
 static void CheckCbrtfConsoleBits() {
     constexpr std::uint32_t cases[][2] = {
         {0x00000123u, 0x28059605u}, {0x00400327u, 0x2a4b334bu}, {0x3f000555u, 0x3f4b32c7u}, {0x3fc00778u, 0x3f9287e5u},
@@ -2190,6 +2254,7 @@ int main() {
     CheckIntegerConversions();
     CheckCbrtfConsoleBits();
     CheckExtendedMath();
+    CheckAcosl();
     CheckTrigfConsoleBits();
     Require(atof_nid_postfix(" -12.5tail") == -12.5);
     char* end = nullptr;

docs/dev/debt/libc.md

@@ -47,8 +47,8 @@
 - [reallocalign](../../../core/libs/prx/libc/src/FileAndHeap.cpp) (libc) - `reallocalign(ptr, 0, alignment)` for a block without that alignment throws, where the console's heap in a payload process on PS5 firmware 13.60 returned a different block with that alignment and 0 usable bytes, `errno` untouched (a 64-aligned block given 256); whether it frees the old block is not known. There, as here, `reallocalign(ptr, 0, 16)` on a 256-aligned block returned it unchanged, `reallocalign(NULL, 0, 256)` a 256-aligned block with 0 usable bytes, and alignment 3 a moved block that kept its bytes
 - [reallocf](../../../core/libs/prx/libc/src/FileAndHeap.cpp) (libSceLibcInternal, implemented in libc) - follows FreeBSD 11 `lib/libc/stdlib/reallocf.c` on top of this libc's `realloc`, as a payload process on a PS5 (firmware 13.60) does with the default heap (null and 0, size 0, 2^40 with `errno` 22 and the block freed, `SIZE_MAX` with `errno` untouched): size 0 does what `realloc` does here, keeping the block with 0 usable bytes with the default heap and freeing it and returning null with an application allocator, so the block is not freed twice; when `realloc` fails, the block is freed through the same allocator and null is returned with the `errno` that `realloc` left. Other failures (a commit failure that is not out of memory, an unknown pointer) throw, as they do for `realloc`
 - [sceLibcMspaceReallocalign](../../../core/libs/prx/libc/src/Mspace.cpp) (libSceLibcInternal, implemented in libc) takes `(void* heap, void* block, size_t alignment, size_t size)`, the order of the console's `libSceLibcInternal.sprx` (SDK 4.5x build of `J03950001`, 0x13880, whose internal routine at 0x12390 treats the third argument as the alignment and, for a null block, passes it to the routine `sceLibcMspaceMemalign` uses). Measured on a PS5 (firmware 13.60) from a homebrew payload: a null block with alignment 4096 and size 100 returns a 4096-aligned block; an alignment of 0 or 3 succeeds and 100 gives a 128-aligned block, as the routine raises it to at least 32 and to the next power of two; an alignment of 2^63 or more returns null, keeps the block and leaves `errno` alone; and a size of 0 returns a block (the same one when it is already suitable) instead of freeing it. A block that is aligned and large enough is kept where it is, without giving back the bytes past the new size as the console's in-place resize does
-- [acosl](../../../core/libs/prx/libc/src/Math.cpp) (libc) - computed by the host's math library, like the other libc math functions here, not by the console's (FreeBSD msun ld80 code in Sony's `libc.prx`): results can differ from the console's in the last bit, and on Linux glibc sets `errno` to `EDOM` for `acosl` outside [-1, 1] where msun only raises the floating-point exception. `long double` is the x87 80-bit format, as Sony's `frexpl` uses it (`fldt 0x8(%rsp)`, result in `st(0)`)
-- [atanl, sinl, cosl, tanl, sinhl, coshl, tanhl, expl, exp2l, logl, log10l and powl](../../../core/libs/prx/libc/src/Math.cpp) (libc) - computed by the host's math library, as `acosl` is, not by the console's FreeBSD msun ld80 code: results can differ from the console's in the last bit (glibc's `log10l(1000)` is not exactly 3), and on Linux glibc sets `errno` to `EDOM` or `ERANGE` where msun only raises the floating-point exception. `ldexpl` and `truncl` are exact operations and give the console's result; `logb` and `remainder` give the console's value too, but on Linux glibc sets `errno` (`ERANGE` for `logb(0)`, `EDOM` for `remainder` by zero or of an infinity) where msun only raises the floating-point exception
+- [acosl](../../../core/libs/prx/libc/src/Math.cpp) (libc) runs at the caller's x87 precision, as the binary does, so on a guest thread that starts at 53 bits (see `asinl`) it returns the binary's 53-bit result, not the console's 64-bit one. It follows the `acosl` of the `libc.prx` bundled with PPSA21564 (0x3e5a0): FreeBSD msun's [`e_acosl.c`](https://github.com/freebsd/freebsd-src/blob/9f0432041b8fb69832f10c32da50a67f34feeb61/lib/msun/src/e_acosl.c) with the ld80 `invtrig.c` coefficients that `asinl` uses, where `_Feraise(FE_INVALID)` raises `FE_INVALID` and sets `errno` to `EDOM` for |x| > 1, the infinities and NaN before `(x - x) / (x - x)`. 20 million random arguments in (-1, 1) give the same bits as that code run natively on Linux at 64-, 53- and 24-bit precision; the host `acosl` differed on 1,603,469 of them at 64 bits. Other SDK versions and the system `libSceLibcInternal` are not compared
+- [atanl, sinl, cosl, tanl, sinhl, coshl, tanhl, expl, exp2l, logl, log10l and powl](../../../core/libs/prx/libc/src/Math.cpp) (libc) - computed by the host's math library, not by the console's FreeBSD msun ld80 code: results can differ from the console's in the last bit (glibc's `log10l(1000)` is not exactly 3), and on Linux glibc sets `errno` to `EDOM` or `ERANGE` where msun only raises the floating-point exception. `ldexpl` and `truncl` are exact operations and give the console's result; `logb` and `remainder` give the console's value too, but on Linux glibc sets `errno` (`ERANGE` for `logb(0)`, `EDOM` for `remainder` by zero or of an infinity) where msun only raises the floating-point exception
 - [Dinkumware atomics](../../../core/libs/prx/libc/src/DinkumwareAtomics.cpp) (libc) - every operation is sequentially consistent whatever memory order it is given, which is never weaker than the console. Matching the `libc.prx` bundled with PPSA08427 (that binary run side by side on every order from 0 to 7), `_Atomic_load_8` returns 0 for a release, acq_rel or invalid order and the 8-byte exchange and fetch operations return 0 without touching memory for an order above `memory_order_seq_cst`; the generic `_Atomic_copy`, `_Atomic_exchange` and `_Atomic_compare_exchange_*` lock the byte flag passed first, as the console does
 - [_Fetch_{and,or,xor}_seq_cst_{1,2,4}](../../../core/libs/prx/libc/src/DinkumwareAtomics.cpp) (libc) follow the two-argument unsigned 8-, 16- and 32-bit helpers in Dinkumware [VC12 xatomic.h](https://github.com/ganboing/crt.vc12/blob/f1d47126bf13d1c9dccaa6c53b72217fd0474968/src/xatomic.h), using sequentially consistent atomic read-modify-write operations and returning the original value. The firmware 3.20 export inventory supplies the names and NIDs; the reference defines inline namespace helpers, so the equivalence to the console exports and their calling convention are inferred, not console-measured. A live, suitably aligned object is required as in the reference; invalid-pointer and misalignment fault behavior on the console is unverified. Windows tests cover raw and patched exports, values, neighboring storage and concurrent updates; a finite ordering litmus cannot prove the memory model.
 - [std_execute_once](../../../core/libs/prx/libc/Export.cpp) (libc) - assumed to be the plain-named twin of `std::_Execute_once` (`_ZSt13_Execute_onceRSt9once_flagPFiPvS1_PS1_ES1_`, the callee of Dinkumware `call_once`): it delegates there and gets the [thread.once.callonce] semantics already verified for that symbol (one invocation per flag, concurrent callers block until it completes, and a failing or throwing invocation leaves the flag unready so the next caller runs it). The name equivalence is not confirmed on a title, and the callback convention (called with (null, arg, null), non-zero return meaning success) is the twin's, unverified for this symbol