[libc][math] Implement acosf function correctly rounded for all rounding modes.

[lntue]

Summary

Implement acosf function correctly rounded for all rounding modes.

We perform range reduction as follows:

• When |x| < 2^(-10), we use cubic Taylor polynomial:

acos(x) = pi/2 - asin(x) ~ pi/2 - x - x^3 / 6.

• When 2^(-10) <= |x| <= 0.5, we use the same approximation that is used for asinf(x) when |x| <= 0.5:

acos(x) = pi/2 - asin(x) ~ pi/2 - x - x^3 * P(x^2).

• When 0.5 < x <= 1, we use the double angle formula: cos(2y) = 1 - 2 * sin^2 (y) to reduce to:

acos(x) = 2 * asin( sqrt( (1 - x)/2 ) )

• When -1 <= x < -0.5, we reduce to the positive case above using the formula:

acos(x) = pi - acos(-x)

Performance benchmark using perf tool from the CORE-MATH project on Ryzen 1700:

$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh acosf
GNU libc version: 2.35
GNU libc release: stable
CORE-MATH reciprocal throughput   : 28.613
System LIBC reciprocal throughput : 29.204
LIBC reciprocal throughput        : 24.271

$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh asinf --latency
GNU libc version: 2.35
GNU libc release: stable
CORE-MATH latency   : 55.554
System LIBC latency : 76.879
LIBC latency        : 62.118

[lntue]

Fix tests and exceptional values.

[zimmermann6]

great work! The reciprocal throughput is indeed slightly better than CORE-MATH, and the latency slightly worse:

# reciprocal throughput
GNU libc version: 2.34
GNU libc release: stable
33.819
37.064
29.462
# latency
GNU libc version: 2.34
GNU libc release: stable
54.951
80.046
62.001

[orex]

libc/src/math/generic/asinf.cpp

17

Please delete this. It is not needed here. My fault.

libc/src/math/generic/inv_trigf_utils.h

95

Don't you think that it is better to put this array to cpp file?

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

39

Do you really need this test until inf? Out of range values can be covered by unittests?

[lntue]

libc/src/math/generic/asinf.cpp

17

