1//===-- Single-precision atan function ------------------------------------===//
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#include "src/math/atanf.h"
10#include "inv_trigf_utils.h"
11#include "src/__support/FPUtil/FPBits.h"
12#include "src/__support/FPUtil/PolyEval.h"
13#include "src/__support/FPUtil/except_value_utils.h"
14#include "src/__support/FPUtil/multiply_add.h"
15#include "src/__support/FPUtil/nearest_integer.h"
16#include "src/__support/FPUtil/rounding_mode.h"
17#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
18
19namespace LIBC_NAMESPACE {
20
21LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
22 using FPBits = typename fputil::FPBits<float>;
23
24 constexpr double FINAL_SIGN[2] = {1.0, -1.0};
25 constexpr double SIGNED_PI_OVER_2[2] = {0x1.921fb54442d18p0,
26 -0x1.921fb54442d18p0};
27
28 FPBits x_bits(x);
29 Sign sign = x_bits.sign();
30 x_bits.set_sign(Sign::POS);
31 uint32_t x_abs = x_bits.uintval();
32
33 // x is inf or nan, |x| < 2^-4 or |x|= > 16.
34 if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U || x_abs >= 0x4180'0000U)) {
35 double x_d = static_cast<double>(x);
36 double const_term = 0.0;
37 if (LIBC_UNLIKELY(x_abs >= 0x4180'0000)) {
38 // atan(+-Inf) = +-pi/2.
39 if (x_bits.is_inf()) {
40 volatile double sign_pi_over_2 = SIGNED_PI_OVER_2[sign.is_neg()];
41 return static_cast<float>(sign_pi_over_2);
42 }
43 if (x_bits.is_nan())
44 return x;
45 // x >= 16
46 x_d = -1.0 / x_d;
47 const_term = SIGNED_PI_OVER_2[sign.is_neg()];
48 }
49 // 0 <= x < 1/16;
50 if (LIBC_UNLIKELY(x_bits.is_zero()))
51 return x;
52 // x <= 2^-12;
53 if (LIBC_UNLIKELY(x_abs < 0x3980'0000)) {
54#if defined(LIBC_TARGET_CPU_HAS_FMA)
55 return fputil::multiply_add(x, y: -0x1.0p-25f, z: x);
56#else
57 double x_d = static_cast<double>(x);
58 return static_cast<float>(fputil::multiply_add(x_d, -0x1.0p-25, x_d));
59#endif // LIBC_TARGET_CPU_HAS_FMA
60 }
61 // Use Taylor polynomial:
62 // atan(x) ~ x * (1 - x^2 / 3 + x^4 / 5 - x^6 / 7 + x^8 / 9 - x^10 / 11).
63 constexpr double ATAN_TAYLOR[6] = {
64 0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,
65 -0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4, -0x1.745d1745d1746p-4,
66 };
67 double x2 = x_d * x_d;
68 double x4 = x2 * x2;
69 double c0 = fputil::multiply_add(x: x2, y: ATAN_TAYLOR[1], z: ATAN_TAYLOR[0]);
70 double c1 = fputil::multiply_add(x: x2, y: ATAN_TAYLOR[3], z: ATAN_TAYLOR[2]);
71 double c2 = fputil::multiply_add(x: x2, y: ATAN_TAYLOR[5], z: ATAN_TAYLOR[4]);
72 double p = fputil::polyeval(x: x4, a0: c0, a: c1, a: c2);
73 double r = fputil::multiply_add(x: x_d, y: p, z: const_term);
74 return static_cast<float>(r);
75 }
76
77 // Range reduction steps:
78 // 1) atan(x) = sign(x) * atan(|x|)
79 // 2) If |x| > 1, atan(|x|) = pi/2 - atan(1/|x|)
80 // 3) For 1/16 < x <= 1, we find k such that: |x - k/16| <= 1/32.
81 // 4) Then we use polynomial approximation:
82 // atan(x) ~ atan((k/16) + (x - (k/16)) * Q(x - k/16)
83 // = P(x - k/16)
84 double x_d, const_term, final_sign;
85 int idx;
86
87 if (x_abs > 0x3f80'0000U) {
88 // |x| > 1, we need to invert x, so we will perform range reduction in
89 // double precision.
90 x_d = 1.0 / static_cast<double>(x_bits.get_val());
91 double k_d = fputil::nearest_integer(x: x_d * 0x1.0p4);
92 x_d = fputil::multiply_add(x: k_d, y: -0x1.0p-4, z: x_d);
93 idx = static_cast<int>(k_d);
94 final_sign = FINAL_SIGN[sign.is_pos()];
95 // Adjust constant term of the polynomial by +- pi/2.
96 const_term = fputil::multiply_add(x: final_sign, y: ATAN_COEFFS[idx][0],
97 z: SIGNED_PI_OVER_2[sign.is_neg()]);
98 } else {
99 // Exceptional value:
100 if (LIBC_UNLIKELY(x_abs == 0x3d8d'6b23U)) { // |x| = 0x1.1ad646p-4
101 return sign.is_pos() ? fputil::round_result_slightly_down(value_rn: 0x1.1a6386p-4f)
102 : fputil::round_result_slightly_up(value_rn: -0x1.1a6386p-4f);
103 }
104 // Perform range reduction in single precision.
105 float x_f = x_bits.get_val();
106 float k_f = fputil::nearest_integer(x: x_f * 0x1.0p4f);
107 x_f = fputil::multiply_add(x: k_f, y: -0x1.0p-4f, z: x_f);
108 x_d = static_cast<double>(x_f);
109 idx = static_cast<int>(k_f);
110 final_sign = FINAL_SIGN[sign.is_neg()];
111 const_term = final_sign * ATAN_COEFFS[idx][0];
112 }
113
114 double p = atan_eval(x: x_d, i: idx);
115 double r = fputil::multiply_add(x: final_sign * x_d, y: p, z: const_term);
116
117 return static_cast<float>(r);
118}
119
120} // namespace LIBC_NAMESPACE
121

source code of libc/src/math/generic/atanf.cpp