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

[lntue]

Summary

Implement asinf function correctly rounded for all rounding modes.

For |x| <= 0.5, we approximate asin(x) by

asin(x) = 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]);

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(2X) = 1 - 2 sin^2(X)

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,

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).

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

$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh asinf
CORE-MATH reciprocal throughput   : 23.418
System LIBC reciprocal throughput : 27.310
LIBC reciprocal throughput        : 22.741

$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh asinf --latency
GNU libc version: 2.35
GNU libc release: stable
CORE-MATH latency   : 58.884
System LIBC latency : 62.055
LIBC latency        : 62.037

[zimmermann6]

the patch fails for me on top of main (revision ea953b9). Is there any other patch to apply first?

[lntue]

Rebase.

[lntue]

In D133400#3773780, @zimmermann6 wrote:

the patch fails for me on top of main (revision ea953b9). Is there any other patch to apply first?

I just uploaded the patch after rebasing. Can you try to see if it works now?

Here is the result of git log --oneline on my local repo:

7b6f9731ca7a (HEAD -> main) [libc][math] Implement asinf function correctly rounded for all rounding modes.
c836ddaf721c (origin/main) [X86][NFC] Refine load/store reg to StackSlot for extensibility
4e5a59a3839f [LLD][COFF] Fix writing a map file when range extension thunks are inserted
ea953b9d9a65 [compiler-rt] [test] Handle missing ld.gold gracefully

[zimmermann6]

it works now, thanks. I confirm it is correctly rounded. I get similar figures on a AMD EPYC 7282 (glibc, core-math, llvm-libc):

# reciprocal throughput
GNU libc version: 2.34
GNU libc release: stable
26.722
31.752
27.841
# latency
GNU libc version: 2.34
GNU libc release: stable
56.310
64.140
61.051

[orex]

Good job, Tue!

I have an optional suggestion. You can use formula
asin(x)=1/2*acos(1-2*x^2)
and, of course,
acos(x)=pi/2-asin(x)
for cases where |x|>=0.5 Probably it can be much faster than sqrt.

[lntue]

In D133400#3774132, @orex wrote:

Good job, Tue!

I have an optional suggestion. You can use formula
asin(x)=1/2*acos(1-2*x^2)
and, of course,
acos(x)=pi/2-asin(x)
for cases where |x|>=0.5 Probably it can be much faster than sqrt.

Thanks @orex for the suggestion! Unfortunately the formula asin(x) = pi/4 - 1/2 * asin(1 - 2*x^2) only works
well when |1 - 2*x^2| <= 0.5, i.e. 0.5 <= |x| <= sqrt(3)/2.

For sqrt(3)/2 < |x| <= 1, we still need some sort of range reduction involving square root.

Nonetheless, I've tried adding another branch for 0.5 <= |x| <= sqrt(3)/2 using the formula
asin(x) = pi/4 - 1/2 * asin(1 - 2*x^2) and it actually makes both the throughput and the latency
slightly worse than 2 branches.

[orex]

In D133400#3774781, @lntue wrote:

In D133400#3774132, @orex wrote:

Good job, Tue!

I have an optional suggestion. You can use formula
asin(x)=1/2*acos(1-2*x^2)
and, of course,
acos(x)=pi/2-asin(x)
for cases where |x|>=0.5 Probably it can be much faster than sqrt.

Thanks @orex for the suggestion! Unfortunately the formula asin(x) = pi/4 - 1/2 * asin(1 - 2*x^2) only works
well when |1 - 2*x^2| <= 0.5, i.e. 0.5 <= |x| <= sqrt(3)/2.

For sqrt(3)/2 < |x| <= 1, we still need some sort of range reduction involving square root.

Nonetheless, I've tried adding another branch for 0.5 <= |x| <= sqrt(3)/2 using the formula
asin(x) = pi/4 - 1/2 * asin(1 - 2*x^2) and it actually makes both the throughput and the latency
slightly worse than 2 branches.

No questions! Thanks for trying the approach.

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.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.src.math.ceill

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.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.src.math.ceill

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.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.src.math.ceill

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.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.src.math.ceill

libc/docs/math.rst

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

 ============== ================ =============== ======================
 
 Accuracy of Higher Math Functions
 =================================
 
 ============== ================ =============== ======================
 <Func>         <Func_f> (float) <Func> (double) <Func_l> (long double)
 ============== ================ =============== ======================
+asin           |check|
 atan           |check|
 atanh          |check|
 cos            |check|          large
 cosh           |check|
 exp            |check|
 exp2           |check|
 expm1          |check|
 fma            |check|          |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 |
 +==============+===========+===================+===========+===================+=====================================+============+=========================+==============+===============+
+| 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           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | coshf        |        23 |                20 |        73 |                49 | :math:`[-10, 10]`                   | 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<"asinf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"atanf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           
           FunctionSpec<"atanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
 
   HeaderSpec StdIO = HeaderSpec<
       "stdio.h",

libc/src/math/CMakeLists.txt

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

libc/src/math/asinf.h

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

libc/src/math/generic/CMakeLists.txt

