148 lines · c
1//===-- Implementation header for acosf -------------------------*- C++ -*-===//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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ACOSF_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSF_H11 12#include "inv_trigf_utils.h"13#include "src/__support/FPUtil/FEnvImpl.h"14#include "src/__support/FPUtil/FPBits.h"15#include "src/__support/FPUtil/except_value_utils.h"16#include "src/__support/FPUtil/multiply_add.h"17#include "src/__support/FPUtil/sqrt.h"18#include "src/__support/macros/config.h"19#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY20 21namespace LIBC_NAMESPACE_DECL {22 23namespace math {24 25namespace acosf_internal {26 27#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS28 29static constexpr size_t N_EXCEPTS = 4;30 31// Exceptional values when |x| <= 0.532static constexpr fputil::ExceptValues<float, N_EXCEPTS> ACOSF_EXCEPTS = {{33 // (inputs, RZ output, RU offset, RD offset, RN offset)34 // x = 0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)35 {0x328885a3, 0x3fc90fda, 1, 0, 1},36 // x = -0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)37 {0xb28885a3, 0x3fc90fda, 1, 0, 1},38 // x = 0x1.04c444p-12, acosf(x) = 0x1.920f68p0 (RZ)39 {0x39826222, 0x3fc907b4, 1, 0, 1},40 // x = -0x1.04c444p-12, acosf(x) = 0x1.923p0 (RZ)41 {0xb9826222, 0x3fc91800, 1, 0, 1},42}};43 44#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS45 46} // namespace acosf_internal47 48LIBC_INLINE static constexpr float acosf(float x) {49 using namespace acosf_internal;50 using namespace inv_trigf_utils_internal;51 using FPBits = typename fputil::FPBits<float>;52 53 FPBits xbits(x);54 uint32_t x_uint = xbits.uintval();55 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;56 uint32_t x_sign = x_uint >> 31;57 58 // |x| <= 0.559 if (LIBC_UNLIKELY(x_abs <= 0x3f00'0000U)) {60 // |x| < 0x1p-1061 if (LIBC_UNLIKELY(x_abs < 0x3a80'0000U)) {62 // When |x| < 2^-10, we use the following approximation:63 // acos(x) = pi/2 - asin(x)64 // ~ pi/2 - x - x^3 / 665 66#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS67 // Check for exceptional values68 if (auto r = ACOSF_EXCEPTS.lookup(x_uint); LIBC_UNLIKELY(r.has_value()))69 return r.value();70#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS71 72 double xd = static_cast<double>(x);73 return static_cast<float>(fputil::multiply_add(74 -0x1.5555555555555p-3 * xd, xd * xd, M_MATH_PI_2 - xd));75 }76 77 // For |x| <= 0.5, we approximate acosf(x) by:78 // acos(x) = pi/2 - asin(x) = pi/2 - x * P(x^2)79 // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating80 // asin(x)/x on [0, 0.5] generated by Sollya with:81 // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],82 // [|1, D...|], [0, 0.5]);83 double xd = static_cast<double>(x);84 double xsq = xd * xd;85 double x3 = xd * xsq;86 double r = asin_eval(xsq);87 return static_cast<float>(fputil::multiply_add(-x3, r, M_MATH_PI_2 - xd));88 }89 90 // |x| >= 1, return 0, 2pi, or NaNs.91 if (LIBC_UNLIKELY(x_abs >= 0x3f80'0000U)) {92 if (x_abs == 0x3f80'0000U)93 return x_sign ? /* x == -1.0f */ fputil::round_result_slightly_down(94 0x1.921fb6p+1f)95 : /* x == 1.0f */ 0.0f;96 97 if (xbits.is_signaling_nan()) {98 fputil::raise_except_if_required(FE_INVALID);99 return FPBits::quiet_nan().get_val();100 }101 102 // |x| <= +/-inf103 if (x_abs <= 0x7f80'0000U) {104 fputil::set_errno_if_required(EDOM);105 fputil::raise_except_if_required(FE_INVALID);106 }107 108 return x + FPBits::quiet_nan().get_val();109 }110 111 // When 0.5 < |x| < 1, we perform range reduction as follow:112 //113 // Assume further that 0.5 < x <= 1, and let:114 // y = acos(x)115 // We use the double angle formula:116 // x = cos(y) = 1 - 2 sin^2(y/2)117 // So:118 // sin(y/2) = sqrt( (1 - x)/2 )119 // And hence:120 // y = 2 * asin( sqrt( (1 - x)/2 ) )121 // Let u = (1 - x)/2, then122 // acos(x) = 2 * asin( sqrt(u) )123 // Moreover, since 0.5 < x <= 1,124 // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,125 // And hence we can reuse the same polynomial approximation of asin(x) when126 // |x| <= 0.5:127 // acos(x) ~ 2 * sqrt(u) * P(u).128 //129 // When -1 < x <= -0.5, we use the identity:130 // acos(x) = pi - acos(-x)131 // which is reduced to the postive case.132 133 xbits.set_sign(Sign::POS);134 double xd = static_cast<double>(xbits.get_val());135 double u = fputil::multiply_add(-0.5, xd, 0.5);136 double cv = 2 * fputil::sqrt<double>(u);137 138 double r3 = asin_eval(u);139 double r = fputil::multiply_add(cv * u, r3, cv);140 return static_cast<float>(x_sign ? M_MATH_PI - r : r);141}142 143} // namespace math144 145} // namespace LIBC_NAMESPACE_DECL146 147#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOS_H148