[libc][math] Added auxiliary function log2_eval for asinhf/acoshf/atanhf.

[orex]

Summary

1. double log2_eval(double) function added with better than float precision is added.
2. Some refactoring done to put all auxiliary functions and corresponding data

to one place to reuse the code.

3. Added tests for new functions.
4. Performance and precision tests of the function shows, that it more precise than exiting log2,

(no exceptional cases), but timing is ~5% higer that on current one.

[lntue]

libc/src/math/generic/supfuncf.cpp

28

Add comments about how these constants are defined and generated. The same for other tables below.

libc/src/math/generic/supfuncf.h

9–10

explogxf.h is a better name for this file and cmake target.

124

Add comments about how the coefficients are generated.

151

Add comment about what the constant is and how it is generated.

158

Add comment about what the constant is and how it is generated.

[orex]

Constants description added.

[orex]

Fix wrong upload.

[orex]

libc/src/math/generic/supfuncf.h

9–10

Please check also D132811 I've added to the file also atan_eval function, so therefore it is named such. If you have any other suggestion about the file name I'll appreciate your suggestion.

[lntue]

libc/src/math/generic/supfuncf.h

9–10

The constants and atan_eval function that you added in https://reviews.llvm.org/D132842 are independent with these exp and log utilities. Maybe you can put them into separate files / target called inv_trig_utils or inv_trig_helpers, as I expect asin and acos will also need some helpers, which could be put in there.

[orex]

Review fixes.

[orex]

The aux functions are split into two files according @lntue suggestion.

libc/config/darwin/arm/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanf
+    libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.coshf
     libc.src.math.cosf

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanf
+    libc.src.math.atanhf    
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.coshf
     libc.src.math.cosf

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanf
+    libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cos
     libc.src.math.coshf

libc/config/windows/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.fenv.fesetenv
     libc.src.fenv.fesetexceptflag
     libc.src.fenv.fesetround
     libc.src.fenv.feraiseexcept
     libc.src.fenv.fetestexcept
     libc.src.fenv.feupdateenv
 
     # math.h entrypoints
+    libc.src.math.atanf    
+    libc.src.math.atanhf
     libc.src.math.copysign
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cos
     libc.src.math.cosf

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
 
           FunctionSpec<"nextafterf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"nextafter", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"nextafterl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"coshf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"sinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          
+          FunctionSpec<"atanf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          
+          FunctionSpec<"atanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
 
   HeaderSpec StdIO = HeaderSpec<
       "stdio.h",
       [
           Macro<"stderr">,
           Macro<"stdout">,

libc/src/__support/FPUtil/FPBits.h

if (likely(ep >= 0)) {
       // Implicit number bit will be removed by mask
       result.set_mantissa(number);
       result.set_unbiased_exponent(ep + 1);
     } else {
       result.set_mantissa(number >> -ep);
     }
     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
 } // namespace __llvm_libc
 
 #ifdef SPECIAL_X86_LONG_DOUBLE
 #include "x86_64/LongDoubleBits.h"
 #endif
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_FP_BITS_H

libc/src/math/CMakeLists.txt

HDRS
     fma.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.fma
   COMPILE_OPTIONS
     -O3
 )
 
+add_math_entrypoint_object(atanf)
+add_math_entrypoint_object(atanhf)
+
 add_math_entrypoint_object(ceil)
 add_math_entrypoint_object(ceilf)
 add_math_entrypoint_object(ceill)
 
 add_math_entrypoint_object(copysign)
 add_math_entrypoint_object(copysignf)
 add_math_entrypoint_object(copysignl)
 

libc/src/math/atanf.h

+//===-- 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/atanhf.h

