1//===-- Implementation header for acoshf ------------------------*- 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-exception
6//
7//===----------------------------------------------------------------------===//
8
9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H
11
12#include "acoshf_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_UNLIKELY
20
21namespace LIBC_NAMESPACE_DECL {
22
23namespace math {
24
25LIBC_INLINE float acoshf(float x) {
26 using namespace acoshf_internal;
27 using FPBits_t = typename fputil::FPBits<float>;
28 FPBits_t xbits(x);
29
30 if (LIBC_UNLIKELY(x <= 1.0f)) {
31 if (x == 1.0f)
32 return 0.0f;
33 // x < 1.
34 fputil::set_errno_if_required(EDOM);
35 fputil::raise_except_if_required(FE_INVALID);
36 return FPBits_t::quiet_nan().get_val();
37 }
38
39 uint32_t x_u = xbits.uintval();
40 double x_d = static_cast<double>(x);
41
42 if (LIBC_UNLIKELY(x_u >= 0x4580'0000U)) {
43 // x >= 2^12.
44 if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
45 if (xbits.is_signaling_nan()) {
46 fputil::raise_except_if_required(FE_INVALID);
47 return FPBits_t::quiet_nan().get_val();
48 }
49 return x;
50 }
51
52 // acosh(x) = log(x + sqrt(x^2 - 1))
53 // For large x:
54 // log(x + sqrt(x^2 - 1)) = log(2x) + log((x + sqrt(x^2 - 1)) / (2x)).
55 // Let U = (x + sqrt(x^2 - 1))/(2x).
56 // Then U = 1 - (x - sqrt(x^2 - 1))/(2x)
57 // = 1 - (1 - sqrt(1 - 1/x^2))/2
58 // = 1 - (1/2) * (1/(2x^2) + 1/(8x^4) + ...)
59 // = 1 - 1/(2x)^2 - 1/(2x)^4 - ...
60 // Hence log(U) = log(1 - 1/(2x)^2 - 1/(2x)^4 - ...)
61 // = -(1/(2x)^2 - 1/(2x)^4 - ...) -
62 // - (1/(2x)^2 - 1/(2x)^4 - ...)^2/2 - ...
63 // ~ -1/(2x)^2 - 1/(2x^4) - ...
64 // For x >= 2^12:
65 // acosh(x) ~ log(2x) - 1/(2x)^2.
66 // > g = log(2*x) + 1/(4 * x^2);
67 // > dirtyinfnorm((acosh(x) - g)/acosh(x), [2^12, 2^20]);
68 // 0x1.54eb81b0c0df3c9bf68c149748e507fa136e2294fp-55
69 //
70 // For x >= 2^26, 1/(2x)^2 <= 2^-54. So we just need log(2x).
71
72 double y = 2.0 * x_d;
73
74 if (x_u <= 0x4c80'0000U) {
75 // x <= 2^26
76#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
77 if (LIBC_UNLIKELY(x_u == 0x45dc'6414U)) // x = 0x1.b8c828p12f
78 return fputil::round_result_slightly_up(value_rn: 0x1.31bcb6p3f);
79#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
80 double y_inv = 0.5 / x_d;
81 return static_cast<float>(
82 fputil::multiply_add(x: y_inv, y: -y_inv, z: log_eval(x: y)));
83
84 } else {
85// x > 2^26
86#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
87 switch (x_u) {
88 case 0x4c803f2c: // x = 0x1.007e58p26f
89 return fputil::round_result_slightly_down(value_rn: 0x1.2b786cp4f);
90 case 0x4f8ffb03: // x = 0x1.1ff606p32f
91 return fputil::round_result_slightly_up(value_rn: 0x1.6fdd34p4f);
92 case 0x5c569e88: // x = 0x1.ad3d1p57f
93 return fputil::round_result_slightly_up(value_rn: 0x1.45c146p5f);
94 case 0x5e68984e: // x = 0x1.d1309cp61f
95 return fputil::round_result_slightly_up(value_rn: 0x1.5c9442p5f);
96 case 0x655890d3: // x = 0x1.b121a6p75f
97 return fputil::round_result_slightly_down(value_rn: 0x1.a9a3f2p5f);
98 case 0x6eb1a8ec: // x = 0x1.6351d8p94f
99 return fputil::round_result_slightly_down(value_rn: 0x1.08b512p6f);
100 case 0x7997f30a: // x = 0x1.2fe614p116f
101 return fputil::round_result_slightly_up(value_rn: 0x1.451436p6f);
102#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
103 case 0x65de7ca6: // x = 0x1.bcf94cp76f
104 return fputil::round_result_slightly_up(value_rn: 0x1.af66cp5f);
105 case 0x7967ec37: // x = 0x1.cfd86ep115f
106 return fputil::round_result_slightly_up(value_rn: 0x1.43ff6ep6f);
107#endif // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
108 }
109#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
110 return static_cast<float>(log_eval(x: y));
111 }
112 }
113
114 // For 1 < x < 2^12, we use the formula:
115 // acosh(x) = log(x + sqrt(x^2 - 1))
116 return static_cast<float>(log_eval(
117 x: x_d + fputil::sqrt<double>(x: fputil::multiply_add(x: x_d, y: x_d, z: -1.0))));
118}
119
120} // namespace math
121
122} // namespace LIBC_NAMESPACE_DECL
123
124#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H
125