aboutsummaryrefslogtreecommitdiff
path: root/libc/src/math/generic
diff options
context:
space:
mode:
Diffstat (limited to 'libc/src/math/generic')
-rw-r--r--libc/src/math/generic/CMakeLists.txt37
-rw-r--r--libc/src/math/generic/asinf16.cpp121
-rw-r--r--libc/src/math/generic/asinhf.cpp106
-rw-r--r--libc/src/math/generic/asinpif16.cpp127
4 files changed, 6 insertions, 385 deletions
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 8116ee2..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(
@@ -3919,25 +3914,6 @@ add_entrypoint_object(
)
add_entrypoint_object(
- asinpif16
- SRCS
- asinpif16.cpp
- HDRS
- ../asinpif16.h
- DEPENDS
- libc.hdr.errno_macros
- libc.hdr.fenv_macros
- libc.src.__support.FPUtil.cast
- libc.src.__support.FPUtil.except_value_utils
- 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
-)
-
-add_entrypoint_object(
atanhf
SRCS
atanhf.cpp
@@ -3987,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
diff --git a/libc/src/math/generic/asinpif16.cpp b/libc/src/math/generic/asinpif16.cpp
deleted file mode 100644
index aabc086..0000000
--- a/libc/src/math/generic/asinpif16.cpp
+++ /dev/null
@@ -1,127 +0,0 @@
-//===-- Half-precision asinpif16(x) function ------------------------------===//
-//
-// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
-// See https://llvm.org/LICENSE.txt for license information.
-// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
-//
-//===----------------------------------------------------------------------===//
-
-#include "src/math/asinpif16.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/except_value_utils.h"
-#include "src/__support/FPUtil/multiply_add.h"
-#include "src/__support/FPUtil/sqrt.h"
-#include "src/__support/macros/optimization.h"
-
-namespace LIBC_NAMESPACE_DECL {
-
-LLVM_LIBC_FUNCTION(float16, asinpif16, (float16 x)) {
- using FPBits = fputil::FPBits<float16>;
-
- FPBits xbits(x);
- bool is_neg = xbits.is_neg();
- double x_abs = fputil::cast<double>(xbits.abs().get_val());
-
- auto signed_result = [is_neg](auto r) -> auto { return is_neg ? -r : r; };
-
- if (LIBC_UNLIKELY(x_abs > 1.0)) {
- // aspinf16(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();
- }
-
- // the coefficients for the polynomial approximation of asin(x)/pi in the
- // range [0, 0.5] extracted using python-sympy
- //
- // Python code to generate the coefficients:
- // > from sympy import *
- // > import math
- // > x = symbols('x')
- // > print(series(asin(x)/math.pi, x, 0, 21))
- //
- // OUTPUT:
- //
- // 0.318309886183791*x + 0.0530516476972984*x**3 + 0.0238732414637843*x**5 +
- // 0.0142102627760621*x**7 + 0.00967087327815336*x**9 +
- // 0.00712127941391293*x**11 + 0.00552355646848375*x**13 +
- // 0.00444514782463692*x**15 + 0.00367705242846804*x**17 +
- // 0.00310721681820837*x**19 + O(x**21)
- //
- // it's very accurate in the range [0, 0.5] and has a maximum error of
- // 0.0000000000000001 in the range [0, 0.5].
- constexpr double POLY_COEFFS[] = {
- 0x1.45f306dc9c889p-2, // x^1
- 0x1.b2995e7b7b5fdp-5, // x^3
- 0x1.8723a1d588a36p-6, // x^5
- 0x1.d1a452f20430dp-7, // x^7
- 0x1.3ce52a3a09f61p-7, // x^9
- 0x1.d2b33e303d375p-8, // x^11
- 0x1.69fde663c674fp-8, // x^13
- 0x1.235134885f19bp-8, // x^15
- };
- // polynomial evaluation using horner's method
- // work only for |x| in [0, 0.5]
- auto asinpi_polyeval = [](double x) -> double {
- return x * fputil::polyeval(x * x, POLY_COEFFS[0], POLY_COEFFS[1],
- POLY_COEFFS[2], POLY_COEFFS[3], POLY_COEFFS[4],
- POLY_COEFFS[5], POLY_COEFFS[6], POLY_COEFFS[7]);
- };
-
- // if |x| <= 0.5:
- if (LIBC_UNLIKELY(x_abs <= 0.5)) {
- // Use polynomial approximation of asin(x)/pi in the range [0, 0.5]
- double result = asinpi_polyeval(fputil::cast<double>(x));
- return fputil::cast<float16>(result);
- }
-
- // If |x| > 0.5, we need to use the range reduction method:
- // y = asin(x) => x = sin(y)
- // because: sin(a) = cos(pi/2 - a)
- // therefore:
- // x = cos(pi/2 - y)
- // let z = pi/2 - y,
- // x = cos(z)
- // because: cos(2a) = 1 - 2 * sin^2(a), z = 2a, a = z/2
- // therefore:
- // cos(z) = 1 - 2 * sin^2(z/2)
- // sin(z/2) = sqrt((1 - cos(z))/2)
- // sin(z/2) = sqrt((1 - x)/2)
- // let u = (1 - x)/2
- // then:
- // sin(z/2) = sqrt(u)
- // z/2 = asin(sqrt(u))
- // z = 2 * asin(sqrt(u))
- // pi/2 - y = 2 * asin(sqrt(u))
- // y = pi/2 - 2 * asin(sqrt(u))
- // y/pi = 1/2 - 2 * asin(sqrt(u))/pi
- //
- // Finally, we can write:
- // asinpi(x) = 1/2 - 2 * asinpi(sqrt(u))
- // where u = (1 - x) /2
- // = 0.5 - 0.5 * x
- // = multiply_add(-0.5, x, 0.5)
-
- double u = fputil::multiply_add(-0.5, x_abs, 0.5);
- double asinpi_sqrt_u = asinpi_polyeval(fputil::sqrt<double>(u));
- double result = fputil::multiply_add(-2.0, asinpi_sqrt_u, 0.5);
-
- return fputil::cast<float16>(signed_result(result));
-}
-
-} // namespace LIBC_NAMESPACE_DECL