|
| 1 | +//===-- Half-precision acospi 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/acospif16.h" |
| 10 | +#include "hdr/errno_macros.h" |
| 11 | +#include "hdr/fenv_macros.h" |
| 12 | +#include "src/__support/FPUtil/FEnvImpl.h" |
| 13 | +#include "src/__support/FPUtil/FPBits.h" |
| 14 | +#include "src/__support/FPUtil/PolyEval.h" |
| 15 | +#include "src/__support/FPUtil/cast.h" |
| 16 | +#include "src/__support/FPUtil/multiply_add.h" |
| 17 | +#include "src/__support/FPUtil/sqrt.h" |
| 18 | +#include "src/__support/macros/optimization.h" |
| 19 | + |
| 20 | +namespace LIBC_NAMESPACE_DECL { |
| 21 | + |
| 22 | +LLVM_LIBC_FUNCTION(float16, acospif16, (float16 x)) { |
| 23 | + using FPBits = fputil::FPBits<float16>; |
| 24 | + FPBits xbits(x); |
| 25 | + |
| 26 | + uint16_t x_u = xbits.uintval(); |
| 27 | + uint16_t x_abs = x_u & 0x7fff; |
| 28 | + uint16_t x_sign = x_u >> 15; |
| 29 | + |
| 30 | + // |x| > 0x1p0, |x| > 1, or x is NaN. |
| 31 | + if (LIBC_UNLIKELY(x_abs > 0x3c00)) { |
| 32 | + // acospif16(NaN) = NaN |
| 33 | + if (xbits.is_nan()) { |
| 34 | + if (xbits.is_signaling_nan()) { |
| 35 | + fputil::raise_except_if_required(FE_INVALID); |
| 36 | + return FPBits::quiet_nan().get_val(); |
| 37 | + } |
| 38 | + |
| 39 | + return x; |
| 40 | + } |
| 41 | + |
| 42 | + // 1 < |x| <= +inf |
| 43 | + fputil::raise_except_if_required(FE_INVALID); |
| 44 | + fputil::set_errno_if_required(EDOM); |
| 45 | + |
| 46 | + return FPBits::quiet_nan().get_val(); |
| 47 | + } |
| 48 | + |
| 49 | + // |x| == 0x1p0, x is 1 or -1 |
| 50 | + // if x is (-)1, return 1 |
| 51 | + // if x is (+)1, return 0 |
| 52 | + if (LIBC_UNLIKELY(x_abs == 0x3c00)) |
| 53 | + return fputil::cast<float16>(x_sign ? 1.0f : 0.0f); |
| 54 | + |
| 55 | + float xf = x; |
| 56 | + float xsq = xf * xf; |
| 57 | + |
| 58 | + // Degree-6 minimax polynomial coefficients of asin(x) generated by Sollya |
| 59 | + // with: > P = fpminimax(asin(x)/(pi * x), [|0, 2, 4, 6, 8|], [|SG...|], [0, |
| 60 | + // 0.5]); |
| 61 | + constexpr float POLY_COEFFS[5] = {0x1.45f308p-2f, 0x1.b2900cp-5f, |
| 62 | + 0x1.897e36p-6f, 0x1.9efafcp-7f, |
| 63 | + 0x1.06d884p-6f}; |
| 64 | + // |x| <= 0x1p-1, |x| <= 0.5 |
| 65 | + if (x_abs <= 0x3800) { |
| 66 | + // if x is 0, return 0.5 |
| 67 | + if (LIBC_UNLIKELY(x_abs == 0)) |
| 68 | + return fputil::cast<float16>(0.5f); |
| 69 | + |
| 70 | + // Note that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x), then |
| 71 | + // acospi(x) = 0.5 - asin(x)/pi |
| 72 | + float interm = |
| 73 | + fputil::polyeval(xsq, POLY_COEFFS[0], POLY_COEFFS[1], POLY_COEFFS[2], |
| 74 | + POLY_COEFFS[3], POLY_COEFFS[4]); |
| 75 | + |
| 76 | + return fputil::cast<float16>(fputil::multiply_add(-xf, interm, 0.5f)); |
| 77 | + } |
| 78 | + |
| 79 | + // When |x| > 0.5, assume that 0.5 < |x| <= 1 |
| 80 | + // |
| 81 | + // Step-by-step range-reduction proof: |
| 82 | + // 1: Let y = asin(x), such that, x = sin(y) |
| 83 | + // 2: From complimentary angle identity: |
| 84 | + // x = sin(y) = cos(pi/2 - y) |
| 85 | + // 3: Let z = pi/2 - y, such that x = cos(z) |
| 86 | + // 4: From double angle formula; cos(2A) = 1 - 2 * sin^2(A): |
| 87 | + // z = 2A, z/2 = A |
| 88 | + // cos(z) = 1 - 2 * sin^2(z/2) |
| 89 | + // 5: Make sin(z/2) subject of the formula: |
| 90 | + // sin(z/2) = sqrt((1 - cos(z))/2) |
| 91 | + // 6: Recall [3]; x = cos(z). Therefore: |
| 92 | + // sin(z/2) = sqrt((1 - x)/2) |
| 93 | + // 7: Let u = (1 - x)/2 |
| 94 | + // 8: Therefore: |
| 95 | + // asin(sqrt(u)) = z/2 |
| 96 | + // 2 * asin(sqrt(u)) = z |
| 97 | + // 9: Recall [3]; z = pi/2 - y. Therefore: |
| 98 | + // y = pi/2 - z |
| 99 | + // y = pi/2 - 2 * asin(sqrt(u)) |
| 100 | + // 10: Recall [1], y = asin(x). Therefore: |
| 101 | + // asin(x) = pi/2 - 2 * asin(sqrt(u)) |
| 102 | + // 11: Recall that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x) |
| 103 | + // Therefore: |
| 104 | + // acos(x) = pi/2 - (pi/2 - 2 * asin(sqrt(u))) |
| 105 | + // acos(x) = 2 * asin(sqrt(u)) |
| 106 | + // acospi(x) = 2 * (asin(sqrt(u)) / pi) |
| 107 | + // |
| 108 | + // THE RANGE REDUCTION, HOW? |
| 109 | + // 12: Recall [7], u = (1 - x)/2 |
| 110 | + // 13: Since 0.5 < x <= 1, therefore: |
| 111 | + // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5 |
| 112 | + // |
| 113 | + // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for |
| 114 | + // Step [11] as `sqrt(u)` is in range. |
| 115 | + // When -1 < x <= -0.5, the identity: |
| 116 | + // acos(x) = pi - acos(-x) |
| 117 | + // acospi(x) = 1 - acos(-x)/pi |
| 118 | + // allows us to compute for the negative x value (lhs) |
| 119 | + // with a positive x value instead (rhs). |
| 120 | + |
| 121 | + float xf_abs = (xf < 0 ? -xf : xf); |
| 122 | + float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f); |
| 123 | + float sqrt_u = fputil::sqrt<float>(u); |
| 124 | + |
| 125 | + float asin_sqrt_u = |
| 126 | + sqrt_u * fputil::polyeval(u, POLY_COEFFS[0], POLY_COEFFS[1], |
| 127 | + POLY_COEFFS[2], POLY_COEFFS[3], POLY_COEFFS[4]); |
| 128 | + |
| 129 | + // Same as acos(x), but devided the expression with pi |
| 130 | + return fputil::cast<float16>( |
| 131 | + x_sign ? fputil::multiply_add(-2.0f, asin_sqrt_u, 1.0f) |
| 132 | + : 2.0f * asin_sqrt_u); |
| 133 | +} |
| 134 | +} // namespace LIBC_NAMESPACE_DECL |
0 commit comments