[libc-commits] [libc] [libc][math] Implement double precision asin correctly rounded for all rounding modes. (PR #134401)

[via libc-commits]

================
@@ -0,0 +1,572 @@
+//===-- Collection of utils for asin/acos -----------------------*- 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_GENERIC_ASIN_UTILS_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_ASIN_UTILS_H
+
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/integer_literals.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace {
+
+using DoubleDouble = fputil::DoubleDouble;
+using Float128 = fputil::DyadicFloat<128>;
+
+constexpr DoubleDouble PI_OVER_TWO = {0x1.1a62633145c07p-54,
+                                      0x1.921fb54442d18p0};
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// When correct rounding is not needed, we use a degree-22 minimax polynomial to
+// approximate asin(x)/x on [0, 0.5] using Sollya with:
+// > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22|],
+//                 [|1, D...|], [0, 0.5]);
+// > dirtyinfnorm(asin(x)/x - P, [0, 0.5]);
+// 0x1.1a71ef0a0f26a9fb7ed7e41dee788b13d1770db3dp-52
+
+constexpr double ASIN_COEFFS[12] = {
+    0x1.0000000000000p0,  0x1.5555555556dcfp-3,  0x1.3333333082e11p-4,
+    0x1.6db6dd14099edp-5, 0x1.f1c69b35bf81fp-6,  0x1.6e97194225a67p-6,
+    0x1.1babddb82ce12p-6, 0x1.d55bd078600d6p-7,  0x1.33328959e63d6p-7,
+    0x1.2b5993bda1d9bp-6, -0x1.806aff270bf25p-7, 0x1.02614e5ed3936p-5,
+};
+
+LIBC_INLINE double asin_eval(double u) {
+  double u2 = u * u;
+  double c0 = fputil::multiply_add(u, ASIN_COEFFS[1], ASIN_COEFFS[0]);
+  double c1 = fputil::multiply_add(u, ASIN_COEFFS[3], ASIN_COEFFS[2]);
+  double c2 = fputil::multiply_add(u, ASIN_COEFFS[5], ASIN_COEFFS[4]);
+  double c3 = fputil::multiply_add(u, ASIN_COEFFS[7], ASIN_COEFFS[6]);
+  double c4 = fputil::multiply_add(u, ASIN_COEFFS[9], ASIN_COEFFS[8]);
+  double c5 = fputil::multiply_add(u, ASIN_COEFFS[11], ASIN_COEFFS[10]);
+
+  double u4 = u2 * u2;
+  double d0 = fputil::multiply_add(u2, c1, c0);
+  double d1 = fputil::multiply_add(u2, c3, c2);
+  double d2 = fputil::multiply_add(u2, c5, c4);
----------------
lntue wrote:

The problem was that the `__builtin_*` functions used for `fma` instructions are not `constexpr`.  We could do some refactoring to make it works in `constexpr` context, but I would leave that for a separate task when needed.

https://github.com/llvm/llvm-project/pull/134401