Skip to content

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

Merged
bassiounix merged 1 commit into
mainfrom
users/bassiounix/spr/08-03-_libc_math_refactor_cbrtf_implementation_to_header-only_in_src___support_math_folder
Aug 8, 2025
Merged

[libc][math] Refactor cbrtf implementation to header-only in src/__support/math folder.#151846
bassiounix merged 1 commit into
mainfrom
users/bassiounix/spr/08-03-_libc_math_refactor_cbrtf_implementation_to_header-only_in_src___support_math_folder

Conversation

@bassiounix

@bassiounix bassiounix commented Aug 3, 2025

Copy link
Copy Markdown
Member

bassiounix commented Aug 3, 2025

Copy link
Copy Markdown
Member Author

This stack of pull requests is managed by Graphite. Learn more about stacking.

@llvmbot

llvmbot commented Aug 3, 2025

Copy link
Copy Markdown
Member

@llvm/pr-subscribers-libc

Author: Muhammad Bassiouni (bassiounix)

Changes

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


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

9 Files Affected:

  • (modified) libc/shared/math.h (+1)
  • (added) libc/shared/math/cbrtf.h (+23)
  • (modified) libc/src/__support/math/CMakeLists.txt (+11)
  • (added) libc/src/__support/math/cbrtf.h (+161)
  • (modified) libc/src/math/generic/CMakeLists.txt (+1-5)
  • (modified) libc/src/math/generic/cbrtf.cpp (+2-145)
  • (modified) libc/test/shared/CMakeLists.txt (+1)
  • (modified) libc/test/shared/shared_math_test.cpp (+1)
  • (modified) utils/bazel/llvm-project-overlay/libc/BUILD.bazel (+13-1)
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 3714f380a27dc..ea645f0afedbc 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -31,6 +31,7 @@
 #include "math/atanhf.h"
 #include "math/atanhf16.h"
 #include "math/cbrt.h"
+#include "math/cbrtf.h"
 #include "math/erff.h"
 #include "math/exp.h"
 #include "math/exp10.h"
diff --git a/libc/shared/math/cbrtf.h b/libc/shared/math/cbrtf.h
new file mode 100644
index 0000000000000..09b86bed3fb7e
--- /dev/null
+++ b/libc/shared/math/cbrtf.h
@@ -0,0 +1,23 @@
+//===-- Shared cbrtf 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 LIBC_SHARED_MATH_CBRTF_H
+#define LIBC_SHARED_MATH_CBRTF_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/cbrtf.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::cbrtf;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LIBC_SHARED_MATH_CBRTF_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index e1076edf1e61c..fe928a8fadd5e 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -346,6 +346,17 @@ add_header_library(
     libc.src.__support.integer_literals
 )
 
