[libc-commits] [libc] [libc][math] Use float-only evaluation in asinbf16 and acosbf16 (PR #228734)
Chen Miao via libc-commits
libc-commits at lists.llvm.org
Sat Oct 3 09:14:05 PDT 2026
https://github.com/ChenMiaoi created https://github.com/llvm/llvm-project/pull/228734
Update `asinbf16` and `acosbf16` to use a shared `float` polynomial instead of the double-precision `asin_eval` helper. Reuse the polynomial coefficients from `asinf16`.
For small `acosbf16` inputs, return pi/2 rounded to `bfloat16` before squaring the input or evaluating the polynomial. This avoids spurious underflow from intermediate products.
Related: #202079
Assisted-by: GPT-6-Astra
>From 9d9b48b7807174d9447db17dafb990e200bcb479 Mon Sep 17 00:00:00 2001
From: Chen Miao <chenmiao.ku at gmail.com>
Date: Sun, 4 Oct 2026 00:01:27 +0800
Subject: [PATCH] [libc][math] Use float-only evaluation in asinbf16 and
acosbf16
Update `asinbf16` and `acosbf16` to use a shared `float` polynomial
instead of the double-precision `asin_eval` helper. Reuse the polynomial
coefficients from `asinf16`.
For small `acosbf16` inputs, return pi/2 rounded to `bfloat16` before
squaring the input or evaluating the polynomial. This avoids spurious
underflow from intermediate products.
Related: #202079
Assisted-by: GPT-6-Astra
---
libc/src/__support/math/acosbf16.h | 21 +++++++++---------
libc/src/__support/math/asin_utils.h | 10 +++++++++
libc/src/__support/math/asinbf16.h | 6 ++----
libc/test/src/math/smoke/acosbf16_test.cpp | 25 ++++++++++++++++++++++
4 files changed, 47 insertions(+), 15 deletions(-)
diff --git a/libc/src/__support/math/acosbf16.h b/libc/src/__support/math/acosbf16.h
index 7b78bd34dd28768..e0fb861accb4755 100644
--- a/libc/src/__support/math/acosbf16.h
+++ b/libc/src/__support/math/acosbf16.h
@@ -35,6 +35,14 @@ LIBC_INLINE bfloat16 acosbf16(bfloat16 x) {
uint16_t x_u = xbits.uintval();
uint16_t x_abs = x_u & 0x7fff;
bool sign = (x_u >> 15);
+
+ // For |x| <= 2^-12, 0x1.92p0 < acos(x) < 0x1.93p0: the lower bound is
+ // a bfloat16 value and the upper bound is the midpoint to its successor.
+ // Thus acos(x) rounds like pi/2 in every rounding mode. Return before
+ // squaring x or evaluating the polynomial to avoid intermediate underflow.
+ if (LIBC_UNLIKELY(x_abs <= 0x3980))
+ return fputil::cast<bfloat16>(PI_2);
+
float xf = x;
float xf_abs = (xf < 0 ? -xf : xf);
@@ -42,12 +50,7 @@ LIBC_INLINE bfloat16 acosbf16(bfloat16 x) {
// case 1: x <= 0.5
if (x_abs <= 0x3F00) {
- // |x| = {0}
- if (LIBC_UNLIKELY(x_abs == 0))
- return fputil::cast<bfloat16>(PI_2);
-
- float xp = fputil::cast<float>(asin_internal::asin_eval(x_sq));
- float result = xf * fputil::multiply_add(x_sq, xp, 1.0f);
+ float result = xf * asin_internal::asin_eval_float(x_sq);
return fputil::cast<bfloat16>(PI_2 - result);
}
@@ -66,11 +69,7 @@ LIBC_INLINE bfloat16 acosbf16(bfloat16 x) {
// and acos(x) = acos(|x|) for x >= 0, pi - acos(|x|) for x < 0
float t = fputil::multiply_add<float>(xf_abs, -0.5f, 0.5f);
float t_sqrt = fputil::sqrt<float>(t);
- // TODO: Use bfloat16 version for inv_trigf_utils_internals after they are
- // available Tracking issue :
- // https://github.com/llvm/llvm-project/issues/202079
- float tp = fputil::cast<float>(asin_internal::asin_eval(t));
- float asin_sqrt_t = t_sqrt * (fputil::multiply_add(t, tp, 1.0f));
+ float asin_sqrt_t = t_sqrt * asin_internal::asin_eval_float(t);
return fputil::cast<bfloat16>(
(sign) ? fputil::multiply_add(asin_sqrt_t, -2.0f, PI)
diff --git a/libc/src/__support/math/asin_utils.h b/libc/src/__support/math/asin_utils.h
index f6f19baa312fb32..f1b15631e262f5f 100644
--- a/libc/src/__support/math/asin_utils.h
+++ b/libc/src/__support/math/asin_utils.h
@@ -21,6 +21,16 @@ namespace LIBC_NAMESPACE_DECL {
namespace asin_internal {
+// Single-precision approximation of asin(x) / x for |x| <= 0.5, evaluated
+// at xsq = x * x. The coefficients are the same as those used by asinf16,
+// generated with Sollya:
+// > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
+// This approximation is intended for 16-bit results, not float results.
+LIBC_INLINE LIBC_CONSTEXPR float asin_eval_float(float xsq) {
+ return fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f,
+ 0x1.43b2d6p-5f, 0x1.a0d73ep-5f);
+}
+
using DoubleDouble = fputil::DoubleDouble;
using DFloat128 = fputil::DyadicFloat<128>;
diff --git a/libc/src/__support/math/asinbf16.h b/libc/src/__support/math/asinbf16.h
index 64c113e68856f54..9083ca9adcce913 100644
--- a/libc/src/__support/math/asinbf16.h
+++ b/libc/src/__support/math/asinbf16.h
@@ -55,8 +55,7 @@ LIBC_INLINE LIBC_CONSTEXPR bfloat16 asinbf16(bfloat16 x) {
return x;
}
- float xp = fputil::cast<float>(asin_internal::asin_eval(x_sq));
- float result = xf * (fputil::multiply_add<float>(x_sq, xp, 1.0f));
+ float result = xf * asin_internal::asin_eval_float(x_sq);
return fputil::cast<bfloat16>(result);
}
@@ -70,8 +69,7 @@ LIBC_INLINE LIBC_CONSTEXPR bfloat16 asinbf16(bfloat16 x) {
float t = fputil::multiply_add<float>(xf_abs, -0.5f, 0.5f);
float t_sqrt = fputil::sqrt<float>(t);
- float tp = fputil::cast<float>(asin_internal::asin_eval(t));
- float asin_sqrt_t = t_sqrt * (fputil::multiply_add<float>(t, tp, 1.0f));
+ float asin_sqrt_t = t_sqrt * asin_internal::asin_eval_float(t);
float result = fputil::multiply_add<float>(-2.0f, asin_sqrt_t, PI_2);
return fputil::cast<bfloat16>(x_sign * result);
}
diff --git a/libc/test/src/math/smoke/acosbf16_test.cpp b/libc/test/src/math/smoke/acosbf16_test.cpp
index 4ca738446b9904c..a477a61c236185b 100644
--- a/libc/test/src/math/smoke/acosbf16_test.cpp
+++ b/libc/test/src/math/smoke/acosbf16_test.cpp
@@ -49,3 +49,28 @@ class LlvmLibcAcosBf16Test : public LIBC_NAMESPACE::testing::FEnvSafeTest {
}
};
TEST_F(LlvmLibcAcosBf16Test, SpecialNumbers) { test_special_numbers(); }
+
+TEST_F(LlvmLibcAcosBf16Test, SmallInputs) {
+ using FPBits = LIBC_NAMESPACE::fputil::FPBits<bfloat16>;
+ auto acos_without_underflow = [](bfloat16 x) {
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT);
+ bfloat16 result = LIBC_NAMESPACE::acosbf16(x);
+ EXPECT_EQ(LIBC_NAMESPACE::fputil::test_except(FE_UNDERFLOW), 0);
+ EXPECT_MATH_ERRNO(0);
+ return result;
+ };
+
+ // Include subnormals, an input that caused an intermediate underflow in
+ // the float polynomial without FMA, and the small-input branch boundary.
+ const uint16_t inputs[] = {0x0001, 0x007f, 0x0080, 0x209d,
+ 0x397f, 0x3980, 0x3981};
+ const bfloat16 pi_2_lo = FPBits(uint16_t(0x3fc9)).get_val();
+ const bfloat16 pi_2_hi = FPBits(uint16_t(0x3fca)).get_val();
+ for (uint16_t bits : inputs) {
+ bfloat16 x = FPBits(bits).get_val();
+ EXPECT_FP_EQ_ALL_ROUNDING(pi_2_lo, pi_2_hi, pi_2_lo, pi_2_lo,
+ acos_without_underflow(x));
+ EXPECT_FP_EQ_ALL_ROUNDING(pi_2_lo, pi_2_hi, pi_2_lo, pi_2_lo,
+ acos_without_underflow(-x));
+ }
+}
More information about the libc-commits
mailing list