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

Code before:
  //
  // 0.318309886183791*x + 0.0530516476972984*x**3 + 0.0238732414637843*x**5 +
  // 0.0142102627760621*x**7 + 0.00967087327815336*x**9 +
  // 0.00712127941391293*x**11 + 0.00552355646848375*x**13 +
  // 0.00444514782463692*x**15 + 0.00367705242846804*x**17 +
  // 0.00310721681820837*x**19 + 0.00267072683660291*x**21 +
  // 0.00232764927854127*x**23 + O(x**25)
  constexpr double ASINPI_POLY_COEFFS[] = {
      0x1.45f306dc9c883p-2, // 0x3fd45f306dc9c883
      0x1.b2995e7b7b604p-5, // 0x3fab2995e7b7b604
      0x1.8723a1d588a37p-6, // 0x3f98723a1d588a37
      0x1.d1a452f204304p-7, // 0x3f8d1a452f204304
      0x1.3ce52a3a09f63p-7, // 0x3f83ce52a3a09f63
      0x1.d2b33e303d376p-8, // 0x3f7d2b33e303d376
      0x1.69fde663c674cp-8, // 0x3f769fde663c674c
      0x1.235134885f19cp-8, // 0x3f7235134885f19c
      0x1.e1f567da15ce4p-9, // 0x3f6e1f567da15ce4
      0x1.9744e53b48525p-9, // 0x3f69744e53b48525
      0x1.5e0eb8d6eaa83p-9, // 0x3f65e0eb8d6eaa83
      0x1.3116f30e47f5dp-9, // 0x3f63116f30e47f5d
      0x1.0cfcc45fca89cp-9  // 0x3f60cfcc45fca89c
  };

  // polynomial evaluation using horner's method
  // work only for |x| in [0, 0.5]
  auto asinpi_polyeval = [&](double v) -> double {
    return v * fputil::polyeval(v * v, 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]);
  };

  // if |x| <= 0.5:
  if (LIBC_UNLIKELY(x_abs <= 0.5)) {
    double result = asinpi_polyeval(fputil::cast<double>(x));
    return fputil::cast<float>(result);
  }

  // If |x| > 0.5, we need to use the range reduction method:
  //    y = asin(x) => x = sin(y)
  //      because: sin(a) = cos(pi/2 - a)
  //      therefore:
  //    x = cos(pi/2 - y)
  //      let z = pi/2 - y,
  //    x = cos(z)
  //      because: cos(2a) = 1 - 2 * sin^2(a), z = 2a, a = z/2
  //      therefore:
  //    cos(z) = 1 - 2 * sin^2(z/2)
  //    sin(z/2) = sqrt((1 - cos(z))/2)
  //    sin(z/2) = sqrt((1 - x)/2)

Review comment (suggestion):
```suggestion
  auto asinpi_polyeval = [&](double v) constexpr -> double {
```

Code after:

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
  auto r = ASINPIF_EXCEPTS.lookup(xbits.uintval());
  if (LIBC_UNLIKELY(r.has_value()))
    return r.value();
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS

  // the coefficients for the polynomial approximation of asin(x)/(pi*x) in the
  // range [0, 0.5] extracted using Sollya.
  //
  // Sollya code:
  // > prec = 200;
  // > display = hexadecimal;
  // > 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
  };
  // 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]);