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

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
```

Reviewed By: orex, zimmermann6

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

1 files changed, 169 insertions, 0 deletions

diff --git a/libc/src/math/generic/asinf.cpp b/libc/src/math/generic/asinf.cpp
new file mode 100644
index 000000000000..262af0fcbd3b
--- /dev/null
+++ b/libc/src/math/generic/asinf.cpp
@@ -0,0 +1,169 @@
+//===-- 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