shylie/llvm-project
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

[libc][math] Refactor asinf implementation to header-only in src/__support/math folder. (#150697)

Browse Source 
Part of #147386

in preparation for:
https://discourse.llvm.org/t/rfc-make-clang-builtin-math-functions-constexpr-with-llvm-libc-to-support-c-23-constexpr-math-functions/86450
This commit is contained in:
Muhammad Bassiouni 2025-07-26 02:05:09 +03:00 committed by GitHub
parent cf6a4bbc42
commit 2eb112d78d
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
9 changed files with 237 additions and 168 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
libc/shared/math.h
View File

@ -18,6 +18,7 @@
#include "math/acoshf16.h"
#include "math/acospif16.h"
#include "math/asin.h"
#include "math/asinf.h"
#include "math/erff.h"
#include "math/exp.h"
#include "math/exp10.h"

23
libc/shared/math/asinf.h Normal file
View File

@ -0,0 +1,23 @@
//===-- Shared asinf function -----------------------------------*- C++ -*-===//
//
// 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
//
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SHARED_MATH_ASINF_H
#define LLVM_LIBC_SHARED_MATH_ASINF_H
#include "shared/libc_common.h"
#include "src/__support/math/asinf.h"
namespace LIBC_NAMESPACE_DECL {
namespace shared {
using math::asinf;
} // namespace shared
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_SHARED_MATH_ASINF_H

16
libc/src/__support/math/CMakeLists.txt
View File

@ -140,6 +140,22 @@ add_header_library(
libc.src.__support.macros.properties.cpu_features
)
add_header_library(
asinf
HDRS
asinf.h
DEPENDS
.inv_trigf_utils
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.config
libc.src.__support.macros.optimization
libc.src.__support.macros.properties.cpu_features
)
add_header_library(
erff
HDRS

175
libc/src/__support/math/asinf.h Normal file
View File

@ -0,0 +1,175 @@
//===-- Implementation header for asinf -------------------------*- C++ -*-===//
//
// 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
//
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H
#include "inv_trigf_utils.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.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/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
namespace LIBC_NAMESPACE_DECL {
namespace math {
LIBC_INLINE static constexpr float asinf(float x) {
using namespace inv_trigf_utils_internal;
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
constexpr size_t N_EXCEPTS = 2;
// Exceptional values when |x| <= 0.5
constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{
// (inputs, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)
{0x3d09bf86, 0x3d09c62c, 1, 0, 1},
// x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)
{0x3de5fa1e, 0x3de6768e, 1, 0, 0},
}};
// Exceptional values when 0.5 < |x| <= 1
constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{
// (inputs, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)
{0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},
// x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)
{0x3f7741b6, 0x3fa7832a, 1, 0, 0},
}};
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
using namespace inv_trigf_utils_internal;
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
uint32_t x_uint = xbits.uintval();
uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
constexpr double SIGN[2] = {1.0, -1.0};
uint32_t x_sign = x_uint >> 31;
// |x| <= 0.5-ish
if (x_abs < 0x3f04'471dU) {
// |x| < 0x1.d12edp-12
if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) {
// When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
// is:
// |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
// = x^2 / 6
// < 2^-25
// < epsilon(1)/2.
// So the correctly rounded values of asin(x) are:
// = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
// or (rounding mode = FE_UPWARD and x is
// negative),
// = x otherwise.
// To simplify the rounding decision and make it more efficient, we use
// fma(x, 2^-25, x) instead.
// An exhaustive test shows that this formula work correctly for all
// rounding modes up to |x| < 0x1.d12edp-12.
// Note: to use the formula x + 2^-25*x to decide the correct rounding, we
// do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
// |x| < 2^-125. For targets without FMA instructions, we simply use
// double for intermediate results as it is more efficient than using an
// emulated version of FMA.
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
return fputil::multiply_add(x, 0x1.0p-25f, x);
#else
double xd = static_cast<double>(x);
return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
}
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Check for exceptional values
if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);
LIBC_UNLIKELY(r.has_value()))
return r.value();
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// For |x| <= 0.5, we approximate asinf(x) by:
// asin(x) = x * P(x^2)
// Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
// asin(x)/x on [0, 0.5] generated by Sollya with:
// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
// [|1, D...|], [0, 0.5]);
// An exhaustive test shows that this approximation works well up to a
// little more than 0.5.
double xd = static_cast<double>(x);
double xsq = xd * xd;
double x3 = xd * xsq;
double r = asin_eval(xsq);
return static_cast<float>(fputil::multiply_add(x3, r, xd));
}
// |x| > 1, return NaNs.
if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
if (xbits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
if (x_abs <= 0x7f80'0000U) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
}
return FPBits::quiet_nan().get_val();
}
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Check for exceptional values
if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);
LIBC_UNLIKELY(r.has_value()))
return r.value();
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// When |x| > 0.5, we perform range reduction as follow:
//
// Assume further that 0.5 < x <= 1, and let:
// y = asin(x)
// We will use the double angle formula:
// cos(2y) = 1 - 2 sin^2(y)
// and the complement angle identity:
// x = sin(y) = cos(pi/2 - y)
// = 1 - 2 sin^2 (pi/4 - y/2)
// So:
// sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
// And hence:
// pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
// Equivalently:
// asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
// Let u = (1 - x)/2, then:
// asin(x) = pi/2 - 2 * asin( sqrt(u) )
// Moreover, since 0.5 < x <= 1:
// 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
// And hence we can reuse the same polynomial approximation of asin(x) when
// |x| <= 0.5:
// asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
xbits.set_sign(Sign::POS);
double sign = SIGN[x_sign];
double xd = static_cast<double>(xbits.get_val());
double u = fputil::multiply_add(-0.5, xd, 0.5);
double c1 = sign * (-2 * fputil::sqrt<double>(u));
double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1);
double c3 = c1 * u;
double r = asin_eval(u);
return static_cast<float>(fputil::multiply_add(c3, r, c2));
}
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H

