[libc][math] Refactor cbrt implementation to header-only in src/__support/math folder. #151837
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This was referenced Aug 2, 2025
Merged
Merged
bassiounix
added
bazel
Graphite App
This was referenced Aug 2, 2025
Merged
Merged
Merged
Merged
|
@llvm/pr-subscribers-libc Author: Muhammad Bassiouni (bassiounix) ChangesPatch is 30.62 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/151837.diff 9 Files Affected:
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 7fb736b78efa5..3714f380a27dc 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -30,6 +30,7 @@
#include "math/atanf16.h"
#include "math/atanhf.h"
#include "math/atanhf16.h"
+#include "math/cbrt.h"
#include "math/erff.h"
#include "math/exp.h"
#include "math/exp10.h"
diff --git a/libc/shared/math/cbrt.h b/libc/shared/math/cbrt.h
new file mode 100644
index 0000000000000..2f49dbd364328
--- /dev/null
+++ b/libc/shared/math/cbrt.h
@@ -0,0 +1,23 @@
+//===-- Shared cbrt 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_CBRT_H
+#define LLVM_LIBC_SHARED_MATH_CBRT_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/cbrt.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::cbrt;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_CBRT_H
\ No newline at end of file
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 9631ab5be7d3b..e1076edf1e61c 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -331,6 +331,21 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ cbrt
+ HDRS
+ cbrt.h
+ DEPENDS
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.dyadic_float
+ 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.macros.optimization
+ libc.src.__support.integer_literals
+)
+
add_header_library(
erff
HDRS
diff --git a/libc/src/__support/math/cbrt.h b/libc/src/__support/math/cbrt.h
new file mode 100644
index 0000000000000..2b9a73c823b14
--- /dev/null
+++ b/libc/src/__support/math/cbrt.h
@@ -0,0 +1,349 @@
+//===-- Implementation header for erff --------------------------*- 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_CBRT_H
+#define LIBC_SRC___SUPPORT_MATH_CBRT_H
+
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/integer_literals.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0)
+#define LIBC_MATH_CBRT_SKIP_ACCURATE_PASS
+#endif
+
+namespace cbrt_internal {
+using namespace fputil;
+
+// Initial approximation of x^(-2/3) for 1 <= x < 2.
+// Polynomial generated by Sollya with:
+// > P = fpminimax(x^(-2/3), 7, [|D...|], [1, 2]);
+// > dirtyinfnorm(P/x^(-2/3) - 1, [1, 2]);
+// 0x1.28...p-21
+LIBC_INLINE static double intial_approximation(double x) {
+ constexpr double COEFFS[8] = {
+ 0x1.bc52aedead5c6p1, -0x1.b52bfebf110b3p2, 0x1.1d8d71d53d126p3,
+ -0x1.de2db9e81cf87p2, 0x1.0154ca06153bdp2, -0x1.5973c66ee6da7p0,
+ 0x1.07bf6ac832552p-2, -0x1.5e53d9ce41cb8p-6,
+ };
+
+ double x_sq = x * x;
+
+ double c0 = fputil::multiply_add(x, COEFFS[1], COEFFS[0]);
+ double c1 = fputil::multiply_add(x, COEFFS[3], COEFFS[2]);
+ double c2 = fputil::multiply_add(x, COEFFS[5], COEFFS[4]);
+ double c3 = fputil::multiply_add(x, COEFFS[7], COEFFS[6]);
+
+ double x_4 = x_sq * x_sq;
+ double d0 = fputil::multiply_add(x_sq, c1, c0);
+ double d1 = fputil::multiply_add(x_sq, c3, c2);
+
+ return fputil::multiply_add(x_4, d1, d0);
+}
+
+// Get the error term for Newton iteration:
+// h(x) = x^3 * a^2 - 1,
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+LIBC_INLINE static double get_error(const DoubleDouble &x_3,
+ const DoubleDouble &a_sq) {
+ return fputil::multiply_add(x_3.hi, a_sq.hi, -1.0) +
+ fputil::multiply_add(x_3.lo, a_sq.hi, x_3.hi * a_sq.lo);
+}
+#else
+LIBC_INLINE static constexpr double get_error(const DoubleDouble &x_3,
+ const DoubleDouble &a_sq) {
+ DoubleDouble x_3_a_sq = fputil::quick_mult(a_sq, x_3);
+ return (x_3_a_sq.hi - 1.0) + x_3_a_sq.lo;
+}
+#endif
+
+} // namespace cbrt_internal
+
+// Correctly rounded cbrt algorithm:
+//
+// === Step 1 - Range reduction ===
+// For x = (-1)^s * 2^e * (1.m), we get 2 reduced arguments x_r and a as:
+// x_r = 1.m
+// a = (-1)^s * 2^(e % 3) * (1.m)
+// Then cbrt(x) = x^(1/3) can be computed as:
+// x^(1/3) = 2^(e / 3) * a^(1/3).
+//
+// In order to avoid division, we compute a^(-2/3) using Newton method and then
+// multiply the results by a:
+// a^(1/3) = a * a^(-2/3).
+//
+// === Step 2 - First approximation to a^(-2/3) ===
+// First, we use a degree-7 minimax polynomial generated by Sollya to
+// approximate x_r^(-2/3) for 1 <= x_r < 2.
+// p = P(x_r) ~ x_r^(-2/3),
+// with relative errors bounded by:
+// | p / x_r^(-2/3) - 1 | < 1.16 * 2^-21.
+//
+// Then we multiply with 2^(e % 3) from a small lookup table to get:
+// x_0 = 2^(-2*(e % 3)/3) * p
+// ~ 2^(-2*(e % 3)/3) * x_r^(-2/3)
+// = a^(-2/3)
+// With relative errors:
+// | x_0 / a^(-2/3) - 1 | < 1.16 * 2^-21.
+// This step is done in double precision.
+//
+// === Step 3 - First Newton iteration ===
+// We follow the method described in:
+// Sibidanov, A. and Zimmermann, P., "Correctly rounded cubic root evaluation
+// in double precision", https://core-math.gitlabpages.inria.fr/cbrt64.pdf
+// to derive multiplicative Newton iterations as below:
+// Let x_n be the nth approximation to a^(-2/3). Define the n^th error as:
+// h_n = x_n^3 * a^2 - 1
+// Then:
+// a^(-2/3) = x_n / (1 + h_n)^(1/3)
+// = x_n * (1 - (1/3) * h_n + (2/9) * h_n^2 - (14/81) * h_n^3 + ...)
+// using the Taylor series expansion of (1 + h_n)^(-1/3).
+//
+// Apply to x_0 above:
+// h_0 = x_0^3 * a^2 - 1
+// = a^2 * (x_0 - a^(-2/3)) * (x_0^2 + x_0 * a^(-2/3) + a^(-4/3)),
+// it's bounded by:
+// |h_0| < 4 * 3 * 1.16 * 2^-21 * 4 < 2^-17.
+// So in the first iteration step, we use:
+// x_1 = x_0 * (1 - (1/3) * h_n + (2/9) * h_n^2 - (14/81) * h_n^3)
+// Its relative error is bounded by:
+// | x_1 / a^(-2/3) - 1 | < 35/242 * |h_0|^4 < 2^-70.
+// Then we perform Ziv's rounding test and check if the answer is exact.
+// This step is done in double-double precision.
+//
+// === Step 4 - Second Newton iteration ===
+// If the Ziv's rounding test from the previous step fails, we define the error
+// term:
+// h_1 = x_1^3 * a^2 - 1,
+// And perform another iteration:
+// x_2 = x_1 * (1 - h_1 / 3)
+// with the relative errors exceed the precision of double-double.
+// We then check the Ziv's accuracy test with relative errors < 2^-102 to
+// compensate for rounding errors.
+//
+// === Step 5 - Final iteration ===
+// If the Ziv's accuracy test from the previous step fails, we perform another
+// iteration in 128-bit precision and check for exact outputs.
+//
+// TODO: It is possible to replace this costly computation step with special
+// exceptional handling, similar to what was done in the CORE-MATH project:
+// https://gitlab.inria.fr/core-math/core-math/-/blob/master/src/binary64/cbrt/cbrt.c
+
+LIBC_INLINE static constexpr double cbrt(double x) {
+ using DoubleDouble = fputil::DoubleDouble;
+ using Float128 = fputil::DyadicFloat<128>;
+ using namespace cbrt_internal;
+ using FPBits = fputil::FPBits<double>;
+
+ uint64_t x_abs = FPBits(x).abs().uintval();
+
+ unsigned exp_bias_correction = 682; // 1023 * 2/3
+
+ if (LIBC_UNLIKELY(x_abs < FPBits::min_normal().uintval() ||
+ x_abs >= FPBits::inf().uintval())) {
+ if (x == 0.0 || x_abs >= FPBits::inf().uintval())
+ // x is 0, Inf, or NaN.
+ // Make sure it works for FTZ/DAZ modes.
+ return static_cast<double>(x + x);
+
+ // x is non-zero denormal number.
+ // Normalize x.
+ x *= 0x1.0p60;
+ exp_bias_correction -= 20;
+ }
+
+ FPBits x_bits(x);
+
+ // 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>(x_bits.get_biased_exponent());
+ unsigned out_e = (x_e / 3 + exp_bias_correction);
+ unsigned shift_e = x_e % 3;
+
+ // Set x_r = 1.mantissa
+ double x_r =
+ FPBits(x_bits.get_mantissa() |
+ (static_cast<uint64_t>(FPBits::EXP_BIAS) << FPBits::FRACTION_LEN))
+ .get_val();
+
+ // Set a = (-1)^x_sign * 2^(x_e % 3) * (1.mantissa)
+ uint64_t a_bits = x_bits.uintval() & 0x800F'FFFF'FFFF'FFFF;
+ a_bits |=
+ (static_cast<uint64_t>(shift_e + static_cast<unsigned>(FPBits::EXP_BIAS))
+ << FPBits::FRACTION_LEN);
+ double a = FPBits(a_bits).get_val();
+
+ // Initial approximation of x_r^(-2/3).
+ double p = intial_approximation(x_r);
+
+ // Look up for 2^(-2*n/3) used for first approximation step.
+ constexpr double EXP2_M2_OVER_3[3] = {1.0, 0x1.428a2f98d728bp-1,
+ 0x1.965fea53d6e3dp-2};
+
+ // x0 is an initial approximation of a^(-2/3) for 1 <= |a| < 8.
+ // Relative error: < 1.16 * 2^(-21).
+ double x0 = static_cast<double>(EXP2_M2_OVER_3[shift_e] * p);
+
+ // First iteration in double precision.
+ DoubleDouble a_sq = fputil::exact_mult(a, a);
+
+ // h0 = x0^3 * a^2 - 1
+ DoubleDouble x0_sq = fputil::exact_mult(x0, x0);
+ DoubleDouble x0_3 = fputil::quick_mult(x0, x0_sq);
+
+ double h0 = get_error(x0_3, a_sq);
+
+#ifdef LIBC_MATH_CBRT_SKIP_ACCURATE_PASS
+ constexpr double REL_ERROR = 0;
+#else
+ constexpr double REL_ERROR = 0x1.0p-51;
+#endif // LIBC_MATH_CBRT_SKIP_ACCURATE_PASS
+
+ // Taylor polynomial of (1 + h)^(-1/3):
+ // (1 + h)^(-1/3) = 1 - h/3 + 2 h^2 / 9 - 14 h^3 / 81 + ...
+ constexpr double ERR_COEFFS[3] = {
+ -0x1.5555555555555p-2 - REL_ERROR, // -1/3 - relative_error
+ 0x1.c71c71c71c71cp-3, // 2/9
+ -0x1.61f9add3c0ca4p-3, // -14/81
+ };
+ // e0 = -14 * h^2 / 81 + 2 * h / 9 - 1/3 - relative_error.
+ double e0 = fputil::polyeval(h0, ERR_COEFFS[0], ERR_COEFFS[1], ERR_COEFFS[2]);
+ double x0_h0 = x0 * h0;
+
+ // x1 = x0 (1 - h0/3 + 2 h0^2 / 9 - 14 h0^3 / 81)
+ // x1 approximate a^(-2/3) with relative errors bounded by:
+ // | x1 / a^(-2/3) - 1 | < (34/243) h0^4 < h0 * REL_ERROR
+ DoubleDouble x1_dd{x0_h0 * e0, x0};
+
+ // r1 = x1 * a ~ a^(-2/3) * a = a^(1/3).
+ DoubleDouble r1 = fputil::quick_mult(a, x1_dd);
+
+ // Lambda function to update the exponent of the result.
+ auto update_exponent = [=](double r) -> double {
+ uint64_t r_m = FPBits(r).uintval() - 0x3FF0'0000'0000'0000;
+ // Adjust exponent and sign.
+ uint64_t r_bits =
+ r_m + (static_cast<uint64_t>(out_e) << FPBits::FRACTION_LEN);
+ return FPBits(r_bits).get_val();
+ };
+
+#ifdef LIBC_MATH_CBRT_SKIP_ACCURATE_PASS
+ // TODO: We probably don't need to use double-double if accurate tests and
+ // passes are skipped.
+ return update_exponent(r1.hi + r1.lo);
+#else
+ // Accurate checks and passes.
+ double r1_lower = r1.hi + r1.lo;
+ double r1_upper =
+ r1.hi + fputil::multiply_add(x0_h0, 2.0 * REL_ERROR * a, r1.lo);
+
+ // Ziv's accuracy test.
+ if (LIBC_LIKELY(r1_upper == r1_lower)) {
+ // Test for exact outputs.
+ // Check if lower (52 - 17 = 35) bits are 0's.
+ if (LIBC_UNLIKELY((FPBits(r1_lower).uintval() & 0x0000'0007'FFFF'FFFF) ==
+ 0)) {
+ double r1_err = (r1_lower - r1.hi) - r1.lo;
+ if (FPBits(r1_err).abs().get_val() < 0x1.0p69)
+ fputil::clear_except_if_required(FE_INEXACT);
+ }
+
+ return update_exponent(r1_lower);
+ }
+
+ // Accuracy test failed, perform another Newton iteration.
+ double x1 = x1_dd.hi + (e0 + REL_ERROR) * x0_h0;
+
+ // Second iteration in double-double precision.
+ // h1 = x1^3 * a^2 - 1.
+ DoubleDouble x1_sq = fputil::exact_mult(x1, x1);
+ DoubleDouble x1_3 = fputil::quick_mult(x1, x1_sq);
+ double h1 = get_error(x1_3, a_sq);
+
+ // e1 = -x1*h1/3.
+ double e1 = h1 * (x1 * -0x1.5555555555555p-2);
+ // x2 = x1*(1 - h1/3) = x1 + e1 ~ a^(-2/3) with relative errors < 2^-101.
+ DoubleDouble x2 = fputil::exact_add(x1, e1);
+ // r2 = a * x2 ~ a * a^(-2/3) = a^(1/3) with relative errors < 2^-100.
+ DoubleDouble r2 = fputil::quick_mult(a, x2);
+
+ double r2_upper = r2.hi + fputil::multiply_add(a, 0x1.0p-102, r2.lo);
+ double r2_lower = r2.hi + fputil::multiply_add(a, -0x1.0p-102, r2.lo);
+
+ // Ziv's accuracy test.
+ if (LIBC_LIKELY(r2_upper == r2_lower))
+ return update_exponent(r2_upper);
+
+ // TODO: Investigate removing float128 and just list exceptional cases.
+ // Apply another Newton iteration with ~126-bit accuracy.
+ Float128 x2_f128 = fputil::quick_add(Float128(x2.hi), Float128(x2.lo));
+ // x2^3
+ Float128 x2_3 =
+ fputil::quick_mul(fputil::quick_mul(x2_f128, x2_f128), x2_f128);
+ // a^2
+ Float128 a_sq_f128 = fputil::quick_mul(Float128(a), Float128(a));
+ // x2^3 * a^2
+ Float128 x2_3_a_sq = fputil::quick_mul(x2_3, a_sq_f128);
+ // h2 = x2^3 * a^2 - 1
+ Float128 h2_f128 = fputil::quick_add(x2_3_a_sq, Float128(-1.0));
+ double h2 = static_cast<double>(h2_f128);
+ // t2 = 1 - h2 / 3
+ Float128 t2 =
+ fputil::quick_add(Float128(1.0), Float128(h2 * (-0x1.5555555555555p-2)));
+ // x3 = x2 * (1 - h2 / 3) ~ a^(-2/3)
+ Float128 x3 = fputil::quick_mul(x2_f128, t2);
+ // r3 = a * x3 ~ a * a^(-2/3) = a^(1/3)
+ Float128 r3 = fputil::quick_mul(Float128(a), x3);
+
+ // Check for exact cases:
+ Float128::MantissaType rounding_bits =
+ r3.mantissa & 0x0000'0000'0000'03FF'FFFF'FFFF'FFFF'FFFF_u128;
+
+ double result = static_cast<double>(r3);
+ if ((rounding_bits < 0x0000'0000'0000'0000'0000'0000'0000'000F_u128) ||
+ (rounding_bits >= 0x0000'0000'0000'03FF'FFFF'FFFF'FFFF'FFF0_u128)) {
+ // Output is exact.
+ r3.mantissa &= 0xFFFF'FFFF'FFFF'FFFF'FFFF'FFFF'FFFF'FFF0_u128;
+
+ if (rounding_bits >= 0x0000'0000'0000'03FF'FFFF'FFFF'FFFF'FFF0_u128) {
+ Float128 tmp{r3.sign, r3.exponent - 123,
+ 0x8000'0000'0000'0000'0000'0000'0000'0000_u128};
+ Float128 r4 = fputil::quick_add(r3, tmp);
+ result = static_cast<double>(r4);
+ } else {
+ result = static_cast<double>(r3);
+ }
+
+ fputil::clear_except_if_required(FE_INEXACT);
+ }
+
+ return update_exponent(result);
+#endif // LIBC_MATH_CBRT_SKIP_ACCURATE_PASS
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LIBC_SRC___SUPPORT_MATH_CBRT_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 9df9973810c77..a86619576cfc6 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4753,15 +4753,7 @@ add_entrypoint_object(
HDRS
../cbrt.h
DEPENDS
- libc.hdr.fenv_macros
- libc.src.__support.FPUtil.double_double
- libc.src.__support.FPUtil.dyadic_float
- 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.macros.optimization
- libc.src.__support.integer_literals
+ libc.src.__support.math.cbrt
)
add_entrypoint_object(
diff --git a/libc/src/math/generic/cbrt.cpp b/libc/src/math/generic/cbrt.cpp
index ce227e6650c84..e9b69bbf35cf6 100644
--- a/libc/src/math/generic/cbrt.cpp
+++ b/libc/src/math/generic/cbrt.cpp
@@ -7,334 +7,10 @@
//===----------------------------------------------------------------------===//
#include "src/math/cbrt.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/double_double.h"
-#include "src/__support/FPUtil/dyadic_float.h"
-#include "src/__support/FPUtil/multiply_add.h"
-#include "src/__support/common.h"
-#include "src/__support/integer_literals.h"
-#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-
-#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0)
-#define LIBC_MATH_CBRT_SKIP_ACCURATE_PASS
-#endif
+#include "src/__support/math/cbrt.h"
namespace LIBC_NAMESPACE_DECL {
-using DoubleDouble = fputil::DoubleDouble;
-using Float128 = fputil::DyadicFloat<128>;
-
-namespace {
-
-// Initial approximation of x^(-2/3) for 1 <= x < 2.
-// Polynomial generated by Sollya with:
-// > P = fpminimax(x^(-2/3), 7, [|D...|], [1, 2]);
-// > dirtyinfnorm(P/x^(-2/3) - 1, [1, 2]);
-// 0x1.28...p-21
-double intial_approximation(double x) {
- constexpr double COEFFS[8] = {
- 0x1.bc52aedead5c6p1, -0x1.b52bfebf110b3p2, 0x1.1d8d71d53d126p3,
- -0x1.de2db9e81cf87p2, 0x1.0154ca06153bdp2, -0x1.5973c66ee6da7p0,
- 0x1.07bf6ac832552p-2, -0x1.5e53d9ce41cb8p-6,
- };
-
- double x_sq = x * x;
-
- double c0 = fputil::multiply_add(x, COEFFS[1], COEFFS[0]);
- double c1 = fputil::multiply_add(x, COEFFS[3], COEFFS[2]);
- double c2 = fputil::multiply_add(x, COEFFS[5], COEFFS[4]);
- double c3 = fputil::multiply_add(x, COEFFS[7], COEFFS[6]);
-
- double x_4 = x_sq * x_sq;
- double d0 = fputil::multiply_add(x_sq, c1, c0);
- double d1 = fputil::multiply_add(x_sq, c3, c2);
-
- return fputil::multiply_add(x_4, d1, d0);
-}
-
-// Get the error term for Newton iteration:
-// h(x) = x^3 * a^2 - 1,
-#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-double get_error(const DoubleDouble &x_3, const DoubleDouble &a_sq) {
- return fputil::multiply_add(x_3.hi, a_sq.hi, -1.0) +
- fputil::multiply_add(x_3.lo, a_sq.hi, x_3.hi * a_sq.lo);
-}
-#else
-double get_error(const DoubleDouble &x_3, const DoubleDouble &a_sq) {
- DoubleDouble x_3_a_sq = fputil::quick_mult(a_sq, x_3);
- return (x_3_a_sq.hi - 1.0) + x_3_a_sq.lo;
-}
-#endif
-
-} // anonymous namespace
-
-// Correctly rounded cbrt algorithm:
-//
-// === Step 1 - Range reduction ===
-// For x = (-1)^s * 2^e * (1.m), we get 2 reduced arguments x_r and a as:
-// x_r = 1.m
-// a = (-1)^s * 2^(e % 3) * (1.m)
-// Then cbrt(x) = x^(1/3) can be computed as:
-// x^(1/3) = 2^(e / 3) * a^(1/3).
-//
-// In order to avoid division, we compute a^(-2/3) using Newton method and then
-// multiply the results by a:
-// a^(1/3) = a * a^(-2/3).
-//
-// === Step 2 - First approximation to a^(-2/3) ===
-// First, we use a degree-7 minimax polynomial generated by Sollya to
-// approximate x_r^(-2/3) for 1 <= x_r < 2.
-// p = P(x_r) ~ x_r^(-2/3),
-// with relative errors bounded by:
-// | p / x_r^(-2/3) - 1 | < 1.16 * 2^-21.
-//
-// Then we multiply with 2^(e % 3) from a small lookup table to get:
-// x_0 = 2^(-2*(e % 3)/3) * p
-// ~ 2^(-2*(e % 3)/3) * x_r^(-2/3)
-// = a^(-2/3)
-// With relative errors:
-// | x_0 / a^(-2/3) - 1 | < 1.16 * 2^-21.
-// This step is done in double precision.
-//
-// === Step 3 - First Newton iteration ===
-// We follow the method described in:
-// Sibidanov, A. and Zimmermann, P., "Correctly rounded cubic root evaluation
-// in double precision", https://core-math.gitlabpages.inria.fr/cbrt64.pdf
-// to derive multiplicative Newton iterations as below:
-// Let x_n be the nth approximation to a^(-2/3). Define the n^th error as:
-// h_n = x_n^3 * a^2 - 1
-// Then:
-// a^(-2/3) = x_n / (1 + h_n)^(1/3)
-// = x_n * (1 - (1/3) * h_n + (2/9) * h_n^2 - (14/81) * h_n^3 + ...)
-// using the Taylor series expansion of (1 + h_n)^(-1/3).
-//
-// Apply to x_0 above:
-// h_0 = x_0^3 * a^2 - 1
-// = a^2 * (x_0 - a^(-2/3)) * (x_0^2 + x_0 * a^(-2/3) + a^(-4/3)),
-// it...
[truncated]
|
bassiounix
requested review from
rupprecht,
keith and
aaronmondal
as code owners
bassiounix
force-pushed
the
users/bassiounix/spr/08-02-_libc_math_refactor_atanhf16_implementation_to_header-only_in_src___support_math_folder
branch
2 times, most recently
from
54e10e5 to
3265e87
Compare
Base automatically changed from
users/bassiounix/spr/08-02-_libc_math_refactor_atanhf16_implementation_to_header-only_in_src___support_math_folder
to
main
August 2, 2025 20:04
bassiounix
force-pushed
the
users/bassiounix/spr/08-02-_libc_math_refactor_cbrt_implementation_to_header-only_in_src___support_math_folder
branch
2 times, most recently
from
1de1830 to
68d35e4
Compare
bassiounix
force-pushed
the
users/bassiounix/spr/08-02-_libc_math_refactor_cbrt_implementation_to_header-only_in_src___support_math_folder
branch
2 times, most recently
from
e2b19e3 to
150b50e
Compare
lntue
approved these changes
Aug 3, 2025
Merge activity
|
…port/math folder.
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
bassiounix
deleted the
users/bassiounix/spr/08-02-_libc_math_refactor_cbrt_implementation_to_header-only_in_src___support_math_folder
branch
|
LLVM Buildbot has detected a new failure on builder Full details are available at: https://lab.llvm.org/buildbot/#/builders/10/builds/10756 Here is the relevant piece of the build log for the reference |
|
LLVM Buildbot has detected a new failure on builder Full details are available at: https://lab.llvm.org/buildbot/#/builders/11/builds/21004 Here is the relevant piece of the build log for the reference |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
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