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
|
//===-- 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/optimization.h" // LIBC_UNLIKELY
#include <errno.h>
#include "inv_trigf_utils.h"
namespace __llvm_libc {
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},
}};
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
// Check for exceptional values
if (auto r = ACOSF_EXCEPTS.lookup(x_uint); LIBC_UNLIKELY(r.has_value()))
return r.value();
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 NaNs.
if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
if (x_abs <= 0x7f80'0000U) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
}
return x + FPBits::build_quiet_nan(0);
}
// 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(false);
double xd = static_cast<double>(xbits.get_val());
double u = fputil::multiply_add(-0.5, xd, 0.5);
double cv = 2 * fputil::sqrt(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 __llvm_libc
|
