Repository: llvm/llvm-project
File: libc/src/__support/math/asinpif.h (C/C++)

Code before:
  // > g = asin(x) / (pi * x);
  // > P = fpminimax(g, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
  // >              [|D...|], [0, 0.5]);
  // > for i from 0 to degree(P) do coeff(P, i);
  // > print("Error:", dirtyinfnorm(P - g, [1e-30; 0.25]));
  // Error: 0x1.45c281e1cf9b58p-50 ~= 2^−49.652
  //
  // Non-zero coefficients (even powers only):
  constexpr double ASINPI_POLY_COEFFS[] = {
      0x1.45f306dc9c881p-2,  // x^0
      0x1.b2995e7b7e756p-5,  // x^2
      0x1.8723a1d12f828p-6,  // x^4
      0x1.d1a45564b9545p-7,  // x^6
      0x1.3ce4ceaa0e1e9p-7,  // x^8
      0x1.d2c305898ea13p-8,  // x^10
      0x1.692212e27a5f9p-8,  // x^12
      0x1.2b22cc744d25bp-8,  // x^14
      0x1.8427b864479ffp-9,  // x^16
      0x1.815522d7a2bf1p-8,  // x^18
      -0x1.f6df98438aef4p-9, // x^20
      0x1.4b50c2eb13708p-7   // x^22
  };
  // polynomial evaluation using horner's method
  // work only for |x| in [0, 0.5]
  // Returns v * P(v2) where P(v2) = c0 + c1*v2 + c2*v2^2 + ...
  auto asinpi_polyeval = [&](double v, double v2) -> double {
    return v * fputil::polyeval(v2, ASINPI_POLY_COEFFS[0],
                                ASINPI_POLY_COEFFS[1], ASINPI_POLY_COEFFS[2],
                                ASINPI_POLY_COEFFS[3], ASINPI_POLY_COEFFS[4],
                                ASINPI_POLY_COEFFS[5], ASINPI_POLY_COEFFS[6],
                                ASINPI_POLY_COEFFS[7], ASINPI_POLY_COEFFS[8],
                                ASINPI_POLY_COEFFS[9], ASINPI_POLY_COEFFS[10],
                                ASINPI_POLY_COEFFS[11]);
  };

  // Returns P(v2) - c0 = c1*v2 + c2*v2^2 + ...
  // This is the "tail" of the polynomial, used to avoid cancellation
  // in the range reduction path.
  auto asinpi_polyeval_tail = [&](double v2) -> double {
    return v2 * fputil::polyeval(v2, ASINPI_POLY_COEFFS[1],
                                 ASINPI_POLY_COEFFS[2], ASINPI_POLY_COEFFS[3],
                                 ASINPI_POLY_COEFFS[4], ASINPI_POLY_COEFFS[5],
                                 ASINPI_POLY_COEFFS[6], ASINPI_POLY_COEFFS[7],
                                 ASINPI_POLY_COEFFS[8], ASINPI_POLY_COEFFS[9],
                                 ASINPI_POLY_COEFFS[10],
                                 ASINPI_POLY_COEFFS[11]);
  };

  // if |x| <= 0.5:
  if (LIBC_UNLIKELY(x_abs <= 0.5)) {
    double x_d = fputil::cast<double>(x);

Review comment (suggestion):
There is only one single call for each `asinpi_polyeval` and `asinpi_polyeval_tail`, and internally, they only differ by `v * (c0 + v2 * P1(v2))` vs `v2 * P1(v2)`.  Let just define a single one `asinpi_polyeval` that returns `P1(v2) = c1 + v2*c2 + ...`.  Then make up the difference at the call sites.  Also technically the tail can be updated to absorb `ONE_OVER_PI_LO + DELTA_C0` into a single `multiply_add`.

Code after:
  // Error: 0x1.45c281e1cf9b58p-50 ~= 2^−49.652
  //
  // Non-zero coefficients (even powers only):
  constexpr double ASINPI_POLY_COEFFS[] = {
      0x1.45f306dc9c881p-2,  // x^0
      0x1.b2995e7b7e756p-5,  // x^2
      0x1.8723a1d12f828p-6,  // x^4
      0x1.d1a45564b9545p-7,  // x^6
      0x1.3ce4ceaa0e1e9p-7,  // x^8
      0x1.d2c305898ea13p-8,  // x^10
      0x1.692212e27a5f9p-8,  // x^12
      0x1.2b22cc744d25bp-8,  // x^14
      0x1.8427b864479ffp-9,  // x^16
      0x1.815522d7a2bf1p-8,  // x^18
      -0x1.f6df98438aef4p-9, // x^20
      0x1.4b50c2eb13708p-7   // x^22
  };
  // Evaluates P1(v2) = c1 + c2*v2 + c3*v2^2 + ... (tail of P without c0)
  auto asinpi_polyeval = [&](double v2) -> double {
    return fputil::polyeval(
        v2, ASINPI_POLY_COEFFS[1], ASINPI_POLY_COEFFS[2], ASINPI_POLY_COEFFS[3],
        ASINPI_POLY_COEFFS[4], ASINPI_POLY_COEFFS[5], ASINPI_POLY_COEFFS[6],
        ASINPI_POLY_COEFFS[7], ASINPI_POLY_COEFFS[8], ASINPI_POLY_COEFFS[9],
        ASINPI_POLY_COEFFS[10], ASINPI_POLY_COEFFS[11]);
  };

  // if |x| <= 0.5:
  //   asinpi(x) = x * (c0 + x^2 * P1(x^2))
  if (LIBC_UNLIKELY(x_abs <= 0.5)) {
    double x_d = fputil::cast<double>(x);
    double v2 = x_d * x_d;
    double result = x_d * fputil::multiply_add(v2, asinpi_polyeval(v2),
                                               ASINPI_POLY_COEFFS[0]);
    return fputil::cast<float>(result);
  }

  // If |x| > 0.5:
  //   asinpi(x) = 0.5 - 2 * sqrt(u) * P(u)
  //             = 0.5 - 2 * sqrt(u) * (c0 + u * P1(u))
  //             = (0.5 - 2*sqrt(u)*ONE_OVER_PI_HI)
  //               - 2*sqrt(u) * (ONE_OVER_PI_LO + DELTA_C0 + u * P1(u))
  //
  // where u = (1 - |x|) / 2, and
  //   ONE_OVER_PI_HI + ONE_OVER_PI_LO = 1/pi to ~106 bits
  //   DELTA_C0 = c0 - ONE_OVER_PI_HI
  //
  // ONE_OVER_PI_LO + DELTA_C0 is a single precomputed constant:
  //   = ONE_OVER_PI_LO + (c0 - ONE_OVER_PI_HI)
  //   = c0 - (ONE_OVER_PI_HI - ONE_OVER_PI_LO)
  //   = c0 - 1/pi  (to ~106 bits)
  constexpr double ONE_OVER_PI_HI = 0x1.45f306dc9c883p-2;