+//===-- Implementation header for atanhf ------------------------*- 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_ATANHF_H
+#define LLVM_LIBC_SRC_MATH_ATANHF_H
+
+namespace __llvm_libc {
+
+float atanhf(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_ATANHF_H

libc/src/math/generic/CMakeLists.txt

DEPENDS
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
+add_entrypoint_object(
+  atanhf
+  SRCS
+    atanhf.cpp
+  HDRS
+    ../atanhf.h
+  DEPENDS
+    .supfuncf
+    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
+)
+
+add_entrypoint_object(
+  atanf
+  SRCS
+    atanf.cpp
+  HDRS
+    ../atanf.h
+  DEPENDS
+    .supfuncf
+    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

+//===-- 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/supfuncf.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_PI_2 : M_PI_2);
+    else
+      return x + 1.0f;
+  }
+
+  if (unlikely(xbits.uintval() == 0x3d8d6b23U)) {
+    if (fputil::get_round() == FE_TONEAREST) {
+      FPBits br(0x3d8d31c3U);
+      br.set_sign(sign);
+      return br.get_val();
+    }
+  }
+
+  if (unlikely(xbits.uintval() == 0x3feefcfbU)) {
+    int rounding_mode = fputil::get_round();
+    if (sign) {
+      if (rounding_mode == FE_DOWNWARD)
+        return FPBits(0xbf8a1f63U).get_val();
+    } else {
+      if (rounding_mode == FE_UPWARD)
+        return FPBits(0x3f8a1f63U).get_val();
+    }
+  }
+
+  return atan_eval(x);
+}
+
+} // namespace __llvm_libc

libc/src/math/generic/atanhf.cpp

+//===-- Single-precision atanh 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/atanhf.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/generic/supfuncf.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(float, atanhf, (float x)) {
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+  bool sign = xbits.get_sign();
+  uint32_t x_abs = xbits.uintval() & FPBits::FloatProp::EXP_MANT_MASK;
+
+  // 0 >= 1.0
+  if (unlikely(x_abs >= 0x3F80'0000U)) {
+    if (xbits.is_nan()) {
+      return x + 1.0f;
+    }
+    if (x_abs == 0x3F80'0000U) {
+      fputil::set_except(FE_DIVBYZERO);
+      return with_errno(FPBits::inf(sign).get_val(), ERANGE);
+    } else {
+      fputil::set_except(FE_INVALID);
+      return with_errno(
+          FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1)),
+          EDOM);
+    }
+  }
+
+  // |x| < ~0.10
+  if (unlikely(x_abs <= 0x3dcc'0000U)) {
+    // |x| <= 2^-26
+    if (unlikely(x_abs <= 0x3280'0000U)) {
+      return unlikely(x_abs == 0) ? x : (x + 0x1.5555555555555p-2 * x * x * x);
+    }
+
+    double xdbl = x;
+    double x2 = xdbl * xdbl;
+    // Pure Taylor series.
+    double pe = fputil::polyeval(x2, 0.0, 0x1.5555555555555p-2,
+                                 0x1.999999999999ap-3, 0x1.2492492492492p-3,
+                                 0x1.c71c71c71c71cp-4, 0x1.745d1745d1746p-4);
+    return fputil::multiply_add(xdbl, pe, xdbl);
+  }
+  double xdbl = x;
+  return 0.5 * log_eval((xdbl + 1.0) / (xdbl - 1.0));
+}
+
+} // namespace __llvm_libc

libc/src/math/generic/supfuncf.h

 //===-- 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
 //
 //===----------------------------------------------------------------------===//
 
 #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"
 #include "src/__support/FPUtil/PolyEval.h"
 #include "src/__support/FPUtil/nearest_integer.h"
 #include "src/__support/common.h"
struct exe_eval_result_t {
   // where
   //   MULT_POWER2 template parameter;
   //   mult_exp = 2^e;
   //   r in range [~-0.3, ~0.41]
   double mult_exp;
   double r;
 };
 
