1//===-- Implementation header for acoshf16 ----------------------*- 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_ACOSHF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "acoshf_utils.h"
17#include "src/__support/FPUtil/FEnvImpl.h"
18#include "src/__support/FPUtil/FPBits.h"
19#include "src/__support/FPUtil/PolyEval.h"
20#include "src/__support/FPUtil/cast.h"
21#include "src/__support/FPUtil/except_value_utils.h"
22#include "src/__support/FPUtil/multiply_add.h"
23#include "src/__support/FPUtil/sqrt.h"
24#include "src/__support/macros/config.h"
25#include "src/__support/macros/optimization.h"
26
27namespace LIBC_NAMESPACE_DECL {
28
29namespace math {
30
31LIBC_INLINE constexpr float16 acoshf16(float16 x) {
32
33 using namespace acoshf_internal;
34 constexpr size_t N_EXCEPTS = 2;
35 constexpr fputil::ExceptValues<float16, N_EXCEPTS> ACOSHF16_EXCEPTS{.values: {
36 // (input, RZ output, RU offset, RD offset, RN offset)
37 // x = 0x1.6dcp+1, acoshf16(x) = 0x1.b6p+0 (RZ)
38 {.input: 0x41B7, .rnd_towardzero_result: 0x3ED8, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
39 // x = 0x1.39p+0, acoshf16(x) = 0x1.4f8p-1 (RZ)
40 {.input: 0x3CE4, .rnd_towardzero_result: 0x393E, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
41 }};
42
43 using FPBits = fputil::FPBits<float16>;
44 FPBits xbits(x);
45 uint16_t x_u = xbits.uintval();
46
47 // Check for NaN input first.
48 if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
49 if (xbits.is_signaling_nan()) {
50 fputil::raise_except_if_required(FE_INVALID);
51 return FPBits::quiet_nan().get_val();
52 }
53 if (xbits.is_neg()) {
54 fputil::set_errno_if_required(EDOM);
55 fputil::raise_except_if_required(FE_INVALID);
56 return FPBits::quiet_nan().get_val();
57 }
58 return x;
59 }
60
61 // Domain error for inputs less than 1.0.
62 if (LIBC_UNLIKELY(x <= 1.0f)) {
63 if (x == 1.0f)
64 return FPBits::zero().get_val();
65 fputil::set_errno_if_required(EDOM);
66 fputil::raise_except_if_required(FE_INVALID);
67 return FPBits::quiet_nan().get_val();
68 }
69
70 if (auto r = ACOSHF16_EXCEPTS.lookup(x_bits: xbits.uintval());
71 LIBC_UNLIKELY(r.has_value()))
72 return r.value();
73
74 float xf = x;
75 // High-precision polynomial approximation for inputs close to 1.0
76 // ([1, 1.25)).
77 //
78 // Brief derivation:
79 // 1. Expand acosh(1 + delta) using Taylor series around delta=0:
80 // acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160
81 // - 5*delta^3/896 + 35*delta^4/18432 + ...]
82 // 2. Truncate the series to fit accurately for delta in [0, 0.25].
83 // 3. Polynomial coefficients (from sollya) used here are:
84 // P(delta) ≈ 1 - 0x1.555556p-4 * delta + 0x1.333334p-6 * delta^2
85 // - 0x1.6db6dcp-8 * delta^3 + 0x1.f1c71cp-10 * delta^4
86 // 4. The Sollya commands used to generate these coefficients were:
87 // > display = hexadecimal;
88 // > round(1/12, SG, RN);
89 // > round(3/160, SG, RN);
90 // > round(5/896, SG, RN);
91 // > round(35/18432, SG, RN);
92 // With hexadecimal display mode enabled, the outputs were:
93 // 0x1.555556p-4
94 // 0x1.333334p-6
95 // 0x1.6db6dcp-8
96 // 0x1.f1c71cp-10
97 // 5. The maximum absolute error, estimated using:
98 // dirtyinfnorm(acosh(1 + x) - sqrt(2*x) * P(x), [0, 0.25])
99 // is:
100 // 0x1.d84281p-22
101 if (LIBC_UNLIKELY(x_u < 0x3D00U)) {
102 float delta = xf - 1.0f;
103 float sqrt_2_delta = fputil::sqrt<float>(x: 2.0 * delta);
104 float pe = fputil::polyeval(x: delta, a0: 0x1p+0f, a: -0x1.555556p-4f, a: 0x1.333334p-6f,
105 a: -0x1.6db6dcp-8f, a: 0x1.f1c71cp-10f);
106 float approx = sqrt_2_delta * pe;
107 return fputil::cast<float16>(x: approx);
108 }
109
110 // acosh(x) = log(x + sqrt(x^2 - 1))
111 float sqrt_term = fputil::sqrt<float>(x: fputil::multiply_add(x: xf, y: xf, z: -1.0f));
112 float result = static_cast<float>(log_eval(x: xf + sqrt_term));
113
114 return fputil::cast<float16>(x: result);
115}
116
117} // namespace math
118
119} // namespace LIBC_NAMESPACE_DECL
120
121#endif // LIBC_TYPES_HAS_FLOAT16
122
123#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF16_H
124