1//===-- Implementation header for atanf16 -----------------------*- 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_ATANF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ATANF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
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 constexpr float16 atanf16(float16 x) {
30 // Generated by Solly using the following command:
31 // > round(pi/2, SG, RN);
32 constexpr float PI_2 = 0x1.921fb6p0;
33
34#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
35 constexpr size_t N_EXCEPTS = 6;
36
37 constexpr fputil::ExceptValues<float16, N_EXCEPTS> ATANF16_EXCEPTS{.values: {
38 // (input, RZ output, RU offset, RD offset, RN offset)
39 {.input: 0x2745, .rnd_towardzero_result: 0x2744, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
40 {.input: 0x3099, .rnd_towardzero_result: 0x3090, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
41 {.input: 0x3c6c, .rnd_towardzero_result: 0x3aae, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
42 {.input: 0x466e, .rnd_towardzero_result: 0x3daa, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
43 {.input: 0x48ae, .rnd_towardzero_result: 0x3ddb, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
44 {.input: 0x5619, .rnd_towardzero_result: 0x3e3d, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
45 }};
46#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
47
48 using FPBits = fputil::FPBits<float16>;
49 FPBits xbits(x);
50
51 uint16_t x_u = xbits.uintval();
52 uint16_t x_abs = x_u & 0x7fff;
53 bool x_sign = x_u >> 15;
54 float sign = (x_sign ? -1.0 : 1.0);
55
56 // |x| >= +/-inf
57 if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
58 if (xbits.is_nan()) {
59 if (xbits.is_signaling_nan()) {
60 fputil::raise_except_if_required(FE_INVALID);
61 return FPBits::quiet_nan().get_val();
62 }
63 return x;
64 }
65
66 // atanf16(+/-inf) = +/-pi/2
67 return fputil::cast<float16>(x: sign * PI_2);
68 }
69
70 float xf = x;
71 float xsq = xf * xf;
72#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
73 // Handle exceptional values
74 if (auto r = ATANF16_EXCEPTS.lookup_odd(x_abs, sign: x_sign);
75 LIBC_UNLIKELY(r.has_value()))
76 return r.value();
77#endif
78
79 // |x| <= 0x1p0, |x| <= 1
80 if (x_abs <= 0x3c00) {
81 // atanf16(+/-0) = +/-0
82 if (LIBC_UNLIKELY(x_abs == 0))
83 return x;
84
85 // Degree-14 minimax odd polynomial of atan(x) generated by Sollya with:
86 // > P = fpminimax(atan(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14|], [|SG...|],
87 // [0, 1]);
88 float result = fputil::polyeval(
89 x: xsq, a0: 0x1.fffffcp-1f, a: -0x1.55519ep-2f, a: 0x1.98f6a8p-3f, a: -0x1.1f0a92p-3f,
90 a: 0x1.95b654p-4f, a: -0x1.e65492p-5f, a: 0x1.8c0c36p-6f, a: -0x1.32316ep-8f);
91 return fputil::cast<float16>(x: xf * result);
92 }
93
94 // If |x| > 1
95 // y = atan(x) = sign(x) * atan(|x|)
96 // atan(|x|) = pi/2 - atan(1/|x|)
97 // Recall, 1/|x| < 1
98 float x_inv_sq = 1.0f / xsq;
99 float x_inv = fputil::sqrt<float>(x: x_inv_sq);
100
101 // Degree-14 minimax odd polynomial of atan(x) generated by Sollya with:
102 // > P = fpminimax(atan(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14|], [|SG...|],
103 // [0, 1]);
104 float interm =
105 fputil::polyeval(x: x_inv_sq, a0: 0x1.fffffcp-1f, a: -0x1.55519ep-2f,
106 a: 0x1.98f6a8p-3f, a: -0x1.1f0a92p-3f, a: 0x1.95b654p-4f,
107 a: -0x1.e65492p-5f, a: 0x1.8c0c36p-6f, a: -0x1.32316ep-8f);
108
109 return fputil::cast<float16>(x: sign *
110 fputil::multiply_add(x: x_inv, y: -interm, z: PI_2));
111}
112
113} // namespace math
114
115} // namespace LIBC_NAMESPACE_DECL
116
117#endif // LIBC_TYPES_HAS_FLOAT16
118
119#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ATANF16_H
120