[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:

Adding `const` here won't change anything with optimization.  It does help with "const correctness" but in this case, it will reduce the readability.  Since these intermediate variables are for improving the readability and order-of-operations of these immediate computations, I opted for readability in all of such computations / polynomial evaluations.  I think I tried `const` for these computations at some point before, and didn't like it much.  So for now, for these computations, I only reserve `const`, `constexpr` for literal constants.

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