diff --git a/core/libs/prx/libc/src/Math.cpp b/core/libs/prx/libc/src/Math.cpp index c53eb092d85d..74a32343507f 100644 --- a/core/libs/prx/libc/src/Math.cpp +++ b/core/libs/prx/libc/src/Math.cpp @@ -352,6 +352,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, @@ -363,6 +368,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; @@ -926,7 +932,46 @@ float APS5_VABI fdimf_nid_postfix(float x, float y) { const float difference = x - y; return difference > 0.0f ? difference : 0.0f; } -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); } diff --git a/core/libs/tests/GuestMath.cpp b/core/libs/tests/GuestMath.cpp index d62bd8577efc..2624db5c303f 100644 --- a/core/libs/tests/GuestMath.cpp +++ b/core/libs/tests/GuestMath.cpp @@ -1697,6 +1697,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((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(test.inputMantissa), test.inputSignExponent); + std::abort(); + } + } + const long double domainErrors[] = { + 1.5L, -2.0L, MakeExtended(0x8000000000000001ull, 0x3fff), MakeExtended(0x8000000000000001ull, 0xbfff), + std::numeric_limits::max(), std::numeric_limits::infinity(), -std::numeric_limits::infinity(), + std::numeric_limits::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(test.inputMantissa), test.inputSignExponent); + std::abort(); + } + } +} + static void CheckCbrtfConsoleBits() { constexpr std::uint32_t cases[][2] = { {0x00000123u, 0x28059605u}, {0x00400327u, 0x2a4b334bu}, {0x3f000555u, 0x3f4b32c7u}, {0x3fc00778u, 0x3f9287e5u}, @@ -3166,6 +3230,7 @@ int main() { CheckCbrtfConsoleBits(); CheckCbrtlConsoleBits(); CheckExtendedMath(); + CheckAcosl(); CheckTrigfConsoleBits(); Require(atof_nid_postfix(" -12.5tail") == -12.5); char* end = nullptr; diff --git a/docs/dev/debt/libc.md b/docs/dev/debt/libc.md index bb79f170717f..9d2ca76e8a53 100644 --- a/docs/dev/debt/libc.md +++ b/docs/dev/debt/libc.md @@ -71,8 +71,8 @@ - [recallocarray](../../../core/libs/prx/libc/src/FileAndHeap.cpp) (libc) - adapted from OpenBSD `lib/libc/stdlib/recallocarray.c` at commit `2cb7cdfc21` (ISC licensed, Otto Moerbeek): console behavior is unmeasured. Delegates to `calloc` on null pointer, checks multiplication overflow on new and old counts setting `ENOMEM` (12) and `EINVAL` (22), zeroes small shrunken tails in place (under half the old size and modeled 16 KiB page size), taking old count on trust as in reference (a caller passing `oldnmemb` larger than the block zeroes past it), and otherwise allocates new storage, copies elements, zeroes new memory or clears the old allocation, and frees old pointer - [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 - [frexp](../../../core/libs/prx/libc/src/Math.cpp) (libc) - follows the two bundled PPSA03766/PPSA07923 functions: exponent-bit decomposition, conditional double-precision subnormal scaling, zero exponent for zero/NaN/infinity, and unchanged signaling NaNs. Masked-exception observations cover rounding, DAZ/FTZ and sticky flags; unmasked traps and other SDK versions are unverified. -- [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}, _Fetch_add_8 and _Fetch_{and,or,xor}_8](../../../core/libs/prx/libc/src/DinkumwareAtomics.cpp) (libc) follow the two-argument unsigned 8-, 16-, 32- and 64-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