+// 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 function correctly calculates exp value with at least float precision
 // in range not narrow than [-log(2^-150), 90]
 template <int MULT_POWER2 = 0>
 inline static exe_eval_result_t exp_eval(double x) {
   double ps_dbl = fputil::nearest_integer(LN2_INV * x);
   // Negative sign due to multiply_add optimization
   double mult_e1, ml;
   {
uint32_t ps = static_cast<int>(kf) + (1 << (EXP_bits_p - 1)) +
                   (fputil::FPBits<double>::EXPONENT_BIAS << EXP_bits_p);
     fputil::FPBits<double> bs;
     bs.set_unbiased_exponent(ps >> EXP_bits_p);
     ml = 1.0 + EXP_2_POW[ps & (EXP_num_p - 1)];
     mult_f = bs.get_val();
   }
 
   // Taylor series coefficients for 2^x
   double pe = fputil::polyeval(
       dx, 1.0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58fp-3, 0x1.c6b08d704a0c0p-5,
       0x1.3b2ab6fba4e77p-7, 0x1.5d87fe78a6731p-10, 0x1.430912f86c787p-13);
 
   return mult_f * ml * pe;
 }
 
 // x should be positive, normal finite value
 inline static double log2_eval(double x) {
inline static double log2_eval(double x) {
   bs.bits &= FPB::FloatProp::MANTISSA_MASK >> LOG_P1_BITS;
   bs.set_unbiased_exponent(FPB::FloatProp::EXPONENT_BIAS);
   double dx = (bs.get_val() - 1.0) * LOG_P1_1_OVER[p1];
 
   // Taylor series for log(2,1+x)
   double c1 = fputil::multiply_add(dx, K_LOG2_ODD[0], K_LOG2_EVEN[0]);
   double c2 = fputil::multiply_add(dx, K_LOG2_ODD[1], K_LOG2_EVEN[1]);
   double c3 = fputil::multiply_add(dx, K_LOG2_ODD[2], K_LOG2_EVEN[2]);
   double c4 = fputil::multiply_add(dx, K_LOG2_ODD[3], K_LOG2_EVEN[3]);
 
   // c0 = dx * (1.0 / ln(2)) + LOG_P1_LOG2[p1]
   double c0 = fputil::multiply_add(dx, 0x1.71547652b82fep+0, LOG_P1_LOG2[p1]);
   result += __llvm_libc::fputil::polyeval(dx * dx, c0, c1, c2, c3, c4);
   return result;
 }
 
 // x should be positive, normal finite value
 inline static double log_eval(double x) {
   // ln(x) = log[2,x] * ln(2)
   return log2_eval(x) * 0x1.62e42fefa39efp-1;
 }
 
+// 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_PI_2) : (M_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_EXPLOGXF_H

libc/src/math/generic/supfuncf.cpp

alignas(64) const double EXP_2_POW[EXP_num_p] = {
     -0x1.53f391822dbc7p-4, -0x1.01b466423250ap-4, -0x1.5b505d5b6f268p-5,
     -0x1.5f134923757f3p-6, 0x0.0000000000000p+0,  0x1.66c34c5615d0fp-6,
     0x1.6ab0d9f3121ecp-5,  0x1.1301d0125b50ap-4,  0x1.72b83c7d517aep-4,
     0x1.d4873168b9aa8p-4,  0x1.1c3d373ab11c3p-3,  0x1.4f4efa8fef709p-3,
     0x1.837f0518db8a9p-3,  0x1.b8d39b9d54e55p-3,  0x1.ef5326091a112p-3,
     0x1.13821818624b4p-2,  0x1.2ff6b54d8a89cp-2,  0x1.4d0ad5a753e07p-2,
     0x1.6ac1f752150a5p-2,  0x1.891fac0e95613p-2};
 
 // N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
 alignas(64) const double LOG_P1_LOG2[LOG_P1_SIZE] = {
     0x0.0000000000000p+0, 0x1.6e79685c2d22ap-6, 0x1.6bad3758efd87p-5,
     0x1.0eb389fa29f9bp-4, 0x1.663f6fac91316p-4, 0x1.bc84240adabbap-4,
     0x1.08c588cda79e4p-3, 0x1.32ae9e278ae1ap-3, 0x1.5c01a39fbd688p-3,
     0x1.84c2bd02f03b3p-3, 0x1.acf5e2db4ec94p-3, 0x1.d49ee4c325970p-3,
     0x1.fbc16b902680ap-3, 0x1.11307dad30b76p-2, 0x1.24407ab0e073ap-2,
     0x1.37124cea4cdedp-2, 0x1.49a784bcd1b8bp-2, 0x1.5c01a39fbd688p-2,
     0x1.6e221cd9d0cdep-2, 0x1.800a563161c54p-2, 0x1.91bba891f1709p-2,

 alignas(64) const
     double K_LOG2_ODD[4] = {0x1.ec709dc3a03fdp-2, 0x1.2776c50ef9bfep-2,
                             0x1.a61762a7aded9p-3, 0x1.484b13d7c02a9p-3};
 
 alignas(64) const
     double K_LOG2_EVEN[4] = {-0x1.71547652b82fep-1, -0x1.71547652b82fep-2,
                              -0x1.ec709dc3a03fdp-3, -0x1.2776c50ef9bfep-3};
 
+// 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/test/src/math/CMakeLists.txt

SUITE
     libc_math_unittests
   SRCS
     tanhf_test.cpp
   DEPENDS
     libc.src.math.tanhf
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  atanhf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    atanhf_test.cpp
+  DEPENDS
+    libc.src.math.atanhf
+    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_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)

libc/test/src/math/atanf_test.cpp

+//===-- Unittests for atanf -----------------------------------------------===//
+//
+// 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/__support/FPUtil/FPBits.h"
+#include "src/math/atanf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+#include <initializer_list>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+TEST(LlvmLibcAtanfTest, SpecialNumbers) {
+  errno = 0;
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanf(aNaN));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(0.0f, __llvm_libc::atanf(0.0f));
+  EXPECT_FP_EXCEPTION(0);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(-0.0f, __llvm_libc::atanf(-0.0f));
+  EXPECT_FP_EXCEPTION(0);
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST(LlvmLibcAtanfTest, InFloatRange) {
+  constexpr uint32_t COUNT = 1000000;
+  const uint32_t STEP = FPBits(inf).uintval() / COUNT;
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = float(FPBits(v));
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atan, x,
+                                   __llvm_libc::atanf(x), 0.5);
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atan, -x,
+                                   __llvm_libc::atanf(-x), 0.5);
+  }
+}
+
+// For small values, tanh(x) is x.
+TEST(LlvmLibcAtanfTest, SpecialValues) {
+  for (uint32_t v : {0x3d8d6b23U, 0x3feefcfbU, 0xbd8d6b23U, 0xbfeefcfbU,
+                     0x7F800000U, 0xFF800000U}) {
+    float x = float(FPBits(v));
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atan, x,
+                                   __llvm_libc::atanf(x), 0.5);
+  }
+}

