Merge pull request #10786 from 00xDisaster/feat/libc-dinkumware-math

feat(libc): implement _Exp and _FExp

core/libs/prx/libc/include/Api.hpp

@@ -625,6 +625,8 @@ float APS5_VABI __powisf2_nid_postfix(float base, int exponent);
 double APS5_VABI __powidf2_nid_postfix(double base, int exponent);
 float APS5_VABI expf_nid_postfix(float x);
 float APS5_VABI exp2f_nid_postfix(float x);
+short APS5_VABI _Exp_nid_postfix(double* value, double y, short offset);
+short APS5_VABI _FExp_nid_postfix(float* value, float y, short offset);
 float APS5_VABI logf_nid_postfix(float x);
 float APS5_VABI log2f_nid_postfix(float x);
 double APS5_VABI log10_nid_postfix(double x);

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

@@ -1,9 +1,13 @@
 #include "prx/libc/include/Api.hpp"
 #include <algorithm>
 #include <cctype>
+#include <cerrno>
+#include <cfenv>
 #include <cmath>
 #include <cstdint>
 #include <cstdlib>
+#include <limits>
+#include <stdexcept>
 #include <string>
 
 #include "prx/libc/include/General.hpp"
@@ -73,6 +77,67 @@ TInteger ParsedInteger(const char* text, char** end, int base, int* error, TInte
     return value;
 }
 
