[libc-commits] [libc] [llvm] [libc][math] Add float-only implementation of tanf (PR #233046)
Prabhu Rajasekaran via libc-commits
libc-commits at lists.llvm.org
Sat Oct 17 00:46:27 PDT 2026
https://github.com/Prabhuk created https://github.com/llvm/llvm-project/pull/233046
Reorganize tanf in the same way as sinf and cosf (#226529):
src/__support/math/tanf.h now selects the implementation, and
tanf_<type>_eval.h implements tanf with <type> as the intermediate
computational type. The existing implementation moves unchanged to
tanf_double_eval.h.
The new float-only implementation in tanf_float_eval.h is selected
under the same conditions as sinf/cosf float_eval:
LIBC_MATH_SKIP_ACCURATE_PASS, LIBC_MATH_SMALL_TABLES or
LIBC_MATH_INTERMEDIATE_COMP_IN_FLOAT, and FMA. It reuses the float range
reduction and the SIN_K_PI_OVER_8 table from sincosf_float_eval.h and
computes:
tan(x) = (sin_k + cos_k * tan(y)) / (cos_k - sin_k * tan(y)),
with x = k * pi/8 + y, |y| <= pi/16, and tan(y) approximated by a degree-9
odd Taylor polynomial (relative truncation error < 2^-30). Cancellation
can only happen when k = 0 or 4 (mod 8), where sin_k or cos_k is exactly
zero.
With this change, tanf no longer pulls in the double-precision
large_range_reduction and the 512-byte SIN_K_PI_OVER_32 table when the
float implementations are selected, so it shares the reduction tables
with sinf/cosf. For x86-64 with FMA and
LIBC_MATH_SKIP_ACCURATE_PASS | LIBC_MATH_SMALL_TABLES |
LIBC_MATH_INTERMEDIATE_COMP_IN_FLOAT, the tanf object goes from 1391 to
1029 bytes standalone (288 bytes of which are tables shared with
sinf/cosf).
Accuracy: the maximum error of the float implementation, checked
exhaustively in round-to-nearest mode against MPFR, is 3.37 ULP, which is
within the 3.5 ULP bound used for sinf/cosf float_eval. The unit and
exhaustive tests now exercise both implementations independently of the
configured selector.
Assisted-by: Gemini
>From e865a2d46f41e4ae74956994e38410fefe56d89d Mon Sep 17 00:00:00 2001
From: prabhukr <prabhukr at google.com>
Date: Sat, 17 Oct 2026 07:44:44 +0000
Subject: [PATCH] [libc][math] Add float-only implementation of tanf
Reorganize tanf in the same way as sinf and cosf (#226529):
src/__support/math/tanf.h now selects the implementation, and
tanf_<type>_eval.h implements tanf with <type> as the intermediate
computational type. The existing implementation moves unchanged to
tanf_double_eval.h.
The new float-only implementation in tanf_float_eval.h is selected
under the same conditions as sinf/cosf float_eval:
LIBC_MATH_SKIP_ACCURATE_PASS, LIBC_MATH_SMALL_TABLES or
LIBC_MATH_INTERMEDIATE_COMP_IN_FLOAT, and FMA. It reuses the float range
reduction and the SIN_K_PI_OVER_8 table from sincosf_float_eval.h and
computes:
tan(x) = (sin_k + cos_k * tan(y)) / (cos_k - sin_k * tan(y)),
with x = k * pi/8 + y, |y| <= pi/16, and tan(y) approximated by a degree-9
odd Taylor polynomial (relative truncation error < 2^-30). Cancellation
can only happen when k = 0 or 4 (mod 8), where sin_k or cos_k is exactly
zero.
With this change, tanf no longer pulls in the double-precision
large_range_reduction and the 512-byte SIN_K_PI_OVER_32 table when the
float implementations are selected, so it shares the reduction tables
with sinf/cosf. For x86-64 with FMA and
LIBC_MATH_SKIP_ACCURATE_PASS | LIBC_MATH_SMALL_TABLES |
LIBC_MATH_INTERMEDIATE_COMP_IN_FLOAT, the tanf object goes from 1391 to
1029 bytes standalone (288 bytes of which are tables shared with
sinf/cosf).
Accuracy: the maximum error of the float implementation, checked
exhaustively in round-to-nearest mode against MPFR, is 3.37 ULP, which is
within the 3.5 ULP bound used for sinf/cosf float_eval. The unit and
exhaustive tests now exercise both implementations independently of the
configured selector.
Assisted-by: Gemini
---
libc/src/__support/math/CMakeLists.txt | 30 ++-
libc/src/__support/math/tanf.h | 167 +++--------------
libc/src/__support/math/tanf_double_eval.h | 172 ++++++++++++++++++
libc/src/__support/math/tanf_float_eval.h | 106 +++++++++++
libc/test/src/math/CMakeLists.txt | 3 +
libc/test/src/math/exhaustive/CMakeLists.txt | 1 +
libc/test/src/math/exhaustive/tanf_test.cpp | 19 ++
libc/test/src/math/tanf_test.cpp | 60 +++++-
.../llvm-project-overlay/libc/BUILD.bazel | 7 +-
9 files changed, 407 insertions(+), 158 deletions(-)
create mode 100644 libc/src/__support/math/tanf_double_eval.h
create mode 100644 libc/src/__support/math/tanf_float_eval.h
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index baac96e76d5435..e339d3422c9ae0 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -5821,9 +5821,9 @@ add_header_library(
)
add_header_library(
- tanf
+ tanf_double_eval
HDRS
- tanf.h
+ tanf_double_eval.h
DEPENDS
.range_reduction
.sincosf_utils
@@ -5836,6 +5836,32 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ tanf_float_eval
+ HDRS
+ tanf_float_eval.h
+ DEPENDS
+ .sincosf_utils
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+)
+
+add_header_library(
+ tanf
+ HDRS
+ tanf.h
+ DEPENDS
+ .tanf_double_eval
+ .tanf_float_eval
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
+)
+
add_header_library(
tanf16
HDRS
diff --git a/libc/src/__support/math/tanf.h b/libc/src/__support/math/tanf.h
index 3a1d8e8bd53b62..9574156801189a 100644
--- a/libc/src/__support/math/tanf.h
+++ b/libc/src/__support/math/tanf.h
@@ -1,165 +1,42 @@
-//===-- Single-precision tan 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
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Implementation selector for tanf.
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANF_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_TANF_H
-#include "sincosf_utils.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/PolyEval.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/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+
+#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) && \
+ (defined(LIBC_MATH_HAS_SMALL_TABLES) || \
+ defined(LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT)) && \
+ defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
+#include "tanf_float_eval.h"
+#define LIBC_MATH_TANF_IMPL float_eval
+#else
+#include "tanf_double_eval.h"
+#define LIBC_MATH_TANF_IMPL double_eval
+#endif
namespace LIBC_NAMESPACE_DECL {
-
namespace math {
-namespace tanf_internal {
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-// Exceptional cases for tanf.
-LIBC_INLINE_VAR constexpr size_t N_EXCEPTS = 6;
-
-LIBC_INLINE_VAR constexpr fputil::ExceptValues<float, N_EXCEPTS> TANF_EXCEPTS{{
- // (inputs, RZ output, RU offset, RD offset, RN offset)
- // x = 0x1.ada6aap27, tan(x) = 0x1.e80304p-3 (RZ)
- {0x4d56d355, 0x3e740182, 1, 0, 0},
- // x = 0x1.862064p33, tan(x) = -0x1.8dee56p-3 (RZ)
- {0x50431032, 0xbe46f72b, 0, 1, 1},
- // x = 0x1.af61dap48, tan(x) = 0x1.60d1c6p-2 (RZ)
- {0x57d7b0ed, 0x3eb068e3, 1, 0, 1},
- // x = 0x1.0088bcp52, tan(x) = 0x1.ca1edp0 (RZ)
- {0x5980445e, 0x3fe50f68, 1, 0, 0},
- // x = 0x1.f90dfcp72, tan(x) = 0x1.597f9cp-1 (RZ)
- {0x63fc86fe, 0x3f2cbfce, 1, 0, 0},
- // x = 0x1.a6ce12p86, tan(x) = -0x1.c5612ep-1 (RZ)
- {0x6ad36709, 0xbf62b097, 0, 1, 0},
-}};
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-} // namespace tanf_internal
-
-LIBC_INLINE float tanf(float x) {
- using namespace sincosf_utils_internal;
- using namespace tanf_internal;
- using FPBits = typename fputil::FPBits<float>;
- FPBits xbits(x);
- uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
-
- // |x| < pi/32
- if (LIBC_UNLIKELY(x_abs <= 0x3dc9'0fdbU)) {
- double xd = static_cast<double>(x);
-
- // |x| < 0x1.0p-12f
- if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) {
- if (LIBC_UNLIKELY(x_abs == 0U)) {
- // For signed zeros.
- return x;
- }
- // When |x| < 2^-12, the relative error of the approximation tan(x) ~ x
- // is:
- // |tan(x) - x| / |tan(x)| < |x^3| / (3|x|)
- // = x^2 / 3
- // < 2^-25
- // < epsilon(1)/2.
- // So the correctly rounded values of tan(x) are:
- // = x + sign(x)*eps(x) if rounding mode = FE_UPWARD and x is positive,
- // or (rounding mode = FE_DOWNWARD and x is
- // negative),
- // = x otherwise.
- // To simplify the rounding decision and make it more efficient, we use
- // fma(x, 2^-25, x) instead.
- // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
- // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
- // |x| < 2^-125. For targets without FMA instructions, we simply use
- // double for intermediate results as it is more efficient than using an
- // emulated version of FMA.
-#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
- return fputil::multiply_add(x, 0x1.0p-25f, x);
-#else
- return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
-#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
- }
-
- // |x| < pi/32
- double xsq = xd * xd;
-
- // Degree-9 minimax odd polynomial of tan(x) generated by Sollya with:
- // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/32]);
- double result =
- fputil::polyeval(xsq, 1.0, 0x1.555555553d022p-2, 0x1.111111ce442c1p-3,
- 0x1.ba180a6bbdecdp-5, 0x1.69c0a88a0b71fp-6);
- return static_cast<float>(xd * result);
- }
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- bool x_sign = xbits.uintval() >> 31;
- // Check for exceptional values
- if (LIBC_UNLIKELY(x_abs == 0x3f8a1f62U)) {
- // |x| = 0x1.143ec4p0
- float sign = x_sign ? -1.0f : 1.0f;
-
- // volatile is used to prevent compiler (gcc) from optimizing the
- // computation, making the results incorrect in different rounding modes.
- volatile float tmp = 0x1.ddf9f4p0f;
- tmp = fputil::multiply_add(sign, tmp, sign * 0x1.1p-24f);
-
- return tmp;
- }
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
- // |x| > 0x1.ada6a8p+27f
- if (LIBC_UNLIKELY(x_abs > 0x4d56'd354U)) {
- // Inf or NaN
- if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
- if (xbits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
- }
-
- if (x_abs == 0x7f80'0000U) {
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- }
- return x + FPBits::quiet_nan().get_val();
- }
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- // Other large exceptional values
- if (auto r = TANF_EXCEPTS.lookup_odd(x_abs, x_sign);
- LIBC_UNLIKELY(r.has_value()))
- return r.value();
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- }
-
- // For |x| >= pi/32, we use the definition of tan(x) function:
- // tan(x) = sin(x) / cos(x)
- // The we follow the same computations of sin(x) and cos(x) as sinf, cosf,
- // and sincosf.
-
- double xd = static_cast<double>(x);
- double sin_k, cos_k, sin_y, cosm1_y;
-
- sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
- // tan(x) = sin(x) / cos(x)
- // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
- using fputil::multiply_add;
- return static_cast<float>(
- multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) /
- multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k)));
-}
+using LIBC_MATH_TANF_IMPL::tanf;
} // namespace math
-
} // namespace LIBC_NAMESPACE_DECL
+#undef LIBC_MATH_TANF_IMPL
+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANF_H
diff --git a/libc/src/__support/math/tanf_double_eval.h b/libc/src/__support/math/tanf_double_eval.h
new file mode 100644
index 00000000000000..42a8a7ee9a62a1
--- /dev/null
+++ b/libc/src/__support/math/tanf_double_eval.h
@@ -0,0 +1,172 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Double-precision implementation of tanf.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANF_DOUBLE_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_TANF_DOUBLE_EVAL_H
+
+#include "sincosf_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.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/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+namespace double_eval {
+
+namespace tanf_internal {
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Exceptional cases for tanf.
+LIBC_INLINE_VAR constexpr size_t N_EXCEPTS = 6;
+
+LIBC_INLINE_VAR constexpr fputil::ExceptValues<float, N_EXCEPTS> TANF_EXCEPTS{{
+ // (inputs, RZ output, RU offset, RD offset, RN offset)
+ // x = 0x1.ada6aap27, tan(x) = 0x1.e80304p-3 (RZ)
+ {0x4d56d355, 0x3e740182, 1, 0, 0},
+ // x = 0x1.862064p33, tan(x) = -0x1.8dee56p-3 (RZ)
+ {0x50431032, 0xbe46f72b, 0, 1, 1},
+ // x = 0x1.af61dap48, tan(x) = 0x1.60d1c6p-2 (RZ)
+ {0x57d7b0ed, 0x3eb068e3, 1, 0, 1},
+ // x = 0x1.0088bcp52, tan(x) = 0x1.ca1edp0 (RZ)
+ {0x5980445e, 0x3fe50f68, 1, 0, 0},
+ // x = 0x1.f90dfcp72, tan(x) = 0x1.597f9cp-1 (RZ)
+ {0x63fc86fe, 0x3f2cbfce, 1, 0, 0},
+ // x = 0x1.a6ce12p86, tan(x) = -0x1.c5612ep-1 (RZ)
+ {0x6ad36709, 0xbf62b097, 0, 1, 0},
+}};
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+} // namespace tanf_internal
+
+LIBC_INLINE float tanf(float x) {
+ using namespace sincosf_utils_internal;
+ using namespace tanf_internal;
+ using FPBits = typename fputil::FPBits<float>;
+ FPBits xbits(x);
+ uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
+
+ // |x| < pi/32
+ if (LIBC_UNLIKELY(x_abs <= 0x3dc9'0fdbU)) {
+ double xd = static_cast<double>(x);
+
+ // |x| < 0x1.0p-12f
+ if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) {
+ if (LIBC_UNLIKELY(x_abs == 0U)) {
+ // For signed zeros.
+ return x;
+ }
+ // When |x| < 2^-12, the relative error of the approximation tan(x) ~ x
+ // is:
+ // |tan(x) - x| / |tan(x)| < |x^3| / (3|x|)
+ // = x^2 / 3
+ // < 2^-25
+ // < epsilon(1)/2.
+ // So the correctly rounded values of tan(x) are:
+ // = x + sign(x)*eps(x) if rounding mode = FE_UPWARD and x is positive,
+ // or (rounding mode = FE_DOWNWARD and x is
+ // negative),
+ // = x otherwise.
+ // To simplify the rounding decision and make it more efficient, we use
+ // fma(x, 2^-25, x) instead.
+ // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
+ // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
+ // |x| < 2^-125. For targets without FMA instructions, we simply use
+ // double for intermediate results as it is more efficient than using an
+ // emulated version of FMA.
+#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
+ return fputil::multiply_add(x, 0x1.0p-25f, x);
+#else
+ return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
+#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
+ }
+
+ // |x| < pi/32
+ double xsq = xd * xd;
+
+ // Degree-9 minimax odd polynomial of tan(x) generated by Sollya with:
+ // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/32]);
+ double result =
+ fputil::polyeval(xsq, 1.0, 0x1.555555553d022p-2, 0x1.111111ce442c1p-3,
+ 0x1.ba180a6bbdecdp-5, 0x1.69c0a88a0b71fp-6);
+ return static_cast<float>(xd * result);
+ }
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ bool x_sign = xbits.uintval() >> 31;
+ // Check for exceptional values
+ if (LIBC_UNLIKELY(x_abs == 0x3f8a1f62U)) {
+ // |x| = 0x1.143ec4p0
+ float sign = x_sign ? -1.0f : 1.0f;
+
+ // volatile is used to prevent compiler (gcc) from optimizing the
+ // computation, making the results incorrect in different rounding modes.
+ volatile float tmp = 0x1.ddf9f4p0f;
+ tmp = fputil::multiply_add(sign, tmp, sign * 0x1.1p-24f);
+
+ return tmp;
+ }
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+ // |x| > 0x1.ada6a8p+27f
+ if (LIBC_UNLIKELY(x_abs > 0x4d56'd354U)) {
+ // Inf or NaN
+ if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ if (x_abs == 0x7f80'0000U) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ }
+ return x + FPBits::quiet_nan().get_val();
+ }
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // Other large exceptional values
+ if (auto r = TANF_EXCEPTS.lookup_odd(x_abs, x_sign);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ }
+
+ // For |x| >= pi/32, we use the definition of tan(x) function:
+ // tan(x) = sin(x) / cos(x)
+ // The we follow the same computations of sin(x) and cos(x) as sinf, cosf,
+ // and sincosf.
+
+ double xd = static_cast<double>(x);
+ double sin_k, cos_k, sin_y, cosm1_y;
+
+ sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
+ // tan(x) = sin(x) / cos(x)
+ // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
+ using fputil::multiply_add;
+ return static_cast<float>(
+ multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) /
+ multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k)));
+}
+
+} // namespace double_eval
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANF_DOUBLE_EVAL_H
diff --git a/libc/src/__support/math/tanf_float_eval.h b/libc/src/__support/math/tanf_float_eval.h
new file mode 100644
index 00000000000000..dadb30ace1b82f
--- /dev/null
+++ b/libc/src/__support/math/tanf_float_eval.h
@@ -0,0 +1,106 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Float-only implementation of tanf.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANF_FLOAT_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_TANF_FLOAT_EVAL_H
+
+#include "sincosf_float_eval.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace float_eval {
+
+// Float-only tanf, requiring FMA. The maximum error, checked exhaustively in
+// round-to-nearest mode against MPFR, is 3.37 ULP.
+LIBC_INLINE float tanf(float x) {
+ using namespace sincosf_float_eval;
+ using FPBits = fputil::FPBits<float>;
+ FPBits xbits(x);
+ uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
+
+ // |x| < 2^-12
+ if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) {
+ // For signed zeros.
+ if (LIBC_UNLIKELY(x_abs == 0U))
+ return x;
+ // When |x| < 2^-12, the relative error of the approximation tan(x) ~ x is:
+ // |tan(x) - x| / |tan(x)| < x^2 / 3 < 2^-25 < epsilon(1)/2.
+ // fma(x, 2^-25, x) gives the correctly rounded result in all rounding modes
+ // without underflowing for subnormal x.
+ return fputil::multiply_add(x, 0x1.0p-25f, x);
+ }
+
+ float y;
+ unsigned k;
+ if (x_abs < 0x4880'0000U) {
+ k = sincosf_range_reduction_small(x, y);
+ } else {
+ // Inf or NaN
+ if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ if (x_abs == 0x7f80'0000U) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ }
+ return x + FPBits::quiet_nan().get_val();
+ }
+
+ k = sincosf_range_reduction_large(x, y);
+ }
+
+ // x = k * pi/8 + y with |y| <= pi/16.
+ float sin_k = SIN_K_PI_OVER_8[k & 15];
+ // cos(k * pi/8) = sin((k + 4) * pi/8).
+ float cos_k = SIN_K_PI_OVER_8[(k + 4) & 15];
+
+ // Taylor polynomial of tan(y) for |y| <= pi/16:
+ // tan(y) ~ y + y^3/3 + 2y^5/15 + 17y^7/315 + 62y^9/2835.
+ // The coefficients of tan's Taylor series decrease by a factor of about
+ // (2/pi)^2 per term, so the relative truncation error is bounded by:
+ // 1.02 * 1382/155925 * (pi/16)^10 < 2^-30.
+ float y_sq = y * y;
+ float y3 = y_sq * y;
+ float p = fputil::multiply_add(y_sq, 0x1.664f48p-6f, 0x1.ba1ba2p-5f);
+ p = fputil::multiply_add(y_sq, p, 0x1.111112p-3f);
+ p = fputil::multiply_add(y_sq, p, 0x1.555556p-2f);
+ // t ~ tan(y)
+ float t = fputil::multiply_add(y3, p, y);
+
+ // tan(x) = sin(x) / cos(x)
+ // = (sin_k * cos(y) + cos_k * sin(y)) /
+ // (cos_k * cos(y) - sin_k * sin(y))
+ // = (sin_k + cos_k * tan(y)) / (cos_k - sin_k * tan(y)).
+ // Cancellation can only happen in the numerator when k = 0 (mod 8) and in
+ // the denominator when k = 4 (mod 8). In those cases, sin_k or cos_k is
+ // exactly 0 respectively, so the remaining term is computed with a small
+ // relative error.
+ float num = fputil::multiply_add(cos_k, t, sin_k);
+ float den = fputil::multiply_add(-sin_k, t, cos_k);
+ return num / den;
+}
+
+} // namespace float_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANF_FLOAT_EVAL_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 921dbf220eb372..bb3e568c2d4c9f 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -206,6 +206,9 @@ add_fp_unittest(
DEPENDS
libc.hdr.errno_macros
libc.src.math.tanf
+ libc.src.__support.math.tanf_double_eval
+ libc.src.__support.math.tanf_float_eval
+ libc.src.__support.macros.properties.cpu_features
libc.src.__support.CPP.array
libc.src.__support.FPUtil.fp_bits
)
diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt
index e1a91fca9b4c81..064a8a81548119 100644
--- a/libc/test/src/math/exhaustive/CMakeLists.txt
+++ b/libc/test/src/math/exhaustive/CMakeLists.txt
@@ -188,6 +188,7 @@ add_fp_unittest(
DEPENDS
.exhaustive_test
libc.src.math.tanf
+ libc.src.__support.math.tanf_float_eval
libc.src.__support.FPUtil.fp_bits
LINK_LIBRARIES
-lpthread
diff --git a/libc/test/src/math/exhaustive/tanf_test.cpp b/libc/test/src/math/exhaustive/tanf_test.cpp
index a3f606afb1c778..2954674ef0ede6 100644
--- a/libc/test/src/math/exhaustive/tanf_test.cpp
+++ b/libc/test/src/math/exhaustive/tanf_test.cpp
@@ -7,6 +7,7 @@
//===----------------------------------------------------------------------===//
#include "exhaustive_test.h"
+#include "src/__support/math/tanf_float_eval.h"
#include "src/math/tanf.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
@@ -31,3 +32,21 @@ static constexpr uint32_t NEG_STOP = 0xff80'0000U;
TEST_F(LlvmLibcTanfExhaustiveTest, NegativeRange) {
test_full_range_all_roundings(NEG_START, NEG_STOP);
}
+
+// Preserve the float implementation's 3.5 ULP MPFR bound in
+// round-to-nearest mode.
+static float tanf_float_eval(float x) {
+ return LIBC_NAMESPACE::math::float_eval::tanf(x);
+}
+
+using LlvmLibcTanfFloatExhaustiveTest =
+ LlvmLibcUnaryOpExhaustiveMathTest<float, mpfr::Operation::Tan,
+ tanf_float_eval, 3>;
+
+TEST_F(LlvmLibcTanfFloatExhaustiveTest, PositiveRange) {
+ test_full_range(mpfr::RoundingMode::Nearest, POS_START, POS_STOP);
+}
+
+TEST_F(LlvmLibcTanfFloatExhaustiveTest, NegativeRange) {
+ test_full_range(mpfr::RoundingMode::Nearest, 0x8000'0000U, NEG_STOP);
+}
diff --git a/libc/test/src/math/tanf_test.cpp b/libc/test/src/math/tanf_test.cpp
index fa04d870912e78..50873769b81d4b 100644
--- a/libc/test/src/math/tanf_test.cpp
+++ b/libc/test/src/math/tanf_test.cpp
@@ -10,6 +10,9 @@
#include "hdr/math_macros.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/tanf_double_eval.h"
+#include "src/__support/math/tanf_float_eval.h"
#include "src/math/tanf.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
@@ -17,19 +20,41 @@
#include "utils/MPFRWrapper/MPFRUtils.h"
#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-#define TOLERANCE 1
+#define TOLERANCE 3
+#define FP_ASSERT ASSERT_MPFR_MATCH
+#define FP_EXPECT EXPECT_MPFR_MATCH
#else
#define TOLERANCE 0
+#define FP_ASSERT ASSERT_MPFR_MATCH_ALL_ROUNDING
+#define FP_EXPECT EXPECT_MPFR_MATCH_ALL_ROUNDING
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
#include "hdr/stdint_proxy.h"
-using LlvmLibcTanfTest = LIBC_NAMESPACE::testing::FPTest<float>;
-
using LIBC_NAMESPACE::testing::SDCOMP26094_VALUES;
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+class LlvmLibcTanfTest : public LIBC_NAMESPACE::testing::FPTest<float> {
+public:
+ void test_eval_in_float_range(float (*func)(float), double tolerance,
+ bool all_rounding) {
+ constexpr uint32_t COUNT = 1'231;
+ constexpr uint32_t STEP = UINT32_MAX / COUNT;
+ for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+ float x = FPBits(v).get_val();
+ if (FPBits(v).is_nan() || FPBits(v).is_inf())
+ continue;
+ if (all_rounding) {
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, x, func(x),
+ tolerance);
+ } else {
+ ASSERT_MPFR_MATCH(mpfr::Operation::Tan, x, func(x), tolerance);
+ }
+ }
+ }
+};
+
TEST_F(LlvmLibcTanfTest, SpecialNumbers) {
EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::tanf(aNaN));
EXPECT_MATH_ERRNO(0);
@@ -54,8 +79,8 @@ TEST_F(LlvmLibcTanfTest, InFloatRange) {
float x = FPBits(v).get_val();
if (FPBits(v).is_nan() || FPBits(v).is_inf())
continue;
- ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, x,
- LIBC_NAMESPACE::tanf(x), 0.5);
+ FP_ASSERT(mpfr::Operation::Tan, x, LIBC_NAMESPACE::tanf(x),
+ TOLERANCE + 0.5);
}
}
@@ -120,10 +145,10 @@ TEST_F(LlvmLibcTanfTest, SpecificBitPatterns) {
for (int i = 0; i < N; ++i) {
float x = FPBits(INPUTS[i]).get_val();
- EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, x,
- LIBC_NAMESPACE::tanf(x), TOLERANCE + 0.5);
- EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, -x,
- LIBC_NAMESPACE::tanf(-x), TOLERANCE + 0.5);
+ FP_EXPECT(mpfr::Operation::Tan, x, LIBC_NAMESPACE::tanf(x),
+ TOLERANCE + 0.5);
+ FP_EXPECT(mpfr::Operation::Tan, -x, LIBC_NAMESPACE::tanf(-x),
+ TOLERANCE + 0.5);
}
}
@@ -132,6 +157,21 @@ TEST_F(LlvmLibcTanfTest, SpecificBitPatterns) {
TEST_F(LlvmLibcTanfTest, SDCOMP_26094) {
for (uint32_t v : SDCOMP26094_VALUES) {
float x = FPBits(v).get_val();
- ASSERT_MPFR_MATCH(mpfr::Operation::Tan, x, LIBC_NAMESPACE::tanf(x), 0.5);
+ ASSERT_MPFR_MATCH(mpfr::Operation::Tan, x, LIBC_NAMESPACE::tanf(x),
+ TOLERANCE + 0.5);
}
}
+
+// Exercise both implementations independently of the configured selector.
+TEST_F(LlvmLibcTanfTest, DoubleEvalInFloatRange) {
+ test_eval_in_float_range(&LIBC_NAMESPACE::math::double_eval::tanf,
+ TOLERANCE + 0.5,
+ /*all_rounding=*/TOLERANCE == 0);
+}
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
+TEST_F(LlvmLibcTanfTest, FloatEvalInFloatRange) {
+ test_eval_in_float_range(&LIBC_NAMESPACE::math::float_eval::tanf, 3.5,
+ /*all_rounding=*/false);
+}
+#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index f2df1ef9466183..1e6ede0039b740 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -11475,8 +11475,13 @@ libc_support_library(
libc_support_library(
name = "__support_math_tanf",
- hdrs = ["src/__support/math/tanf.h"],
+ hdrs = [
+ "src/__support/math/tanf.h",
+ "src/__support/math/tanf_double_eval.h",
+ "src/__support/math/tanf_float_eval.h",
+ ],
deps = [
+ ":__support_common",
":__support_fputil_except_value_utils",
":__support_fputil_fenv_impl",
":__support_fputil_fma",
More information about the libc-commits
mailing list