libc/test/src/math/atanhf_test.cpp

+//===-- Unittests for atanhf ----------------------------------------------===//
+//
+// 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/__support/FPUtil/FPBits.h"
+#include "src/math/atanhf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+TEST(LlvmLibcAtanhfTest, SpecialNumbers) {
+  errno = 0;
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(aNaN));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(0.0f, __llvm_libc::atanhf(0.0f));
+  EXPECT_FP_EXCEPTION(0);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(-0.0f, __llvm_libc::atanhf(-0.0f));
+  EXPECT_FP_EXCEPTION(0);
+  EXPECT_MATH_ERRNO(0);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(inf, __llvm_libc::atanhf(1.0f));
+  EXPECT_FP_EXCEPTION(FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(neg_inf, __llvm_libc::atanhf(-1.0f));
+  EXPECT_FP_EXCEPTION(FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  auto bt = FPBits(1.0f);
+  bt.bits += 1;
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(bt.get_val()));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  bt.set_sign(true);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(bt.get_val()));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(2.0f));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(-2.0f));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(inf));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+
+  bt.set_sign(true);
+  EXPECT_FP_EQ(aNaN, __llvm_libc::atanhf(neg_inf));
+  EXPECT_FP_EXCEPTION(FE_INVALID);
+  EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST(LlvmLibcAtanhfTest, InFloatRange) {
+  constexpr uint32_t COUNT = 1000000;
+  const uint32_t STEP = FPBits(1.0f).uintval() / COUNT;
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = float(FPBits(v));
+    ASSERT_MPFR_MATCH(mpfr::Operation::Atanh, x, __llvm_libc::atanhf(x), 0.5);
+    ASSERT_MPFR_MATCH(mpfr::Operation::Atanh, -x, __llvm_libc::atanhf(-x), 0.5);
+  }
+}
+
+// For small values, tanh(x) is x.
+TEST(LlvmLibcAtanhfTest, SmallValues) {
+  float x = float(FPBits(uint32_t(0x17800000)));
+  float result = __llvm_libc::atanhf(x);
+  EXPECT_MPFR_MATCH(mpfr::Operation::Atanh, x, result, 0.5);
+  EXPECT_FP_EQ(x, result);
+
+  x = float(FPBits(uint32_t(0x00400000)));
+  result = __llvm_libc::atanhf(x);
+  EXPECT_MPFR_MATCH(mpfr::Operation::Atanh, x, result, 0.5);
+  EXPECT_FP_EQ(x, result);
+}

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

