[libc][math][c23] Add tanf16 function (#121018)

- Implementation of tan for 16-bit floating point inputs.
- Exhaustive tests across the 16-bit input range

Author: wldfngrs

libc/config/linux/x86_64/entrypoints.txt

@@ -721,6 +721,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.sinhf16
     libc.src.math.sinpif16
     libc.src.math.sqrtf16
+    libc.src.math.tanf16
     libc.src.math.tanhf16
     libc.src.math.tanpif16
     libc.src.math.totalorderf16

libc/docs/headers/math/index.rst

@@ -346,7 +346,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | sqrt      | |check|          | |check|         | |check|                | |check|              | |check|                | 7.12.7.10              | F.10.4.10                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| tan       | |check|          | |check|         |                        |                      |                        | 7.12.4.7               | F.10.1.7                   |
+| tan       | |check|          | |check|         |                        | |check|              |                        | 7.12.4.7               | F.10.1.7                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | tanh      | |check|          |                 |                        | |check|              |                        | 7.12.5.6               | F.10.2.6                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

libc/include/math.yaml

@@ -2418,6 +2418,13 @@ functions:
     return_type: float
     arguments:
       - type: float
+  - name: tanf16
+    standards:
+      - stdc
+    return_type: _Float16
+    arguments:
+      - type: _Float16
+    guard: LIBC_TYPES_HAS_FLOAT16
   - name: tanhf
     standards:
       - stdc

libc/src/math/CMakeLists.txt

@@ -501,6 +501,7 @@ add_math_entrypoint_object(sqrtf128)
 
 add_math_entrypoint_object(tan)
 add_math_entrypoint_object(tanf)
+add_math_entrypoint_object(tanf16)
 
 add_math_entrypoint_object(tanh)
 add_math_entrypoint_object(tanhf)

libc/src/math/generic/CMakeLists.txt

@@ -661,6 +661,25 @@ add_entrypoint_object(
     ${libc_opt_high_flag}
 )
 
+add_entrypoint_object(
+  tanf16
+  SRCS
+    tanf16.cpp
+  HDRS
+    ../tanf16.h
+  DEPENDS
+    .sincosf16_utils
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.cast
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.except_value_utils
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.macros.optimization
+    libc.src.__support.macros.properties.types
+)
+
 add_entrypoint_object(
   tanpif16
   SRCS

libc/src/math/generic/sincosf16_utils.h

@@ -47,24 +47,23 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
 
 // Recall, range reduction:
 //   k = round(x * 32/pi)
-//   y = x * 32/pi - k
 //
-// The constant 0x1.45f306dc9c883p3 is 32/pi rounded to double-precision.
-// 32/pi is generated by Sollya with the following commands:
-// > display = hexadecimal;
-// > round(32/pi, D, RN);
-//
-// The precision choice of 'double' is to minimize rounding errors
-// in this initial scaling step, preserving enough bits so errors accumulated
-// while computing the subtraction: y = x * 32/pi - round(x * 32/pi)
+// The precision choice of 'double' in the following function is to minimize
+// rounding errors in this initial scaling step,
+// preserving enough bits so errors accumulated while computing the subtraction:
+// y = x * 32/pi - round(x * 32/pi)
 // are beyond the least-significant bit of single-precision used during
 // further intermediate computation.
 LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) {
-  double prod = x * 0x1.45f306dc9c883p3;
-  double kf = fputil::nearest_integer(prod);
-  y = static_cast<float>(prod - kf);
+  // Generated by Sollya with:
+  // > D(32/pi);
+  constexpr double THIRTYTWO_OVER_PI = 0x1.45f306dc9c883p3;
 
-  return static_cast<int32_t>(kf);
+  double prod = x * THIRTYTWO_OVER_PI;
+  double kd = fputil::nearest_integer(prod);
+  y = static_cast<float>(prod - kd);
+
+  return static_cast<int32_t>(kd);
 }
 
 static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k,

libc/src/math/generic/tanf16.cpp

@@ -0,0 +1,115 @@
+//===-- Half-precision tan(x) 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/tanf16.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "sincosf16_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+constexpr size_t N_EXCEPTS = 9;
+
+constexpr fputil::ExceptValues<float16, N_EXCEPTS> TANF16_EXCEPTS{{
+    // (input, RZ output, RU offset, RD offset, RN offset)
+    {0x2894, 0x2894, 1, 0, 1},
+    {0x3091, 0x3099, 1, 0, 0},
+    {0x3098, 0x30a0, 1, 0, 0},
+    {0x55ed, 0x3911, 1, 0, 0},
+    {0x607b, 0xc638, 0, 1, 1},
+    {0x674e, 0x3b7d, 1, 0, 0},
+    {0x6807, 0x4014, 1, 0, 1},
+    {0x6f4d, 0xbe19, 0, 1, 1},
+    {0x7330, 0xcb62, 0, 1, 0},
+}};
+
+LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits xbits(x);
+
+  uint16_t x_u = xbits.uintval();
+  uint16_t x_abs = x_u & 0x7fff;
+  bool x_sign = x_u >> 15;
+  float xf = x;
+
+  // Handle exceptional values
+  if (auto r = TANF16_EXCEPTS.lookup_odd(x_abs, x_sign);
+      LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+
+  // |x| <= 0x1.d1p-5
+  if (LIBC_UNLIKELY(x_abs <= 0x2b44)) {
+    // |x| <= 0x1.398p-11
+    if (LIBC_UNLIKELY(x_abs <= 0x10e6)) {
+      // tan(+/-0) = +/-0
+      if (LIBC_UNLIKELY(x_abs == 0))
+        return x;
+
+      int rounding = fputil::quick_get_round();
+
+      // Exhaustive tests show that, when:
+      // x > 0, and rounding upward or
+      // x < 0, and rounding downward then,
+      // tan(x) = x * 2^-11 + x
+      if ((xbits.is_pos() && rounding == FE_UPWARD) ||
+          (xbits.is_neg() && rounding == FE_DOWNWARD))
+        return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
+      return x;
+    }
+
+    float xsq = xf * xf;
+
+    // Degree-6 minimax odd polynomial of tan(x) generated by Sollya with:
+    // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6|], [|1, SG...|], [0, pi/32]);
+    float result = fputil::polyeval(xsq, 0x1p0f, 0x1.555556p-2f, 0x1.110ee4p-3f,
+                                    0x1.be80f6p-5f);
+
+    return fputil::cast<float16>(xf * result);
+  }
+
+  // tan(+/-inf) = NaN, and tan(NaN) = NaN
+  if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
+    // x = +/-inf
+    if (x_abs == 0x7c00) {
+      fputil::set_errno_if_required(EDOM);
+      fputil::raise_except_if_required(FE_INVALID);
+    }
+
+    return x + FPBits::quiet_nan().get_val();
+  }
+
+  // Range reduction:
+  // For |x| > pi/32, we perform range reduction as follows:
+  // Find k and y such that:
+  //   x = (k + y) * pi/32;
+  //   k is an integer, |y| < 0.5
+  //
+  // This is done by performing:
+  //   k = round(x * 32/pi)
+  //   y = x * 32/pi - k
+  //
+  // Once k and y are computed, we then deduce the answer by the formula:
+  // tan(x) = sin(x) / cos(x)
+  // 	    = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
+  float sin_k, cos_k, sin_y, cosm1_y;
+  sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
+
+  // Note that, cosm1_y = cos_y - 1:
+  using fputil::multiply_add;
+  return fputil::cast<float16>(
+      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 LIBC_NAMESPACE_DECL

libc/src/math/generic/tanpif16.cpp

@@ -79,7 +79,7 @@ LLVM_LIBC_FUNCTION(float16, tanpif16, (float16 x)) {
   //   k = round(x * 32)
   //   y = x * 32 - k
   //
-  // Once k and y are computed, we then deduce the answer by tthe formula:
+  // Once k and y are computed, we then deduce the answer by the formula:
   // tan(x) = sin(x) / cos(x)
   //        = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
   float xf = x;

libc/src/math/tanf16.h

@@ -0,0 +1,21 @@
+//===-- Implementation header for tanf16 ------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_TANF16_H
+#define LLVM_LIBC_SRC_MATH_TANF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 tanf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_TANF16_H

libc/test/src/math/CMakeLists.txt

@@ -190,6 +190,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  tanf16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    tanf16_test.cpp
+  DEPENDS
+    libc.src.math.tanf16  
+)
+
 add_fp_unittest(
   tanpif16_test
   NEED_MPFR

libc/test/src/math/cosf16_test.cpp

@@ -17,7 +17,7 @@ namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
 
 // Range: [0, Inf]
 static constexpr uint16_t POS_START = 0x0000U;
-static constexpr uint16_t POS_STOP = 0x7c00u;
+static constexpr uint16_t POS_STOP = 0x7c00U;
 
 // Range: [-Inf, 0]
 static constexpr uint16_t NEG_START = 0x8000U;

libc/test/src/math/smoke/CMakeLists.txt

@@ -121,6 +121,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  tanf16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    tanf16_test.cpp
+  DEPENDS
+    libc.src.errno.errno
+    libc.src.math.tanf16 
+)
+
 add_fp_unittest(
   tanpif16_test
   SUITE

libc/test/src/math/smoke/tanf16_test.cpp

@@ -0,0 +1,34 @@
+//===-- Unittests for tanf16 ----------------------------------------------===//
+//
+//
+// 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/errno/libc_errno.h"
+#include "src/math/tanf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcTanf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcTanf16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::tanf16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(zero, LIBC_NAMESPACE::tanf16(zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::tanf16(neg_zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::tanf16(inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::tanf16(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+}

libc/test/src/math/tanf16_test.cpp

@@ -0,0 +1,40 @@
+//===-- Exhaustive test for tanf16 ----------------------------------------===//
+//
+// 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/tanf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcTanf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf]
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0]
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcTanf16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, x,
+                                   LIBC_NAMESPACE::tanf16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcTanf16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, x,
+                                   LIBC_NAMESPACE::tanf16(x), 0.5);
+  }
+}