8
libc/src/math/generic/CMakeLists.txt
View File

@ -3958,13 +3958,7 @@ add_entrypoint_object(
HDRS
../asinf.h
DEPENDS
libc.src.__support.FPUtil.except_value_utils
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.inv_trigf_utils
libc.src.__support.math.asinf
)
add_entrypoint_object(

155
libc/src/math/generic/asinf.cpp
View File

@ -7,161 +7,10 @@
//===----------------------------------------------------------------------===//
#include "src/math/asinf.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.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/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
#include "src/__support/math/inv_trigf_utils.h"
#include "src/__support/math/asinf.h"
namespace LIBC_NAMESPACE_DECL {
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
static constexpr size_t N_EXCEPTS = 2;
// Exceptional values when |x| <= 0.5
static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{
// (inputs, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)
{0x3d09bf86, 0x3d09c62c, 1, 0, 1},
// x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)
{0x3de5fa1e, 0x3de6768e, 1, 0, 0},
}};
// Exceptional values when 0.5 < |x| <= 1
static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{
// (inputs, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)
{0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},
// x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)
{0x3f7741b6, 0x3fa7832a, 1, 0, 0},
}};
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
LLVM_LIBC_FUNCTION(float, asinf, (float x)) {
using namespace inv_trigf_utils_internal;
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
uint32_t x_uint = xbits.uintval();
uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
constexpr double SIGN[2] = {1.0, -1.0};
uint32_t x_sign = x_uint >> 31;
// |x| <= 0.5-ish
if (x_abs < 0x3f04'471dU) {
// |x| < 0x1.d12edp-12
if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) {
// When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
// is:
// |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
// = x^2 / 6
// < 2^-25
// < epsilon(1)/2.
// So the correctly rounded values of asin(x) are:
// = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
// or (rounding mode = FE_UPWARD and x is
// negative),
// = x otherwise.
// To simplify the rounding decision and make it more efficient, we use
// fma(x, 2^-25, x) instead.
// An exhaustive test shows that this formula work correctly for all
// rounding modes up to |x| < 0x1.d12edp-12.
// Note: to use the formula x + 2^-25*x to decide the correct rounding, we
// do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
// |x| < 2^-125. For targets without FMA instructions, we simply use
// double for intermediate results as it is more efficient than using an
// emulated version of FMA.
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
return fputil::multiply_add(x, 0x1.0p-25f, x);
#else
double xd = static_cast<double>(x);
return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
}
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Check for exceptional values
if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);
LIBC_UNLIKELY(r.has_value()))
return r.value();
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// For |x| <= 0.5, we approximate asinf(x) by:
// asin(x) = x * P(x^2)
// Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
// asin(x)/x on [0, 0.5] generated by Sollya with:
// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
// [|1, D...|], [0, 0.5]);
// An exhaustive test shows that this approximation works well up to a
// little more than 0.5.
double xd = static_cast<double>(x);
double xsq = xd * xd;
double x3 = xd * xsq;
double r = asin_eval(xsq);
return static_cast<float>(fputil::multiply_add(x3, r, xd));
}
// |x| > 1, return NaNs.
if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
if (xbits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
if (x_abs <= 0x7f80'0000U) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
}
return FPBits::quiet_nan().get_val();
}
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Check for exceptional values
if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);
LIBC_UNLIKELY(r.has_value()))
return r.value();
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// When |x| > 0.5, we perform range reduction as follow:
//
// Assume further that 0.5 < x <= 1, and let:
// y = asin(x)
// We will use the double angle formula:
// cos(2y) = 1 - 2 sin^2(y)
// and the complement angle identity:
// x = sin(y) = cos(pi/2 - y)
// = 1 - 2 sin^2 (pi/4 - y/2)
// So:
// sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
// And hence:
// pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
// Equivalently:
// asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
// Let u = (1 - x)/2, then:
// asin(x) = pi/2 - 2 * asin( sqrt(u) )
// Moreover, since 0.5 < x <= 1:
// 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
// And hence we can reuse the same polynomial approximation of asin(x) when
// |x| <= 0.5:
// asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
xbits.set_sign(Sign::POS);
double sign = SIGN[x_sign];
double xd = static_cast<double>(xbits.get_val());
double u = fputil::multiply_add(-0.5, xd, 0.5);
double c1 = sign * (-2 * fputil::sqrt<double>(u));
double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1);
double c3 = c1 * u;
double r = asin_eval(u);
return static_cast<float>(fputil::multiply_add(c3, r, c2));
}
LLVM_LIBC_FUNCTION(float, asinf, (float x)) { return math::asinf(x); }
} // namespace LIBC_NAMESPACE_DECL