DEPENDS
     .exhaustive_test
     libc.include.math
     libc.src.math.tanhf
     libc.src.__support.FPUtil.fputil
   LINK_LIBRARIES
     -lpthread
 )
 
+add_fp_unittest(
+  atanhf_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    atanhf_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.include.math
+    libc.src.math.atanhf
+    libc.src.__support.FPUtil.fputil
+  LINK_LIBRARIES
+    -lpthread
+)
+
+add_fp_unittest(
+  atanf_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    atanf_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.include.math
+    libc.src.math.atanf
+    libc.src.__support.FPUtil.fputil
+  LINK_LIBRARIES
+    -lpthread
+)
+

libc/test/src/math/exhaustive/atanf_test.cpp

+//===-- Exhaustive test for atanf -----------------------------------------===//
+//
+// 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 "exhaustive_test.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/atanf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <thread>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+struct LlvmLibcAtanfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
+  bool check(uint32_t start, uint32_t stop,
+             mpfr::RoundingMode rounding) override {
+    mpfr::ForceRoundingMode r(rounding);
+    uint32_t bits = start;
+    bool result = true;
+    do {
+      FPBits xbits(bits);
+      float x = float(xbits);
+      result &= EXPECT_MPFR_MATCH(mpfr::Operation::Atan, x,
+                                  __llvm_libc::atanf(x), 0.5, rounding);
+    } while (bits++ < stop);
+    return result;
+  }
+};
+
+static const int NUM_THREADS = std::thread::hardware_concurrency();
+
+// Range: [0, 1.0];
+static const uint32_t POS_START = 0x0000'0000U;
+static const uint32_t POS_STOP = FPBits::inf().uintval();
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, PostiveRangeRoundUp) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, PostiveRangeRoundDown) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, PostiveRangeRoundTowardZero) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero);
+}
+
+// Range: [-1.0, 0];
+static const uint32_t NEG_START = 0x8000'0000U;
+static const uint32_t NEG_STOP = FPBits::neg_inf().uintval();
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, NegativeRangeRoundUp) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, NegativeRangeRoundDown) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAtanfExhaustiveTest, NegativeRangeRoundTowardZero) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
+}

libc/test/src/math/exhaustive/atanhf_test.cpp

+//===-- Exhaustive test for atanhf ----------------------------------------===//
+//
+// 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 "exhaustive_test.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/atanhf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <thread>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+struct LlvmLibcAtanhfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
+  bool check(uint32_t start, uint32_t stop,
+             mpfr::RoundingMode rounding) override {
+    mpfr::ForceRoundingMode r(rounding);
+    uint32_t bits = start;
+    bool result = true;
+    do {
+      FPBits xbits(bits);
+      float x = float(xbits);
+      result &= EXPECT_MPFR_MATCH(mpfr::Operation::Atanh, x,
+                                  __llvm_libc::atanhf(x), 0.5, rounding);
+    } while (bits++ < stop);
+    return result;
+  }
+};
+
+static const int NUM_THREADS = std::thread::hardware_concurrency();
+
+// Range: [0, 1.0];
+static const uint32_t POS_START = 0x0000'0000U;
+static const uint32_t POS_STOP = FPBits(1.0f).uintval();
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, PostiveRangeRoundUp) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, PostiveRangeRoundDown) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, PostiveRangeRoundTowardZero) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero);
+}
+
+// Range: [-1.0, 0];
+static const uint32_t NEG_START = 0x8000'0000U;
+static const uint32_t NEG_STOP = FPBits(-1.0f).uintval();
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, NegativeRangeRoundUp) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, NegativeRangeRoundDown) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAtanhfExhaustiveTest, NegativeRangeRoundTowardZero) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
+}

