165 lines · c
1//===-- Implementation header for acosf16 -----------------------*- 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_ACOSF16_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSF16_H11 12#include "include/llvm-libc-macros/float16-macros.h"13 14#ifdef LIBC_TYPES_HAS_FLOAT1615 16#include "src/__support/FPUtil/FEnvImpl.h"17#include "src/__support/FPUtil/FPBits.h"18#include "src/__support/FPUtil/PolyEval.h"19#include "src/__support/FPUtil/cast.h"20#include "src/__support/FPUtil/except_value_utils.h"21#include "src/__support/FPUtil/multiply_add.h"22#include "src/__support/FPUtil/sqrt.h"23#include "src/__support/macros/optimization.h"24 25namespace LIBC_NAMESPACE_DECL {26 27namespace math {28 29LIBC_INLINE static constexpr float16 acosf16(float16 x) {30 31 // Generated by Sollya using the following command:32 // > round(pi/2, SG, RN);33 // > round(pi, SG, RN);34 constexpr float PI_OVER_2 = 0x1.921fb6p0f;35 constexpr float PI = 0x1.921fb6p1f;36 37#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS38 constexpr size_t N_EXCEPTS = 2;39 40 constexpr fputil::ExceptValues<float16, N_EXCEPTS> ACOSF16_EXCEPTS{{41 // (input, RZ output, RU offset, RD offset, RN offset)42 {0xacaf, 0x3e93, 1, 0, 0},43 {0xb874, 0x4052, 1, 0, 1},44 }};45#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS46 47 using FPBits = fputil::FPBits<float16>;48 FPBits xbits(x);49 50 uint16_t x_u = xbits.uintval();51 uint16_t x_abs = x_u & 0x7fff;52 uint16_t x_sign = x_u >> 15;53 54 // |x| > 0x1p0, |x| > 1, or x is NaN.55 if (LIBC_UNLIKELY(x_abs > 0x3c00)) {56 // acosf16(NaN) = NaN57 if (xbits.is_nan()) {58 if (xbits.is_signaling_nan()) {59 fputil::raise_except_if_required(FE_INVALID);60 return FPBits::quiet_nan().get_val();61 }62 63 return x;64 }65 66 // 1 < |x| <= +/-inf67 fputil::raise_except_if_required(FE_INVALID);68 fputil::set_errno_if_required(EDOM);69 70 return FPBits::quiet_nan().get_val();71 }72 73 float xf = x;74 75#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS76 // Handle exceptional values77 if (auto r = ACOSF16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))78 return r.value();79#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS80 81 // |x| == 0x1p0, x is 1 or -182 // if x is (-)1, return pi, else83 // if x is (+)1, return 084 if (LIBC_UNLIKELY(x_abs == 0x3c00))85 return fputil::cast<float16>(x_sign ? PI : 0.0f);86 87 float xsq = xf * xf;88 89 // |x| <= 0x1p-1, |x| <= 0.590 if (x_abs <= 0x3800) {91 // if x is 0, return pi/292 if (LIBC_UNLIKELY(x_abs == 0))93 return fputil::cast<float16>(PI_OVER_2);94 95 // Note that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x)96 // Degree-6 minimax polynomial of asin(x) generated by Sollya with:97 // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);98 float interm =99 fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f,100 0x1.43b2d6p-5f, 0x1.a0d73ep-5f);101 return fputil::cast<float16>(fputil::multiply_add(-xf, interm, PI_OVER_2));102 }103 104 // When |x| > 0.5, assume that 0.5 < |x| <= 1105 //106 // Step-by-step range-reduction proof:107 // 1: Let y = asin(x), such that, x = sin(y)108 // 2: From complimentary angle identity:109 // x = sin(y) = cos(pi/2 - y)110 // 3: Let z = pi/2 - y, such that x = cos(z)111 // 4: From double angle formula; cos(2A) = 1 - 2 * sin^2(A):112 // z = 2A, z/2 = A113 // cos(z) = 1 - 2 * sin^2(z/2)114 // 5: Make sin(z/2) subject of the formula:115 // sin(z/2) = sqrt((1 - cos(z))/2)116 // 6: Recall [3]; x = cos(z). Therefore:117 // sin(z/2) = sqrt((1 - x)/2)118 // 7: Let u = (1 - x)/2119 // 8: Therefore:120 // asin(sqrt(u)) = z/2121 // 2 * asin(sqrt(u)) = z122 // 9: Recall [3]; z = pi/2 - y. Therefore:123 // y = pi/2 - z124 // y = pi/2 - 2 * asin(sqrt(u))125 // 10: Recall [1], y = asin(x). Therefore:126 // asin(x) = pi/2 - 2 * asin(sqrt(u))127 // 11: Recall that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x)128 // Therefore:129 // acos(x) = pi/2 - (pi/2 - 2 * asin(sqrt(u)))130 // acos(x) = 2 * asin(sqrt(u))131 //132 // THE RANGE REDUCTION, HOW?133 // 12: Recall [7], u = (1 - x)/2134 // 13: Since 0.5 < x <= 1, therefore:135 // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5136 //137 // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for138 // Step [11] as `sqrt(u)` is in range.139 // When -1 < x <= -0.5, the identity:140 // acos(x) = pi - acos(-x)141 // allows us to compute for the negative x value (lhs)142 // with a positive x value instead (rhs).143 144 float xf_abs = (xf < 0 ? -xf : xf);145 float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f);146 float sqrt_u = fputil::sqrt<float>(u);147 148 // Degree-6 minimax polynomial of asin(x) generated by Sollya with:149 // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);150 float asin_sqrt_u =151 sqrt_u * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f,152 0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f);153 154 return fputil::cast<float16>(155 x_sign ? fputil::multiply_add(-2.0f, asin_sqrt_u, PI) : 2 * asin_sqrt_u);156}157 158} // namespace math159 160} // namespace LIBC_NAMESPACE_DECL161 162#endif // LIBC_TYPES_HAS_FLOAT16163 164#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOS_H165