We do need this to set EDOM for out of range inputs. errno.h header is generated separately (it's kind of complete I think) and does not depend on math.h, so it's safe to include and use in math entrypoints. That's not true for math.h constants though, as we are building that header's implementation here.

libc/src/math/generic/inv_trigf_utils.h

95

This is a small table that we do want to inline, and we don't take its address anywhere, so by leaving its definition in the header, it does improve the performance. I got the reciprocal throughput of 24 for leaving its definition here, vs 26 for putting its definition in the cpp file and providing external linkage in the header.

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

39

I think it's better to leave the whole range in the committed tests to safe guard against future changes. Ideally (and soon) we will run these tests automatically with the CI's. For manual testing, we can always restrict the the range that we are interested. The main reason that I didn't make other exhaustive tests running full range is that the ulp function comparing to mpfr outputs did not handle NaN properly. It was fixed in https://reviews.llvm.org/D133400, so I think it's better to just test everything in these exhaustive tests.

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.acosf
     libc.src.math.asinf
     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/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.acosf
     libc.src.math.asinf
     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/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.acosf
     libc.src.math.asinf
     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/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.acosf
     libc.src.math.asinf
     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/docs/math.rst

 ============== ================ =============== ======================
 
 Higher Math Functions
 ---------------------
 
 ============== ================ =============== ======================
 <Func>         <Func_f> (float) <Func> (double) <Func_l> (long double)
 ============== ================ =============== ======================
-acos
+acos           |check|
 acosh
 asin           |check|
 asinh
 atan           |check|
 atan2
 atanh          |check|
 cbrt
 cos            |check|          |check|

 ============== ================ =============== ======================
 
 Accuracy of Higher Math Functions
 =================================
 
 ============== ================ =============== ======================
 <Func>         <Func_f> (float) <Func> (double) <Func_l> (long double)
 ============== ================ =============== ======================
+acos           |check|
 asin           |check|
 atan           |check|
 atanh          |check|
 cos            |check|          large
 cosh           |check|
 exp            |check|
 exp2           |check|
 expm1          |check|
* We also use the *perf* tool from the `CORE-MATH <https://core-math.gitlabpages.inria.fr/>`_
   table below, using the system library as reference (such as the `GNU C library <https://www.gnu.org/software/libc/>`_
   on Linux). Fmod performance results obtained with "differential_testing".
 
 +--------------+-------------------------------+-------------------------------+-------------------------------------+---------------------------------------------------------------------+
 | <Func>       | Reciprocal throughput (ns)    | Latency (ns)                  | Testing ranges                      | Testing configuration                                               |
 |              +-----------+-------------------+-----------+-------------------+                                     +------------+-------------------------+--------------+---------------+
 |              | LLVM libc | Reference (glibc) | LLVM libc | Reference (glibc) |                                     | CPU        | OS                      | Compiler     | Special flags |
 +==============+===========+===================+===========+===================+=====================================+============+=========================+==============+===============+
+| acosf        |        24 |                29 |        62 |                77 | :math:`[-1, 1]`                     | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
++--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | asinf        |        23 |                27 |        62 |                62 | :math:`[-1, 1]`                     | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | atanf        |        27 |                29 |        79 |                68 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | atanhf       |        20 |                66 |        71 |               133 | :math:`[-1, 1]`                     | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | cosf         |        13 |                32 |        53 |                59 | :math:`[0, 2\pi]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+

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<"acosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"asinf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"atanf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           
           FunctionSpec<"atanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
 
   HeaderSpec StdIO = HeaderSpec<

libc/src/math/CMakeLists.txt

add_entrypoint_object(
   HDRS
     fma.h
   DEPENDS
     libc.src.__support.FPUtil.fma
   COMPILE_OPTIONS
     -O3
 )
 
+add_math_entrypoint_object(acosf)
 add_math_entrypoint_object(asinf)
 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)

libc/src/math/acosf.h

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

libc/src/math/generic/CMakeLists.txt

add_entrypoint_object(
   HDRS
     ../asinf.h
   DEPENDS
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.polyeval
     libc.src.__support.FPUtil.sqrt
+    .inv_trigf_utils
+  COMPILE_OPTIONS
+    -O3
+)
+
+add_entrypoint_object(
+  acosf
+  SRCS
+    acosf.cpp
+  HDRS
+    ../acosf.h
+  DEPENDS
+    libc.src.__support.FPUtil.except_value_utils
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
+    libc.src.__support.FPUtil.sqrt
+    .inv_trigf_utils
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   atanf
   SRCS
     atanf.cpp

libc/src/math/generic/acosf.cpp

+//===-- Single-precision acos 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/acosf.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/sqrt.h"
+
+#include <errno.h>
+
+#include "inv_trigf_utils.h"
+
+namespace __llvm_libc {
+
+static constexpr size_t N_EXCEPTS = 2;
+
+// Exceptional values when |x| <= 0.5
+static constexpr fputil::ExceptValues<float, N_EXCEPTS> ACOSF_EXCEPTS = {{
+    // (inputs, RZ output, RU offset, RD offset, RN offset)
+    // x = 0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)
+    {0x328885a3, 0x3d09c62c, 1, 0, 1},
+    // x = 0x1.04c444p-12, acosf(x) = 0x1.920f68p0 (RZ)
+    {0x39826222, 0x3fc907b4, 1, 0, 1},
+}};
+
+LLVM_LIBC_FUNCTION(float, acosf, (float x)) {
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+  uint32_t x_uint = xbits.uintval();
+  uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
+  uint32_t x_sign = x_uint >> 31;
+
+  // |x| <= 0.5
+  if (unlikely(x_abs <= 0x3f00'0000U)) {
+    // |x| < 0x1p-10
+    if (unlikely(x_abs < 0x3a80'0000U)) {
+      // When |x| < 2^-10, we use the following approximation:
+      //   acos(x) = pi/2 - asin(x)
+      //           ~ pi/2 - x - x^3 / 6
+
+      // Check for exceptional values
+      if (auto r = ACOSF_EXCEPTS.lookup(x_abs); unlikely(r.has_value()))
+        return r.value();
+
+      double xd = static_cast<double>(x);
+      return fputil::multiply_add(-0x1.5555555555555p-3 * xd, xd * xd,
+                                  M_MATH_PI_2 - xd);
+    }
+
+    // For |x| <= 0.5, we approximate acosf(x) by:
+    //   acos(x) = pi/2 - asin(x) = pi/2 - x * P(x^2)
+    // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
+    // asin(x)/x on [0, 0.5] generated by Sollya with:
+    // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
+    //                 [|1, D...|], [0, 0.5]);
+    double xd = static_cast<double>(x);
+    double xsq = xd * xd;
+    double x3 = xd * xsq;
+    double r = asin_eval(xsq);
+    return fputil::multiply_add(-x3, r, M_MATH_PI_2 - xd);
+  }
+
+  // |x| > 1, return NaNs.
+  if (unlikely(x_abs > 0x3f80'0000U)) {
+    if (x_abs <= 0x7f80'0000U) {
+      errno = EDOM;
+      fputil::set_except(FE_INVALID);
+    }
+    return x +
+           FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
+  }
+
+  // When 0.5 < |x| <= 1, we perform range reduction as follow:
+  //
+  // Assume further that 0.5 < x <= 1, and let:
+  //   y = acos(x)
+  // We use the double angle formula:
+  //   x = cos(y) = 1 - 2 sin^2(y/2)
+  // So:
+  //   sin(y/2) = sqrt( (1 - x)/2 )
+  // And hence:
+  //   y = 2 * asin( sqrt( (1 - x)/2 ) )
+  // Let u = (1 - x)/2, then
+  //   acos(x) = 2 * asin( sqrt(u) )
+  // Moreover, since 0.5 < x <= 1,
+  //   0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
+  // And hence we can reuse the same polynomial approximation of asin(x) when
+  // |x| <= 0.5:
+  //   acos(x) ~ 2 * sqrt(u) * P(u).
+  //
+  // When -1 <= x <= -0.5, we use the identity:
+  //   acos(x) = pi - acos(-x)
+  // which is reduced to the postive case.
+
+  xbits.set_sign(false);
+  double xd = static_cast<double>(xbits.get_val());
+  double u = fputil::multiply_add(-0.5, xd, 0.5);
+  double cv = 2 * fputil::sqrt(u);
+
+  double r3 = asin_eval(u);
+  double r = fputil::multiply_add(cv * u, r3, cv);
+  return x_sign ? M_MATH_PI - r : r;
+}
+
+} // namespace __llvm_libc

libc/src/math/generic/asinf.cpp

 //===-- Single-precision asin 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/asinf.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/sqrt.h"
 
 #include <errno.h>
 
-namespace __llvm_libc {
+#include "inv_trigf_utils.h"
 
-// PI / 2
-constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
+namespace __llvm_libc {
 
 static constexpr size_t N_EXCEPTS = 2;
 
 // Exceptional values when |x| <= 0.5
 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{
     // (inputs, RZ output, RU offset, RD offset, RN offset)
     // x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)
     {0x3d09bf86, 0x3d09c62c, 1, 0, 1},
     // x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)
     {0x3de5fa1e, 0x3de6768e, 1, 0, 0},
 }};
 
 // Exceptional values when 0.5 < |x| <= 1
 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{
     // (inputs, RZ output, RU offset, RD offset, RN offset)
     // x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)
     {0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},
     // x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)
     {0x3f7741b6, 0x3fa7832a, 1, 0, 0},
 }};
 
-// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
-//                 [|1, D...|], [0, 0.5]);
-static constexpr double COEFFS[10] = {
-    0x1.5555555540fa1p-3, 0x1.333333512edc2p-4, 0x1.6db6cc1541b31p-5,
-    0x1.f1caff324770ep-6, 0x1.6e43899f5f4f4p-6, 0x1.1f847cf652577p-6,
-    0x1.9b60f47f87146p-7, 0x1.259e2634c494fp-6, -0x1.df946fa875ddp-8,
-    0x1.02311ecf99c28p-5};
-
 LLVM_LIBC_FUNCTION(float, asinf, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
   uint32_t x_uint = xbits.uintval();
   uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
   constexpr double SIGN[2] = {1.0, -1.0};
   uint32_t x_sign = x_uint >> 31;
 
#endif // LIBC_TARGET_HAS_FMA
     // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
     // asin(x)/x on [0, 0.5] generated by Sollya with:
     // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
     //                 [|1, D...|], [0, 0.5]);
     // An exhaustive test shows that this approximation works well up to a
     // little more than 0.5.
     double xd = static_cast<double>(x);
     double xsq = xd * xd;
-    double x4 = xsq * xsq;
-    double r1 = fputil::polyeval(x4, COEFFS[0], COEFFS[2], COEFFS[4], COEFFS[6],
-                                 COEFFS[8]);
-    double r2 = fputil::polyeval(x4, COEFFS[1], COEFFS[3], COEFFS[5], COEFFS[7],
-                                 COEFFS[9]);
-    double r3 = fputil::multiply_add(xsq, r2, r1);
     double x3 = xd * xsq;
-    return fputil::multiply_add(x3, r3, xd);
+    double r = asin_eval(xsq);
+    return fputil::multiply_add(x3, r, xd);
   }
 
   // |x| > 1, return NaNs.
   if (unlikely(x_abs > 0x3f80'0000U)) {
     if (x_abs <= 0x7f80'0000U) {
       errno = EDOM;
       fputil::set_except(FE_INVALID);
     }
     return x +
            FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
   }
 
   // Check for exceptional values
   if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);
       unlikely(r.has_value()))
     return r.value();
 
   // When |x| > 0.5, we perform range reduction as follow:
+  //
   // Assume further that 0.5 < x <= 1, and let:
   //   y = asin(x)
   // We will use the double angle formula:
   //   cos(2y) = 1 - 2 sin^2(y)
   // and the complement angle identity:
   //   x = sin(y) = cos(pi/2 - y)
   //              = 1 - 2 sin^2 (pi/4 - y/2)
   // So:
   //   sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
   // And hence:
   //   pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
   // Equivalently:
   //   asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
-  // Let u = (1 - x)/2, then
-  //   asin(x) = pi/2 - 2 * asin(u)
-  // Moreover, since 0.5 < x <= 1,
+  // Let u = (1 - x)/2, then:
+  //   asin(x) = pi/2 - 2 * asin( sqrt(u) )
+  // Moreover, since 0.5 < x <= 1:
   //   0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
   // And hence we can reuse the same polynomial approximation of asin(x) when
   // |x| <= 0.5:
-  //   asin(x) = pi/2 - 2 * u * P(u^2),
+  //   asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
 
   xbits.set_sign(false);
+  double sign = SIGN[x_sign];
   double xd = static_cast<double>(xbits.get_val());
   double u = fputil::multiply_add(-0.5, xd, 0.5);
-  double cv = -2 * fputil::sqrt(u);
+  double c1 = sign * (-2 * fputil::sqrt(u));
+  double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1);
+  double c3 = c1 * u;
 
-  double usq = u * u;
-  double r1 = fputil::polyeval(usq, COEFFS[0], COEFFS[2], COEFFS[4], COEFFS[6],
-                               COEFFS[8]);
-  double r2 = fputil::polyeval(usq, COEFFS[1], COEFFS[3], COEFFS[5], COEFFS[7],
-                               COEFFS[9]);
-  double r3 = fputil::multiply_add(u, r2, r1);
-  double r = fputil::multiply_add(cv * u, r3, M_MATH_PI_2 + cv);
-  return SIGN[x_sign] * r;
+  double r = asin_eval(u);
+  return fputil::multiply_add(c3, r, c2);
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/inv_trigf_utils.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
+// PI and PI / 2
+constexpr double M_MATH_PI = 0x1.921fb54442d18p+1;
 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) {
+static inline 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 =
static inline double atan_eval(double x) {
   double result;
   if (one_over_x)
     result = M_MATH_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;
 }
 
+// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
+//                 [|1, D...|], [0, 0.5]);
+constexpr double ASIN_COEFFS[10] = {0x1.5555555540fa1p-3, 0x1.333333512edc2p-4,
+                                    0x1.6db6cc1541b31p-5, 0x1.f1caff324770ep-6,
+                                    0x1.6e43899f5f4f4p-6, 0x1.1f847cf652577p-6,
+                                    0x1.9b60f47f87146p-7, 0x1.259e2634c494fp-6,
+                                    -0x1.df946fa875ddp-8, 0x1.02311ecf99c28p-5};
+
+// Evaluate P(x^2) - 1, where P(x^2) ~ asin(x)/x
+static inline double asin_eval(double xsq) {
+  double x4 = xsq * xsq;
+  double r1 = fputil::polyeval(x4, ASIN_COEFFS[0], ASIN_COEFFS[2],
+                               ASIN_COEFFS[4], ASIN_COEFFS[6], ASIN_COEFFS[8]);
+  double r2 = fputil::polyeval(x4, ASIN_COEFFS[1], ASIN_COEFFS[3],
+                               ASIN_COEFFS[5], ASIN_COEFFS[7], ASIN_COEFFS[9]);
+  return fputil::multiply_add(xsq, r2, r1);
+}
+
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H

libc/test/src/math/CMakeLists.txt

add_fp_unittest(
   DEPENDS
     libc.include.errno
     libc.src.errno.errno
     libc.src.math.asinf
     libc.src.__support.FPUtil.fp_bits
 )
 
 add_fp_unittest(
+  acosf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    acosf_test.cpp
+  DEPENDS
+    libc.include.errno
+    libc.src.errno.errno
+    libc.src.math.acosf
+    libc.src.__support.FPUtil.fp_bits
+)
+
+add_fp_unittest(
   atanf_test
   NEED_MPFR
   SUITE
     libc_math_unittests
   SRCS
     atanf_test.cpp
   DEPENDS
     libc.include.errno

libc/test/src/math/acosf_test.cpp

+//===-- Unittests for acosf
+//------------------------------------------------===//
+//
+// 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/acosf.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(LlvmLibcAcosfTest, SpecialNumbers) {
+  errno = 0;
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::acosf(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(0.0f, __llvm_libc::acosf(0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(-0.0f, __llvm_libc::acosf(-0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::acosf(inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::acosf(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST(LlvmLibcAcosfTest, InFloatRange) {
+  constexpr uint32_t COUNT = 1000000;
+  constexpr uint32_t STEP = UINT32_MAX / COUNT;
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = float(FPBits(v));
+    if (isnan(x) || isinf(x))
+      continue;
+    ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acos, x,
+                                   __llvm_libc::acosf(x), 0.5);
+  }
+}
+
+TEST(LlvmLibcAcosfTest, SpecificBitPatterns) {
+  constexpr int N = 11;
+  constexpr uint32_t INPUTS[N] = {
+      0x3f000000, // x = 0.5f
+      0x3f3504f3, // x = sqrt(2)/2, FE_DOWNWARD
+      0x3f3504f4, // x = sqrt(2)/2, FE_UPWARD
+      0x3f5db3d7, // x = sqrt(3)/2, FE_DOWNWARD
+      0x3f5db3d8, // x = sqrt(3)/2, FE_UPWARD
+      0x3f800000, // x = 1.0f
+      0x40000000, // x = 2.0f
+      0x3d09bf86, // x = 0x1.137f0cp-5f
+      0x3de5fa1e, // x = 0x1.cbf43cp-4f
+      0x3f083a1a, // x = 0x1.107434p-1f
+      0x3f7741b6, // x = 0x1.ee836cp-1f
+  };
+
+  for (int i = 0; i < N; ++i) {
+    float x = float(FPBits(INPUTS[i]));
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acos, x,
+                                   __llvm_libc::acosf(x), 0.5);
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acos, -x,
+                                   __llvm_libc::acosf(-x), 0.5);
+  }
+}

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

add_fp_unittest(
   DEPENDS
     .exhaustive_test
     libc.include.math
     libc.src.math.asinf
     libc.src.__support.FPUtil.fp_bits
   LINK_LIBRARIES
     -lpthread
 )
+
+add_fp_unittest(
+  acosf_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    acosf_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.include.math
+    libc.src.math.acosf
+    libc.src.__support.FPUtil.fp_bits
+  LINK_LIBRARIES
+    -lpthread
+)

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

-//===-- Exhaustive test for asinf -----------------------------------------===//
+//===-- Exhaustive test for acosf -----------------------------------------===//
 //
 // 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/asinf.h"
+#include "src/math/acosf.h"
 #include "utils/MPFRWrapper/MPFRUtils.h"
 
 #include <thread>
 
 using FPBits = __llvm_libc::fputil::FPBits<float>;
 
 namespace mpfr = __llvm_libc::testing::mpfr;
 
-struct LlvmLibcAsinfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
+struct LlvmLibcAcosfExhaustiveTest : 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;
+    int count = 0;
     do {
       FPBits xbits(bits);
       float x = float(xbits);
-      result &= EXPECT_MPFR_MATCH(mpfr::Operation::Asin, x,
-                                  __llvm_libc::asinf(x), 0.5, rounding);
-      // if (!result) break;
+      result &= EXPECT_MPFR_MATCH(mpfr::Operation::Acos, x,
+                                  __llvm_libc::acosf(x), 0.5, rounding);
     } while (bits++ < stop);
     return result;
   }
 };
 
 static const int NUM_THREADS = std::thread::hardware_concurrency();
 
 // Range: [0, Inf];
 static const uint32_t POS_START = 0x0000'0000U;
-static const uint32_t POS_STOP = 0x7f80'0000U;
-/
+static const uint32_t POS_STOP = 0x3f80'0000U; // 0x7f80'0000U;
 
-    TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
   test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
 }
 
-TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundUp) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundUp) {
   test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
 }
 
-TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundDown) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundDown) {
   test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
 }
 
-TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundTowardZero) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundTowardZero) {
   test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero);
 }
 
 // Range: [-Inf, 0];
-static const uint32_t NEG_START = 0x8000'0000U;
-static const uint32_t NEG_STOP = 0xff80'0000U;
+static const uint32_t NEG_START = 0xb000'0000U;
+static const uint32_t NEG_STOP = 0xbf80'0000U; // 0xff80'0000U;
 
-TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
   test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
 }
 
-TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundUp) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundUp) {
   test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
 }
 
-TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundDown) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundDown) {
   test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
 }
 
-TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundTowardZero) {
+TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundTowardZero) {
   test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
 }

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

struct LlvmLibcAsinfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
   }
 };
 
 static const int NUM_THREADS = std::thread::hardware_concurrency();
 
 // Range: [0, Inf];
 static const uint32_t POS_START = 0x0000'0000U;
 static const uint32_t POS_STOP = 0x7f80'0000U;
-/
 
-    TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
   test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
 }
 
 TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundUp) {
   test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
 }
 
 TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundDown) {

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,
+  Acos,
   Asin,
   Atan,
   Atanh,
   Ceil,
   Cos,
   Cosh,
   Exp,
   Exp2,

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   bool is_nan() const { return mpfr_nan_p(value); }
 
   MPFRNumber abs() const {
     MPFRNumber result(*this);
     mpfr_abs(result.value, value, mpfr_rounding);
     return result;
   }
 
+  MPFRNumber acos() const {
+    MPFRNumber result(*this);
+    mpfr_acos(result.value, value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber asin() const {
     MPFRNumber result(*this);
     mpfr_asin(result.value, value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber atan() 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::Acos:
+    return mpfrInput.acos();
   case Operation::Asin:
     return mpfrInput.asin();
   case Operation::Atan:
     return mpfrInput.atan();
   case Operation::Atanh:
     return mpfrInput.atanh();
   case Operation::Ceil:
     return mpfrInput.ceil();