libc/test/src/math/supfuncf_test.cpp

#define CHECK_DATA(start, stop, mfp_op, f, f_check, count, prec)               \
   {                                                                            \
     uint64_t ustart = FPBits(start).uintval();                                 \
     uint64_t ustop = FPBits(stop).uintval();                                   \
     for (uint64_t i = 0;; ++i) {                                               \
       uint64_t v = ustart + (ustop - ustart) * i / count;                      \
       if (v > ustop)                                                           \
         break;                                                                 \
       float x = FPBits(uint32_t(v)).get_val();                                 \
-      if ((f_check)(x))                                                        \
+      if ((f_check)(x)) {                                                      \
         EXPECT_MPFR_MATCH_ALL_ROUNDING(mfp_op, x, static_cast<float>((f)(x)),  \
                                        (prec));                                \
-    }                                                                          \
+      }                                                                        \
+    }                                                                          \
   }
 
 TEST(LlvmLibcExpxfTest, InFloatRange) {
   auto fx = [](float x) -> float {
     auto result = __llvm_libc::exp_eval<-1>(x);
     return static_cast<float>(2 * result.mult_exp * result.r +
                               2 * result.mult_exp);
   };
TEST(LlvmLibcLog2xfTest, InFloatRange) {
   CHECK_DATA(0.0f, inf, mpfr::Operation::Log2, __llvm_libc::log2_eval, isnormal,
              def_count, def_prec);
 }
 
 TEST(LlvmLibcLogxfTest, InFloatRange) {
   CHECK_DATA(0.0f, inf, mpfr::Operation::Log, __llvm_libc::log_eval, isnormal,
              def_count, def_prec);
 }
+
+TEST(LlvmLibcAtanfPosTest, InFloatRange) {
+  CHECK_DATA(0.0f, inf, mpfr::Operation::Atan, __llvm_libc::atan_eval, isfinite,
+             def_count, def_prec);
+}
+
+TEST(LlvmLibcAtanfNegTest, InFloatRange) {
+  CHECK_DATA(-0.0f, neg_inf, mpfr::Operation::Atan, __llvm_libc::atan_eval,
+             isfinite, def_count, def_prec);
+}

libc/utils/MPFRWrapper/MPFRUtils.h

 namespace mpfr {
 
 enum class Operation : int {
   // Operations with take a single floating point number as input
   // and produce a single floating point number as output. The input
   // and output floating point numbers are of the same kind.
   BeginUnaryOperationsSingleOutput,
   Abs,
+  Atan,
+  Atanh,
   Ceil,
   Cos,
   Cosh,
   Exp,
   Exp2,
   Expm1,
   Floor,
   Log,

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber abs() const {
     MPFRNumber result(*this);
     mpfr_abs(result.value, value, mpfr_rounding);
     return result;
   }
 
+  MPFRNumber atan() const {
+    MPFRNumber result(*this);
+    mpfr_atan(result.value, value, mpfr_rounding);
+    return result;
+  }
+
+  MPFRNumber atanh() const {
+    MPFRNumber result(*this);
+    mpfr_atanh(result.value, value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber ceil() const {
     MPFRNumber result(*this);
     mpfr_ceil(result.value, value);
     return result;
   }
 
   MPFRNumber cos() const {
     MPFRNumber result(*this);

 template <typename InputType>
 cpp::enable_if_t<cpp::is_floating_point_v<InputType>, MPFRNumber>
 unary_operation(Operation op, InputType input, unsigned int precision,
                 RoundingMode rounding) {
   MPFRNumber mpfrInput(input, precision, rounding);
   switch (op) {
   case Operation::Abs:
     return mpfrInput.abs();
+  case Operation::Atan:
+    return mpfrInput.atan();
+  case Operation::Atanh:
+    return mpfrInput.atanh();
   case Operation::Ceil:
     return mpfrInput.ceil();
   case Operation::Cos:
     return mpfrInput.cos();
   case Operation::Cosh:
     return mpfrInput.cosh();
   case Operation::Exp:
     return mpfrInput.exp();