Skip to content

[libc][math] Refactor asinf16 implementation to header-only in src/__support/math folder. #150800

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jul 27, 2025
Merged

[libc][math] Refactor asinf16 implementation to header-only in src/__support/math folder. #150800

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions libc/shared/math.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include "math/acospif16.h"
#include "math/asin.h"
#include "math/asinf.h"
#include "math/asinf16.h"
#include "math/erff.h"
#include "math/exp.h"
#include "math/exp10.h"
Expand Down
28 changes: 28 additions & 0 deletions libc/shared/math/asinf16.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
//===-- Shared asinf16 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_ASINF16_H
#define LLVM_LIBC_SHARED_MATH_ASINF16_H

#include "shared/libc_common.h"

#ifdef LIBC_TYPES_HAS_FLOAT16

#include "src/__support/math/asinf16.h"

namespace LIBC_NAMESPACE_DECL {
namespace shared {

using math::asinf16;

} // namespace shared
} // namespace LIBC_NAMESPACE_DECL

#endif // LIBC_TYPES_HAS_FLOAT16

#endif // LLVM_LIBC_SHARED_MATH_ASINF_H
14 changes: 14 additions & 0 deletions libc/src/__support/math/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,20 @@ add_header_library(
libc.src.__support.macros.properties.cpu_features
)

add_header_library(
asinf16
HDRS
asinf16.h
DEPENDS
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.polyeval
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.optimization
)

add_header_library(
erff
HDRS
Expand Down
146 changes: 146 additions & 0 deletions libc/src/__support/math/asinf16.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
//===-- Implementation header for asinf16 -----------------------*- 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_ASINF16_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF16_H

#include "include/llvm-libc-macros/float16-macros.h"

#ifdef LIBC_TYPES_HAS_FLOAT16

#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"

namespace LIBC_NAMESPACE_DECL {

namespace math {

LIBC_INLINE static constexpr float16 asinf16(float16 x) {

// Generated by Sollya using the following command:
// > round(pi/2, D, RN);
constexpr float PI_2 = 0x1.921fb54442d18p0f;

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));
}

} // namespace math

} // namespace LIBC_NAMESPACE_DECL

#endif // LIBC_TYPES_HAS_FLOAT16

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF16_H
11 changes: 1 addition & 10 deletions libc/src/math/generic/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3968,16 +3968,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(
Expand Down
121 changes: 2 additions & 119 deletions libc/src/math/generic/asinf16.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
1 change: 1 addition & 0 deletions libc/test/shared/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ add_fp_unittest(
libc.src.__support.math.acospif16
libc.src.__support.math.asin
libc.src.__support.math.asinf
libc.src.__support.math.asinf16
libc.src.__support.math.erff
libc.src.__support.math.exp
libc.src.__support.math.exp10
Expand Down
1 change: 1 addition & 0 deletions libc/test/shared/shared_math_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) {

EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::acoshf16(1.0f16));
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::acospif16(1.0f16));
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::asinf16(0.0f16));

EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::exp10f16(0.0f16));

Expand Down
Loading
Loading