[libc-commits] [libc] [libc][math] Use float-only evaluation in asinbf16 and acosbf16 (PR #228734)

Sat Oct 3 09:14:59 PDT 2026

llvmorg-github-actions[bot] wrote:




@llvm/pr-subscribers-libc

Author: Chen Miao (ChenMiaoi)

<details>
<summary>Changes</summary>

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

---
Full diff: https://github.com/llvm/llvm-project/pull/228734.diff


4 Files Affected:

- (modified) libc/src/__support/math/acosbf16.h (+10-11) 
- (modified) libc/src/__support/math/asin_utils.h (+10) 
- (modified) libc/src/__support/math/asinbf16.h (+2-4) 
- (modified) libc/test/src/math/smoke/acosbf16_test.cpp (+25) 


``````````diff

diff --git a/libc/src/__support/math/acosbf16.h b/libc/src/__support/math/acosbf16.h
index 7b78bd34dd2876..e0fb861accb475 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 f6f19baa312fb3..f1b15631e262f5 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 64c113e68856f5..9083ca9adcce91 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 4ca738446b9904..a477a61c236185 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));
+  }
+}

``````````

</details>


https://github.com/llvm/llvm-project/pull/228734
