[libc][math] Implement double precision cbrt correctly rounded to all rounding modes. #99262
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
lntue
requested review from
SchrodingerZhu,
jhuber6,
michaelrj-google and
overmighty
|
@llvm/pr-subscribers-libc Author: None (lntue) ChangesDivision-less Newton iterations algorithm for cube roots.
For Then In order to avoid division, we compute
First, we use a degree-7 minimax polynomial generated by Sollya to with relative errors bounded by: Then we multiply with with relative errors: This step is done in double precision.
We follow the method described in: Then: using the Taylor series expansion of Apply to it's bounded by: So in the first iteration step, we use: Its relative error is bounded by: Then we perform Ziv's rounding test and check if the answer is exact.
If the Ziv's rounding test from the previous step fails, we define the error And perform another iteration: with the relative errors exceed the precision of double-double.
If the Ziv's accuracy test from the previous step fails, we perform another Patch is 24.49 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/99262.diff 15 Files Affected:
diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt
index 383118dc781e5..32a08f20b328f 100644
--- a/libc/config/darwin/arm/entrypoints.txt
+++ b/libc/config/darwin/arm/entrypoints.txt
@@ -123,6 +123,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
+ libc.src.math.cbrt
libc.src.math.cbrtf
libc.src.math.copysign
libc.src.math.copysignf
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index dee6ac673643e..9b718c3f81151 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -345,6 +345,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
+ libc.src.math.cbrt
libc.src.math.cbrtf
libc.src.math.ceil
libc.src.math.ceilf
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index b0ee0e989b5ed..a72f8668808a5 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -216,6 +216,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
+ libc.src.math.cbrt
libc.src.math.cbrtf
libc.src.math.ceil
libc.src.math.ceilf
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 516a4b6ce3433..266c94d54a9df 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -347,6 +347,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
+ libc.src.math.cbrt
libc.src.math.cbrtf
libc.src.math.ceil
libc.src.math.ceilf
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index b6c55e7aa3033..4d19a28f4a2b3 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -370,6 +370,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.canonicalize
libc.src.math.canonicalizef
libc.src.math.canonicalizel
+ libc.src.math.cbrt
libc.src.math.cbrtf
libc.src.math.ceil
libc.src.math.ceilf
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 499c6bfe3a229..afc9ca87ff094 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -121,6 +121,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
+ libc.src.math.cbrt
libc.src.math.cbrtf
libc.src.math.copysign
libc.src.math.copysignf
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index aa56152aee141..a4c6b40b98388 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -382,6 +382,7 @@ def StdC : StandardSpec<"stdc"> {
],
[], // Enumerations
[
+ FunctionSpec<"cbrt", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"cbrtf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"copysign", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 6462afbc54a4f..dc2339896f2bb 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -65,6 +65,7 @@ add_math_entrypoint_object(canonicalizel)
add_math_entrypoint_object(canonicalizef16)
add_math_entrypoint_object(canonicalizef128)
+add_math_entrypoint_object(cbrt)
add_math_entrypoint_object(cbrtf)
add_math_entrypoint_object(ceil)
diff --git a/libc/src/math/cbrt.h b/libc/src/math/cbrt.h
new file mode 100644
index 0000000000000..a7d5fe80e57b3
--- /dev/null
+++ b/libc/src/math/cbrt.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for cbrt --------------------------*- 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_MATH_CBRT_H
+#define LLVM_LIBC_SRC_MATH_CBRT_H
+
+namespace LIBC_NAMESPACE {
+
+double cbrt(double x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_CBRT_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index c2f58fb1a4f71..318728d6e315c 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4180,3 +4180,19 @@ add_entrypoint_object(
libc.src.__support.FPUtil.multiply_add
libc.src.__support.macros.optimization
)
+
+add_entrypoint_object(
+ cbrt
+ SRCS
+ cbrt.cpp
+ HDRS
+ ../cbrt.h
+ COMPILE_OPTIONS
+ -O3
+ 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
+)
diff --git a/libc/src/math/generic/cbrt.cpp b/libc/src/math/generic/cbrt.cpp
new file mode 100644
index 0000000000000..a60e2ea44b6e9
--- /dev/null
+++ b/libc/src/math/generic/cbrt.cpp
@@ -0,0 +1,340 @@
+//===-- Implementation of cbrt function -----------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#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
+
+namespace LIBC_NAMESPACE_DECL {
+
+using DoubleDouble = fputil::DoubleDouble;
+using Float128 = typename 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
+constexpr 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
+constexpr 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
+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
+
+} // 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'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
+
+LLVM_LIBC_FUNCTION(double, cbrt, (double x)) {
+ using FPBits = typename fputil::FPBits<double>;
+
+ uint64_t x_u = FPBits(x).uintval();
+ uint64_t x_abs = x_u & 0x7fff'ffff'ffff'ffff;
+
+ unsigned exp_bias_correction = 682; // 1023 * 2/3
+
+ if (LIBC_UNLIKELY(x_abs < FPBits::min_normal().uintval() ||
+ x_abs >= FPBits::inf().uintval())) {
+ if (x_abs == 0 || x_abs >= FPBits::inf().uintval())
+ // x is 0, Inf, or NaN.
+ return 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() & 0x800F'FFFF'FFFF'FFFF;
+ // 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 LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 0dc7ae6aae2df..64b4d2c58fb6a 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2225,6 +2225,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ cbrt_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ cbrt_test.cpp
+ DEPENDS
+ libc.src.math.cbrt
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_subdirectory(generic)
add_subdirectory(smoke)
diff --git a/libc/test/src/math/cbrt_test.cpp b/libc/test/src/math/cbrt_test.cpp
new file mode 100644
index 0000000000000..123351496118b
--- /dev/null
+++ b/libc/test/src/math/cbrt_test.cpp
@@ -0,0 +1,104 @@
+//===-- Unittests for cbrt ------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "hdr/math_macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/cbrt.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcCbrtTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LIBC_NAMESPACE::testing::tlog;
+
+TEST_F(LlvmLibcCbrtTest, InDoubleRange) {
+ constexpr uint64_t COUNT = 123'451;
+ uint64_t START = LIBC_NAMESPACE::fputil::FPBits<double>(1.0).uintval();
+ uint64_t STOP = LIBC_NAMESPACE::fputil::FPBits<double>(8.0).uintval();
+ uint64_t STEP = (STOP - START) / COUNT;
+
+ auto test = [&](mpfr::RoundingMode rounding_mode) {
+ mpfr::ForceRoundingMode force_rounding(rounding_mode);
+ if (!force_rounding.success)
+ return;
+
+ uint64_t fails = 0;
+ uint64_t tested = 0;
+ uint64_t total = 0;
+ double worst_input, worst_output = 0.0;
+ double ulp = 0.5;
+
+ for (uint64_t i = 0, v = START; i <= COUNT; ++i, v += STEP) {
+ double x = FPBits(v).get_val();
+ if (isnan(x) || isinf(x))
+ continue;
+
+ double result = LIB...
[truncated]
|
|
all tests are ok on my side |
overmighty
reviewed
Jul 17, 2024
Comment on lines
+39
to
+43
| constexpr double COEFFS[8] = { | ||
| 0x1.bc52aedead5c6p1, -0x1.b52bfebf110b3p2, 0x1.1d8d71d53d126p3, | ||
| -0x1.de2db9e81cf87p2, 0x1.0154ca06153bdp2, -0x1.5973c66ee6da7p0, | ||
| 0x1.07bf6ac832552p-2, -0x1.5e53d9ce41cb8p-6, | ||
| }; |
There was a problem hiding this comment.
Do we have a policy on when to use C-style arrays, and when to use cpp::array?
There was a problem hiding this comment.
I don't think we have one.
overmighty
approved these changes
Jul 17, 2024
SchrodingerZhu
approved these changes
Jul 17, 2024
|
LLVM Buildbot has detected a new failure on builder Full details are available at: https://lab.llvm.org/buildbot/#/builders/104/builds/2394 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/71/builds/2383 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/131/builds/2342 Here is the relevant piece of the build log for the reference: |
Closed
|
LLVM Buildbot has detected a new failure on builder Full details are available at: https://lab.llvm.org/buildbot/#/builders/182/builds/1172 Here is the relevant piece of the build log for the reference: |
petrhosek
added a commit
to petrhosek/llvm-project
that referenced
this pull request
Jul 22, 2024
This is a follow up to llvm#99262.
petrhosek
added a commit
that referenced
this pull request
Jul 22, 2024
This is a follow up to #99262.
yuxuanchen1997
pushed a commit
that referenced
this pull request
Jul 25, 2024
… rounding modes. (#99262) Division-less Newton iterations algorithm for cube roots. 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). ``` 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. 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. 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. 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.
yuxuanchen1997
pushed a commit
that referenced
this pull request
Jul 25, 2024
This is a follow up to #99262.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
6 participants
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.
Division-less Newton iterations algorithm for cube roots.
For
x = (-1)^s * 2^e * (1.m), we get 2 reduced argumentsx_randaas:Then
cbrt(x) = x^(1/3)can be computed as:In order to avoid division, we compute
a^(-2/3)using Newton method and thenmultiply the results by a:
First, we use a degree-7 minimax polynomial generated by Sollya to
approximate
x_r^(-2/3)for1 <= x_r < 2.with relative errors bounded by:
Then we multiply with
2^(e % 3)from a small lookup table to get:with relative errors:
This step is done in double precision.
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_nbe the nth approximation toa^(-2/3). Define the n^th error as:Then:
using the Taylor series expansion of
(1 + h_n)^(-1/3).Apply to
x_0above:it's bounded by:
So in the first iteration step, we use:
Its relative error is bounded by:
Then we perform Ziv's rounding test and check if the answer is exact.
This step is done in double-double precision.
If the Ziv's rounding test from the previous step fails, we define the error
term:
And perform another iteration:
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.
If the Ziv's accuracy test from the previous step fails, we perform another
iteration in 128-bit precision and check for exact outputs.