DEPENDS
     libc.src.__support.FPUtil.nearest_integer_operations
     libc.src.__support.FPUtil.polyeval
     libc.include.errno
     libc.src.errno.errno
     libc.include.math
 )
 
 add_entrypoint_object(
+  asinf
+  SRCS
+    asinf.cpp
+  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
+  COMPILE_OPTIONS
+    -O3
+)
+
+add_entrypoint_object(
   atanf
   SRCS
     atanf.cpp
   HDRS
     ../atanf.h
   DEPENDS
     .inv_trigf_utils
     .math_utils
     libc.src.__support.FPUtil.fp_bits
   COMPILE_OPTIONS
     -O3
 )
 

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 {
+
+// PI / 2
+constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
+
+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;
+
+  // |x| <= 0.5-ish
+  if (x_abs < 0x3f04'471dU) {
+    // |x| < 0x1.d12edp-12
+    if (unlikely(x_abs < 0x39e8'9768U)) {
+      // When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
+      // is:
+      //   |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
+      //                             = x^2 / 6
+      //                             < 2^-25
+      //                             < epsilon(1)/2.
+      // So the correctly rounded values of asin(x) are:
+      //   = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
+      //                        or (rounding mode = FE_UPWARD and x is
+      //                        negative),
+      //   = x otherwise.
+      // To simplify the rounding decision and make it more efficient, we use
+      //   fma(x, 2^-25, x) instead.
+      // An exhaustive test shows that this formula work correctly for all
+      // rounding modes up to |x| < 0x1.d12edp-12.
+      // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
+      // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
+      // |x| < 2^-125. For targets without FMA instructions, we simply use
+      // double for intermediate results as it is more efficient than using an
+      // emulated version of FMA.
+#if defined(LIBC_TARGET_HAS_FMA)
+      return fputil::multiply_add(x, 0x1.0p-25f, x);
+#else
+      double xd = static_cast<double>(x);
+      return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
+#endif // LIBC_TARGET_HAS_FMA
+    }
+
+    // Check for exceptional values
+    if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);
+        unlikely(r.has_value()))
+      return r.value();
+
+    // For |x| <= 0.5, we approximate asinf(x) by:
+    //   asin(x) = 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]);
+    // 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);
+  }
+
+  // |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,
+  //   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),
+
+  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 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;
+}
+
+} // namespace __llvm_libc

libc/test/src/math/CMakeLists.txt

add_fp_unittest(
   DEPENDS
     libc.include.errno
     libc.src.errno.errno
     libc.src.math.atanhf
     libc.src.__support.FPUtil.fp_bits
 )
 
 add_fp_unittest(
+  asinf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    asinf_test.cpp
+  DEPENDS
+    libc.include.errno
+    libc.src.errno.errno
+    libc.src.math.asinf
+    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/asinf_test.cpp

+//===-- Unittests for asinf
+//------------------------------------------------===//
+//
+// 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/asinf.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(LlvmLibcAsinfTest, SpecialNumbers) {
+  errno = 0;
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::asinf(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(0.0f, __llvm_libc::asinf(0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(-0.0f, __llvm_libc::asinf(-0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::asinf(inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::asinf(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST(LlvmLibcAsinfTest, 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::Asin, x,
+                                   __llvm_libc::asinf(x), 0.5);
+  }
+}
+
+TEST(LlvmLibcAsinfTest, 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::Asin, x,
+                                   __llvm_libc::asinf(x), 0.5);
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, -x,
+                                   __llvm_libc::asinf(-x), 0.5);
+  }
+}

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

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

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

+//===-- Exhaustive test for asinf -----------------------------------------===//
+//
+// 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 "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> {
+  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::Asin, x,
+                                  __llvm_libc::asinf(x), 0.5, rounding);
+      // if (!result) break;
+    } 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;
+/
+
+    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) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, 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;
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundUp) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundDown) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundTowardZero) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
+}

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

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
 
   MPFRNumber &operator=(const MPFRNumber &rhs) {
     mpfr_precision = rhs.mpfr_precision;
     mpfr_rounding = rhs.mpfr_rounding;
     mpfr_set(value, rhs.value, mpfr_rounding);
     return *this;
   }
 
+  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 asin() const {
+    MPFRNumber result(*this);
+    mpfr_asin(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);
public:
   //    of N between this number and [input].
   // 4. A values of +0.0 and -0.0 are treated as equal.
   template <typename T>
   cpp::enable_if_t<cpp::is_floating_point_v<T>, double> ulp(T input) {
     T thisAsT = as<T>();
     if (thisAsT == input)
       return T(0.0);
 
+    if (is_nan()) {
+      if (fputil::FPBits<T>(input).is_nan())
+        return T(0.0);
+      return T(fputil::FPBits<T>::inf());
+    }
+
     int thisExponent = fputil::FPBits<T>(thisAsT).get_exponent();
     int inputExponent = fputil::FPBits<T>(input).get_exponent();
     // Adjust the exponents for denormal numbers.
     if (fputil::FPBits<T>(thisAsT).get_unbiased_exponent() == 0)
       ++thisExponent;
     if (fputil::FPBits<T>(input).get_unbiased_exponent() == 0)
       ++inputExponent;
 

 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::Asin:
+    return mpfrInput.asin();
   case Operation::Atan:
     return mpfrInput.atan();
   case Operation::Atanh:
     return mpfrInput.atanh();
   case Operation::Ceil:
     return mpfrInput.ceil();
   case Operation::Cos:
     return mpfrInput.cos();