diff options
Diffstat (limited to 'libc/src/math/generic/asinf16.cpp')
| -rw-r--r-- | libc/src/math/generic/asinf16.cpp | 121 |
1 files changed, 2 insertions, 119 deletions
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp index 518c384..af8dbfe 100644 --- a/libc/src/math/generic/asinf16.cpp +++ b/libc/src/math/generic/asinf16.cpp @@ -7,127 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/asinf16.h" -#include "hdr/errno_macros.h" -#include "hdr/fenv_macros.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/FPUtil/cast.h" -#include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/sqrt.h" -#include "src/__support/macros/optimization.h" +#include "src/__support/math/asinf16.h" namespace LIBC_NAMESPACE_DECL { -// Generated by Sollya using the following command: -// > round(pi/2, D, RN); -static constexpr float PI_2 = 0x1.921fb54442d18p0f; - -LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) { - using FPBits = fputil::FPBits<float16>; - FPBits xbits(x); - - uint16_t x_u = xbits.uintval(); - uint16_t x_abs = x_u & 0x7fff; - float xf = x; - - // |x| > 0x1p0, |x| > 1, or x is NaN. - if (LIBC_UNLIKELY(x_abs > 0x3c00)) { - // asinf16(NaN) = NaN - if (xbits.is_nan()) { - if (xbits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - - return x; - } - - // 1 < |x| <= +/-inf - fputil::raise_except_if_required(FE_INVALID); - fputil::set_errno_if_required(EDOM); - - return FPBits::quiet_nan().get_val(); - } - - float xsq = xf * xf; - - // |x| <= 0x1p-1, |x| <= 0.5 - if (x_abs <= 0x3800) { - // asinf16(+/-0) = +/-0 - if (LIBC_UNLIKELY(x_abs == 0)) - return x; - - // Exhaustive tests show that, - // for |x| <= 0x1.878p-9, when: - // x > 0, and rounding upward, or - // x < 0, and rounding downward, then, - // asin(x) = x * 2^-11 + x - // else, in other rounding modes, - // asin(x) = x - if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) { - int rounding = fputil::quick_get_round(); - - if ((xbits.is_pos() && rounding == FE_UPWARD) || - (xbits.is_neg() && rounding == FE_DOWNWARD)) - return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf)); - return x; - } - - // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: - // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); - float result = - fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f, - 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); - return fputil::cast<float16>(xf * result); - } - - // When |x| > 0.5, assume that 0.5 < |x| <= 1, - // - // Step-by-step range-reduction proof: - // 1: Let y = asin(x), such that, x = sin(y) - // 2: From complimentary angle identity: - // x = sin(y) = cos(pi/2 - y) - // 3: Let z = pi/2 - y, such that x = cos(z) - // 4: From double angle formula; cos(2A) = 1 - sin^2(A): - // z = 2A, z/2 = A - // cos(z) = 1 - 2 * sin^2(z/2) - // 5: Make sin(z/2) subject of the formula: - // sin(z/2) = sqrt((1 - cos(z))/2) - // 6: Recall [3]; x = cos(z). Therefore: - // sin(z/2) = sqrt((1 - x)/2) - // 7: Let u = (1 - x)/2 - // 8: Therefore: - // asin(sqrt(u)) = z/2 - // 2 * asin(sqrt(u)) = z - // 9: Recall [3], z = pi/2 - y. Therefore: - // y = pi/2 - z - // y = pi/2 - 2 * asin(sqrt(u)) - // 10: Recall [1], y = asin(x). Therefore: - // asin(x) = pi/2 - 2 * asin(sqrt(u)) - // - // WHY? - // 11: Recall [7], u = (1 - x)/2 - // 12: Since 0.5 < x <= 1, therefore: - // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5 - // - // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for - // Step [10] as `sqrt(u)` is in range. - - // 0x1p-1 < |x| <= 0x1p0, 0.5 < |x| <= 1.0 - float xf_abs = (xf < 0 ? -xf : xf); - float sign = (xbits.uintval() >> 15 == 1 ? -1.0 : 1.0); - float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f); - float u_sqrt = fputil::sqrt<float>(u); - - // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: - // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); - float asin_sqrt_u = - u_sqrt * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f, - 0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); - - return fputil::cast<float16>(sign * - fputil::multiply_add(-2.0f, asin_sqrt_u, PI_2)); -} +LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) { return math::asinf16(x); } } // namespace LIBC_NAMESPACE_DECL |