+constexpr short Denormal = -2, Finite = -1, Zero = 0, Infinite = 1, NotANumber = 2;
+
+void RaiseLibcError(int exceptions) {
+    if ((exceptions & (FE_OVERFLOW | FE_UNDERFLOW)) != 0) exceptions |= FE_INEXACT;
+    std::feraiseexcept(exceptions);
+    const bool range = (exceptions & FE_INVALID) == 0 &&
+        (exceptions & (FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW)) != 0;
+    errno = range ? ERANGE : EDOM;
+}
+
+template<typename TFloat>
+short ClassifyValue(TFloat value) {
+    switch (std::fpclassify(value)) {
+    case FP_NAN:
+        return NotANumber;
+    case FP_INFINITE:
+        return Infinite;
+    case FP_ZERO:
+        return Zero;
+    case FP_SUBNORMAL:
+        return Denormal;
+    default:
+        return Finite;
+    }
+}
+
+template<typename TFloat>
+short ScaleProduct(TFloat& value, TFloat factor, TFloat multiplier, long exponent) {
+    std::fexcept_t saved;
+    std::fegetexceptflag(&saved, FE_ALL_EXCEPT);
+    const int savedErrno = errno;
+    TFloat result = factor * multiplier;
+    if (!std::isfinite(result) || result == static_cast<TFloat>(0)) {
+        int factorExponent = 0;
+        int multiplierExponent = 0;
+        const TFloat factorMantissa = std::frexp(factor, &factorExponent);
+        const TFloat multiplierMantissa = std::frexp(multiplier, &multiplierExponent);
+        const long long combinedExponent = static_cast<long long>(factorExponent) + multiplierExponent + exponent;
+        result = std::scalbn(factorMantissa * multiplierMantissa, static_cast<int>(combinedExponent));
+    } else {
+        result = std::scalbn(result, static_cast<int>(exponent));
+    }
+    value = result;
+    const short kind = ClassifyValue(result);
+    const int raised = std::fetestexcept(FE_ALL_EXCEPT) & ~static_cast<int>(saved) & FE_ALL_EXCEPT;
+    std::fesetexceptflag(&saved, FE_ALL_EXCEPT);
+    errno = savedErrno;
+    if (kind == Infinite) {
+        RaiseLibcError(FE_OVERFLOW);
+        return Infinite;
+    }
+    if (kind == Zero) {
+        RaiseLibcError(FE_UNDERFLOW);
+        return Zero;
+    }
+    if (kind == Denormal) return Finite;
+    if (kind == NotANumber) return NotANumber;
+    if ((raised & FE_INEXACT) != 0) std::feraiseexcept(FE_INEXACT);
+    return Finite;
+}
+
 }
 
 extern "C" {
@@ -131,6 +196,200 @@ long long APS5_VABI _Stoll_nid_postfix(const char* text, char** end, int base) {
     return _Stollx_nid_postfix(text, end, base, nullptr);
 }
 
+short APS5_VABI _Exp_nid_postfix(double* value, double y, short offset) {
+    if (value == nullptr) throw std::invalid_argument("_Exp: null value");
+    std::fexcept_t entryFlags;
+    std::fegetexceptflag(&entryFlags, FE_ALL_EXCEPT);
+    const int entryErrno = errno;
+    const double x = *value;
+    constexpr double HugeExp = 1842.0;
+    constexpr double Eps = std::numeric_limits<double>::epsilon();
+    constexpr double InverseLn2 = 1.4426950408889634073599246810018921;
+    constexpr double Ln2High = 22713.0 / 32768.0;
+    constexpr double Ln2Low = 1.4286068203094172321214581765680755e-6;
+    constexpr double P1 = 420.30235984910635;
+    constexpr double P2 = 15132.70094680474802;
+    constexpr double Q0 = 30.01511290683317;
+    constexpr double Q1 = 3362.72154416553028;
+    constexpr double Q2 = 30265.40189360949691;
+    constexpr double E = 0x1.5bf0a8b145769p+1;
+    const short xKind = ClassifyValue(x);
+    const short yKind = ClassifyValue(y);
+    if (xKind == NotANumber) return NotANumber;
+    if (yKind == NotANumber) return NotANumber;
+    if (xKind == Infinite) {
+        if (x > 0.0) {
+            if (yKind == Zero) {
+                *value = std::numeric_limits<double>::quiet_NaN();
+                RaiseLibcError(FE_INVALID);
+                return NotANumber;
+            }
+            if (yKind == Infinite) {
+                *value = y;
+            } else {
+                *value = std::copysign(std::numeric_limits<double>::infinity(), y);
+            }
+            return Infinite;
+        }
+        if (yKind == Infinite) {
+            *value = std::numeric_limits<double>::quiet_NaN();
+            RaiseLibcError(FE_INVALID);
+            return NotANumber;
+        }
+        *value = std::copysign(0.0, y);
+        return Zero;
+    }
+    if (yKind == Zero) {
+        *value = y;
+        return Zero;
+    }
+    if (yKind == Infinite) {
+        *value = y;
+        return Infinite;
+    }
+    if (xKind == Zero) {
+        if (offset == 0) {
+            *value = y;
+            return Finite;
+        }
+        const short result = ScaleProduct(*value, 1.0, y, offset);
+        if (result == Finite && std::fpclassify(*value) == FP_SUBNORMAL) {
+            std::fesetexceptflag(&entryFlags, FE_ALL_EXCEPT);
+            errno = entryErrno;
+        }
+        return result;
+    }
+    if (offset == 0 && x == 1.0 && y == 1.0) {
+        *value = E;
+        return Finite;
+    }
+    if (x < -HugeExp) {
+        *value = std::copysign(0.0, y);
+        return Zero;
+    }
+    if (x > HugeExp) {
+        *value = std::copysign(std::numeric_limits<double>::infinity(), y);
+        RaiseLibcError(FE_OVERFLOW);
+        return Infinite;
+    }
+    double g = x * InverseLn2;
+    long exponent = static_cast<long>(g + (g < 0.0 ? -0.5 : 0.5));
+    g = static_cast<double>(exponent);
+    g = (x - g * Ln2High) - g * Ln2Low;
+    double factor;
+    if (-Eps < g && g < Eps) {
+        factor = 1.0;
+    } else {
+        const double z = g * g;
+        const double w = (Q0 * z + Q1) * z + Q2;
+        g *= (z + P1) * z + P2;
+        factor = (w + g) / (w - g) * 2.0;
+        --exponent;
+    }
+    const short result = ScaleProduct(*value, factor, y, exponent + offset);
+    if (result == Finite && std::fpclassify(*value) == FP_SUBNORMAL) {
+        std::fesetexceptflag(&entryFlags, FE_ALL_EXCEPT);
+        errno = entryErrno;
+    }
+    return result;
+}
+
+short APS5_VABI _FExp_nid_postfix(float* value, float y, short offset) {
+    if (value == nullptr) throw std::invalid_argument("_FExp: null value");
+    std::fexcept_t entryFlags;
+    std::fegetexceptflag(&entryFlags, FE_ALL_EXCEPT);
+    const int entryErrno = errno;
+    const float x = *value;
+    constexpr float HugeExp = 229.0f;
+    constexpr float Eps = std::numeric_limits<float>::epsilon();
+    constexpr float InverseLn2 = 1.4426950216f;
+    constexpr float Ln2High = 0.69287109375f;
+    constexpr float Ln2Low = 0.00027608682f;
+    constexpr float Q = 0.09997485578f;
+    constexpr float P0 = 0.5f;
+    constexpr float P1 = 0.0083208261f;
+    constexpr float E = 0x1.5bf0a8p+1f;
+    const short xKind = ClassifyValue(x);
+    const short yKind = ClassifyValue(y);
+    if (xKind == NotANumber) return NotANumber;
+    if (yKind == NotANumber) return NotANumber;
+    if (xKind == Infinite) {
+        if (x > 0.0f) {
+            if (yKind == Zero) {
+                *value = std::numeric_limits<float>::quiet_NaN();
+                RaiseLibcError(FE_INVALID);
+                return NotANumber;
+            }
+            if (yKind == Infinite) {
+                *value = y;
+            } else {
+                *value = std::copysign(std::numeric_limits<float>::infinity(), y);
+            }
+            return Infinite;
+        }
+        if (yKind == Infinite) {
+            *value = std::numeric_limits<float>::quiet_NaN();
+            RaiseLibcError(FE_INVALID);
+            return NotANumber;
+        }
+        *value = std::copysign(0.0f, y);
+        return Zero;
+    }
+    if (yKind == Zero) {
+        *value = y;
+        return Zero;
+    }
+    if (yKind == Infinite) {
+        *value = y;
+        return Infinite;
+    }
+    if (xKind == Zero) {
+        if (offset == 0) {
+            *value = y;
+            return Finite;
+        }
+        const short result = ScaleProduct(*value, 1.0f, y, offset);
+        if (result == Finite && std::fpclassify(*value) == FP_SUBNORMAL) {
+            std::fesetexceptflag(&entryFlags, FE_ALL_EXCEPT);
+            errno = entryErrno;
+        }
+        return result;
+    }
+    if (offset == 0 && x == 1.0f && y == 1.0f) {
+        *value = E;
+        return Finite;
+    }
+    if (x < -HugeExp) {
+        *value = std::copysign(0.0f, y);
+        return Zero;
+    }
+    if (x > HugeExp) {
+        *value = std::copysign(std::numeric_limits<float>::infinity(), y);
+        RaiseLibcError(FE_OVERFLOW);
+        return Infinite;
+    }
+    float g = x * InverseLn2;
+    long exponent = static_cast<long>(g + (g < 0.0f ? -0.5f : 0.5f));
+    g = static_cast<float>(exponent);
+    g = (x - g * Ln2High) - g * Ln2Low;
+    float factor;
+    if (-Eps < g && g < Eps) {
+        factor = 1.0f;
+    } else {
+        const float z = g * g;
+        const float w = Q * z + 1.0f;
+        g *= P0 + P1 * z;
+        factor = (w + g) / (w - g) * 2.0f;
+        --exponent;
+    }
+    const short result = ScaleProduct(*value, factor, y, exponent + offset);
+    if (result == Finite && std::fpclassify(*value) == FP_SUBNORMAL) {
+        std::fesetexceptflag(&entryFlags, FE_ALL_EXCEPT);
+        errno = entryErrno;
+    }
+    return result;
+}
+
 double APS5_VABI _Sin_nid_postfix(double x, unsigned int quadrant) {
     return quadrant != 0 ? cos_nid_postfix(x) : sin_nid_postfix(x);
 }

core/libs/tests/GuestDinkumwareInternals.cpp

@@ -1,16 +1,25 @@
 #include "prx/libc/include/Api.hpp"
 #include "prx/libc/include/general/VabiMacros.hpp"
+#include <cerrno>
+#include <cfenv>
 #include <cmath>
 #include <cstdint>
 #include <cstdio>
 #include <cstdlib>
 #include <cstring>
 #include <initializer_list>
+#include <limits>
 
 namespace {
 
 constexpr int RangeError = 34;
 constexpr int InvalidArgument = 22;
+constexpr short Denormal = -2;
+constexpr short Finite = -1;
+constexpr short Zero = 0;
+constexpr short Infinite = 1;
+constexpr short NotANumber = 2;
+constexpr int FloatingExceptions = FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW | FE_INEXACT;
 
 void Require(bool condition, const char* message) {
     if (!condition) {
@@ -27,6 +36,24 @@ void CheckStod(const char* text, std::int64_t pten, double expected, long consum
     Require(same && end - text == consumed && *__error_nid_postfix() == error, message);
 }
 
+void CheckExp(double x, double y, short offset, double expected, short expectedKind, int expectedErrno, int expectedFlags, const char* message) {
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    double result = x;
+    const short kind = _Exp_nid_postfix(&result, y, offset);
+    const bool same = std::isnan(expected) ? std::isnan(result) : result == expected && std::signbit(result) == std::signbit(expected);
+    Require(kind == expectedKind && same && errno == expectedErrno && (std::fetestexcept(FE_ALL_EXCEPT) & FloatingExceptions) == expectedFlags, message);
+}
+
+void CheckFExp(float x, float y, short offset, float expected, short expectedKind, int expectedErrno, int expectedFlags, const char* message) {
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    float result = x;
+    const short kind = _FExp_nid_postfix(&result, y, offset);
+    const bool same = std::isnan(expected) ? std::isnan(result) : result == expected && std::signbit(result) == std::signbit(expected);
+    Require(kind == expectedKind && same && errno == expectedErrno && (std::fetestexcept(FE_ALL_EXCEPT) & FloatingExceptions) == expectedFlags, message);
+}
+
 void TestStod() {
     CheckStod("1.5", 0, 1.5, 3, 0, "unscaled decimal");
     CheckStod("1.5", 1, 15.0, 3, 0, "decimal scaled up");
@@ -182,6 +209,182 @@ void TestIntegerErrors() {
     Require(_Stoullx_nid_postfix("18446744073709551616", &end, 10, &error) == 18446744073709551615ull && error == RangeError && *__error_nid_postfix() == RangeError, "_Stoullx overflow");
 }
 
+void TestDoubleHelpers() {
+    const double e = 0x1.5bf0a8b145769p+1;
+    CheckExp(1.0, 1.0, 0, e, Finite, 0, 0, "_Exp e");
+    CheckExp(0.0, 1.0, 0, 1.0, Finite, 0, 0, "_Exp zero exponent");
+    CheckExp(0.0, 5.0, 0, 5.0, Finite, 0, 0, "_Exp multiplier at zero exponent");
+    CheckExp(1.0, -0.0, 0, -0.0, Zero, 0, 0, "_Exp keeps a negative zero multiplier");
+    CheckExp(1.0, 0.0, 0, 0.0, Zero, 0, 0, "_Exp zero multiplier");
+    CheckExp(-1.0, -0.0, 0, -0.0, Zero, 0, 0, "_Exp negative zero multiplier");
+    CheckExp(INFINITY, -0.0, 0, std::numeric_limits<double>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_Exp positive infinity with zero multiplier is invalid");
+    CheckExp(INFINITY, 0.0, 0, std::numeric_limits<double>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_Exp positive infinity with positive zero is invalid");
+    CheckExp(INFINITY, 1.0, 0, INFINITY, Infinite, 0, 0, "_Exp positive infinity");
+    CheckExp(INFINITY, -1.0, 0, -INFINITY, Infinite, 0, 0, "_Exp positive infinity preserves the multiplier sign");
+    CheckExp(INFINITY, INFINITY, 0, INFINITY, Infinite, 0, 0, "_Exp positive infinity with infinite multiplier");
+    CheckExp(INFINITY, -INFINITY, 0, -INFINITY, Infinite, 0, 0, "_Exp positive infinity with negative infinite multiplier");
+    CheckExp(-INFINITY, 1.0, 0, 0.0, Zero, 0, 0, "_Exp negative infinity");
+    CheckExp(-INFINITY, -1.0, 0, -0.0, Zero, 0, 0, "_Exp negative infinity preserves the multiplier sign");
+    CheckExp(-INFINITY, 0.0, 0, 0.0, Zero, 0, 0, "_Exp negative infinity with positive zero");
+    CheckExp(-INFINITY, -0.0, 0, -0.0, Zero, 0, 0, "_Exp negative infinity with negative zero");
+    CheckExp(-INFINITY, INFINITY, 0, std::numeric_limits<double>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_Exp negative infinity with infinite multiplier is invalid");
+    CheckExp(-INFINITY, -INFINITY, 0, std::numeric_limits<double>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_Exp negative infinity with negative infinite multiplier is invalid");
+    CheckExp(0.0, INFINITY, 0, INFINITY, Infinite, 0, 0, "_Exp zero exponent with infinite multiplier");
+    CheckExp(0.0, -INFINITY, 0, -INFINITY, Infinite, 0, 0, "_Exp zero exponent with negative infinite multiplier");
+    CheckExp(1.0, INFINITY, 0, INFINITY, Infinite, 0, 0, "_Exp finite exponent with infinite multiplier");
+    CheckExp(1.0, -INFINITY, 0, -INFINITY, Infinite, 0, 0, "_Exp finite exponent with negative infinite multiplier");
+    CheckExp(1.0, std::numeric_limits<double>::quiet_NaN(), 0, 1.0, NotANumber, 0, 0, "_Exp NaN multiplier leaves the exponent");
+    CheckExp(INFINITY, std::numeric_limits<double>::quiet_NaN(), 0, INFINITY, NotANumber, 0, 0, "_Exp NaN multiplier leaves an infinite exponent");
+    CheckExp(std::numeric_limits<double>::quiet_NaN(), 0.0, 0, std::numeric_limits<double>::quiet_NaN(), NotANumber, 0, 0, "_Exp NaN exponent with zero multiplier");
+    CheckExp(std::numeric_limits<double>::quiet_NaN(), 1.0, 0, std::numeric_limits<double>::quiet_NaN(), NotANumber, 0, 0, "_Exp NaN exponent");
+    CheckExp(std::numeric_limits<double>::quiet_NaN(), INFINITY, 0, std::numeric_limits<double>::quiet_NaN(), NotANumber, 0, 0, "_Exp NaN exponent with infinite multiplier");
+    CheckExp(1843.0, 1.0, 0, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_Exp positive cutoff");
+    CheckExp(1843.0, -1.0, 0, -INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_Exp positive cutoff preserves the multiplier sign");
+    CheckExp(1843.0, 1.0, -2660, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_Exp positive cutoff ignores a compensating offset");
+    CheckExp(-1843.0, 1.0, 0, 0.0, Zero, 0, 0, "_Exp negative cutoff");
+    CheckExp(-1843.0, -1.0, 0, -0.0, Zero, 0, 0, "_Exp negative cutoff preserves the multiplier sign");
+    CheckExp(-1843.0, 1.0, 2660, 0.0, Zero, 0, 0, "_Exp negative cutoff ignores a compensating offset");
+    CheckExp(1843.0, 0.0, 0, 0.0, Zero, 0, 0, "_Exp positive cutoff with zero multiplier");
+    CheckExp(1843.0, std::numeric_limits<double>::quiet_NaN(), 0, 1843.0, NotANumber, 0, 0, "_Exp positive cutoff with NaN multiplier");
+    CheckExp(0.0, 1.0, -1074, std::numeric_limits<double>::denorm_min(), Finite, 0, 0, "_Exp exact minimum subnormal");
+    CheckExp(0.0, 1.0, -1075, 0.0, Zero, ERANGE, FE_UNDERFLOW | FE_INEXACT, "_Exp scaled underflow");
+    CheckExp(0.0, -1.0, -1075, -0.0, Zero, ERANGE, FE_UNDERFLOW | FE_INEXACT, "_Exp scaled underflow preserves the multiplier sign");
+    CheckExp(0.0, 1.0, 100, std::scalbn(1.0, 100), Finite, 0, 0, "_Exp scaled finite result");
+    CheckExp(0.0, 1.0, 2000, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_Exp scaled overflow");
+    CheckExp(1842.0, 1.0, 0, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_Exp inside the positive cutoff overflows");
+    CheckExp(-1842.0, 1.0, 0, 0.0, Zero, ERANGE, FE_UNDERFLOW | FE_INEXACT, "_Exp inside the negative cutoff underflows");
+
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    double result = 1.0;
+    Require(_Exp_nid_postfix(&result, 3.0, 0) == Finite && std::fabs(result - 3.0 * e) <= 3.0 * e * 1e-14, "_Exp multiplier");
+    result = 1.0;
+    Require(_Exp_nid_postfix(&result, -2.0, 0) == Finite && std::fabs(result + 2.0 * e) <= 2.0 * e * 1e-14, "_Exp negative multiplier");
+    result = 1.0;
+    Require(_Exp_nid_postfix(&result, 1.0, 10) == Finite && std::fabs(result - e * 1024.0) <= e * 1024.0 * 1e-14, "_Exp scale");
+    for (int i = -200; i <= 200; ++i) {
+        const double x = static_cast<double>(i) / 16.0;
+        result = x;
+        Require(_Exp_nid_postfix(&result, 1.0, 0) == Finite && std::fabs(result - std::exp(x)) <= std::fabs(std::exp(x)) * 1e-13, "_Exp finite sweep");
+    }
+
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    result = 1.0;
+    Require(_Exp_nid_postfix(&result, std::numeric_limits<double>::denorm_min(), 0) == Finite && std::fpclassify(result) == FP_SUBNORMAL && errno == 0 && std::fetestexcept(FE_ALL_EXCEPT) == 0, "_Exp subnormal result classification");
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    result = -740.0;
+    Require(_Exp_nid_postfix(&result, 1.0, 0) == Finite && std::fpclassify(result) == FP_SUBNORMAL && errno == 0 && std::fetestexcept(FE_ALL_EXCEPT) == 0, "_Exp negative exponent subnormal result");
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    result = 1.0;
+    Require(_Exp_nid_postfix(&result, std::numeric_limits<double>::max(), 0) == Infinite && std::isinf(result) && errno == ERANGE && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_INEXACT)) == (FE_OVERFLOW | FE_INEXACT), "_Exp finite multiplier overflow");
+
+    const double largeExpected = std::scalbn(std::exp(710.0 - 1024.0 * std::log(2.0)) * 1e-100, 1024);
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    result = 710.0;
+    Require(_Exp_nid_postfix(&result, 1e-100, 0) == Finite && std::isfinite(result) && std::fabs(result - largeExpected) <= largeExpected * 1e-13 && errno == 0 && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_UNDERFLOW)) == 0, "_Exp large exponent stays finite");
+    const double smallExpected = std::scalbn(std::exp(-710.0 + 1024.0 * std::log(2.0)) * 1e300, -1024);
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    result = -710.0;
+    Require(_Exp_nid_postfix(&result, 1e300, 0) == Finite && std::isfinite(result) && std::fabs(result - smallExpected) <= smallExpected * 1e-13 && errno == 0 && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_UNDERFLOW)) == 0, "_Exp small exponent stays finite");
+    const double multiplierExpected = std::exp(-0.5) * std::numeric_limits<double>::max();
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    result = -0.5;
+    Require(_Exp_nid_postfix(&result, std::numeric_limits<double>::max(), 0) == Finite && std::isfinite(result) && std::fabs(result - multiplierExpected) <= multiplierExpected * 1e-13 && errno == 0 && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_UNDERFLOW)) == 0, "_Exp large multiplier stays finite");
+
+    const float ef = 0x1.5bf0a8p+1f;
+    CheckFExp(1.0f, 1.0f, 0, ef, Finite, 0, 0, "_FExp e");
+    CheckFExp(0.0f, 1.0f, 0, 1.0f, Finite, 0, 0, "_FExp zero exponent");
+    CheckFExp(0.0f, 5.0f, 0, 5.0f, Finite, 0, 0, "_FExp multiplier at zero exponent");
+    CheckFExp(1.0f, -0.0f, 0, -0.0f, Zero, 0, 0, "_FExp keeps a negative zero multiplier");
+    CheckFExp(1.0f, 0.0f, 0, 0.0f, Zero, 0, 0, "_FExp zero multiplier");
+    CheckFExp(-1.0f, -0.0f, 0, -0.0f, Zero, 0, 0, "_FExp negative zero multiplier");
+    CheckFExp(INFINITY, -0.0f, 0, std::numeric_limits<float>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_FExp positive infinity with zero multiplier is invalid");
+    CheckFExp(INFINITY, 0.0f, 0, std::numeric_limits<float>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_FExp positive infinity with positive zero is invalid");
+    CheckFExp(INFINITY, 1.0f, 0, INFINITY, Infinite, 0, 0, "_FExp positive infinity");
+    CheckFExp(INFINITY, -1.0f, 0, -INFINITY, Infinite, 0, 0, "_FExp positive infinity preserves the multiplier sign");
+    CheckFExp(INFINITY, INFINITY, 0, INFINITY, Infinite, 0, 0, "_FExp positive infinity with infinite multiplier");
+    CheckFExp(INFINITY, -INFINITY, 0, -INFINITY, Infinite, 0, 0, "_FExp positive infinity with negative infinite multiplier");
+    CheckFExp(-INFINITY, 1.0f, 0, 0.0f, Zero, 0, 0, "_FExp negative infinity");
+    CheckFExp(-INFINITY, -1.0f, 0, -0.0f, Zero, 0, 0, "_FExp negative infinity preserves the multiplier sign");
+    CheckFExp(-INFINITY, 0.0f, 0, 0.0f, Zero, 0, 0, "_FExp negative infinity with positive zero");
+    CheckFExp(-INFINITY, -0.0f, 0, -0.0f, Zero, 0, 0, "_FExp negative infinity with negative zero");
+    CheckFExp(-INFINITY, INFINITY, 0, std::numeric_limits<float>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_FExp negative infinity with infinite multiplier is invalid");
+    CheckFExp(-INFINITY, -INFINITY, 0, std::numeric_limits<float>::quiet_NaN(), NotANumber, EDOM, FE_INVALID, "_FExp negative infinity with negative infinite multiplier is invalid");
+    CheckFExp(0.0f, INFINITY, 0, INFINITY, Infinite, 0, 0, "_FExp zero exponent with infinite multiplier");
+    CheckFExp(0.0f, -INFINITY, 0, -INFINITY, Infinite, 0, 0, "_FExp zero exponent with negative infinite multiplier");
+    CheckFExp(1.0f, INFINITY, 0, INFINITY, Infinite, 0, 0, "_FExp finite exponent with infinite multiplier");
+    CheckFExp(1.0f, -INFINITY, 0, -INFINITY, Infinite, 0, 0, "_FExp finite exponent with negative infinite multiplier");
+    CheckFExp(1.0f, std::numeric_limits<float>::quiet_NaN(), 0, 1.0f, NotANumber, 0, 0, "_FExp NaN multiplier leaves the exponent");
+    CheckFExp(INFINITY, std::numeric_limits<float>::quiet_NaN(), 0, INFINITY, NotANumber, 0, 0, "_FExp NaN multiplier leaves an infinite exponent");
+    CheckFExp(std::numeric_limits<float>::quiet_NaN(), 0.0f, 0, std::numeric_limits<float>::quiet_NaN(), NotANumber, 0, 0, "_FExp NaN exponent with zero multiplier");
+    CheckFExp(std::numeric_limits<float>::quiet_NaN(), 1.0f, 0, std::numeric_limits<float>::quiet_NaN(), NotANumber, 0, 0, "_FExp NaN exponent");
+    CheckFExp(std::numeric_limits<float>::quiet_NaN(), INFINITY, 0, std::numeric_limits<float>::quiet_NaN(), NotANumber, 0, 0, "_FExp NaN exponent with infinite multiplier");
+    CheckFExp(230.0f, 1.0f, 0, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_FExp positive cutoff");
+    CheckFExp(230.0f, -1.0f, 0, -INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_FExp positive cutoff preserves the multiplier sign");
+    CheckFExp(230.0f, 1.0f, -332, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_FExp positive cutoff ignores the offset");
+    CheckFExp(-230.0f, 1.0f, 0, 0.0f, Zero, 0, 0, "_FExp negative cutoff");
+    CheckFExp(-230.0f, -1.0f, 0, -0.0f, Zero, 0, 0, "_FExp negative cutoff preserves the multiplier sign");
+    CheckFExp(-230.0f, 1.0f, 332, 0.0f, Zero, 0, 0, "_FExp negative cutoff ignores the offset");
+    CheckFExp(230.0f, 0.0f, 0, 0.0f, Zero, 0, 0, "_FExp positive cutoff with zero multiplier");
+    CheckFExp(230.0f, std::numeric_limits<float>::quiet_NaN(), 0, 230.0f, NotANumber, 0, 0, "_FExp positive cutoff with NaN multiplier");
+    CheckFExp(0.0f, 1.0f, -149, std::numeric_limits<float>::denorm_min(), Finite, 0, 0, "_FExp exact minimum subnormal");
+    CheckFExp(0.0f, 1.0f, -150, 0.0f, Zero, ERANGE, FE_UNDERFLOW | FE_INEXACT, "_FExp scaled underflow");
+    CheckFExp(0.0f, -1.0f, -150, -0.0f, Zero, ERANGE, FE_UNDERFLOW | FE_INEXACT, "_FExp scaled underflow preserves the multiplier sign");
+    CheckFExp(0.0f, 1.0f, 100, std::scalbn(1.0f, 100), Finite, 0, 0, "_FExp scaled finite result");
+    CheckFExp(0.0f, 1.0f, 200, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_FExp scaled overflow");
+    CheckFExp(229.0f, 1.0f, 0, INFINITY, Infinite, ERANGE, FE_OVERFLOW | FE_INEXACT, "_FExp inside the positive cutoff overflows");
+    CheckFExp(-229.0f, 1.0f, 0, 0.0f, Zero, ERANGE, FE_UNDERFLOW | FE_INEXACT, "_FExp inside the negative cutoff underflows");
+
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    float resultf = 1.0f;
+    Require(_FExp_nid_postfix(&resultf, 3.0f, 0) == Finite && std::fabs(resultf - 3.0f * ef) <= 3.0f * ef * 1e-6f, "_FExp multiplier");
+    resultf = 1.0f;
+    Require(_FExp_nid_postfix(&resultf, -2.0f, 0) == Finite && std::fabs(resultf + 2.0f * ef) <= 2.0f * ef * 1e-6f, "_FExp negative multiplier");
+    resultf = 1.0f;
+    Require(_FExp_nid_postfix(&resultf, 1.0f, 10) == Finite && std::fabs(resultf - ef * 1024.0f) <= ef * 1024.0f * 1e-6f, "_FExp scale");
+    for (int i = -120; i <= 120; ++i) {
+        const float x = static_cast<float>(i) / 16.0f;
+        resultf = x;
+        Require(_FExp_nid_postfix(&resultf, 1.0f, 0) == Finite && std::fabs(resultf - std::exp(x)) <= std::fabs(std::exp(x)) * 1e-6f, "_FExp finite sweep");
+    }
+
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    resultf = 1.0f;
+    Require(_FExp_nid_postfix(&resultf, std::numeric_limits<float>::denorm_min(), 0) == Finite && std::fpclassify(resultf) == FP_SUBNORMAL && errno == 0 && std::fetestexcept(FE_ALL_EXCEPT) == 0, "_FExp subnormal result classification");
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    resultf = -100.0f;
+    Require(_FExp_nid_postfix(&resultf, 1.0f, 0) == Finite && std::fpclassify(resultf) == FP_SUBNORMAL && errno == 0 && std::fetestexcept(FE_ALL_EXCEPT) == 0, "_FExp negative exponent subnormal result");
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    resultf = 1.0f;
+    Require(_FExp_nid_postfix(&resultf, std::numeric_limits<float>::max(), 0) == Infinite && std::isinf(resultf) && errno == ERANGE && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_INEXACT)) == (FE_OVERFLOW | FE_INEXACT), "_FExp finite multiplier overflow");
+
+    const float largeExpectedf = static_cast<float>(std::exp(100.0) * 1e-35);
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    resultf = 100.0f;
+    Require(_FExp_nid_postfix(&resultf, 1e-35f, 0) == Finite && std::isfinite(resultf) && std::fabs(resultf - largeExpectedf) <= largeExpectedf * 1e-6f && errno == 0 && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_UNDERFLOW)) == 0, "_FExp large exponent stays finite");
+    const float smallExpectedf = static_cast<float>(std::exp(-100.0) * 1e30);
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    resultf = -100.0f;
+    Require(_FExp_nid_postfix(&resultf, 1e30f, 0) == Finite && std::isfinite(resultf) && std::fabs(resultf - smallExpectedf) <= smallExpectedf * 1e-6f && errno == 0 && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_UNDERFLOW)) == 0, "_FExp small exponent stays finite");
+    const float multiplierExpectedf = static_cast<float>(std::exp(-0.5) * std::numeric_limits<float>::max());
+    errno = 0;
+    std::feclearexcept(FE_ALL_EXCEPT);
+    resultf = -0.5f;
+    Require(_FExp_nid_postfix(&resultf, std::numeric_limits<float>::max(), 0) == Finite && std::isfinite(resultf) && std::fabs(resultf - multiplierExpectedf) <= multiplierExpectedf * 1e-6f && errno == 0 && (std::fetestexcept(FE_ALL_EXCEPT) & (FE_OVERFLOW | FE_UNDERFLOW)) == 0, "_FExp large multiplier stays finite");
+}
+
 void TestDispatch() {
     Require(_Sin_nid_postfix(0.5, 0) == std::sin(0.5) && _Sin_nid_postfix(0.5, 1) == std::cos(0.5) && _Sin_nid_postfix(0.5, 7) == std::cos(0.5), "_Sin");
     Require(_FSin_nid_postfix(0.5f, 0) == std::sin(0.5f) && _FSin_nid_postfix(0.5f, 2) == std::cos(0.5f), "_FSin");
@@ -206,6 +409,7 @@ void TestLongDoubleHelpers() {
 int main() {
     TestStod();
     TestStoll();
+    TestDoubleHelpers();
     TestFloatErrors();
     TestIntegerErrors();
     TestArithmetic();

docs/dev/debt/libc.md

@@ -173,6 +173,7 @@
 - [std::_Throw_Cpp_error](../../../core/libs/prx/libc/src/exception/StandardExceptions.cpp) (libc) follows Dinkumware `thread0.cpp` ([VC12](https://github.com/ganboing/crt.vc12/blob/f1d47126/src/thread0.cpp), [STL](https://github.com/microsoft/STL/blob/da52dc0f/stl/src/thread0.cpp)), without console verification. Codes outside 0..6 throw `std::out_of_range`, the deleting destructor throws, and STREAMS errnos absent from FreeBSD 11.4 [`errno.h`](https://github.com/freebsd/freebsd-src/blob/releng/11.4/sys/sys/errno.h) (`ENODATA`, `ENOSR`, `ENOSTR`, `ETIME`) are omitted. Codes 0..6 throw a 32-byte [VC12 `system_error`](https://github.com/ganboing/crt.vc12/blob/f1d47126/src/system_error) with FreeBSD errnos (16, 22, 3, 12, 1, 11, 35), a static `generic_category` (not `_Generic_object`), and unverified 40-byte [GuestLocale::String](../../../core/libs/prx/libc/include/GuestLocale.hpp) messages from [`syserror.cpp`](https://github.com/ganboing/crt.vc12/blob/f1d47126/src/syserror.cpp)
 - [Call tracing](../../../core/libs/prx/libc/src/specifics/x86_64/CallTrace.cpp) (`--trace-imports`, `APS5_TRACE_CALLS`) is missing on macOS and records only the entry of a call, not its arguments or result. One process-wide ring keeps the last 65536 calls, so a busy thread pushes the others out; the reporting thread's own calls are listed only while they are still in that ring. Calls are counted per traced target; an import bound to a target already traced under another name is counted under the first name. The thunk saves the integer argument registers and `xmm0` to `xmm7`; wider vector arguments survive only while the recording code avoids AVX, which is untested. Linux has no crash report, so it prints the trace only on `SIGABRT` (guest or host `abort()`), not on `SIGSEGV`; a caller in a module loaded after the last traced import prints without a module name, and a thread aborting while another writes the report waits for it at most 2 seconds
 - [_Dtest](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows Dinkumware double classification matching _FDtest: finite values return -1, subnormals -2, zeros 0, infinities 1, and NaNs 2.
+- [_Exp and _FExp](../../../core/libs/prx/libc/src/DinkumwareInternals.cpp) (libc) follow PRAGMATA PPSA02530's bundled `libc.prx` (static disassembly): both operands are classified first and only finite/denormal pairs take the finite path. The strict x-only cutoffs at ±229/±1842 store `0.0*y` without an error below and `y*inf` with overflow above, so an offset does not compensate them. A zero multiplier keeps its sign, a NaN multiplier leaves `*px` unchanged, `+inf` with zero or `-inf` with infinity is invalid, and finite multipliers with `+inf`/`-inf` return signed infinity/zero without flags. Inside the cutoffs subnormals are `Finite` without flags, a zero from scaling raises underflow, and the negative cutoff itself raises nothing. `_FExp` uses the measured `c1`, `c2`, `invln2`, `Q` and `P`; `_Exp` finite constants, exact finite bits, signaling NaNs, directed rounding and flush-to-zero modes remain unverified.
 - [exp2f](../../../core/libs/prx/libc/src/Math.cpp) (libc) is FreeBSD msun's [`s_exp2f.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_exp2f.c), with `huge` and `twom100` read at run time and `x + redux` rounded to float, as msun's `volatile` reads and `STRICT_ASSIGN` do. As in the `libc.prx` that titles ship, `errno` is set to `ERANGE` only in msun's overflow (`x >= 128`) and underflow (`x <= -150`) branches, so a subnormal result leaves it unchanged. Floating-point exceptions are not compared. It gives the same bits and `errno` as the `exp2f` of the `libc.prx` bundled with PPSA07923, PPSA05916 and PPSA13934 for all 2^32 arguments, run on a Windows host.
 - [app0.zar](../../../core/libs/prx/libc/src/AppArchive.cpp) (libc) - a [ZArchive](https://github.com/Exzap/ZArchive) beside `app0/` is a read-only layer under it: a path missing from `app0/` is looked up in the archive, case-insensitively. It is not a console format; titles ship their files in a PFS image instead. Resolving an archive-only path through any call other than the archive-aware libkernel ones throws, so `fopen`, `stat`, `access`, `opendir` and the rest do not read the archive yet, and `opendir` of a directory that exists in both places, `app0/` itself included, and `fdopen` of an archived descriptor throw instead of missing its entries or reading the null device. A damaged compressed block throws instead of returning `EIO`
 - [fopen_s](../../../core/libs/prx/libc/src/FileAndHeap.cpp) (libc) - returns the errno of the failed open, as C11 K.3.5.2.1 leaves the nonzero value open and Microsoft's CRT returns errno; the console's Dinkumware `fopen_s` is not verified