[libc-commits] [libc] [libc][math] Implement double precision asin correctly rounded for all rounding modes. (PR #134401)
via libc-commits
libc-commits at lists.llvm.org
Mon Apr 14 16:11:42 PDT 2025
================
@@ -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);
+
+ return fputil::polyeval(u4, d0, d1, d2);
+}
+
+#else
+
+// The Taylor expansion of asin(x) around 0 is:
+// asin(x) = x + x^3/6 + 3x^5/40 + ...
+// ~ x * P(x^2).
+// Let u = x^2, then P(x^2) = P(u), and |x| = sqrt(u). Note that when
+// |x| <= 0.5, we have |u| <= 0.25.
+// We approximate P(u) by breaking it down by performing range reduction mod
+// 2^-5 = 1/32.
+// So for:
+// k = round(u * 32),
+// y = u - k/32,
+// we have that:
+// x = sqrt(u) = sqrt(k/32 + y),
+// |y| <= 2^-5 = 1/32,
+// and:
+// P(u) = P(k/32 + y) = Q_k(y).
+// Hence :
+// asin(x) = sqrt(k/32 + y) * Q_k(y),
+// Or equivalently:
+// Q_k(y) = asin(sqrt(k/32 + y)) / sqrt(k/32 + y).
+// We generate the coefficients of Q_k by Sollya as following:
+// > procedure ASIN_APPROX(N, Deg) {
+// abs_error = 0;
+// rel_error = 0;
+// deg = [||];
+// for i from 2 to Deg do deg = deg :. i;
+// for i from 1 to N/4 do {
+// F = asin(sqrt(i/N + x))/sqrt(i/N + x);
+// T = taylor(F, 1, 0);
+// T_DD = roundcoefficients(T, [|DD...|]);
+// I = [-1/(2*N), 1/(2*N)];
+// Q = fpminimax(F, deg, [|D...|], I, T_DD);
+// abs_err = dirtyinfnorm(F - Q, I);
+// rel_err = dirtyinfnorm((F - Q)/x^2, I);
+// if (abs_err > abs_error) then abs_error = abs_err;
+// if (rel_err > rel_error) then rel_error = rel_err;
+// d0 = D(coeff(Q, 0));
+// d1 = coeff(Q, 0) - d0;
+// write("{", d0, ", ", d1);
+// d0 = D(coeff(Q, 1)); d1 = coeff(Q, 1) - d0; write(", ", d0, ", ", d1);
+// for j from 2 to Deg do {
+// write(", ", coeff(Q, j));
+// };
+// print("},");
+// };
+// print("Absolute Errors:", D(abs_error));
+// print("Relative Errors:", D(rel_error));
+// };
+// > ASIN_APPROX(32, 9);
+// Absolute Errors: 0x1.69837b5183654p-72
+// Relative Errors: 0x1.4d7f82835bf64p-55
+
+// For k = 0, we use the degree-18 Taylor polynomial of asin(x)/x:
+//
+// > P = 1 + x^2 * DD(1/6) + x^4 * D(3/40) + x^6 * D(5/112) + x^8 * D(35/1152) +
+// x^10 * D(63/2816) + x^12 * D(231/13312) + x^14 * D(143/10240) +
+// x^16 * D(6435/557056) + x^18 * D(12155/1245184);
+// > dirtyinfnorm(asin(x)/x - P, [-1/64, 1/64]);
+// 0x1.999075402cafp-83
+
+constexpr double ASIN_COEFFS[9][12] = {
+ {1.0, 0.0, 0x1.5555555555555p-3, 0x1.5555555555555p-57,
+ 0x1.3333333333333p-4, 0x1.6db6db6db6db7p-5, 0x1.f1c71c71c71c7p-6,
+ 0x1.6e8ba2e8ba2e9p-6, 0x1.1c4ec4ec4ec4fp-6, 0x1.c99999999999ap-7,
+ 0x1.7a87878787878p-7, 0x1.3fde50d79435ep-7},
+ {0x1.015a397cf0f1cp0, -0x1.eebd6ccfe3ee3p-55, 0x1.5f3581be7b08bp-3,
+ -0x1.5df80d0e7237dp-57, 0x1.4519ddf1ae53p-4, 0x1.8eb4b6eeb1696p-5,
+ 0x1.17bc85420fec8p-5, 0x1.a8e39b5dcad81p-6, 0x1.53f8df127539bp-6,
+ 0x1.1a485a0b0130ap-6, 0x1.e20e6e493002p-7, 0x1.a466a7030f4c9p-7},
+ {0x1.02be9ce0b87cdp0, 0x1.e5d09da2e0f04p-56, 0x1.69ab5325bc359p-3,
+ -0x1.92f480cfede2dp-57, 0x1.58a4c3097aab1p-4, 0x1.b3db36068dd8p-5,
+ 0x1.3b9482184625p-5, 0x1.eedc823765d21p-6, 0x1.98e35d756be6bp-6,
+ 0x1.5ea4f1b32731ap-6, 0x1.355115764148ep-6, 0x1.16a5853847c91p-6},
+ {0x1.042dc6a65ffbfp0, -0x1.c7ea28dce95d1p-55, 0x1.74c4bd7412f9dp-3,
+ 0x1.447024c0a3c87p-58, 0x1.6e09c6d2b72b9p-4, 0x1.ddd9dcdae5315p-5,
+ 0x1.656f1f64058b8p-5, 0x1.21a42e4437101p-5, 0x1.eed0350b7edb2p-6,
+ 0x1.b6bc877e58c52p-6, 0x1.903a0872eb2a4p-6, 0x1.74da839ddd6d8p-6},
+ {0x1.05a8621feb16bp0, -0x1.e5b33b1407c5fp-56, 0x1.809186c2e57ddp-3,
+ -0x1.3dcb4d6069407p-60, 0x1.8587d99442dc5p-4, 0x1.06c23d1e75be3p-4,
+ 0x1.969024051c67dp-5, 0x1.54e4f934aacfdp-5, 0x1.2d60a732dbc9cp-5,
+ 0x1.149f0c046eac7p-5, 0x1.053a56dba1fbap-5, 0x1.f7face3343992p-6},
+ {0x1.072f2b6f1e601p0, -0x1.2dcbb0541997p-54, 0x1.8d2397127aebap-3,
+ 0x1.ead0c497955fbp-57, 0x1.9f68df88da518p-4, 0x1.21ee26a5900d7p-4,
+ 0x1.d08e7081b53a9p-5, 0x1.938dd661713f7p-5, 0x1.71b9f299b72e6p-5,
+ 0x1.5fbc7d2450527p-5, 0x1.58573247ec325p-5, 0x1.585a174a6a4cep-5},
+ {0x1.08c2f1d638e4cp0, 0x1.b47c159534a3dp-56, 0x1.9a8f592078624p-3,
+ -0x1.ea339145b65cdp-57, 0x1.bc04165b57aabp-4, 0x1.410df5f58441dp-4,
+ 0x1.0ab6bdf5f8f7p-4, 0x1.e0b92eea1fce1p-5, 0x1.c9094e443a971p-5,
+ 0x1.c34651d64bc74p-5, 0x1.caa008d1af08p-5, 0x1.dc165bc0c4fc5p-5},
+ {0x1.0a649a73e61f2p0, 0x1.74ac0d817e9c7p-55, 0x1.a8ec30dc9389p-3,
+ -0x1.8ab1c0eef300cp-59, 0x1.dbc11ea95061bp-4, 0x1.64e371d661328p-4,
+ 0x1.33e0023b3d895p-4, 0x1.2042269c243cep-4, 0x1.1cce74bda223p-4,
+ 0x1.244d425572ce9p-4, 0x1.34d475c7f1e3ep-4, 0x1.4d4e653082ad3p-4},
+ {0x1.0c152382d7366p0, -0x1.ee6913347c2a6p-54, 0x1.b8550d62bfb6dp-3,
+ -0x1.d10aec3f116d5p-57, 0x1.ff1bde0fa3cap-4, 0x1.8e5f3ab69f6a4p-4,
+ 0x1.656be8b6527cep-4, 0x1.5c39755dc041ap-4, 0x1.661e6ebd40599p-4,
+ 0x1.7ea3dddee2a4fp-4, 0x1.a4f439abb4869p-4, 0x1.d9181c0fda658p-4},
+};
+
+// We calculate the lower part of the approximation P(u).
+LIBC_INLINE DoubleDouble asin_eval(const DoubleDouble &u, unsigned &idx,
+ double &err) {
+ using fputil::multiply_add;
+ // k = round(u * 64).
+ double k = fputil::nearest_integer(u.hi * 0x1.0p5);
+ idx = static_cast<unsigned>(k);
+ // y = u - k/64.
+ double y_hi = multiply_add(k, -0x1.0p-5, u.hi); // Exact
+ DoubleDouble y = fputil::exact_add(y_hi, u.lo);
+ double y2 = y.hi * y.hi;
+ err *= y2 + 0x1.0p-30;
+ DoubleDouble c0 = fputil::quick_mult(
+ y, DoubleDouble{ASIN_COEFFS[idx][3], ASIN_COEFFS[idx][2]});
+ double c1 = multiply_add(y.hi, ASIN_COEFFS[idx][5], ASIN_COEFFS[idx][4]);
+ double c2 = multiply_add(y.hi, ASIN_COEFFS[idx][7], ASIN_COEFFS[idx][6]);
+ double c3 = multiply_add(y.hi, ASIN_COEFFS[idx][9], ASIN_COEFFS[idx][8]);
+ double c4 = multiply_add(y.hi, ASIN_COEFFS[idx][11], ASIN_COEFFS[idx][10]);
+
+ double y4 = y2 * y2;
+ double d0 = multiply_add(y2, c2, c1);
+ double d1 = multiply_add(y2, c4, c3);
+
+ DoubleDouble r = fputil::exact_add(ASIN_COEFFS[idx][0], c0.hi);
+
+ double e1 = multiply_add(y4, d1, d0);
+
+ r.lo = multiply_add(y2, e1, ASIN_COEFFS[idx][1] + c0.lo + r.lo);
+
+ return r;
+}
+
+// Follow the discussion above, we generate the coefficients of Q_k by Sollya as
+// following:
+// > procedure PRINTF128(a) {
+// write("{");
+// if (a < 0)
+// then write("Sign::NEG, ") else write("Sign::POS, ");
+// a_exp = floor(log2(a)) + 1;
+// write((a + 2 ^ a_exp) * 2 ^ -128);
+// print("},");
+// };
+// > verbosity = 0;
+// > procedure ASIN_APPROX(N, Deg) {
+// abs_error = 0;
+// rel_error = 0;
+// for
+// i from 1 to N / 4 do {
+// Q = fpminimax(asin(sqrt(i / N + x)) / sqrt(i / N + x), Deg,
+// [| 128... | ], [ -1 / (2 * N), 1 / (2 * N) ]);
+// abs_err = dirtyinfnorm(asin(sqrt(i / N + x)) - sqrt(i / N + x) * Q,
+// [ -1 / (2 * N), 1 / (2 * N) ]);
+// rel_err = dirtyinfnorm(asin(sqrt(i / N + x)) / sqrt(i / N + x) - Q,
+// [ -1 / (2 * N), 1 / (2 * N) ]);
+// if (abs_err > abs_error)
+// then abs_error = abs_err;
+// if (rel_err > rel_error)
+// then rel_error = rel_err;
+// write("{");
+// for
+// j from 0 to Deg do PRINTF128(coeff(Q, j));
+// print("},");
+// };
----------------
lntue wrote:
Done.
https://github.com/llvm/llvm-project/pull/134401
More information about the libc-commits
mailing list