libc/src/math/generic/asinf.cpp

   1 //===-- Single-precision asin function ------------------------------------===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 #include "src/math/asinf.h"
  10 #include "src/__support/FPUtil/FEnvImpl.h"
  11 #include "src/__support/FPUtil/FPBits.h"
  12 #include "src/__support/FPUtil/PolyEval.h"
  13 #include "src/__support/FPUtil/except_value_utils.h"
  14 #include "src/__support/FPUtil/multiply_add.h"
  15 #include "src/__support/FPUtil/sqrt.h"
  16 #include "src/__support/macros/optimization.h"            // LIBC_UNLIKELY
  17 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
  18
  19 #include <errno.h>
  20
  21 #include "inv_trigf_utils.h"
  22
  23 namespace __llvm_libc {
  24
  25 static constexpr size_t N_EXCEPTS = 2;
  26
  27 // Exceptional values when |x| <= 0.5
  28 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{
  29     // (inputs, RZ output, RU offset, RD offset, RN offset)
  30     // x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)
  31     {0x3d09bf86, 0x3d09c62c, 1, 0, 1},
  32     // x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)
  33     {0x3de5fa1e, 0x3de6768e, 1, 0, 0},
  34 }};
  35
  36 // Exceptional values when 0.5 < |x| <= 1
  37 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{
  38     // (inputs, RZ output, RU offset, RD offset, RN offset)
  39     // x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)
  40     {0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},
  41     // x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)
  42     {0x3f7741b6, 0x3fa7832a, 1, 0, 0},
  43 }};
  44
  45 LLVM_LIBC_FUNCTION(float, asinf, (float x)) {
  46   using FPBits = typename fputil::FPBits<float>;
  47   FPBits xbits(x);
  48   uint32_t x_uint = xbits.uintval();
  49   uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
  50   constexpr double SIGN[2] = {1.0, -1.0};
  51   uint32_t x_sign = x_uint >> 31;
  52
  53   // |x| <= 0.5-ish
  54   if (x_abs < 0x3f04'471dU) {
  55     // |x| < 0x1.d12edp-12
  56     if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) {
  57       // When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
  58       // is:
  59       //   |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
  60       //                             = x^2 / 6
  61       //                             < 2^-25
  62       //                             < epsilon(1)/2.
  63       // So the correctly rounded values of asin(x) are:
  64       //   = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
  65       //                        or (rounding mode = FE_UPWARD and x is
  66       //                        negative),
  67       //   = x otherwise.
  68       // To simplify the rounding decision and make it more efficient, we use
  69       //   fma(x, 2^-25, x) instead.
  70       // An exhaustive test shows that this formula work correctly for all
  71       // rounding modes up to |x| < 0x1.d12edp-12.
  72       // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
  73       // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
  74       // |x| < 2^-125. For targets without FMA instructions, we simply use
  75       // double for intermediate results as it is more efficient than using an
  76       // emulated version of FMA.
  77 #if defined(LIBC_TARGET_CPU_HAS_FMA)
  78       return fputil::multiply_add(x, 0x1.0p-25f, x);
  79 #else
  80       double xd = static_cast<double>(x);
  81       return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
  82 #endif // LIBC_TARGET_CPU_HAS_FMA
  83     }
  84
  85     // Check for exceptional values
  86     if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);
  87         LIBC_UNLIKELY(r.has_value()))
  88       return r.value();
  89
  90     // For |x| <= 0.5, we approximate asinf(x) by:
  91     //   asin(x) = x * P(x^2)
  92     // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
  93     // asin(x)/x on [0, 0.5] generated by Sollya with:
  94     // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
  95     //                 [|1, D...|], [0, 0.5]);
  96     // An exhaustive test shows that this approximation works well up to a
  97     // little more than 0.5.
  98     double xd = static_cast<double>(x);
  99     double xsq = xd * xd;
 100     double x3 = xd * xsq;
 101     double r = asin_eval(xsq);
 102     return static_cast<float>(fputil::multiply_add(x3, r, xd));
 103   }
 104
 105   // |x| > 1, return NaNs.
 106   if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
 107     if (x_abs <= 0x7f80'0000U) {
 108       fputil::set_errno_if_required(EDOM);
 109       fputil::raise_except_if_required(FE_INVALID);
 110     }
 111     return x +
 112            FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
 113   }
 114
 115   // Check for exceptional values
 116   if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);
 117       LIBC_UNLIKELY(r.has_value()))
 118     return r.value();
 119
 120   // When |x| > 0.5, we perform range reduction as follow:
 121   //
 122   // Assume further that 0.5 < x <= 1, and let:
 123   //   y = asin(x)
 124   // We will use the double angle formula:
 125   //   cos(2y) = 1 - 2 sin^2(y)
 126   // and the complement angle identity:
 127   //   x = sin(y) = cos(pi/2 - y)
 128   //              = 1 - 2 sin^2 (pi/4 - y/2)
 129   // So:
 130   //   sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
 131   // And hence:
 132   //   pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
 133   // Equivalently:
 134   //   asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
 135   // Let u = (1 - x)/2, then:
 136   //   asin(x) = pi/2 - 2 * asin( sqrt(u) )
 137   // Moreover, since 0.5 < x <= 1:
 138   //   0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
 139   // And hence we can reuse the same polynomial approximation of asin(x) when
 140   // |x| <= 0.5:
 141   //   asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
 142
 143   xbits.set_sign(false);
 144   double sign = SIGN[x_sign];
 145   double xd = static_cast<double>(xbits.get_val());
 146   double u = fputil::multiply_add(-0.5, xd, 0.5);
 147   double c1 = sign * (-2 * fputil::sqrt(u));
 148   double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1);
 149   double c3 = c1 * u;
 150
 151   double r = asin_eval(u);
 152   return static_cast<float>(fputil::multiply_add(c3, r, c2));
 153 }
 154
 155 } // namespace __llvm_libc