diff options
Diffstat (limited to 'libc/src/math/generic')
| -rw-r--r-- | libc/src/math/generic/CMakeLists.txt | 18 | ||||
| -rw-r--r-- | libc/src/math/generic/asinf16.cpp | 121 | ||||
| -rw-r--r-- | libc/src/math/generic/asinhf.cpp | 106 |
3 files changed, 6 insertions, 239 deletions
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index ecf0967..f91feacb 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3889,12 +3889,7 @@ add_entrypoint_object( HDRS ../asinhf.h DEPENDS - .explogxf - libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.multiply_add - libc.src.__support.FPUtil.polyeval - libc.src.__support.FPUtil.sqrt - libc.src.__support.macros.optimization + libc.src.__support.math.asinhf ) add_entrypoint_object( @@ -3968,16 +3963,7 @@ add_entrypoint_object( HDRS ../asinf16.h DEPENDS - libc.hdr.errno_macros - libc.hdr.fenv_macros - libc.src.__support.FPUtil.cast - libc.src.__support.FPUtil.fenv_impl - libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.multiply_add - libc.src.__support.FPUtil.polyeval - libc.src.__support.FPUtil.sqrt - libc.src.__support.macros.optimization - libc.src.__support.macros.properties.types + libc.src.__support.math.asinf16 ) add_entrypoint_object( diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp index 518c384..af8dbfe 100644 --- a/libc/src/math/generic/asinf16.cpp +++ b/libc/src/math/generic/asinf16.cpp @@ -7,127 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/asinf16.h" -#include "hdr/errno_macros.h" -#include "hdr/fenv_macros.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/FPUtil/cast.h" -#include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/sqrt.h" -#include "src/__support/macros/optimization.h" +#include "src/__support/math/asinf16.h" namespace LIBC_NAMESPACE_DECL { -// Generated by Sollya using the following command: -// > round(pi/2, D, RN); -static constexpr float PI_2 = 0x1.921fb54442d18p0f; - -LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) { - using FPBits = fputil::FPBits<float16>; - FPBits xbits(x); - - uint16_t x_u = xbits.uintval(); - uint16_t x_abs = x_u & 0x7fff; - float xf = x; - - // |x| > 0x1p0, |x| > 1, or x is NaN. - if (LIBC_UNLIKELY(x_abs > 0x3c00)) { - // asinf16(NaN) = NaN - if (xbits.is_nan()) { - if (xbits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - - return x; - } - - // 1 < |x| <= +/-inf - fputil::raise_except_if_required(FE_INVALID); - fputil::set_errno_if_required(EDOM); - - return FPBits::quiet_nan().get_val(); - } - - float xsq = xf * xf; - - // |x| <= 0x1p-1, |x| <= 0.5 - if (x_abs <= 0x3800) { - // asinf16(+/-0) = +/-0 - if (LIBC_UNLIKELY(x_abs == 0)) - return x; - - // Exhaustive tests show that, - // for |x| <= 0x1.878p-9, when: - // x > 0, and rounding upward, or - // x < 0, and rounding downward, then, - // asin(x) = x * 2^-11 + x - // else, in other rounding modes, - // asin(x) = x - if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) { - int rounding = fputil::quick_get_round(); - - if ((xbits.is_pos() && rounding == FE_UPWARD) || - (xbits.is_neg() && rounding == FE_DOWNWARD)) - return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf)); - return x; - } - - // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: - // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); - float result = - fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f, - 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); - return fputil::cast<float16>(xf * result); - } - - // When |x| > 0.5, assume that 0.5 < |x| <= 1, - // - // Step-by-step range-reduction proof: - // 1: Let y = asin(x), such that, x = sin(y) - // 2: From complimentary angle identity: - // x = sin(y) = cos(pi/2 - y) - // 3: Let z = pi/2 - y, such that x = cos(z) - // 4: From double angle formula; cos(2A) = 1 - sin^2(A): - // z = 2A, z/2 = A - // cos(z) = 1 - 2 * sin^2(z/2) - // 5: Make sin(z/2) subject of the formula: - // sin(z/2) = sqrt((1 - cos(z))/2) - // 6: Recall [3]; x = cos(z). Therefore: - // sin(z/2) = sqrt((1 - x)/2) - // 7: Let u = (1 - x)/2 - // 8: Therefore: - // asin(sqrt(u)) = z/2 - // 2 * asin(sqrt(u)) = z - // 9: Recall [3], z = pi/2 - y. Therefore: - // y = pi/2 - z - // y = pi/2 - 2 * asin(sqrt(u)) - // 10: Recall [1], y = asin(x). Therefore: - // asin(x) = pi/2 - 2 * asin(sqrt(u)) - // - // WHY? - // 11: Recall [7], u = (1 - x)/2 - // 12: Since 0.5 < x <= 1, therefore: - // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5 - // - // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for - // Step [10] as `sqrt(u)` is in range. - - // 0x1p-1 < |x| <= 0x1p0, 0.5 < |x| <= 1.0 - float xf_abs = (xf < 0 ? -xf : xf); - float sign = (xbits.uintval() >> 15 == 1 ? -1.0 : 1.0); - float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f); - float u_sqrt = fputil::sqrt<float>(u); - - // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: - // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); - float asin_sqrt_u = - u_sqrt * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f, - 0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); - - return fputil::cast<float16>(sign * - fputil::multiply_add(-2.0f, asin_sqrt_u, PI_2)); -} +LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) { return math::asinf16(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/asinhf.cpp b/libc/src/math/generic/asinhf.cpp index 3aed3bc..45023c8 100644 --- a/libc/src/math/generic/asinhf.cpp +++ b/libc/src/math/generic/asinhf.cpp @@ -7,112 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/asinhf.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/sqrt.h" -#include "src/__support/macros/config.h" -#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY -#include "src/math/generic/common_constants.h" -#include "src/math/generic/explogxf.h" +#include "src/__support/math/asinhf.h" namespace LIBC_NAMESPACE_DECL { -LLVM_LIBC_FUNCTION(float, asinhf, (float x)) { - using namespace acoshf_internal; - using FPBits_t = typename fputil::FPBits<float>; - FPBits_t xbits(x); - uint32_t x_u = xbits.uintval(); - uint32_t x_abs = xbits.abs().uintval(); - - // |x| <= 2^-3 - if (LIBC_UNLIKELY(x_abs <= 0x3e80'0000U)) { - // |x| <= 2^-26 - if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) { - return static_cast<float>(LIBC_UNLIKELY(x_abs == 0) - ? x - : (x - 0x1.5555555555555p-3 * x * x * x)); - } - - double x_d = x; - double x_sq = x_d * x_d; - // Generated by Sollya with: - // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16|], [|D...|], - // [0, 2^-2]); - double p = fputil::polyeval( - x_sq, 0.0, -0x1.555555555551ep-3, 0x1.3333333325495p-4, - -0x1.6db6db5a7622bp-5, 0x1.f1c70f82928c6p-6, -0x1.6e893934266b7p-6, - 0x1.1c0b41d3fbe78p-6, -0x1.c0f47810b3c4fp-7, 0x1.2c8602690143dp-7); - return static_cast<float>(fputil::multiply_add(x_d, p, x_d)); - } - - const double SIGN[2] = {1.0, -1.0}; - double x_sign = SIGN[x_u >> 31]; - double x_d = x; - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - // Helper functions to set results for exceptional cases. - auto round_result_slightly_down = [x_sign](float r) -> float { - return fputil::multiply_add(static_cast<float>(x_sign), r, - static_cast<float>(x_sign) * (-0x1.0p-24f)); - }; - auto round_result_slightly_up = [x_sign](float r) -> float { - return fputil::multiply_add(static_cast<float>(x_sign), r, - static_cast<float>(x_sign) * 0x1.0p-24f); - }; - - if (LIBC_UNLIKELY(x_abs >= 0x4bdd'65a5U)) { - if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) { - if (xbits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits_t::quiet_nan().get_val(); - } - - return x; - } - - // Exceptional cases when x > 2^24. - switch (x_abs) { - case 0x4bdd65a5: // |x| = 0x1.bacb4ap24f - return round_result_slightly_down(0x1.1e0696p4f); - case 0x4c803f2c: // |x| = 0x1.007e58p26f - return round_result_slightly_down(0x1.2b786cp4f); - case 0x4f8ffb03: // |x| = 0x1.1ff606p32f - return round_result_slightly_up(0x1.6fdd34p4f); - case 0x5c569e88: // |x| = 0x1.ad3d1p57f - return round_result_slightly_up(0x1.45c146p5f); - case 0x5e68984e: // |x| = 0x1.d1309cp61f - return round_result_slightly_up(0x1.5c9442p5f); - case 0x655890d3: // |x| = 0x1.b121a6p75f - return round_result_slightly_down(0x1.a9a3f2p5f); - case 0x65de7ca6: // |x| = 0x1.bcf94cp76f - return round_result_slightly_up(0x1.af66cp5f); - case 0x6eb1a8ec: // |x| = 0x1.6351d8p94f - return round_result_slightly_down(0x1.08b512p6f); - case 0x7997f30a: // |x| = 0x1.2fe614p116f - return round_result_slightly_up(0x1.451436p6f); - } - } else { - // Exceptional cases when x < 2^24. - if (LIBC_UNLIKELY(x_abs == 0x45abaf26)) { - // |x| = 0x1.575e4cp12f - return round_result_slightly_down(0x1.29becap3f); - } - if (LIBC_UNLIKELY(x_abs == 0x49d29048)) { - // |x| = 0x1.a5209p20f - return round_result_slightly_down(0x1.e1b92p3f); - } - } -#else - if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) - return x; -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - // asinh(x) = log(x + sqrt(x^2 + 1)) - return static_cast<float>( - x_sign * log_eval(fputil::multiply_add( - x_d, x_sign, - fputil::sqrt<double>(fputil::multiply_add(x_d, x_d, 1.0))))); -} +LLVM_LIBC_FUNCTION(float, asinhf, (float x)) { return math::asinhf(x); } } // namespace LIBC_NAMESPACE_DECL |