diff --git a/libc/shared/math.h b/libc/shared/math.h index a331abcc033c7..5126e46c9772e 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -85,5 +85,6 @@ #include "math/rsqrtf16.h" #include "math/sin.h" #include "math/tan.h" +#include "math/tanf.h" #endif // LLVM_LIBC_SHARED_MATH_H diff --git a/libc/shared/math/tanf.h b/libc/shared/math/tanf.h new file mode 100644 index 0000000000000..ce0f11db1ee6a --- /dev/null +++ b/libc/shared/math/tanf.h @@ -0,0 +1,23 @@ +//===-- Shared tanf function ------------------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SHARED_MATH_TANF_H +#define LLVM_LIBC_SHARED_MATH_TANF_H + +#include "shared/libc_common.h" +#include "src/__support/math/tanf.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::tanf; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SHARED_MATH_TANF_H diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index a03952d5c5ed0..d18ebc03ba7b8 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -1270,3 +1270,19 @@ add_header_library( libc.src.__support.FPUtil.multiply_add libc.src.__support.macros.optimization ) + +add_header_library( + tanf + HDRS + tanf.h + DEPENDS + .range_reduction + .sincosf_utils + libc.src.__support.FPUtil.except_value_utils + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fma + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.nearest_integer + libc.src.__support.FPUtil.polyeval + libc.src.__support.macros.optimization +) diff --git a/libc/src/__support/math/tanf.h b/libc/src/__support/math/tanf.h new file mode 100644 index 0000000000000..980bb98e1cb0e --- /dev/null +++ b/libc/src/__support/math/tanf.h @@ -0,0 +1,164 @@ +//===-- 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 +// +//===----------------------------------------------------------------------===// + +#ifndef LIBC_SRC___SUPPORT_MATH_TANF_H +#define 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 + +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 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 static float tanf(float x) { + using namespace tanf_internal; + using FPBits = typename fputil::FPBits; + 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(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(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(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(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( + 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 math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LIBC_SRC___SUPPORT_MATH_TANF_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 04ad1aa896613..3addfab6c6603 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -496,17 +496,8 @@ add_entrypoint_object( HDRS ../tanf.h DEPENDS - libc.src.__support.math.range_reduction - libc.src.__support.math.sincosf_utils + libc.src.__support.math.tanf libc.src.errno.errno - libc.src.__support.FPUtil.fenv_impl - libc.src.__support.FPUtil.fenv_impl - libc.src.__support.FPUtil.except_value_utils - libc.src.__support.FPUtil.nearest_integer - libc.src.__support.FPUtil.fma - libc.src.__support.FPUtil.multiply_add - libc.src.__support.FPUtil.polyeval - libc.src.__support.macros.optimization ) add_entrypoint_object( diff --git a/libc/src/math/generic/tanf.cpp b/libc/src/math/generic/tanf.cpp index a8c557b14808c..aee5a6cab2b56 100644 --- a/libc/src/math/generic/tanf.cpp +++ b/libc/src/math/generic/tanf.cpp @@ -7,146 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/tanf.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/common.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/math/sincosf_utils.h" +#include "src/__support/math/tanf.h" namespace LIBC_NAMESPACE_DECL { -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS -// Exceptional cases for tanf. -constexpr size_t N_EXCEPTS = 6; - -constexpr fputil::ExceptValues 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 - -LLVM_LIBC_FUNCTION(float, tanf, (float x)) { - using FPBits = typename fputil::FPBits; - 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(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(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(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(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( - 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))); -} +LLVM_LIBC_FUNCTION(float, tanf, (float x)) { return math::tanf(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt index 4f82d00cc6c84..dabf5f6b168cc 100644 --- a/libc/test/shared/CMakeLists.txt +++ b/libc/test/shared/CMakeLists.txt @@ -81,4 +81,5 @@ add_fp_unittest( libc.src.__support.math.rsqrtf16 libc.src.__support.math.sin libc.src.__support.math.tan + libc.src.__support.math.tanf ) diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp index 5ec8a7b23081e..cd389c41cb7b2 100644 --- a/libc/test/shared/shared_math_test.cpp +++ b/libc/test/shared/shared_math_test.cpp @@ -89,6 +89,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat) { EXPECT_EQ(long(0), LIBC_NAMESPACE::shared::llogbf(1.0f)); EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::logbf(1.0f)); EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::rsqrtf(1.0f)); + EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::shared::tanf(0.0f)); } TEST(LlvmLibcSharedMathTest, AllDouble) { diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 853182ce18b57..346864d2276c0 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -3492,6 +3492,22 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_tanf", + hdrs = ["src/__support/math/tanf.h"], + deps = [ + ":__support_fputil_except_value_utils", + ":__support_fputil_fma", + ":__support_fputil_multiply_add", + ":__support_fputil_nearest_integer", + ":__support_fputil_polyeval", + ":__support_macros_optimization", + ":__support_macros_properties_cpu_features", + ":__support_range_reduction", + ":__support_sincosf_utils", + ], +) + ############################### complex targets ################################ libc_function( @@ -5166,14 +5182,7 @@ libc_math_function( libc_math_function( name = "tanf", additional_deps = [ - ":__support_fputil_fma", - ":__support_fputil_multiply_add", - ":__support_fputil_nearest_integer", - ":__support_fputil_polyeval", - ":__support_macros_optimization", - ":__support_macros_properties_cpu_features", - ":__support_range_reduction", - ":__support_sincosf_utils", + ":__support_math_tanf", ], )