File size: 4,683 Bytes
00df61d | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 | //===-- Single-precision acos function ------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#include "src/math/acosf.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/except_value_utils.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/sqrt.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "inv_trigf_utils.h"
namespace LIBC_NAMESPACE_DECL {
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
static constexpr size_t N_EXCEPTS = 4;
// Exceptional values when |x| <= 0.5
static constexpr fputil::ExceptValues<float, N_EXCEPTS> ACOSF_EXCEPTS = {{
// (inputs, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)
{0x328885a3, 0x3fc90fda, 1, 0, 1},
// x = -0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)
{0xb28885a3, 0x3fc90fda, 1, 0, 1},
// x = 0x1.04c444p-12, acosf(x) = 0x1.920f68p0 (RZ)
{0x39826222, 0x3fc907b4, 1, 0, 1},
// x = -0x1.04c444p-12, acosf(x) = 0x1.923p0 (RZ)
{0xb9826222, 0x3fc91800, 1, 0, 1},
}};
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
LLVM_LIBC_FUNCTION(float, acosf, (float x)) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
uint32_t x_uint = xbits.uintval();
uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
uint32_t x_sign = x_uint >> 31;
// |x| <= 0.5
if (LIBC_UNLIKELY(x_abs <= 0x3f00'0000U)) {
// |x| < 0x1p-10
if (LIBC_UNLIKELY(x_abs < 0x3a80'0000U)) {
// When |x| < 2^-10, we use the following approximation:
// acos(x) = pi/2 - asin(x)
// ~ pi/2 - x - x^3 / 6
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Check for exceptional values
if (auto r = ACOSF_EXCEPTS.lookup(x_uint); LIBC_UNLIKELY(r.has_value()))
return r.value();
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
double xd = static_cast<double>(x);
return static_cast<float>(fputil::multiply_add(
-0x1.5555555555555p-3 * xd, xd * xd, M_MATH_PI_2 - xd));
}
// For |x| <= 0.5, we approximate acosf(x) by:
// acos(x) = pi/2 - asin(x) = pi/2 - x * P(x^2)
// Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
// asin(x)/x on [0, 0.5] generated by Sollya with:
// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
// [|1, D...|], [0, 0.5]);
double xd = static_cast<double>(x);
double xsq = xd * xd;
double x3 = xd * xsq;
double r = asin_eval(xsq);
return static_cast<float>(fputil::multiply_add(-x3, r, M_MATH_PI_2 - xd));
}
// |x| >= 1, return 0, 2pi, or NaNs.
if (LIBC_UNLIKELY(x_abs >= 0x3f80'0000U)) {
if (x_abs == 0x3f80'0000U)
return x_sign ? /* x == -1.0f */ fputil::round_result_slightly_down(
0x1.921fb6p+1f)
: /* x == 1.0f */ 0.0f;
if (xbits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
// |x| <= +/-inf
if (x_abs <= 0x7f80'0000U) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
}
return x + FPBits::quiet_nan().get_val();
}
// When 0.5 < |x| < 1, we perform range reduction as follow:
//
// Assume further that 0.5 < x <= 1, and let:
// y = acos(x)
// We use the double angle formula:
// x = cos(y) = 1 - 2 sin^2(y/2)
// So:
// sin(y/2) = sqrt( (1 - x)/2 )
// And hence:
// y = 2 * asin( sqrt( (1 - x)/2 ) )
// Let u = (1 - x)/2, then
// acos(x) = 2 * asin( sqrt(u) )
// Moreover, since 0.5 < x <= 1,
// 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
// And hence we can reuse the same polynomial approximation of asin(x) when
// |x| <= 0.5:
// acos(x) ~ 2 * sqrt(u) * P(u).
//
// When -1 < x <= -0.5, we use the identity:
// acos(x) = pi - acos(-x)
// which is reduced to the postive case.
xbits.set_sign(Sign::POS);
double xd = static_cast<double>(xbits.get_val());
double u = fputil::multiply_add(-0.5, xd, 0.5);
double cv = 2 * fputil::sqrt<double>(u);
double r3 = asin_eval(u);
double r = fputil::multiply_add(cv * u, r3, cv);
return static_cast<float>(x_sign ? M_MATH_PI - r : r);
}
} // namespace LIBC_NAMESPACE_DECL
|