1
libc/test/shared/CMakeLists.txt
View File

@ -14,6 +14,7 @@ add_fp_unittest(
libc.src.__support.math.acoshf16
libc.src.__support.math.acospif16
libc.src.__support.math.asin
libc.src.__support.math.asinf
libc.src.__support.math.erff
libc.src.__support.math.exp
libc.src.__support.math.exp10

1
libc/test/shared/shared_math_test.cpp
View File

@ -40,6 +40,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat) {
EXPECT_FP_EQ(0x1.921fb6p+0, LIBC_NAMESPACE::shared::acosf(0.0f));
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::acoshf(1.0f));
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::asinf(0.0f));
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::erff(0.0f));
EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::exp10f(0.0f));
EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::expf(0.0f));

25
utils/bazel/llvm-project-overlay/libc/BUILD.bazel
View File

@ -2231,6 +2231,22 @@ libc_support_library(
],
)
libc_support_library(
name = "__support_math_asinf",
hdrs = ["src/__support/math/asinf.h"],
deps = [
":__support_math_inv_trigf_utils",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
":__support_fputil_except_value_utils",
":__support_fputil_multiply_add",
":__support_fputil_sqrt",
":__support_macros_config",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
],
)
libc_support_library(
name = "__support_math_erff",
hdrs = ["src/__support/math/erff.h"],
@ -2786,14 +2802,7 @@ libc_math_function(
libc_math_function(
name = "asinf",
additional_deps = [
":__support_fputil_fma",
":__support_fputil_multiply_add",
":__support_fputil_nearest_integer",
":__support_fputil_polyeval",
":__support_fputil_sqrt",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
":__support_math_inv_trigf_utils",
":__support_math_asinf",
],
)