[libc][math] Added atanf function.

Performance by core-math (core-math/glibc 2.31/current llvm-14):
28.879/20.843/20.15

Differential Revision: https://reviews.llvm.org/D132842

libc/config/darwin/arm/entrypoints.txt

@@ -105,6 +105,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanf
     libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf

libc/config/linux/aarch64/entrypoints.txt

@@ -149,6 +149,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanf
     libc.src.math.atanhf    
     libc.src.math.copysign
     libc.src.math.copysignf

libc/config/linux/x86_64/entrypoints.txt

@@ -149,7 +149,8 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
-    libc.src.math.atanhf    
+    libc.src.math.atanf
+    libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl

libc/config/windows/entrypoints.txt

@@ -106,6 +106,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanf    
     libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf

libc/spec/stdc.td

@@ -480,7 +480,9 @@ def StdC : StandardSpec<"stdc"> {
           FunctionSpec<"sinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
-          FunctionSpec<"atanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,          
+          FunctionSpec<"atanf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+
+          FunctionSpec<"atanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
 

libc/src/__support/FPUtil/FPBits.h

@@ -196,6 +196,15 @@ template <typename T> struct FPBits {
     }
     return result;
   }
+
+  inline static FPBits<T> create_value(bool sign, UIntType unbiased_exp,
+                                       UIntType mantissa) {
+    FPBits<T> result;
+    result.set_sign(sign);
+    result.set_unbiased_exponent(unbiased_exp);
+    result.set_mantissa(mantissa);
+    return result;
+  }
 };
 
 } // namespace fputil

libc/src/__support/FPUtil/x86_64/LongDoubleBits.h

@@ -184,6 +184,15 @@ template <> struct FPBits<long double> {
     bits.set_mantissa(v);
     return bits;
   }
