[libc-commits] [libc] [libc][math] Implement atan2f correctly rounded to all rounding modes. (PR #86716)
via libc-commits
libc-commits at lists.llvm.org
Fri Mar 29 08:46:29 PDT 2024
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/86716
>From 93d41b7e2b54f37e042904386ba31f025188f795 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue at google.com>
Date: Tue, 26 Mar 2024 15:08:31 -0400
Subject: [PATCH 1/5] [libc][math] Implement atan2f correctly rounded to all
rounding modes.
---
libc/config/baremetal/arm/entrypoints.txt | 1 +
libc/config/baremetal/riscv/entrypoints.txt | 1 +
libc/config/darwin/arm/entrypoints.txt | 1 +
libc/config/linux/aarch64/entrypoints.txt | 1 +
libc/config/linux/arm/entrypoints.txt | 1 +
libc/config/linux/riscv/entrypoints.txt | 1 +
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/config/windows/entrypoints.txt | 1 +
libc/docs/math/index.rst | 3 +-
libc/spec/stdc.td | 4 +
libc/src/math/generic/CMakeLists.txt | 19 ++
libc/src/math/generic/atan2f.cpp | 297 ++++++++++++++++++++
libc/src/math/generic/atanf.cpp | 12 +-
libc/src/math/generic/inv_trigf_utils.cpp | 107 +++----
libc/src/math/generic/inv_trigf_utils.h | 5 +-
libc/test/src/math/CMakeLists.txt | 13 +
libc/test/src/math/atan2f_test.cpp | 132 +++++++++
libc/test/src/math/atanf_test.cpp | 5 +-
libc/test/src/math/smoke/CMakeLists.txt | 12 +
libc/test/src/math/smoke/atan2f_test.cpp | 50 ++++
libc/utils/MPFRWrapper/MPFRUtils.cpp | 8 +
libc/utils/MPFRWrapper/MPFRUtils.h | 1 +
22 files changed, 611 insertions(+), 65 deletions(-)
create mode 100644 libc/src/math/generic/atan2f.cpp
create mode 100644 libc/test/src/math/atan2f_test.cpp
create mode 100644 libc/test/src/math/smoke/atan2f_test.cpp
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index 17ce56e228a6ac..9e21f5c20d9207 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -207,6 +207,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.ceil
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index 39756e1ee29f54..7664937da0f6e0 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -207,6 +207,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.ceil
diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt
index 02a09256606956..6b89ce55d72b65 100644
--- a/libc/config/darwin/arm/entrypoints.txt
+++ b/libc/config/darwin/arm/entrypoints.txt
@@ -118,6 +118,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 78da7f0b334b1f..4ba2d8387b07eb 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -330,6 +330,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index 6e63e270280e7a..04baa4c1cf93ab 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -198,6 +198,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.ceil
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 5aae4e246cfb3c..25745513b920ba 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -338,6 +338,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5b428e51aee620..56ac08ab23c609 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -348,6 +348,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.canonicalize
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index f4456f561ec017..c38125a6462272 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -116,6 +116,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.atan2f
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index 080b6a4427f511..b7f1b8739648ca 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -422,7 +422,7 @@ Higher Math Functions
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
| atan2 | | | | | | | | | | | | |
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
-| atan2f | | | | | | | | | | | | |
+| atan2f | |check| | |check| | |check| | |check| | |check| | | | |check| | |check| | |check| | | |
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
| atan2l | | | | | | | | | | | | |
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
@@ -591,6 +591,7 @@ acosh |check|
asin |check|
asinh |check|
atan |check|
+atan2 |check|
atanh |check|
cos |check| large
cosh |check|
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index ac6e1d1801ba55..719bb9aa18cb0a 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -620,10 +620,14 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"acosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+
FunctionSpec<"asinf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"asin", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
+
FunctionSpec<"atanf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+ FunctionSpec<"atan2f", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
+
FunctionSpec<"acoshf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"asinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"atanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 4d9b91150d0200..7e9c9ce7bb2748 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2854,6 +2854,25 @@ add_entrypoint_object(
-O3
)
+add_entrypoint_object(
+ atan2f
+ SRCS
+ atan2f.cpp
+ HDRS
+ ../atan2f.h
+ COMPILE_OPTIONS
+ -O3
+ DEPENDS
+ .inv_trigf_utils
+ libc.src.__support.FPUtil.except_value_utils
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.rounding_mode
+ libc.src.__support.macros.optimization
+)
+
add_entrypoint_object(
scalbn
SRCS
diff --git a/libc/src/math/generic/atan2f.cpp b/libc/src/math/generic/atan2f.cpp
new file mode 100644
index 00000000000000..089e6c61984bea
--- /dev/null
+++ b/libc/src/math/generic/atan2f.cpp
@@ -0,0 +1,297 @@
+//===-- Single-precision atan2f 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/atan2f.h"
+#include "inv_trigf_utils.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE {
+
+namespace {
+
+// Look up tables for accurate pass:
+
+// atan(i/16) with i = 0..16, generated by Sollya with:
+// > for i from 0 to 16 do {
+// a = round(atan(i/16), D, RN);
+// b = round(atan(i/16) - a, D, RN);
+// print("{", b, ",", a, "},");
+// };
+constexpr fputil::DoubleDouble ATAN_I[17] = {
+ {0.0, 0.0},
+ {-0x1.c934d86d23f1dp-60, 0x1.ff55bb72cfdeap-5},
+ {-0x1.cd37686760c17p-59, 0x1.fd5ba9aac2f6ep-4},
+ {0x1.347b0b4f881cap-58, 0x1.7b97b4bce5b02p-3},
+ {0x1.8ab6e3cf7afbdp-57, 0x1.f5b75f92c80ddp-3},
+ {-0x1.963a544b672d8p-57, 0x1.362773707ebccp-2},
+ {-0x1.c63aae6f6e918p-56, 0x1.6f61941e4def1p-2},
+ {-0x1.24dec1b50b7ffp-56, 0x1.a64eec3cc23fdp-2},
+ {0x1.a2b7f222f65e2p-56, 0x1.dac670561bb4fp-2},
+ {-0x1.d5b495f6349e6p-56, 0x1.0657e94db30dp-1},
+ {-0x1.928df287a668fp-58, 0x1.1e00babdefeb4p-1},
+ {0x1.1021137c71102p-55, 0x1.345f01cce37bbp-1},
+ {0x1.2419a87f2a458p-56, 0x1.4978fa3269ee1p-1},
+ {0x1.0028e4bc5e7cap-57, 0x1.5d58987169b18p-1},
+ {-0x1.8c34d25aadef6p-56, 0x1.700a7c5784634p-1},
+ {-0x1.bf76229d3b917p-56, 0x1.819d0b7158a4dp-1},
+ {0x1.1a62633145c07p-55, 0x1.921fb54442d18p-1},
+};
+
+// Taylor polynomial, generated by Sollya with:
+// > for i from 0 to 8 do {
+// j = (-1)^(i + 1)/(2*i + 1);
+// a = round(j, D, RN);
+// b = round(j - a, D, RN);
+// print("{", b, ",", a, "},");
+// };
+constexpr fputil::DoubleDouble COEFFS[9] = {
+ {0.0, 1.0}, // 1
+ {-0x1.5555555555555p-56, -0x1.5555555555555p-2}, // -1/3
+ {-0x1.999999999999ap-57, 0x1.999999999999ap-3}, // 1/5
+ {-0x1.2492492492492p-57, -0x1.2492492492492p-3}, // -1/7
+ {0x1.c71c71c71c71cp-58, 0x1.c71c71c71c71cp-4}, // 1/9
+ {0x1.745d1745d1746p-59, -0x1.745d1745d1746p-4}, // -1/11
+ {-0x1.3b13b13b13b14p-58, 0x1.3b13b13b13b14p-4}, // 1/13
+ {-0x1.1111111111111p-60, -0x1.1111111111111p-4}, // -1/15
+ {0x1.e1e1e1e1e1e1ep-61, 0x1.e1e1e1e1e1e1ep-5}, // 1/17
+};
+
+// Veltkamp's splitting of a double precision into hi + lo, where the hi part is
+// slightly smaller than an even split, so that the product of
+// hi * s * k is exact,
+// where:
+// s is single precsion,
+// 0 < k < 16 is an integer.
+// This is used when FMA instruction is not available.
+[[maybe_unused]] LIBC_INLINE constexpr fputil::DoubleDouble split_d(double a) {
+ fputil::DoubleDouble r{0.0, 0.0};
+ constexpr double C = 0x1.0p33 + 1.0;
+ double t1 = C * a;
+ double t2 = a - t1;
+ r.hi = t1 + t2;
+ r.lo = a - r.hi;
+ return r;
+}
+
+} // anonymous namespace
+
+// There are several range reduction steps we can take for atan2(y, x) as
+// follow:
+
+// * Range reduction 1: signness
+// atan2(y, x) will return a number between -PI and PI representing the angle
+// forming by the 0x axis and the vector (x, y) on the 0xy-plane.
+// In particular, we have that:
+// atan2(y, x) = atan( y/x ) if x >= 0 and y >= 0 (I-quadrant)
+// = pi + atan( y/x ) if x < 0 and y >= 0 (II-quadrant)
+// = -pi + atan( y/x ) if x < 0 and y < 0 (III-quadrant)
+// = atan( y/x ) if x >= 0 and y < 0 (IV-quadrant)
+// Since atan function is odd, we can use the formula:
+// atan(-u) = -atan(u)
+// to adjust the above conditions a bit further:
+// atan2(y, x) = atan( |y|/|x| ) if x >= 0 and y >= 0 (I-quadrant)
+// = pi - atan( |y|/|x| ) if x < 0 and y >= 0 (II-quadrant)
+// = -pi + atan( |y|/|x| ) if x < 0 and y < 0 (III-quadrant)
+// = -atan( |y|/|x| ) if x >= 0 and y < 0 (IV-quadrant)
+// Which can be simplified to:
+// atan2(y, x) = sign(y) * atan( |y|/|x| ) if x >= 0
+// = sign(y) * (pi - atan( |y|/|x| )) if x < 0
+
+// * Range reduction 2: reciprocal
+// Now that the argument inside atan is positive, we can use the formula:
+// atan(1/x) = pi/2 - atan(x)
+// to make the argument inside atan <= 1 as follow:
+// atan2(y, x) = sign(y) * atan( |y|/|x|) if 0 <= |y| <= x
+// = sign(y) * (pi/2 - atan( |x|/|y| ) if 0 <= x < |y|
+// = sign(y) * (pi - atan( |y|/|x| )) if 0 <= |y| <= -x
+// = sign(y) * (pi/2 + atan( |x|/|y| )) if 0 <= -x < |y|
+
+// * Range reduction 3: look up table.
+// After the previous two range reduction steps, we reduce the problem to
+// compute atan(u) with 0 <= u <= 1, or to be precise:
+// atan( n / d ) where n = min(|x|, |y|) and d = max(|x|, |y|).
+// An accurate polynomial approximation for the whole [0, 1] input range will
+// require a very large degree. To make it more efficient, we reduce the input
+// range further by finding an integer idx such that:
+// | n/d - idx/16 | <= 1/32.
+// In particular,
+// idx := 2^-4 * round(2^4 * n/d)
+// Then for the fast pass, we find a polynomial approximation for:
+// atan( n/d ) ~ atan( idx/16 ) + (n/d - idx/16) * Q(n/d - idx/16)
+// For the accurate pass, we use the addition formula:
+// atan( n/d ) - atan( idx/16 ) = atan( (n/d - idx/16)/(1 + (n*idx)/(16*d)) )
+// = atan( (n - d * idx/16)/(d + n * idx/16) )
+// And finally we use Taylor polynomial to compute the RHS in the accurate pass:
+// atan(u) ~ P(u) = u - u^3/3 + u^5/5 - u^7/7 + u^9/9 - u^11/11 + u^13/13 -
+// - u^15/15 + u^17/17
+// It's error in double-double precision is estimated in Sollya to be:
+// > P = x - x^3/3 + x^5/5 -x^7/7 + x^9/9 - x^11/11 + x^13/13 - x^15/15
+// + x^17/17;
+// > dirtyinfnorm(atan(x) - P, [-2^-5, 2^-5]);
+// 0x1.aec6f...p-100
+// which is about rounding errors of double-double (2^-104).
+
+LLVM_LIBC_FUNCTION(float, atan2f, (float y, float x)) {
+ using FPBits = typename fputil::FPBits<float>;
+ constexpr double IS_NEG[2] = {1.0, -1.0};
+ constexpr double PI = 0x1.921fb54442d18p1;
+ constexpr double PI_LO = 0x1.1a62633145c07p-53;
+ constexpr double PI_OVER_4 = 0x1.921fb54442d18p-1;
+ constexpr double PI_OVER_2 = 0x1.921fb54442d18p0;
+ constexpr double THREE_PI_OVER_4 = 0x1.2d97c7f3321d2p+1;
+ // Adjustment for constant term:
+ // CONST_ADJ[x_sign][y_sign][recip]
+ constexpr fputil::DoubleDouble CONST_ADJ[2][2][2] = {
+ {{{0.0, 0.0}, {-PI_LO / 2, -PI_OVER_2}},
+ {{-0.0, -0.0}, {-PI_LO / 2, -PI_OVER_2}}},
+ {{{-PI_LO, -PI}, {PI_LO / 2, PI_OVER_2}},
+ {{-PI_LO, -PI}, {PI_LO / 2, PI_OVER_2}}}};
+
+ FPBits x_bits(x), y_bits(y);
+ bool x_sign = x_bits.sign().is_neg();
+ bool y_sign = y_bits.sign().is_neg();
+ x_bits.set_sign(Sign::POS);
+ y_bits.set_sign(Sign::POS);
+ uint32_t x_abs = x_bits.uintval();
+ uint32_t y_abs = y_bits.uintval();
+ uint32_t max_abs = x_abs > y_abs ? x_abs : y_abs;
+ uint32_t min_abs = x_abs <= y_abs ? x_abs : y_abs;
+ bool recip = x_abs < y_abs;
+
+ if (LIBC_UNLIKELY(max_abs >= 0x7f80'0000U || min_abs == 0U)) {
+ if (x_bits.is_nan() || y_bits.is_nan())
+ return FPBits::quiet_nan().get_val();
+ size_t x_except = x_abs == 0 ? 0 : (x_abs == 0x7f80'0000 ? 2 : 1);
+ size_t y_except = y_abs == 0 ? 0 : (y_abs == 0x7f80'0000 ? 2 : 1);
+
+ // Exceptional cases:
+ // EXCEPT[y_except][x_except][x_is_neg]
+ // with x_except & y_except:
+ // 0: zero
+ // 1: finite, non-zero
+ // 2: infinity
+ constexpr double EXCEPTS[3][3][2] = {
+ {{0.0, PI}, {0.0, PI}, {0.0, PI}},
+ {{PI_OVER_2, PI_OVER_2}, {0.0, 0.0}, {0.0, PI}},
+ {{PI_OVER_2, PI_OVER_2},
+ {PI_OVER_2, PI_OVER_2},
+ {PI_OVER_4, THREE_PI_OVER_4}},
+ };
+
+ double r = IS_NEG[y_sign] * EXCEPTS[y_except][x_except][x_sign];
+
+ return static_cast<float>(r);
+ }
+
+ double final_sign = IS_NEG[(x_sign != y_sign) != recip];
+ fputil::DoubleDouble const_term = CONST_ADJ[x_sign][y_sign][recip];
+ double num_d = static_cast<double>(FPBits(min_abs).get_val());
+ double den_d = static_cast<double>(FPBits(max_abs).get_val());
+ double q_d = num_d / den_d;
+ int idx;
+
+ double k_d = fputil::nearest_integer(q_d * 0x1.0p4f);
+ idx = static_cast<int>(k_d);
+ q_d = fputil::multiply_add(k_d, -0x1.0p-4, q_d);
+
+ double p, r;
+
+ if (LIBC_UNLIKELY(idx == 0)) {
+ // For |x| < 1/16, we use Taylor polynomial:
+ // atan(x) ~ x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - x^11/11 + x^13/13
+ double q2 = q_d * q_d;
+ double c0 = fputil::multiply_add(q2, ATAN_COEFFS[0][2], ATAN_COEFFS[0][1]);
+ double c1 = fputil::multiply_add(q2, ATAN_COEFFS[0][4], ATAN_COEFFS[0][3]);
+ double c2 = fputil::multiply_add(q2, ATAN_COEFFS[0][6], ATAN_COEFFS[0][5]);
+ double q4 = q2 * q2;
+ double q3 = q_d * q2;
+ double d0 = fputil::polyeval(q4, c0, c1, c2);
+ p = fputil::multiply_add(q3, d0, q_d);
+ r = final_sign * (p + const_term.hi);
+ } else {
+ p = atan_eval(q_d, idx);
+ r = final_sign *
+ fputil::multiply_add(q_d, p, const_term.hi + ATAN_COEFFS[idx][0]);
+ }
+
+ constexpr uint32_t LOWER_ERR = 4;
+ // Mask sticky bits in double precision before rounding to single precision.
+ constexpr uint32_t MASK =
+ mask_trailing_ones<uint32_t, fputil::FPBits<double>::SIG_LEN -
+ FPBits::SIG_LEN - 1>();
+ constexpr uint32_t UPPER_ERR = MASK - LOWER_ERR;
+
+ uint32_t r_bits = static_cast<uint32_t>(cpp::bit_cast<uint64_t>(r)) & MASK;
+
+ // Ziv's rounding test.
+ if (LIBC_LIKELY(r_bits > LOWER_ERR && r_bits < UPPER_ERR))
+ return static_cast<float>(r);
+
+ // Use double-double.
+ fputil::DoubleDouble q;
+ double num_r, den_r;
+
+ if (idx != 0) {
+ // The following range reduction is accurate even without fma for
+ // 1/16 <= n/d <= 1.
+ // atan(n/d) - atan(idx/16) = atan((n/d - idx/16) / (1 + (n/d) * (idx/16)))
+ // = atan((n - d*(idx/16)) / (d + n*idx/16))
+ k_d *= 0x1.0p-4;
+ num_r = fputil::multiply_add(k_d, -den_d, num_d); // Exact
+ den_r = fputil::multiply_add(k_d, num_d, den_d); // Exact
+ q.hi = num_r / den_r;
+ } else {
+ // For 0 < n/d < 1/16, we just need to calculate the lower part of their
+ // quotient.
+ q.hi = q_d;
+ num_r = num_d;
+ den_r = den_d;
+ }
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ q.lo = fputil::multiply_add(q.hi, -den_r, num_r) / den_r;
+#else
+ fputil::DoubleDouble q_hi_dd = split_d(q.hi);
+ double t1 = fputil::multiply_add(q_hi_dd.hi, -den_r, num_r); // Exact
+ double t2 = fputil::multiply_add(q_hi_dd.lo, -den_r, t1);
+ q.lo = t2 / den_r;
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+ fputil::DoubleDouble q2 = fputil::quick_mult(q, q);
+ fputil::DoubleDouble p_dd =
+ fputil::polyeval(q2, COEFFS[0], COEFFS[1], COEFFS[2], COEFFS[3],
+ COEFFS[4], COEFFS[5], COEFFS[6], COEFFS[7], COEFFS[8]);
+ fputil::DoubleDouble r_dd =
+ fputil::add(const_term, fputil::multiply_add(q, p_dd, ATAN_I[idx]));
+ r_dd.hi *= final_sign;
+ r_dd.lo *= final_sign;
+
+ // Make sure the sum is normalized:
+ fputil::DoubleDouble rr = fputil::exact_add(r_dd.hi, r_dd.lo);
+ // Round to odd.
+ uint64_t rr_bits = cpp::bit_cast<uint64_t>(rr.hi);
+ if (LIBC_UNLIKELY(((rr_bits & 0xfff'ffff) == 0) && (rr.lo != 0.0))) {
+ Sign hi_sign = fputil::FPBits<double>(rr.hi).sign();
+ Sign lo_sign = fputil::FPBits<double>(rr.lo).sign();
+ if (hi_sign == lo_sign) {
+ ++rr_bits;
+ } else if ((rr_bits & fputil::FPBits<double>::FRACTION_MASK) > 0) {
+ --rr_bits;
+ }
+ }
+
+ return static_cast<float>(cpp::bit_cast<double>(rr_bits));
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/math/generic/atanf.cpp b/libc/src/math/generic/atanf.cpp
index 4adda429cc041c..a58bbdd5033f87 100644
--- a/libc/src/math/generic/atanf.cpp
+++ b/libc/src/math/generic/atanf.cpp
@@ -81,11 +81,6 @@ LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
int idx;
if (x_abs > 0x3f80'0000U) {
- // Exceptional value:
- if (LIBC_UNLIKELY(x_abs == 0x3ffe'2ec1U)) { // |x| = 0x1.fc5d82p+0
- return sign.is_pos() ? fputil::round_result_slightly_up(0x1.1ab2fp0f)
- : fputil::round_result_slightly_down(-0x1.1ab2fp0f);
- }
// |x| > 1, we need to invert x, so we will perform range reduction in
// double precision.
x_d = 1.0 / static_cast<double>(x_bits.get_val());
@@ -98,10 +93,9 @@ LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
SIGNED_PI_OVER_2[sign.is_neg()]);
} else {
// Exceptional value:
- if (LIBC_UNLIKELY(x_abs == 0x3dbb'6ac7U)) { // |x| = 0x1.76d58ep-4
- return sign.is_pos()
- ? fputil::round_result_slightly_up(0x1.75cb06p-4f)
- : fputil::round_result_slightly_down(-0x1.75cb06p-4f);
+ if (LIBC_UNLIKELY(x_abs == 0x3d8d'6b23U)) { // |x| = 0x1.1ad646p-4
+ return sign.is_pos() ? fputil::round_result_slightly_down(0x1.1a6386p-4f)
+ : fputil::round_result_slightly_up(-0x1.1a6386p-4f);
}
// Perform range reduction in single precision.
float x_f = x_bits.get_val();
diff --git a/libc/src/math/generic/inv_trigf_utils.cpp b/libc/src/math/generic/inv_trigf_utils.cpp
index 93d5bcbf7b567d..06da93db96350b 100644
--- a/libc/src/math/generic/inv_trigf_utils.cpp
+++ b/libc/src/math/generic/inv_trigf_utils.cpp
@@ -16,65 +16,68 @@ namespace LIBC_NAMESPACE {
// Generated by Sollya with:
// > for i from 1 to 16 do {
// mid_point = i/16;
-// P = fpminimax(atan(mid_point + x), 7, [|D...|], [-1/32, 1/32]);
+// P = fpminimax(atan(mid_point + x), 8, [|D...|], [-1/32, 1/32]);
// print("{", coeff(P, 0), ",", coeff(P, 1), ",", coeff(P, 2), ",",
// coeff(P, 3), ",", coeff(P, 4), ",", coeff(P, 5), ",", coeff(P, 6),
-// ",", coeff(P, 7), "},");
+// ",", coeff(P, 7), ",", coeff(P, 8), "},");
// };
// For i = 0, ATAN_COEFFS[0][j] = (-1)^j * (1/(2*j + 1)) is the odd coefficients
// of the Taylor polynomial of atan(x).
-double ATAN_COEFFS[17][8] = {
+// Notice that degree-7 is good enough for atanf, but degree-8 helps reduce the
+// error bounds for atan2f's fast pass 16 times, and it does not affect the
+// performance of atanf much.
+double ATAN_COEFFS[17][9] = {
{0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,
-0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4, -0x1.745d1745d1746p-4,
- 0x1.3b13b13b13b14p-4, -0x1.1111111111111p-4},
- {0x1.ff55bb72cfdb1p-5, 0x1.fe01fe01fe1bp-1, -0x1.fc05f80821d1ap-5,
- -0x1.4d6930419fc5fp-2, 0x1.f61b9f6d69313p-5, 0x1.8208a32f4346cp-3,
- -0x1.ecb8fc53d04efp-5, -0x1.060710cb59cbcp-3},
- {0x1.fd5ba9aac2f3cp-4, 0x1.f81f81f81f96ap-1, -0x1.f05e09cf4c1b2p-4,
- -0x1.368c3aac7543ep-2, 0x1.d9b14bddfac55p-4, 0x1.4048e55ec725ep-3,
- -0x1.b98ca3c1594b5p-4, -0x1.664eabaabbc16p-4},
- {0x1.7b97b4bce5ae7p-3, 0x1.ee9c7f8458f06p-1, -0x1.665c226c8dc69p-3,
- -0x1.1344bb77961b7p-2, 0x1.42ac97745d3ccp-3, 0x1.c32e142047ec1p-4,
- -0x1.137ae41ab96cbp-3, -0x1.1a6ae8c09a4b6p-5},
- {0x1.f5b75f92c80c6p-3, 0x1.e1e1e1e1e1ed4p-1, -0x1.c5894d101ad4p-3,
- -0x1.ce6de02b38c38p-3, 0x1.78a3920c336b9p-3, 0x1.dd5ff94a9d499p-5,
- -0x1.1ac2d3f9d072ep-3, 0x1.0af9735dff373p-6},
- {0x1.362773707ebc5p-2, 0x1.d272ca3fc5b8bp-1, -0x1.0997e8ae90cb6p-2,
- -0x1.6cf6667146798p-3, 0x1.8dd1dff17f3d3p-3, 0x1.24860eced656fp-7,
- -0x1.f4220e8f18ed5p-4, 0x1.b700aed7cdc34p-5},
- {0x1.6f61941e4deeep-2, 0x1.c0e070381c115p-1, -0x1.2726dd1347c7ep-2,
- -0x1.09f37b3ad010dp-3, 0x1.85eaca5196f5cp-3, -0x1.04d640117852ap-5,
- -0x1.802c2956871c7p-4, 0x1.2992b45df0ee7p-4},
- {0x1.a64eec3cc23fep-2, 0x1.adbe87f94906bp-1, -0x1.3b9d8eab5eae5p-2,
- -0x1.57c09646faabbp-4, 0x1.6795330e73aep-3, -0x1.f2d89a702a652p-5,
- -0x1.f3afd90a9d4d7p-5, 0x1.3261723d3f153p-4},
- {0x1.dac670561bb53p-2, 0x1.999999999998fp-1, -0x1.47ae147afd8cap-2,
- -0x1.5d867c3dfd72ap-5, 0x1.3a92a76cba833p-3, -0x1.3ec460286928ap-4,
- -0x1.ed02ff86892acp-6, 0x1.0a674c8f05727p-4},
- {0x1.0657e94db30d2p-1, 0x1.84f00c27805ffp-1, -0x1.4c62cb564f677p-2,
- -0x1.e6495b262dfe7p-8, 0x1.063c34eca262bp-3, -0x1.58b78dc79b5aep-4,
- -0x1.4623815233be1p-8, 0x1.93afe94328089p-5},
- {0x1.1e00babdefeb6p-1, 0x1.702e05c0b8159p-1, -0x1.4af2b78236bd6p-2,
- 0x1.5d0b7ea46ed08p-6, 0x1.a124870236935p-4, -0x1.519e1ec133a88p-4,
- 0x1.a54632a3f48c7p-7, 0x1.099ca0945096dp-5},
- {0x1.345f01cce37bdp-1, 0x1.5babcc647fa7ep-1, -0x1.449db09443a67p-2,
- 0x1.655caac78a0fcp-5, 0x1.3bbbdb0d09efap-4, -0x1.34a306c27e021p-4,
- 0x1.83fe749c7966p-6, 0x1.2057cc96d9edcp-6},
- {0x1.4978fa3269ee2p-1, 0x1.47ae147ae146bp-1, -0x1.3a92a305652e1p-2,
- 0x1.ec21b5172657fp-5, 0x1.c2f8c45d2f4eep-5, -0x1.0ba99c4aeb8acp-4,
- 0x1.d716a4af4d1d6p-6, 0x1.97fba0a9696dep-8},
- {0x1.5d58987169b19p-1, 0x1.34679ace0133cp-1, -0x1.2ddfb03920e2fp-2,
- 0x1.2491307c0fa0bp-4, 0x1.29c7eca0136fp-5, -0x1.bca792caa6f1cp-5,
- 0x1.e5d92545576bcp-6, -0x1.8ca76fcf5ccd2p-10},
- {0x1.700a7c5784634p-1, 0x1.21fb78121fb71p-1, -0x1.1f6a8499ea541p-2,
- 0x1.41b15e5e77bcfp-4, 0x1.59bc9bf54fb02p-6, -0x1.63b54ff058e0fp-5,
- 0x1.c8da01221306fp-6, -0x1.906b2c274c39cp-8},
- {0x1.819d0b7158a4dp-1, 0x1.107fbbe01107cp-1, -0x1.0feeb40897d4ep-2,
- 0x1.50e5afb95f5d6p-4, 0x1.2a7c2f0c7495dp-7, -0x1.12bd2bb5062cdp-5,
- 0x1.93e8ceb89afebp-6, -0x1.10da9b8c6b731p-7},
- {0x1.921fb54442d18p-1, 0x1.fffffffffffebp-2, -0x1.fffffffffcbbcp-3,
- 0x1.555555564e2fep-4, -0x1.20b17d5dd89dcp-30, -0x1.9999c5ad71711p-6,
- 0x1.5558b76e7aaf9p-6, -0x1.236e803c6c1f6p-7},
+ 0x1.3b13b13b13b14p-4, -0x1.1111111111111p-4, 0x1.e1e1e1e1e1e1ep-5},
+ {0x1.ff55bb72cfde9p-5, 0x1.fe01fe01fe007p-1, -0x1.fc05f809ed8dap-5,
+ -0x1.4d69303afe04ep-2, 0x1.f61bc3e8349cp-5, 0x1.820839278756bp-3,
+ -0x1.eda4de1c6bf3fp-5, -0x1.0514d42d64a63p-3, 0x1.db3746a442dcbp-5},
+ {0x1.fd5ba9aac2f6ep-4, 0x1.f81f81f81f813p-1, -0x1.f05e09d0dc378p-4,
+ -0x1.368c3aa719215p-2, 0x1.d9b16b33ff9c9p-4, 0x1.40488f9c6262ap-3,
+ -0x1.ba55933e62ea5p-4, -0x1.64c6a15cd9116p-4, 0x1.9273d5939a75ap-4},
+ {0x1.7b97b4bce5b02p-3, 0x1.ee9c7f8458e05p-1, -0x1.665c226d6961p-3,
+ -0x1.1344bb7391703p-2, 0x1.42aca8b0081b9p-3, 0x1.c32d9381d7c03p-4,
+ -0x1.13e970672e246p-3, -0x1.181ed934dd733p-5, 0x1.bad81ea190c08p-4},
+ {0x1.f5b75f92c80ddp-3, 0x1.e1e1e1e1e1e2cp-1, -0x1.c5894d10d363dp-3,
+ -0x1.ce6de025f9f5ep-3, 0x1.78a3a07c8dd7fp-3, 0x1.dd5f5180f386ep-5,
+ -0x1.1b1f513c4536bp-3, 0x1.0df852e58c43cp-6, 0x1.722e7a7e42505p-4},
+ {0x1.362773707ebccp-2, 0x1.d272ca3fc5b2ep-1, -0x1.0997e8aeca8fbp-2,
+ -0x1.6cf6666e5e693p-3, 0x1.8dd1e907e88adp-3, 0x1.24849ac0caa5dp-7,
+ -0x1.f496be486229dp-4, 0x1.b7d54b8e759ecp-5, 0x1.d39c0d39c3922p-5},
+ {0x1.6f61941e4def1p-2, 0x1.c0e070381c0f2p-1, -0x1.2726dd135d9eep-2,
+ -0x1.09f37b39b70e4p-3, 0x1.85eacdaadd712p-3, -0x1.04d66340d5b9p-5,
+ -0x1.8056b15a22b98p-4, 0x1.29baf494ad3ddp-4, 0x1.52d5881322a7ap-6},
+ {0x1.a64eec3cc23fdp-2, 0x1.adbe87f94906ap-1, -0x1.3b9d8eab55addp-2,
+ -0x1.57c09646eb7p-4, 0x1.6795319e3b8dfp-3, -0x1.f2d89b5ef31bep-5,
+ -0x1.f38aac26203cap-5, 0x1.3262802235e3fp-4, -0x1.2afd6b9a57d66p-7},
+ {0x1.dac670561bb4fp-2, 0x1.99999999999ap-1, -0x1.47ae147adff11p-2,
+ -0x1.5d867c40188b7p-5, 0x1.3a92a2df85e7ap-3, -0x1.3ec457c46e851p-4,
+ -0x1.ec1b9777e2e5bp-6, 0x1.0a542992a821ep-4, -0x1.ccffbe2f0d945p-6},
+ {0x1.0657e94db30dp-1, 0x1.84f00c2780615p-1, -0x1.4c62cb562defap-2,
+ -0x1.e6495b3c14e03p-8, 0x1.063c2fa617bfcp-3, -0x1.58b782d9907aap-4,
+ -0x1.41e6ff524b7fp-8, 0x1.937dfff3205a7p-5, -0x1.0fb1fd1c729dp-5},
+ {0x1.1e00babdefeb4p-1, 0x1.702e05c0b816ep-1, -0x1.4af2b78215fbep-2,
+ 0x1.5d0b7e9f36997p-6, 0x1.a1247cb978debp-4, -0x1.519e1457734cap-4,
+ 0x1.a755cf86b5bfbp-7, 0x1.096d174284564p-5, -0x1.081adf539ad58p-5},
+ {0x1.345f01cce37bbp-1, 0x1.5babcc647fa8ep-1, -0x1.449db09426a6dp-2,
+ 0x1.655caac5896dap-5, 0x1.3bbbd22d05a61p-4, -0x1.34a2febee042fp-4,
+ 0x1.84df9c8269e34p-6, 0x1.200e8176c899ap-6, -0x1.c00b23c3ce222p-6},
+ {0x1.4978fa3269ee1p-1, 0x1.47ae147ae1477p-1, -0x1.3a92a3055231ap-2,
+ 0x1.ec21b515a4a2p-5, 0x1.c2f8b81f9a0d2p-5, -0x1.0ba9964125453p-4,
+ 0x1.d7b5614777a05p-6, 0x1.971e91ed73595p-8, -0x1.3fc375a78dc74p-6},
+ {0x1.5d58987169b18p-1, 0x1.34679ace01343p-1, -0x1.2ddfb039136e5p-2,
+ 0x1.2491307b9fb73p-4, 0x1.29c7e4886dc22p-5, -0x1.bca78bcca83ap-5,
+ 0x1.e63efd7cbe1ddp-6, -0x1.8ea6c4f03b42dp-10, -0x1.9385b5c3a6997p-7},
+ {0x1.700a7c5784634p-1, 0x1.21fb78121fb76p-1, -0x1.1f6a8499e5d1ap-2,
+ 0x1.41b15e5e29423p-4, 0x1.59bc953163345p-6, -0x1.63b54b13184ddp-5,
+ 0x1.c9086666d213p-6, -0x1.90c3b4ad8d4bcp-8, -0x1.80f08ed9f6f57p-8},
+ {0x1.819d0b7158a4dp-1, 0x1.107fbbe01107ep-1, -0x1.0feeb4089670ep-2,
+ 0x1.50e5afb93f5cbp-4, 0x1.2a7c2adffeffbp-7, -0x1.12bd29b4f1b43p-5,
+ 0x1.93f71f0eb00eap-6, -0x1.10ece5ad30e28p-7, -0x1.db1a76bcd2b9cp-10},
+ {0x1.921fb54442d18p-1, 0x1.ffffffffffffep-2, -0x1.fffffffffc51cp-3,
+ 0x1.555555557002ep-4, -0x1.a88260c338e75p-30, -0x1.99999f9a7614fp-6,
+ 0x1.555e31a1e15e9p-6, -0x1.245240d65e629p-7, -0x1.fa9ba66478903p-11},
};
} // namespace LIBC_NAMESPACE
diff --git a/libc/src/math/generic/inv_trigf_utils.h b/libc/src/math/generic/inv_trigf_utils.h
index c621f27e101460..8ee478c5bfa4d8 100644
--- a/libc/src/math/generic/inv_trigf_utils.h
+++ b/libc/src/math/generic/inv_trigf_utils.h
@@ -19,7 +19,7 @@ namespace LIBC_NAMESPACE {
constexpr double M_MATH_PI = 0x1.921fb54442d18p+1;
constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
-extern double ATAN_COEFFS[17][8];
+extern double ATAN_COEFFS[17][9];
// For |x| <= 1/32 and 1 <= i <= 16, return Q(x) such that:
// Q(x) ~ (atan(x + i/16) - atan(i/16)) / x.
@@ -29,10 +29,11 @@ LIBC_INLINE double atan_eval(double x, int i) {
double c0 = fputil::multiply_add(x, ATAN_COEFFS[i][2], ATAN_COEFFS[i][1]);
double c1 = fputil::multiply_add(x, ATAN_COEFFS[i][4], ATAN_COEFFS[i][3]);
double c2 = fputil::multiply_add(x, ATAN_COEFFS[i][6], ATAN_COEFFS[i][5]);
+ double c3 = fputil::multiply_add(x, ATAN_COEFFS[i][8], ATAN_COEFFS[i][7]);
double x4 = x2 * x2;
double d1 = fputil::multiply_add(x2, c1, c0);
- double d2 = fputil::multiply_add(x2, ATAN_COEFFS[i][7], c2);
+ double d2 = fputil::multiply_add(x2, c3, c2);
double p = fputil::multiply_add(x4, d2, d1);
return p;
}
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 026bcd12928bdf..f8f0f8ba7b6f63 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1719,6 +1719,19 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ atan2f_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ atan2f_test.cpp
+ DEPENDS
+ libc.include.math
+ libc.src.math.atan2f
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_subdirectory(generic)
add_subdirectory(smoke)
diff --git a/libc/test/src/math/atan2f_test.cpp b/libc/test/src/math/atan2f_test.cpp
new file mode 100644
index 00000000000000..71b36195c96306
--- /dev/null
+++ b/libc/test/src/math/atan2f_test.cpp
@@ -0,0 +1,132 @@
+//===-- Unittests for atan2f ----------------------------------------------===//
+//
+// 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 "include/llvm-libc-macros/math-macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/atan2f.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcAtan2fTest = LIBC_NAMESPACE::testing::FPTest<float>;
+using LIBC_NAMESPACE::testing::tlog;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+TEST_F(LlvmLibcAtan2fTest, TrickyInputs) {
+ constexpr int N = 17;
+ mpfr::BinaryInput<float> INPUTS[N] = {
+ {0x1.0cb3a4p+20f, 0x1.4ebacp+22f}, {0x1.12215p+1f, 0x1.4fabfcp+22f},
+ {-0x1.13baaep+41f, 0x1.5bd22ep+23f}, {0x1.1ff7dcp+41f, 0x1.aec0a6p+23f},
+ {0x1.2bc794p+23f, 0x1.0bc0c6p+23f}, {0x1.2fba3ap+42f, 0x1.f99456p+23f},
+ {0x1.5ea1f8p+27f, 0x1.f2a1aep+23f}, {0x1.7a931p+44f, 0x1.352ac4p+22f},
+ {0x1.8802bcp+21f, 0x1.8f130ap+23f}, {0x1.658ef8p+17f, 0x1.3c00f4p+22f},
+ {0x1.69fb0cp+21f, 0x1.39e4c4p+23f}, {0x1.8eb24cp+11f, 0x1.36518p+23f},
+ {0x1.9e7ebp+30f, 0x1.d80522p+23f}, {0x1.b4bdeep+19f, 0x1.c19b4p+23f},
+ {0x1.bc201p+43f, 0x1.617346p+23f}, {0x1.c96c3cp+20f, 0x1.c01d1ep+23f},
+ {0x1.781fcp+28f, 0x1.dcb3cap+23f},
+ };
+
+ for (int i = 0; i < N; ++i) {
+ float x = INPUTS[i].x;
+ float y = INPUTS[i].y;
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atan2, INPUTS[i],
+ LIBC_NAMESPACE::atan2f(x, y), 0.5);
+ INPUTS[i].x = -INPUTS[i].x;
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atan2, INPUTS[i],
+ LIBC_NAMESPACE::atan2f(-x, y), 0.5);
+ INPUTS[i].y = -INPUTS[i].y;
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atan2, INPUTS[i],
+ LIBC_NAMESPACE::atan2f(-x, -y), 0.5);
+ INPUTS[i].x = -INPUTS[i].x;
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atan2, INPUTS[i],
+ LIBC_NAMESPACE::atan2f(x, -y), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcAtan2fTest, InFloatRange) {
+ constexpr uint32_t X_COUNT = 1'23;
+ constexpr uint32_t X_START = FPBits(0.25f).uintval();
+ constexpr uint32_t X_STOP = FPBits(4.0f).uintval();
+ constexpr uint32_t X_STEP = (X_STOP - X_START) / X_COUNT;
+
+ constexpr uint32_t Y_COUNT = 1'37;
+ constexpr uint32_t Y_START = FPBits(0.25f).uintval();
+ constexpr uint32_t Y_STOP = FPBits(4.0f).uintval();
+ constexpr uint32_t Y_STEP = (Y_STOP - Y_START) / Y_COUNT;
+
+ auto test = [&](mpfr::RoundingMode rounding_mode) {
+ mpfr::ForceRoundingMode __r(rounding_mode);
+ if (!__r.success)
+ return;
+
+ uint64_t fails = 0;
+ uint64_t count = 0;
+ uint64_t cc = 0;
+ float mx, my, mr = 0.0;
+ double tol = 0.5;
+
+ for (uint32_t i = 0, v = X_START; i <= X_COUNT; ++i, v += X_STEP) {
+ float x = FPBits(v).get_val();
+ if (isnan(x) || isinf(x) || x < 0.0)
+ continue;
+
+ for (uint32_t j = 0, w = Y_START; j <= Y_COUNT; ++j, w += Y_STEP) {
+ float y = FPBits(w).get_val();
+ if (isnan(y) || isinf(y))
+ continue;
+
+ LIBC_NAMESPACE::libc_errno = 0;
+ float result = LIBC_NAMESPACE::atan2f(x, y);
+ ++cc;
+ if (isnan(result) || isinf(result))
+ continue;
+
+ ++count;
+ mpfr::BinaryInput<float> inputs{x, y};
+
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Atan2, inputs,
+ result, 0.5, rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
+ mpfr::Operation::Atan2, inputs, result, tol, rounding_mode)) {
+ mx = x;
+ my = y;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ }
+ if (fails || (count < cc)) {
+ tlog << " Atan2f failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n"
+ << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
+ }
+ if (fails) {
+ mpfr::BinaryInput<float> inputs{mx, my};
+ EXPECT_MPFR_MATCH(mpfr::Operation::Atan2, inputs, mr, 0.5, rounding_mode);
+ }
+ };
+
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
+
+ tlog << " Test Rounding Downward...\n";
+ test(mpfr::RoundingMode::Downward);
+
+ tlog << " Test Rounding Upward...\n";
+ test(mpfr::RoundingMode::Upward);
+
+ tlog << " Test Rounding Toward Zero...\n";
+ test(mpfr::RoundingMode::TowardZero);
+}
diff --git a/libc/test/src/math/atanf_test.cpp b/libc/test/src/math/atanf_test.cpp
index e51932fc495525..58b0eadd63f8d6 100644
--- a/libc/test/src/math/atanf_test.cpp
+++ b/libc/test/src/math/atanf_test.cpp
@@ -55,12 +55,15 @@ TEST_F(LlvmLibcAtanfTest, InFloatRange) {
TEST_F(LlvmLibcAtanfTest, SpecialValues) {
uint32_t val_arr[] = {
0x3d8d6b23U, // x = 0x1.1ad646p-4f
+ 0x3dbb6ac7U, // x = 0x1.76d58ep-4f
0x3feefcfbU, // x = 0x1.ddf9f6p+0f
+ 0x3ffe2ec1U, // x = 0x1.fc5d82p+0f
0xbd8d6b23U, // x = -0x1.1ad646p-4f
+ 0xbdbb6ac7U, // x = -0x1.76d58ep-4f
0xbfeefcfbU, // x = -0x1.ddf9f6p+0f
+ 0xbffe2ec1U, // x = -0x1.fc5d82p+0
0x7F800000U, // x = +Inf
0xFF800000U, // x = -Inf
- 0xbffe2ec1U, // x = -0x1.fc5d82p+0f
};
for (uint32_t v : val_arr) {
float x = FPBits(v).get_val();
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 3b756127fe21ea..28141d46c9334f 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -2751,6 +2751,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ atan2f_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ atan2f_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.math.atan2f
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
scalbn_test
SUITE
diff --git a/libc/test/src/math/smoke/atan2f_test.cpp b/libc/test/src/math/smoke/atan2f_test.cpp
new file mode 100644
index 00000000000000..ecac36b3a8c01f
--- /dev/null
+++ b/libc/test/src/math/smoke/atan2f_test.cpp
@@ -0,0 +1,50 @@
+//===-- Unittests for atan2f ----------------------------------------------===//
+//
+// 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 "include/llvm-libc-macros/math-macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/atan2f.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcAtan2fTest = LIBC_NAMESPACE::testing::FPTest<float>;
+
+TEST_F(LlvmLibcAtan2fTest, SpecialNumbers) {
+ LIBC_NAMESPACE::libc_errno = 0;
+
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT);
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::atan2f(aNaN, zero));
+ EXPECT_FP_EXCEPTION(0);
+ EXPECT_MATH_ERRNO(0);
+
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT);
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::atan2f(1.0f, aNaN));
+ EXPECT_FP_EXCEPTION(0);
+ EXPECT_MATH_ERRNO(0);
+
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT);
+ EXPECT_FP_EQ_ALL_ROUNDING(0.0f, LIBC_NAMESPACE::atan2f(zero, zero));
+ EXPECT_FP_EXCEPTION(0);
+ EXPECT_MATH_ERRNO(0);
+
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT);
+ EXPECT_FP_EQ_ALL_ROUNDING(-0.0f, LIBC_NAMESPACE::atan2f(-0.0f, zero));
+ EXPECT_FP_EXCEPTION(0);
+ EXPECT_MATH_ERRNO(0);
+
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT);
+ EXPECT_FP_EQ_ALL_ROUNDING(0.0f, LIBC_NAMESPACE::atan2f(1.0f, inf));
+ EXPECT_FP_EXCEPTION(0);
+ EXPECT_MATH_ERRNO(0);
+
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT);
+ EXPECT_FP_EQ_ALL_ROUNDING(-0.0f, LIBC_NAMESPACE::atan2f(-1.0f, inf));
+ EXPECT_FP_EXCEPTION(0);
+ EXPECT_MATH_ERRNO(0);
+}
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 2e1c44e6fd5da9..a83f7a7ceb922c 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -181,6 +181,12 @@ class MPFRNumber {
return result;
}
+ MPFRNumber atan2(const MPFRNumber &b) {
+ MPFRNumber result(*this);
+ mpfr_atan2(result.value, value, b.value, mpfr_rounding);
+ return result;
+ }
+
MPFRNumber atanh() const {
MPFRNumber result(*this);
mpfr_atanh(result.value, value, mpfr_rounding);
@@ -623,6 +629,8 @@ binary_operation_one_output(Operation op, InputType x, InputType y,
MPFRNumber inputX(x, precision, rounding);
MPFRNumber inputY(y, precision, rounding);
switch (op) {
+ case Operation::Atan2:
+ return inputX.atan2(inputY);
case Operation::Fmod:
return inputX.fmod(inputY);
case Operation::Hypot:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 6164d78fa5adc4..d5ff590cd7bb69 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -68,6 +68,7 @@ enum class Operation : int {
// input and produce a single floating point number of the same type as
// output.
BeginBinaryOperationsSingleOutput,
+ Atan2,
Fmod,
Hypot,
Pow,
>From f70a641958bb7a62c2743726f9cc8f2c5d4a67f8 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue at google.com>
Date: Wed, 27 Mar 2024 00:23:11 -0400
Subject: [PATCH 2/5] Address comments.
---
libc/src/math/generic/CMakeLists.txt | 1 -
libc/src/math/generic/atan2f.cpp | 132 +++++++++++++++------------
libc/test/src/math/atan2f_test.cpp | 29 +++---
3 files changed, 90 insertions(+), 72 deletions(-)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 7e9c9ce7bb2748..b164d33e204b1a 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2864,7 +2864,6 @@ add_entrypoint_object(
-O3
DEPENDS
.inv_trigf_utils
- libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
diff --git a/libc/src/math/generic/atan2f.cpp b/libc/src/math/generic/atan2f.cpp
index 089e6c61984bea..31f3dcc89ea789 100644
--- a/libc/src/math/generic/atan2f.cpp
+++ b/libc/src/math/generic/atan2f.cpp
@@ -11,7 +11,6 @@
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/double_double.h"
-#include "src/__support/FPUtil/except_value_utils.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/FPUtil/rounding_mode.h"
@@ -85,6 +84,77 @@ constexpr fputil::DoubleDouble COEFFS[9] = {
return r;
}
+// Compute atan( num_d / den_d ) in double-double precision.
+// num_d = min(|x|, |y|)
+// den_d = max(|x|, |y|)
+// q_d = num_d / den_d
+// idx, k_d = round( 2^4 * num_d / den_d )
+// final_sign = sign of the final result
+// const_term = the constant term in the final expression.
+LIBC_INLINE float atan2f_double_double(double num_d, double den_d, double q_d,
+ int idx, double k_d, double final_sign,
+ const fputil::DoubleDouble &const_term) {
+ fputil::DoubleDouble q;
+ double num_r, den_r;
+
+ if (idx != 0) {
+ // The following range reduction is accurate even without fma for
+ // 1/16 <= n/d <= 1.
+ // atan(n/d) - atan(idx/16) = atan((n/d - idx/16) / (1 + (n/d) * (idx/16)))
+ // = atan((n - d*(idx/16)) / (d + n*idx/16))
+ k_d *= 0x1.0p-4;
+ num_r = fputil::multiply_add(k_d, -den_d, num_d); // Exact
+ den_r = fputil::multiply_add(k_d, num_d, den_d); // Exact
+ q.hi = num_r / den_r;
+ } else {
+ // For 0 < n/d < 1/16, we just need to calculate the lower part of their
+ // quotient.
+ q.hi = q_d;
+ num_r = num_d;
+ den_r = den_d;
+ }
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ q.lo = fputil::multiply_add(q.hi, -den_r, num_r) / den_r;
+#else
+ // Compute `(num_r - q.hi * den_r) / den_r` accurately without FMA
+ // instructions.
+ fputil::DoubleDouble q_hi_dd = split_d(q.hi);
+ double t1 = fputil::multiply_add(q_hi_dd.hi, -den_r, num_r); // Exact
+ double t2 = fputil::multiply_add(q_hi_dd.lo, -den_r, t1);
+ q.lo = t2 / den_r;
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+ // Taylor polynomial, evaluating using Horner's scheme:
+ // P = x - x^3/3 + x^5/5 -x^7/7 + x^9/9 - x^11/11 + x^13/13 - x^15/15
+ // + x^17/17
+ // = x*(1 + x^2*(-1/3 + x^2*(1/5 + x^2*(-1/7 + x^2*(1/9 + x^2*
+ // *(-1/11 + x^2*(1/13 + x^2*(-1/15 + x^2 * 1/17))))))))
+ fputil::DoubleDouble q2 = fputil::quick_mult(q, q);
+ fputil::DoubleDouble p_dd =
+ fputil::polyeval(q2, COEFFS[0], COEFFS[1], COEFFS[2], COEFFS[3],
+ COEFFS[4], COEFFS[5], COEFFS[6], COEFFS[7], COEFFS[8]);
+ fputil::DoubleDouble r_dd =
+ fputil::add(const_term, fputil::multiply_add(q, p_dd, ATAN_I[idx]));
+ r_dd.hi *= final_sign;
+ r_dd.lo *= final_sign;
+
+ // Make sure the sum is normalized:
+ fputil::DoubleDouble rr = fputil::exact_add(r_dd.hi, r_dd.lo);
+ // Round to odd.
+ uint64_t rr_bits = cpp::bit_cast<uint64_t>(rr.hi);
+ if (LIBC_UNLIKELY(((rr_bits & 0xfff'ffff) == 0) && (rr.lo != 0.0))) {
+ Sign hi_sign = fputil::FPBits<double>(rr.hi).sign();
+ Sign lo_sign = fputil::FPBits<double>(rr.lo).sign();
+ if (hi_sign == lo_sign) {
+ ++rr_bits;
+ } else if ((rr_bits & fputil::FPBits<double>::FRACTION_MASK) > 0) {
+ --rr_bits;
+ }
+ }
+
+ return static_cast<float>(cpp::bit_cast<double>(rr_bits));
+}
+
} // anonymous namespace
// There are several range reduction steps we can take for atan2(y, x) as
@@ -168,7 +238,6 @@ LLVM_LIBC_FUNCTION(float, atan2f, (float y, float x)) {
uint32_t y_abs = y_bits.uintval();
uint32_t max_abs = x_abs > y_abs ? x_abs : y_abs;
uint32_t min_abs = x_abs <= y_abs ? x_abs : y_abs;
- bool recip = x_abs < y_abs;
if (LIBC_UNLIKELY(max_abs >= 0x7f80'0000U || min_abs == 0U)) {
if (x_bits.is_nan() || y_bits.is_nan())
@@ -195,15 +264,15 @@ LLVM_LIBC_FUNCTION(float, atan2f, (float y, float x)) {
return static_cast<float>(r);
}
+ bool recip = x_abs < y_abs;
double final_sign = IS_NEG[(x_sign != y_sign) != recip];
fputil::DoubleDouble const_term = CONST_ADJ[x_sign][y_sign][recip];
double num_d = static_cast<double>(FPBits(min_abs).get_val());
double den_d = static_cast<double>(FPBits(max_abs).get_val());
double q_d = num_d / den_d;
- int idx;
double k_d = fputil::nearest_integer(q_d * 0x1.0p4f);
- idx = static_cast<int>(k_d);
+ int idx = static_cast<int>(k_d);
q_d = fputil::multiply_add(k_d, -0x1.0p-4, q_d);
double p, r;
@@ -239,59 +308,8 @@ LLVM_LIBC_FUNCTION(float, atan2f, (float y, float x)) {
if (LIBC_LIKELY(r_bits > LOWER_ERR && r_bits < UPPER_ERR))
return static_cast<float>(r);
- // Use double-double.
- fputil::DoubleDouble q;
- double num_r, den_r;
-
- if (idx != 0) {
- // The following range reduction is accurate even without fma for
- // 1/16 <= n/d <= 1.
- // atan(n/d) - atan(idx/16) = atan((n/d - idx/16) / (1 + (n/d) * (idx/16)))
- // = atan((n - d*(idx/16)) / (d + n*idx/16))
- k_d *= 0x1.0p-4;
- num_r = fputil::multiply_add(k_d, -den_d, num_d); // Exact
- den_r = fputil::multiply_add(k_d, num_d, den_d); // Exact
- q.hi = num_r / den_r;
- } else {
- // For 0 < n/d < 1/16, we just need to calculate the lower part of their
- // quotient.
- q.hi = q_d;
- num_r = num_d;
- den_r = den_d;
- }
-#ifdef LIBC_TARGET_CPU_HAS_FMA
- q.lo = fputil::multiply_add(q.hi, -den_r, num_r) / den_r;
-#else
- fputil::DoubleDouble q_hi_dd = split_d(q.hi);
- double t1 = fputil::multiply_add(q_hi_dd.hi, -den_r, num_r); // Exact
- double t2 = fputil::multiply_add(q_hi_dd.lo, -den_r, t1);
- q.lo = t2 / den_r;
-#endif // LIBC_TARGET_CPU_HAS_FMA
-
- fputil::DoubleDouble q2 = fputil::quick_mult(q, q);
- fputil::DoubleDouble p_dd =
- fputil::polyeval(q2, COEFFS[0], COEFFS[1], COEFFS[2], COEFFS[3],
- COEFFS[4], COEFFS[5], COEFFS[6], COEFFS[7], COEFFS[8]);
- fputil::DoubleDouble r_dd =
- fputil::add(const_term, fputil::multiply_add(q, p_dd, ATAN_I[idx]));
- r_dd.hi *= final_sign;
- r_dd.lo *= final_sign;
-
- // Make sure the sum is normalized:
- fputil::DoubleDouble rr = fputil::exact_add(r_dd.hi, r_dd.lo);
- // Round to odd.
- uint64_t rr_bits = cpp::bit_cast<uint64_t>(rr.hi);
- if (LIBC_UNLIKELY(((rr_bits & 0xfff'ffff) == 0) && (rr.lo != 0.0))) {
- Sign hi_sign = fputil::FPBits<double>(rr.hi).sign();
- Sign lo_sign = fputil::FPBits<double>(rr.lo).sign();
- if (hi_sign == lo_sign) {
- ++rr_bits;
- } else if ((rr_bits & fputil::FPBits<double>::FRACTION_MASK) > 0) {
- --rr_bits;
- }
- }
-
- return static_cast<float>(cpp::bit_cast<double>(rr_bits));
+ return atan2f_double_double(num_d, den_d, q_d, idx, k_d, final_sign,
+ const_term);
}
} // namespace LIBC_NAMESPACE
diff --git a/libc/test/src/math/atan2f_test.cpp b/libc/test/src/math/atan2f_test.cpp
index 71b36195c96306..22562328d17b5c 100644
--- a/libc/test/src/math/atan2f_test.cpp
+++ b/libc/test/src/math/atan2f_test.cpp
@@ -62,13 +62,13 @@ TEST_F(LlvmLibcAtan2fTest, InFloatRange) {
auto test = [&](mpfr::RoundingMode rounding_mode) {
mpfr::ForceRoundingMode __r(rounding_mode);
- if (!__r.success)
+ if (!__r.sutotal_countess)
return;
uint64_t fails = 0;
- uint64_t count = 0;
- uint64_t cc = 0;
- float mx, my, mr = 0.0;
+ uint64_t finite_count = 0;
+ uint64_t total_count = 0;
+ float failed_x, failed_y, failed_r = 0.0;
double tol = 0.5;
for (uint32_t i = 0, v = X_START; i <= X_COUNT; ++i, v += X_STEP) {
@@ -83,11 +83,11 @@ TEST_F(LlvmLibcAtan2fTest, InFloatRange) {
LIBC_NAMESPACE::libc_errno = 0;
float result = LIBC_NAMESPACE::atan2f(x, y);
- ++cc;
+ ++total_count;
if (isnan(result) || isinf(result))
continue;
- ++count;
+ ++finite_count;
mpfr::BinaryInput<float> inputs{x, y};
if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Atan2, inputs,
@@ -95,9 +95,9 @@ TEST_F(LlvmLibcAtan2fTest, InFloatRange) {
++fails;
while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
mpfr::Operation::Atan2, inputs, result, tol, rounding_mode)) {
- mx = x;
- my = y;
- mr = result;
+ failed_x = x;
+ failed_y = y;
+ failed_r = result;
if (tol > 1000.0)
break;
@@ -107,14 +107,15 @@ TEST_F(LlvmLibcAtan2fTest, InFloatRange) {
}
}
}
- if (fails || (count < cc)) {
- tlog << " Atan2f failed: " << fails << "/" << count << "/" << cc
- << " tests.\n"
+ if (fails || (finite_count < total_count)) {
+ tlog << " Atan2f failed: " << fails << "/" << finite_count << "/"
+ << total_count << " tests.\n"
<< " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
}
if (fails) {
- mpfr::BinaryInput<float> inputs{mx, my};
- EXPECT_MPFR_MATCH(mpfr::Operation::Atan2, inputs, mr, 0.5, rounding_mode);
+ mpfr::BinaryInput<float> inputs{failed_x, failed_y};
+ EXPECT_MPFR_MATCH(mpfr::Operation::Atan2, inputs, failed_r, 0.5,
+ rounding_mode);
}
};
>From 3601695f421d6be218082e2f9924d49a91329100 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue at google.com>
Date: Wed, 27 Mar 2024 00:28:13 -0400
Subject: [PATCH 3/5] Fix test.
---
libc/test/src/math/atan2f_test.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/test/src/math/atan2f_test.cpp b/libc/test/src/math/atan2f_test.cpp
index 22562328d17b5c..343e7601b0392b 100644
--- a/libc/test/src/math/atan2f_test.cpp
+++ b/libc/test/src/math/atan2f_test.cpp
@@ -62,7 +62,7 @@ TEST_F(LlvmLibcAtan2fTest, InFloatRange) {
auto test = [&](mpfr::RoundingMode rounding_mode) {
mpfr::ForceRoundingMode __r(rounding_mode);
- if (!__r.sutotal_countess)
+ if (!__r.success)
return;
uint64_t fails = 0;
>From e73c7510fc45d6eb5260bbf02756d2be91eec9e6 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Thu, 28 Mar 2024 17:50:59 +0000
Subject: [PATCH 4/5] Simplify the implementation by using different polynomial
for small cases.
---
libc/src/math/generic/atan2f.cpp | 22 +++-------------------
libc/src/math/generic/atanf.cpp | 10 +++++++---
libc/src/math/generic/inv_trigf_utils.cpp | 12 +++++++-----
libc/src/math/generic/inv_trigf_utils.h | 2 +-
4 files changed, 18 insertions(+), 28 deletions(-)
diff --git a/libc/src/math/generic/atan2f.cpp b/libc/src/math/generic/atan2f.cpp
index 31f3dcc89ea789..b999f8ec7578e7 100644
--- a/libc/src/math/generic/atan2f.cpp
+++ b/libc/src/math/generic/atan2f.cpp
@@ -275,25 +275,9 @@ LLVM_LIBC_FUNCTION(float, atan2f, (float y, float x)) {
int idx = static_cast<int>(k_d);
q_d = fputil::multiply_add(k_d, -0x1.0p-4, q_d);
- double p, r;
-
- if (LIBC_UNLIKELY(idx == 0)) {
- // For |x| < 1/16, we use Taylor polynomial:
- // atan(x) ~ x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - x^11/11 + x^13/13
- double q2 = q_d * q_d;
- double c0 = fputil::multiply_add(q2, ATAN_COEFFS[0][2], ATAN_COEFFS[0][1]);
- double c1 = fputil::multiply_add(q2, ATAN_COEFFS[0][4], ATAN_COEFFS[0][3]);
- double c2 = fputil::multiply_add(q2, ATAN_COEFFS[0][6], ATAN_COEFFS[0][5]);
- double q4 = q2 * q2;
- double q3 = q_d * q2;
- double d0 = fputil::polyeval(q4, c0, c1, c2);
- p = fputil::multiply_add(q3, d0, q_d);
- r = final_sign * (p + const_term.hi);
- } else {
- p = atan_eval(q_d, idx);
- r = final_sign *
- fputil::multiply_add(q_d, p, const_term.hi + ATAN_COEFFS[idx][0]);
- }
+ double p = atan_eval(q_d, idx);
+ double r = final_sign *
+ fputil::multiply_add(q_d, p, const_term.hi + ATAN_COEFFS[idx][0]);
constexpr uint32_t LOWER_ERR = 4;
// Mask sticky bits in double precision before rounding to single precision.
diff --git a/libc/src/math/generic/atanf.cpp b/libc/src/math/generic/atanf.cpp
index a58bbdd5033f87..9fa1a331c9c038 100644
--- a/libc/src/math/generic/atanf.cpp
+++ b/libc/src/math/generic/atanf.cpp
@@ -60,11 +60,15 @@ LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
}
// Use Taylor polynomial:
// atan(x) ~ x * (1 - x^2 / 3 + x^4 / 5 - x^6 / 7 + x^8 / 9 - x^10 / 11).
+ constexpr double ATAN_TAYLOR[6] = {
+ 0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,
+ -0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4, -0x1.745d1745d1746p-4,
+ };
double x2 = x_d * x_d;
double x4 = x2 * x2;
- double c0 = fputil::multiply_add(x2, ATAN_COEFFS[0][1], ATAN_COEFFS[0][0]);
- double c1 = fputil::multiply_add(x2, ATAN_COEFFS[0][3], ATAN_COEFFS[0][2]);
- double c2 = fputil::multiply_add(x2, ATAN_COEFFS[0][5], ATAN_COEFFS[0][4]);
+ double c0 = fputil::multiply_add(x2, ATAN_TAYLOR[1], ATAN_TAYLOR[0]);
+ double c1 = fputil::multiply_add(x2, ATAN_TAYLOR[3], ATAN_TAYLOR[2]);
+ double c2 = fputil::multiply_add(x2, ATAN_TAYLOR[5], ATAN_TAYLOR[4]);
double p = fputil::polyeval(x4, c0, c1, c2);
double r = fputil::multiply_add(x_d, p, const_term);
return static_cast<float>(r);
diff --git a/libc/src/math/generic/inv_trigf_utils.cpp b/libc/src/math/generic/inv_trigf_utils.cpp
index 06da93db96350b..19c8b997dc4ed5 100644
--- a/libc/src/math/generic/inv_trigf_utils.cpp
+++ b/libc/src/math/generic/inv_trigf_utils.cpp
@@ -21,15 +21,17 @@ namespace LIBC_NAMESPACE {
// coeff(P, 3), ",", coeff(P, 4), ",", coeff(P, 5), ",", coeff(P, 6),
// ",", coeff(P, 7), ",", coeff(P, 8), "},");
// };
-// For i = 0, ATAN_COEFFS[0][j] = (-1)^j * (1/(2*j + 1)) is the odd coefficients
-// of the Taylor polynomial of atan(x).
+// For i = 0, the polynomial is generated by:
+// > P = fpminimax(atan(x)/x, 7, [|1, D...|], [0, 1/32]);
+// > dirtyinfnorm((atan(x) - x*P)/x, [0, 1/32]);
+// 0x1.feb2fcdba66447ccbe28a1a0f935b51678a718fb1p-59
// Notice that degree-7 is good enough for atanf, but degree-8 helps reduce the
// error bounds for atan2f's fast pass 16 times, and it does not affect the
// performance of atanf much.
double ATAN_COEFFS[17][9] = {
- {0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,
- -0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4, -0x1.745d1745d1746p-4,
- 0x1.3b13b13b13b14p-4, -0x1.1111111111111p-4, 0x1.e1e1e1e1e1e1ep-5},
+ {0.0, 1.0, 0x1.3f8d76d26d61bp-47, -0x1.5555555574cd8p-2,
+ 0x1.0dde5d06878eap-29, 0x1.99997738acc77p-3, 0x1.2c43eac9797cap-16,
+ -0x1.25fb020007dbdp-3, 0x1.c1b6c31d7b0aep-7},
{0x1.ff55bb72cfde9p-5, 0x1.fe01fe01fe007p-1, -0x1.fc05f809ed8dap-5,
-0x1.4d69303afe04ep-2, 0x1.f61bc3e8349cp-5, 0x1.820839278756bp-3,
-0x1.eda4de1c6bf3fp-5, -0x1.0514d42d64a63p-3, 0x1.db3746a442dcbp-5},
diff --git a/libc/src/math/generic/inv_trigf_utils.h b/libc/src/math/generic/inv_trigf_utils.h
index 8ee478c5bfa4d8..e60c367d7b46b8 100644
--- a/libc/src/math/generic/inv_trigf_utils.h
+++ b/libc/src/math/generic/inv_trigf_utils.h
@@ -21,7 +21,7 @@ constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
extern double ATAN_COEFFS[17][9];
-// For |x| <= 1/32 and 1 <= i <= 16, return Q(x) such that:
+// For |x| <= 1/32 and 0 <= i <= 16, return Q(x) such that:
// Q(x) ~ (atan(x + i/16) - atan(i/16)) / x.
LIBC_INLINE double atan_eval(double x, int i) {
double x2 = x * x;
>From 667096e1433b2bc12a7dbc161d02acca3ce3d449 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 29 Mar 2024 15:45:26 +0000
Subject: [PATCH 5/5] Update Veltkamp splitting constant comments.
---
libc/src/math/generic/atan2f.cpp | 13 ++++++++++---
1 file changed, 10 insertions(+), 3 deletions(-)
diff --git a/libc/src/math/generic/atan2f.cpp b/libc/src/math/generic/atan2f.cpp
index b999f8ec7578e7..d50ac660883684 100644
--- a/libc/src/math/generic/atan2f.cpp
+++ b/libc/src/math/generic/atan2f.cpp
@@ -69,10 +69,17 @@ constexpr fputil::DoubleDouble COEFFS[9] = {
// Veltkamp's splitting of a double precision into hi + lo, where the hi part is
// slightly smaller than an even split, so that the product of
-// hi * s * k is exact,
+// hi * (s1 * k + s2) is exact,
// where:
-// s is single precsion,
-// 0 < k < 16 is an integer.
+// s1, s2 are single precsion,
+// 1/16 <= s1/s2 <= 1
+// 1/16 <= k <= 1 is an integer.
+// So the maximal precision of (s1 * k + s2) is:
+// prec(s1 * k + s2) = 2 + log2(msb(s2)) - log2(lsb(k_d * s1))
+// = 2 + log2(msb(s1)) + 4 - log2(lsb(k_d)) - log2(lsb(s1))
+// = 2 + log2(lsb(s1)) + 23 + 4 - (-4) - log2(lsb(s1))
+// = 33.
+// Thus, the Veltkamp splitting constant is C = 2^33 + 1.
// This is used when FMA instruction is not available.
[[maybe_unused]] LIBC_INLINE constexpr fputil::DoubleDouble split_d(double a) {
fputil::DoubleDouble r{0.0, 0.0};
More information about the libc-commits
mailing list