1//===-- Implementation header for atan2f16 ----------------------*- 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_ATAN2F16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ATAN2F16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "inv_trigf_utils.h"
17#include "src/__support/FPUtil/FEnvImpl.h"
18#include "src/__support/FPUtil/FPBits.h"
19#include "src/__support/FPUtil/cast.h"
20#include "src/__support/FPUtil/multiply_add.h"
21#include "src/__support/FPUtil/nearest_integer.h"
22#include "src/__support/macros/optimization.h"
23
24namespace LIBC_NAMESPACE_DECL {
25
26namespace math {
27
28LIBC_INLINE float16 atan2f16(float16 y, float16 x) {
29 using namespace inv_trigf_utils_internal;
30 using FPBits = fputil::FPBits<float16>;
31
32 constexpr double IS_NEG[2] = {1.0, -1.0};
33 constexpr double PI = 0x1.921fb54442d18p+1;
34 constexpr double PI_OVER_2 = 0x1.921fb54442d18p+0;
35 constexpr double PI_OVER_4 = 0x1.921fb54442d18p-1;
36 constexpr double THREE_PI_OVER_4 = 0x1.2d97c7f3321d2p+1;
37
38 // const_term[x_sign][recip]; recip = (|x| < |y|)
39 constexpr double CONST_TERM[2][2] = {
40 {0.0, -PI_OVER_2},
41 {-PI, PI_OVER_2},
42 };
43
44 FPBits x_bits(x), y_bits(y);
45 bool x_sign = x_bits.sign().is_neg();
46 bool y_sign = y_bits.sign().is_neg();
47 x_bits.set_sign(Sign::POS);
48 y_bits.set_sign(Sign::POS);
49 uint16_t x_abs = x_bits.uintval();
50 uint16_t y_abs = y_bits.uintval();
51 uint16_t max_abs = x_abs > y_abs ? x_abs : y_abs;
52 uint16_t min_abs = x_abs <= y_abs ? x_abs : y_abs;
53
54 if (LIBC_UNLIKELY(max_abs >= 0x7c00U || min_abs == 0)) {
55 if (x_bits.is_nan() || y_bits.is_nan()) {
56 if (x_bits.is_signaling_nan() || y_bits.is_signaling_nan())
57 fputil::raise_except_if_required(FE_INVALID);
58 return FPBits::quiet_nan().get_val();
59 }
60 size_t x_except = (x_abs == 0) ? 0 : (x_abs == 0x7c00U ? 2 : 1);
61 size_t y_except = (y_abs == 0) ? 0 : (y_abs == 0x7c00U ? 2 : 1);
62 constexpr double EXCEPTS[3][3][2] = {
63 {{0.0, PI}, {0.0, PI}, {0.0, PI}},
64 {{PI_OVER_2, PI_OVER_2}, {0.0, 0.0}, {0.0, PI}},
65 {{PI_OVER_2, PI_OVER_2},
66 {PI_OVER_2, PI_OVER_2},
67 {PI_OVER_4, THREE_PI_OVER_4}},
68 };
69 double r = IS_NEG[y_sign] * EXCEPTS[y_except][x_except][x_sign];
70 return fputil::cast<float16>(x: r);
71 }
72
73 bool recip = x_abs < y_abs;
74 double final_sign = IS_NEG[x_sign ^ y_sign ^ recip];
75 double const_term = CONST_TERM[x_sign][recip];
76
77 // atan2(y,x) = final_sign * (const_term + atan(n/d)),
78 // where n = min(|x|,|y|), d = max(|x|,|y|), so 0 <= n/d <= 1.
79 //
80 // To compute atan(n/d), we use a lookup table with 16 equally-spaced knots:
81 // idx = round(16 * n/d), so |n/d - idx/16| <= 1/32.
82 // q_d = n/d - idx/16.
83 // Then by the atan addition formula:
84 // atan(n/d) = atan(idx/16) + q_d * Q(q_d)
85 // where Q(q_d) approximates
86 // [atan(idx/16 + q_d) - atan(idx/16)] / q_d
87 // via atan_eval(q_d, idx) from inv_trigf_utils.
88 double n = static_cast<double>(FPBits(min_abs).get_val());
89 double d = static_cast<double>(FPBits(max_abs).get_val());
90 double q_d = n / d;
91
92 double k_d = fputil::nearest_integer(x: q_d * 0x1.0p4);
93 int idx = static_cast<int>(k_d);
94 q_d = fputil::multiply_add(x: k_d, y: -0x1.0p-4, z: q_d);
95
96 double p = atan_eval(x: q_d, i: static_cast<unsigned>(idx));
97 double r = final_sign *
98 fputil::multiply_add(x: q_d, y: p, z: const_term + ATAN_K_OVER_16[idx]);
99 return fputil::cast<float16>(x: r);
100}
101
102} // namespace math
103
104} // namespace LIBC_NAMESPACE_DECL
105
106#endif // LIBC_TYPES_HAS_FLOAT16
107
108#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ATAN2F16_H
109