+
+  inline static FPBits<long double>
+  create_value(bool sign, UIntType unbiased_exp, UIntType mantissa) {
+    FPBits<long double> result;
+    result.set_sign(sign);
+    result.set_unbiased_exponent(unbiased_exp);
+    result.set_mantissa(mantissa);
+    return result;
+  }
 };
 
 static_assert(

libc/src/math/CMakeLists.txt

@@ -66,6 +66,7 @@ add_entrypoint_object(
     -O3
 )
 
+add_math_entrypoint_object(atanf)
 add_math_entrypoint_object(atanhf)
 
 add_math_entrypoint_object(ceil)

libc/src/math/atanf.h

@@ -0,0 +1,18 @@
+//===-- Implementation header for atanf -------------------------*- 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_ATANF_H
+#define LLVM_LIBC_SRC_MATH_ATANF_H
+
+namespace __llvm_libc {
+
+float atanf(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_ATANF_H

libc/src/math/generic/CMakeLists.txt

@@ -1242,3 +1242,28 @@ add_entrypoint_object(
     -O3
 )
 
+add_object_library(
+  inv_trigf_utils
+  HDRS
+    inv_trigf_utils.h
+  SRCS
+    inv_trigf_utils.cpp
+)
+
+add_entrypoint_object(
+  atanf
+  SRCS
+    atanf.cpp
+  HDRS
+    ../atanf.h
+  DEPENDS
+    .inv_trigf_utils
+    libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.nearest_integer
+    libc.src.__support.FPUtil.polyeval
+    libc.include.math
+  COMPILE_OPTIONS
+    -O3
+)
+

libc/src/math/generic/atanf.cpp

@@ -0,0 +1,57 @@
+//===-- Single-precision atan 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/atanf.h"
+#include "math_utils.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/generic/inv_trigf_utils.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+  bool sign = xbits.get_sign();
+  xbits.set_sign(false);
+
+  if (unlikely(xbits.is_inf_or_nan())) {
+    if (xbits.is_inf())
+      return opt_barrier(sign ? -M_MATH_PI_2 : M_MATH_PI_2);
+    else
+      return x + 1.0f;
+  }
+  // |x| == 0.06905200332403183
+  if (unlikely(xbits.uintval() == 0x3d8d6b23U)) {
+    if (fputil::get_round() == FE_TONEAREST) {
+      // 0.06894256919622421
+      FPBits br(0x3d8d31c3U);
+      br.set_sign(sign);
+      return br.get_val();
+    }
+  }
+
+  // |x| == 1.8670953512191772
+  if (unlikely(xbits.uintval() == 0x3feefcfbU)) {
+    int rounding_mode = fputil::get_round();
+    if (sign) {
+      if (rounding_mode == FE_DOWNWARD) {
+        // -1.0790828466415405
+        return FPBits(0xbf8a1f63U).get_val();
+      }
+    } else {
+      if (rounding_mode == FE_UPWARD) {
+        // 1.0790828466415405
+        return FPBits(0x3f8a1f63U).get_val();
+      }
+    }
+  }
+
+  return atan_eval(x);
+}
+
+} // namespace __llvm_libc

libc/src/math/generic/explogxf.h

@@ -9,7 +9,6 @@
 #ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H
 #define LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H
 
-#include "common_constants.h" // Lookup tables EXP_M
 #include "math_utils.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"

libc/src/math/generic/inv_trigf_utils.cpp

@@ -0,0 +1,28 @@
+//===-- Single-precision general exp/log functions ------------------------===//
+//
+// 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 "inv_trigf_utils.h"
+
+namespace __llvm_libc {
+
+// N[Table[ArcTan[x], {x, 1/16, 16/16, 1/16}], 40]
+alignas(64) const double ATAN_T[ATAN_T_SIZE] = {
+    0x1.ff55bb72cfdeap-5, 0x1.fd5ba9aac2f6ep-4, 0x1.7b97b4bce5b02p-3,
+    0x1.f5b75f92c80ddp-3, 0x1.362773707ebccp-2, 0x1.6f61941e4def1p-2,
+    0x1.a64eec3cc23fdp-2, 0x1.dac670561bb4fp-2, 0x1.0657e94db30d0p-1,
+    0x1.1e00babdefeb4p-1, 0x1.345f01cce37bbp-1, 0x1.4978fa3269ee1p-1,
+    0x1.5d58987169b18p-1, 0x1.700a7c5784634p-1, 0x1.819d0b7158a4dp-1,
+    0x1.921fb54442d18p-1};
+
+// for(int i = 0; i < 5; i++)
+//     printf("%.13a,\n", (-2 * (i % 2) + 1) * 1.0 / (2 * i + 1));
+alignas(64) const double ATAN_K[5] = {
+    0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,
+    -0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4};
+
+} // namespace __llvm_libc

libc/src/math/generic/inv_trigf_utils.h

@@ -0,0 +1,94 @@
+//===-- Single-precision general inverse trigonometric functions ----------===//
+//
+// 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_GENERIC_INV_TRIGF_UTILS_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H
+
+#include "math_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include <src/__support/FPUtil/NearestIntegerOperations.h>
+
+#include <errno.h>
+
+namespace __llvm_libc {
+
+// PI / 2
+constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
+
+// atan table size
+constexpr int ATAN_T_BITS = 4;
+constexpr int ATAN_T_SIZE = 1 << ATAN_T_BITS;
+
+// N[Table[ArcTan[x], {x, 1/8, 8/8, 1/8}], 40]
+extern const double ATAN_T[ATAN_T_SIZE];
+extern const double ATAN_K[5];
+
+// The main idea of the function is to use formula
+// atan(u) + atan(v) = atan((u+v)/(1-uv))
+
+// x should be positive, normal finite value
+inline static double atan_eval(double x) {
+  using FPB = fputil::FPBits<double>;
+  // Added some small value to umin and umax mantissa to avoid possible rounding
+  // errors.
+  FPB::UIntType umin =
+      FPB::create_value(false, FPB::EXPONENT_BIAS - ATAN_T_BITS - 1,
+                        0x100000000000UL)
+          .uintval();
+  FPB::UIntType umax =
+      FPB::create_value(false, FPB::EXPONENT_BIAS + ATAN_T_BITS,
+                        0xF000000000000UL)
+          .uintval();
+
+  FPB bs(x);
+  bool sign = bs.get_sign();
+  auto x_abs = bs.uintval() & FPB::FloatProp::EXP_MANT_MASK;
+
+  if (x_abs <= umin) {
+    double pe = __llvm_libc::fputil::polyeval(x * x, 0.0, ATAN_K[1], ATAN_K[2],
+                                              ATAN_K[3], ATAN_K[4]);
+    return fputil::multiply_add(pe, x, x);
+  }
+
+  if (x_abs >= umax) {
+    double one_over_x_m = -1.0 / x;
+    double one_over_x2 = one_over_x_m * one_over_x_m;
+    double pe = __llvm_libc::fputil::polyeval(one_over_x2, ATAN_K[0], ATAN_K[1],
+                                              ATAN_K[2], ATAN_K[3]);
+    return fputil::multiply_add(pe, one_over_x_m, sign ? (-M_MATH_PI_2) : (M_MATH_PI_2));
+  }
+
+  double pos_x = FPB(x_abs).get_val();
+  bool one_over_x = pos_x > 1.0;
+  if (one_over_x) {
+    pos_x = 1.0 / pos_x;
+  }
+
+  double near_x = fputil::nearest_integer(pos_x * ATAN_T_SIZE);
+  int val = static_cast<int>(near_x);
+  near_x *= 1.0 / ATAN_T_SIZE;
+
+  double v = (pos_x - near_x) / fputil::multiply_add(near_x, pos_x, 1.0);
+  double v2 = v * v;
+  double pe = __llvm_libc::fputil::polyeval(v2, ATAN_K[0], ATAN_K[1], ATAN_K[2],
+                                            ATAN_K[3], ATAN_K[4]);
+  double result;
+  if (one_over_x)
+    result = M_PI_2 - fputil::multiply_add(pe, v, ATAN_T[val - 1]);
+  else
+    result = fputil::multiply_add(pe, v, ATAN_T[val - 1]);
+  return sign ? -result : result;
+}
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H

libc/test/src/math/CMakeLists.txt

@@ -1388,6 +1388,32 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  atanf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    atanf_test.cpp
+  DEPENDS
+    libc.src.math.atanf
+    libc.src.__support.FPUtil.fputil
+)
+
+add_fp_unittest(
+  inv_trigf_utils_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  HDRS
+    in_float_range_test_helper.h
+  SRCS
+    inv_trigf_utils_test.cpp
+  DEPENDS
+    libc.src.math.generic.inv_trigf_utils
+    libc.src.__support.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)