+add_header_library(
+  cbrtf
+  HDRS
+    cbrtf.h
+  DEPENDS
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.macros.optimization
+)
+
 add_header_library(
   erff
   HDRS
diff --git a/libc/src/__support/math/cbrtf.h b/libc/src/__support/math/cbrtf.h
new file mode 100644
index 0000000000000..f82892bbbe61b
--- /dev/null
+++ b/libc/src/__support/math/cbrtf.h
@@ -0,0 +1,161 @@
+//===-- Implementation header for cbrtf -------------------------*- 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 LIBC_SRC___SUPPORT_MATH_CBRTF_H
+#define LIBC_SRC___SUPPORT_MATH_CBRTF_H
+
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+LIBC_INLINE static constexpr float cbrtf(float x) {
+  // Look up table for 2^(i/3) for i = 0, 1, 2.
+  constexpr double CBRT2[3] = {1.0, 0x1.428a2f98d728bp0, 0x1.965fea53d6e3dp0};
+
+  // Degree-7 polynomials approximation of ((1 + x)^(1/3) - 1)/x for 0 <= x <= 1
+  // generated by Sollya with:
+  // > for i from 0 to 15 do {
+  //     P = fpminimax(((1 + x)^(1/3) - 1)/x, 6, [|D...|], [i/16, (i + 1)/16]);
+  //     print("{", coeff(P, 0), ",", coeff(P, 1), ",", coeff(P, 2), ",",
+  //           coeff(P, 3), ",", coeff(P, 4), ",", coeff(P, 5), ",",
+  //           coeff(P, 6), "},");
+  // };
+  // Then (1 + x)^(1/3) ~ 1 + x * P(x).
+  constexpr double COEFFS[16][7] = {
+      {0x1.55555555554ebp-2, -0x1.c71c71c678c0cp-4, 0x1.f9add2776de81p-5,
+       -0x1.511e10aa964a7p-5, 0x1.ee44165937fa2p-6, -0x1.7c5c9e059345dp-6,
+       0x1.047f75e0aff14p-6},
+      {0x1.5555554d1149ap-2, -0x1.c71c676fcb5bp-4, 0x1.f9ab127dc57ebp-5,
+       -0x1.50ea8fd1d4c15p-5, 0x1.e9d68f28ced43p-6, -0x1.60e0e1e661311p-6,
+       0x1.716eca1d6e3bcp-7},
+      {0x1.5555546377d45p-2, -0x1.c71bc1c6d49d2p-4, 0x1.f9924cc0ed24dp-5,
+       -0x1.4fea3beb53b3bp-5, 0x1.de028a9a07b1bp-6, -0x1.3b090d2233524p-6,
+       0x1.0aeca34893785p-7},
+      {0x1.55554dce9f649p-2, -0x1.c7188b34b98f8p-4, 0x1.f93e1af34af49p-5,
+       -0x1.4d9a06be75c63p-5, 0x1.cb943f4f68992p-6, -0x1.139a685a5e3c4p-6,
+       0x1.88410674c6a5dp-8},
+      {0x1.5555347d211c3p-2, -0x1.c70f2a4b1a5fap-4, 0x1.f88420e8602c3p-5,
+       -0x1.49becfa4ed3ep-5, 0x1.b475cd9013162p-6, -0x1.dcfee1dd2f8efp-7,
+       0x1.249bb51a1c498p-8},
+      {0x1.5554f01b33dbap-2, -0x1.c6facb929dbf1p-4, 0x1.f73fb7861252ep-5,
+       -0x1.4459a4a0071fap-5, 0x1.9a8df2b504fc2p-6, -0x1.9a7ce3006d06ep-7,
+       0x1.ba9230918fa2ep-9},
+      {0x1.55545c695db5fp-2, -0x1.c6d6089f20275p-4, 0x1.f556e0ea80efp-5,
+       -0x1.3d91372d083f4p-5, 0x1.7f66cff331f4p-6, -0x1.606a562491737p-7,
+       0x1.52e3e17c71069p-9},
+      {0x1.55534a879232ap-2, -0x1.c69b836998b84p-4, 0x1.f2bb26dac0e4cp-5,
+       -0x1.359eed43716d7p-5, 0x1.64218cd824fbcp-6, -0x1.2e703e2e091e8p-7,
+       0x1.0677d9af6aad4p-9},
+      {0x1.5551836bb5494p-2, -0x1.c64658c15353bp-4, 0x1.ef68517451a6ep-5,
+       -0x1.2cc20a980dceep-5, 0x1.49843e0fad93ap-6, -0x1.03c59ccb68e54p-7,
+       0x1.9ad325dc7adcbp-10},
+      {0x1.554ecacb0d035p-2, -0x1.c5d2664026ffcp-4, 0x1.eb624796ba809p-5,
+       -0x1.233803d19a535p-5, 0x1.300decb1c3c28p-6, -0x1.befe18031ec3dp-8,
+       0x1.449f5ee175c69p-10},
+      {0x1.554ae1f5ae815p-2, -0x1.c53c6b14ff6b2p-4, 0x1.e6b2d5127bb5bp-5,
+       -0x1.19387336788a3p-5, 0x1.180955a6ab255p-6, -0x1.81696703ba369p-8,
+       0x1.02cb36389bd79p-10},
+      {0x1.55458a59f356ep-2, -0x1.c4820dd631ae9p-4, 0x1.e167af818bd15p-5,
+       -0x1.0ef35f6f72e52p-5, 0x1.019c33b65e4ebp-6, -0x1.4d25bdd52d3a5p-8,
+       0x1.a008ae91f5936p-11},
+      {0x1.553e878eafee1p-2, -0x1.c3a1d0b2a3db2p-4, 0x1.db90d8ed9f89bp-5,
+       -0x1.0490e20f1ae91p-5, 0x1.d9a5d1fc42fe3p-7, -0x1.20bf8227c2abfp-8,
+       0x1.50f8174cdb6e9p-11},
+      {0x1.5535a0dedf1b1p-2, -0x1.c29afb8bd01a1p-4, 0x1.d53f6371c1e27p-5,
+       -0x1.f463209b433e2p-6, 0x1.b35222a17e44p-7, -0x1.f5efbf505e133p-9,
+       0x1.12e0e94e8586dp-11},
+      {0x1.552aa25e57bfdp-2, -0x1.c16d811e4acadp-4, 0x1.ce8489b47aa51p-5,
+       -0x1.dfde7ff758ea8p-6, 0x1.901f43aac38c8p-7, -0x1.b581d07df5ad5p-9,
+       0x1.c3726535f1fc6p-12},
+      {0x1.551d5d9b204d3p-2, -0x1.c019e328f8db1p-4, 0x1.c7710f44fc3cep-5,
+       -0x1.cbbbe25ea8ba4p-6, 0x1.6fe270088623dp-7, -0x1.7e6fc79733761p-9,
+       0x1.75077abf18d84p-12},
+  };
+
+  using FloatBits = typename fputil::FPBits<float>;
+  using DoubleBits = typename fputil::FPBits<double>;
+
+  FloatBits x_bits(x);
+
+  uint32_t x_abs = x_bits.uintval() & 0x7fff'ffff;
+  uint32_t sign_bit = (x_bits.uintval() >> 31) << DoubleBits::EXP_LEN;
+
+  if (LIBC_UNLIKELY(x == 0.0f || x_abs >= 0x7f80'0000)) {
+    // x is 0, Inf, or NaN.
+    // Make sure it works for FTZ/DAZ modes.
+    return x + x;
+  }
+
+  double xd = static_cast<double>(x);
+  DoubleBits xd_bits(xd);
+
+  // When using biased exponent of x in double precision,
+  //   x_e = real_exponent_of_x + 1023
+  // Then:
+  //   x_e / 3 = real_exponent_of_x / 3 + 1023/3
+  //           = real_exponent_of_x / 3 + 341
+  // So to make it the correct biased exponent of x^(1/3), we add
+  //   1023 - 341 = 682
+  // to the quotient x_e / 3.
+  unsigned x_e = static_cast<unsigned>(xd_bits.get_biased_exponent());
+  unsigned out_e = (x_e / 3 + 682) | sign_bit;
+  unsigned shift_e = x_e % 3;
+
+  // Set x_m = 2^(x_e % 3) * (1.mantissa)
+  uint64_t x_m = xd_bits.get_mantissa();
+  // Use the leading 4 bits for look up table
+  unsigned idx = static_cast<unsigned>(x_m >> (DoubleBits::FRACTION_LEN - 4));
+
+  x_m |= static_cast<uint64_t>(DoubleBits::EXP_BIAS)
+         << DoubleBits::FRACTION_LEN;
+
+  double x_reduced = DoubleBits(x_m).get_val();
+  double dx = x_reduced - 1.0;
+
+  double dx_sq = dx * dx;
+  double c0 = fputil::multiply_add(dx, COEFFS[idx][0], 1.0);
+  double c1 = fputil::multiply_add(dx, COEFFS[idx][2], COEFFS[idx][1]);
+  double c2 = fputil::multiply_add(dx, COEFFS[idx][4], COEFFS[idx][3]);
+  double c3 = fputil::multiply_add(dx, COEFFS[idx][6], COEFFS[idx][5]);
+
+  double dx_4 = dx_sq * dx_sq;
+  double p0 = fputil::multiply_add(dx_sq, c1, c0);
+  double p1 = fputil::multiply_add(dx_sq, c3, c2);
+
+  double r = fputil::multiply_add(dx_4, p1, p0) * CBRT2[shift_e];
+
+  uint64_t r_m = DoubleBits(r).get_mantissa();
+  // Check if the output is exact.  To be exact, the smallest 1-bit of the
+  // output has to be at least 2^-7 or higher.  So we check the lowest 44 bits
+  // to see if they are within 2^(-52 + 3) errors from all zeros, then the
+  // result cube root is exact.
+  if (LIBC_UNLIKELY(((r_m + 8) & 0xfffffffffff) <= 16)) {
+    if ((r_m & 0xfffffffffff) <= 8)
+      r_m &= 0xffff'ffff'ffff'ffe0;
+    else
+      r_m = (r_m & 0xffff'ffff'ffff'ffe0) + 0x20;
+    fputil::clear_except_if_required(FE_INEXACT);
+  }
+  // Adjust exponent and sign.
+  uint64_t r_bits =
+      r_m | (static_cast<uint64_t>(out_e) << DoubleBits::FRACTION_LEN);
+
+  return static_cast<float>(DoubleBits(r_bits).get_val());
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LIBC_SRC___SUPPORT_MATH_CBRTF_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index a86619576cfc6..44c17cf5b0770 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4739,11 +4739,7 @@ add_entrypoint_object(
   HDRS
     ../cbrtf.h
   DEPENDS
-    libc.hdr.fenv_macros
-    libc.src.__support.FPUtil.fenv_impl
-    libc.src.__support.FPUtil.fp_bits
-    libc.src.__support.FPUtil.multiply_add
-    libc.src.__support.macros.optimization
+    libc.src.__support.math.cbrtf
 )
 
 add_entrypoint_object(
diff --git a/libc/src/math/generic/cbrtf.cpp b/libc/src/math/generic/cbrtf.cpp
index 71b23c4a8c742..0bd8f71365300 100644
--- a/libc/src/math/generic/cbrtf.cpp
+++ b/libc/src/math/generic/cbrtf.cpp
@@ -7,153 +7,10 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/cbrtf.h"
-#include "hdr/fenv_macros.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/multiply_add.h"
-#include "src/__support/common.h"
-#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/math/cbrtf.h"
 
 namespace LIBC_NAMESPACE_DECL {
 
-namespace {
-
-// Look up table for 2^(i/3) for i = 0, 1, 2.
-constexpr double CBRT2[3] = {1.0, 0x1.428a2f98d728bp0, 0x1.965fea53d6e3dp0};
-
-// Degree-7 polynomials approximation of ((1 + x)^(1/3) - 1)/x for 0 <= x <= 1
-// generated by Sollya with:
-// > for i from 0 to 15 do {
-//     P = fpminimax(((1 + x)^(1/3) - 1)/x, 6, [|D...|], [i/16, (i + 1)/16]);
-//     print("{", coeff(P, 0), ",", coeff(P, 1), ",", coeff(P, 2), ",",
-//           coeff(P, 3), ",", coeff(P, 4), ",", coeff(P, 5), ",",
-//           coeff(P, 6), "},");
-// };
-// Then (1 + x)^(1/3) ~ 1 + x * P(x).
-constexpr double COEFFS[16][7] = {
-    {0x1.55555555554ebp-2, -0x1.c71c71c678c0cp-4, 0x1.f9add2776de81p-5,
-     -0x1.511e10aa964a7p-5, 0x1.ee44165937fa2p-6, -0x1.7c5c9e059345dp-6,
-     0x1.047f75e0aff14p-6},
-    {0x1.5555554d1149ap-2, -0x1.c71c676fcb5bp-4, 0x1.f9ab127dc57ebp-5,
-     -0x1.50ea8fd1d4c15p-5, 0x1.e9d68f28ced43p-6, -0x1.60e0e1e661311p-6,
-     0x1.716eca1d6e3bcp-7},
-    {0x1.5555546377d45p-2, -0x1.c71bc1c6d49d2p-4, 0x1.f9924cc0ed24dp-5,
-     -0x1.4fea3beb53b3bp-5, 0x1.de028a9a07b1bp-6, -0x1.3b090d2233524p-6,
-     0x1.0aeca34893785p-7},
-    {0x1.55554dce9f649p-2, -0x1.c7188b34b98f8p-4, 0x1.f93e1af34af49p-5,
-     -0x1.4d9a06be75c63p-5, 0x1.cb943f4f68992p-6, -0x1.139a685a5e3c4p-6,
-     0x1.88410674c6a5dp-8},
-    {0x1.5555347d211c3p-2, -0x1.c70f2a4b1a5fap-4, 0x1.f88420e8602c3p-5,
-     -0x1.49becfa4ed3ep-5, 0x1.b475cd9013162p-6, -0x1.dcfee1dd2f8efp-7,
-     0x1.249bb51a1c498p-8},
-    {0x1.5554f01b33dbap-2, -0x1.c6facb929dbf1p-4, 0x1.f73fb7861252ep-5,
-     -0x1.4459a4a0071fap-5, 0x1.9a8df2b504fc2p-6, -0x1.9a7ce3006d06ep-7,
-     0x1.ba9230918fa2ep-9},
-    {0x1.55545c695db5fp-2, -0x1.c6d6089f20275p-4, 0x1.f556e0ea80efp-5,
-     -0x1.3d91372d083f4p-5, 0x1.7f66cff331f4p-6, -0x1.606a562491737p-7,
-     0x1.52e3e17c71069p-9},
-    {0x1.55534a879232ap-2, -0x1.c69b836998b84p-4, 0x1.f2bb26dac0e4cp-5,
-     -0x1.359eed43716d7p-5, 0x1.64218cd824fbcp-6, -0x1.2e703e2e091e8p-7,
-     0x1.0677d9af6aad4p-9},
-    {0x1.5551836bb5494p-2, -0x1.c64658c15353bp-4, 0x1.ef68517451a6ep-5,
-     -0x1.2cc20a980dceep-5, 0x1.49843e0fad93ap-6, -0x1.03c59ccb68e54p-7,
-     0x1.9ad325dc7adcbp-10},
-    {0x1.554ecacb0d035p-2, -0x1.c5d2664026ffcp-4, 0x1.eb624796ba809p-5,
-     -0x1.233803d19a535p-5, 0x1.300decb1c3c28p-6, -0x1.befe18031ec3dp-8,
-     0x1.449f5ee175c69p-10},
-    {0x1.554ae1f5ae815p-2, -0x1.c53c6b14ff6b2p-4, 0x1.e6b2d5127bb5bp-5,
-     -0x1.19387336788a3p-5, 0x1.180955a6ab255p-6, -0x1.81696703ba369p-8,
-     0x1.02cb36389bd79p-10},
-    {0x1.55458a59f356ep-2, -0x1.c4820dd631ae9p-4, 0x1.e167af818bd15p-5,
-     -0x1.0ef35f6f72e52p-5, 0x1.019c33b65e4ebp-6, -0x1.4d25bdd52d3a5p-8,
-     0x1.a008ae91f5936p-11},
-    {0x1.553e878eafee1p-2, -0x1.c3a1d0b2a3db2p-4, 0x1.db90d8ed9f89bp-5,
-     -0x1.0490e20f1ae91p-5, 0x1.d9a5d1fc42fe3p-7, -0x1.20bf8227c2abfp-8,
-     0x1.50f8174cdb6e9p-11},
-    {0x1.5535a0dedf1b1p-2, -0x1.c29afb8bd01a1p-4, 0x1.d53f6371c1e27p-5,
-     -0x1.f463209b433e2p-6, 0x1.b35222a17e44p-7, -0x1.f5efbf505e133p-9,
-     0x1.12e0e94e8586dp-11},
-    {0x1.552aa25e57bfdp-2, -0x1.c16d811e4acadp-4, 0x1.ce8489b47aa51p-5,
-     -0x1.dfde7ff758ea8p-6, 0x1.901f43aac38c8p-7, -0x1.b581d07df5ad5p-9,
-     0x1.c3726535f1fc6p-12},
-    {0x1.551d5d9b204d3p-2, -0x1.c019e328f8db1p-4, 0x1.c7710f44fc3cep-5,
-     -0x1.cbbbe25ea8ba4p-6, 0x1.6fe270088623dp-7, -0x1.7e6fc79733761p-9,
-     0x1.75077abf18d84p-12},
-};
-
-} // anonymous namespace
-
-LLVM_LIBC_FUNCTION(float, cbrtf, (float x)) {
-  using FloatBits = typename fputil::FPBits<float>;
-  using DoubleBits = typename fputil::FPBits<double>;
-
-  FloatBits x_bits(x);
-
-  uint32_t x_abs = x_bits.uintval() & 0x7fff'ffff;
-  uint32_t sign_bit = (x_bits.uintval() >> 31) << DoubleBits::EXP_LEN;
-
-  if (LIBC_UNLIKELY(x == 0.0f || x_abs >= 0x7f80'0000)) {
-    // x is 0, Inf, or NaN.
-    // Make sure it works for FTZ/DAZ modes.
-    return x + x;
-  }
-
-  double xd = static_cast<double>(x);
-  DoubleBits xd_bits(xd);
-
-  // When using biased exponent of x in double precision,
-  //   x_e = real_exponent_of_x + 1023
-  // Then:
-  //   x_e / 3 = real_exponent_of_x / 3 + 1023/3
-  //           = real_exponent_of_x / 3 + 341
-  // So to make it the correct biased exponent of x^(1/3), we add
-  //   1023 - 341 = 682
-  // to the quotient x_e / 3.
-  unsigned x_e = static_cast<unsigned>(xd_bits.get_biased_exponent());
-  unsigned out_e = (x_e / 3 + 682) | sign_bit;
-  unsigned shift_e = x_e % 3;
-
-  // Set x_m = 2^(x_e % 3) * (1.mantissa)
-  uint64_t x_m = xd_bits.get_mantissa();
-  // Use the leading 4 bits for look up table
-  unsigned idx = static_cast<unsigned>(x_m >> (DoubleBits::FRACTION_LEN - 4));
-
-  x_m |= static_cast<uint64_t>(DoubleBits::EXP_BIAS)
-         << DoubleBits::FRACTION_LEN;
-
-  double x_reduced = DoubleBits(x_m).get_val();
-  double dx = x_reduced - 1.0;
-
-  double dx_sq = dx * dx;
-  double c0 = fputil::multiply_add(dx, COEFFS[idx][0], 1.0);
-  double c1 = fputil::multiply_add(dx, COEFFS[idx][2], COEFFS[idx][1]);
-  double c2 = fputil::multiply_add(dx, COEFFS[idx][4], COEFFS[idx][3]);
-  double c3 = fputil::multiply_add(dx, COEFFS[idx][6], COEFFS[idx][5]);
-
-  double dx_4 = dx_sq * dx_sq;
-  double p0 = fputil::multiply_add(dx_sq, c1, c0);
-  double p1 = fputil::multiply_add(dx_sq, c3, c2);
-
-  double r = fputil::multiply_add(dx_4, p1, p0) * CBRT2[shift_e];
-
-  uint64_t r_m = DoubleBits(r).get_mantissa();
-  // Check if the output is exact.  To be exact, the smallest 1-bit of the
-  // output has to be at least 2^-7 or higher.  So we check the lowest 44 bits
-  // to see if they are within 2^(-52 + 3) errors from all zeros, then the
-  // result cube root is exact.
-  if (LIBC_UNLIKELY(((r_m + 8) & 0xfffffffffff) <= 16)) {
-    if ((r_m & 0xfffffffffff) <= 8)
-      r_m &= 0xffff'ffff'ffff'ffe0;
-    else
-      r_m = (r_m & 0xffff'ffff'ffff'ffe0) + 0x20;
-    fputil::clear_except_if_required(FE_INEXACT);
-  }
-  // Adjust exponent and sign.
-  uint64_t r_bits =
-      r_m | (static_cast<uint64_t>(out_e) << DoubleBits::FRACTION_LEN);
-
-  return static_cast<float>(DoubleBits(r_bits).get_val());
-}
+LLVM_LIBC_FUNCTION(float, cbrtf, (float x)) { return math::cbrtf(x); }
 
 } // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index f5ea510180366..9685aeaee7d8d 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -27,6 +27,7 @@ add_fp_unittest(
     libc.src.__support.math.atanhf
     libc.src.__support.math.atanhf16
     libc.src.__support.math.cbrt
+    libc.src.__support.math.cbrtf
     libc.src.__support.math.erff
     libc.src.__support.math.exp
     libc.src.__support.math.exp10
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 3d64e5e68882b..5e57c49ed23d9 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -49,6 +49,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat) {
   EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::atan2f(0.0f, 0.0f));
   EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::atanf(0.0f));
   EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::atanhf(0.0f));
+  EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::cbrtf(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));
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index b46f334512979..d30c725524d17 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -2412,6 +2412,18 @@ libc_support_library(
     ],
 )
 
+libc_support_library(
+    name = "__support_math_cbrtf",
+    hdrs = ["src/__support/math/cbrtf.h"],
+    deps = [
+        ":__support_fputil_double_double",
+        ":__support_fputil_polyeval",
+        ":__support_fputil_fenv_impl",
+        ":__support_fputil_dyadic_float",
+        ":__support_integer_literals",
+    ],
+)
+
 libc_support_library(
     name = "__support_math_erff",
     hdrs = ["src/__support/math/erff.h"],
@@ -3088,7 +3100,7 @@ libc_math_function(
 libc_math_function(
     name = "cbrtf",
     additional_deps = [
-        ":__support_fputil_polyeval",
+        ":__support_math_cbrtf",
     ],
 )
 

@bassiounix bassiounix force-pushed the users/bassiounix/spr/08-02-_libc_math_refactor_cbrt_implementation_to_header-only_in_src___support_math_folder branch from 150b50e to b200607 Compare August 3, 2025 17:54
Base automatically changed from users/bassiounix/spr/08-02-_libc_math_refactor_cbrt_implementation_to_header-only_in_src___support_math_folder to main August 3, 2025 18:02
@bassiounix bassiounix force-pushed the users/bassiounix/spr/08-03-_libc_math_refactor_cbrtf_implementation_to_header-only_in_src___support_math_folder branch 2 times, most recently from 2539a02 to 4aa6063 Compare August 3, 2025 19:25
@bassiounix bassiounix force-pushed the users/bassiounix/spr/08-03-_libc_math_refactor_cbrtf_implementation_to_header-only_in_src___support_math_folder branch from 4aa6063 to ae312dc Compare August 5, 2025 03:31
@bassiounix bassiounix force-pushed the users/bassiounix/spr/08-03-_libc_math_refactor_cbrtf_implementation_to_header-only_in_src___support_math_folder branch 7 times, most recently from 847602c to fd28f6c Compare August 8, 2025 15:13
@bassiounix bassiounix force-pushed the users/bassiounix/spr/08-03-_libc_math_refactor_cbrtf_implementation_to_header-only_in_src___support_math_folder branch from fd28f6c to 50a8445 Compare August 8, 2025 15:22
@bassiounix bassiounix merged commit 66734f4 into main Aug 8, 2025
19 checks passed
@bassiounix bassiounix deleted the users/bassiounix/spr/08-03-_libc_math_refactor_cbrtf_implementation_to_header-only_in_src___support_math_folder branch August 8, 2025 15:28
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bazel "Peripheral" support tier build system: utils/bazel libc

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants