[libc][math] Implement double precision cbrt correctly rounded to all rounding modes. (#99262)

lntue · yuxuanchen1997 · commit b878d8ae17e1 · 2024-07-25T12:51:29.000-07:00
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.
diff --git 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/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
@@ -245,6 +245,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.atanf
     libc.src.math.atanh
     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/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
@@ -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
@@ -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
@@ -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
@@ -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/docs/math/index.rst b/libc/docs/math/index.rst
@@ -266,7 +266,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | atanpi    |                  |                 |                        |                      |                        | 7.12.4.10              | F.10.1.10                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| cbrt      | |check|          |                 |                        |                      |                        | 7.12.7.1               | F.10.4.1                   |
+| cbrt      | |check|          | |check|         |                        |                      |                        | 7.12.7.1               | F.10.4.1                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | compoundn |                  |                 |                        |                      |                        | 7.12.7.2               | F.10.4.2                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git 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
@@ -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
@@ -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
@@ -4180,3 +4180,23 @@ 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.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
+)
diff --git a/libc/src/math/generic/cbrt.cpp b/libc/src/math/generic/cbrt.cpp
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
diff --git a/libc/test/src/math/cbrt_test.cpp b/libc/test/src/math/cbrt_test.cpp
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
diff --git a/libc/test/src/math/smoke/cbrt_test.cpp b/libc/test/src/math/smoke/cbrt_test.cpp

Original file line number	Diff line number	Diff line change
`@@ -382,6 +382,7 @@ def StdC : StandardSpec<"stdc"> {`
`382`	`382`	`],`
`383`	`383`	`[], // Enumerations`
`384`	`384`	`[`
	`385`	`+ FunctionSpec<"cbrt", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,`
`385`	`386`	`FunctionSpec<"cbrtf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,`
`386`	`387`
`387`	`388`	`FunctionSpec<"copysign", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,`