1//===-- Half-precision tan(x) 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/tanf16.h"
10#include "hdr/errno_macros.h"
11#include "hdr/fenv_macros.h"
12#include "sincosf16_utils.h"
13#include "src/__support/FPUtil/FEnvImpl.h"
14#include "src/__support/FPUtil/FPBits.h"
15#include "src/__support/FPUtil/cast.h"
16#include "src/__support/FPUtil/except_value_utils.h"
17#include "src/__support/FPUtil/multiply_add.h"
18#include "src/__support/macros/optimization.h"
19
20namespace LIBC_NAMESPACE_DECL {
21
22#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
23constexpr size_t N_EXCEPTS = 9;
24
25constexpr fputil::ExceptValues<float16, N_EXCEPTS> TANF16_EXCEPTS{{
26 // (input, RZ output, RU offset, RD offset, RN offset)
27 {0x2894, 0x2894, 1, 0, 1},
28 {0x3091, 0x3099, 1, 0, 0},
29 {0x3098, 0x30a0, 1, 0, 0},
30 {0x55ed, 0x3911, 1, 0, 0},
31 {0x607b, 0xc638, 0, 1, 1},
32 {0x674e, 0x3b7d, 1, 0, 0},
33 {0x6807, 0x4014, 1, 0, 1},
34 {0x6f4d, 0xbe19, 0, 1, 1},
35 {0x7330, 0xcb62, 0, 1, 0},
36}};
37#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
38
39LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) {
40 using FPBits = fputil::FPBits<float16>;
41 FPBits xbits(x);
42
43 uint16_t x_u = xbits.uintval();
44 uint16_t x_abs = x_u & 0x7fff;
45 float xf = x;
46
47#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
48 bool x_sign = x_u >> 15;
49 // Handle exceptional values
50 if (auto r = TANF16_EXCEPTS.lookup_odd(x_abs, x_sign);
51 LIBC_UNLIKELY(r.has_value()))
52 return r.value();
53#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
54
55 // |x| <= 0x1.d1p-5
56 if (LIBC_UNLIKELY(x_abs <= 0x2b44)) {
57 // |x| <= 0x1.398p-11
58 if (LIBC_UNLIKELY(x_abs <= 0x10e6)) {
59 // tan(+/-0) = +/-0
60 if (LIBC_UNLIKELY(x_abs == 0))
61 return x;
62
63 int rounding = fputil::quick_get_round();
64
65 // Exhaustive tests show that, when:
66 // x > 0, and rounding upward or
67 // x < 0, and rounding downward then,
68 // tan(x) = x * 2^-11 + x
69 if ((xbits.is_pos() && rounding == FE_UPWARD) ||
70 (xbits.is_neg() && rounding == FE_DOWNWARD))
71 return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
72 return x;
73 }
74
75 float xsq = xf * xf;
76
77 // Degree-6 minimax odd polynomial of tan(x) generated by Sollya with:
78 // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6|], [|1, SG...|], [0, pi/32]);
79 float result = fputil::polyeval(xsq, 0x1p0f, 0x1.555556p-2f, 0x1.110ee4p-3f,
80 0x1.be80f6p-5f);
81
82 return fputil::cast<float16>(xf * result);
83 }
84
85 // tan(+/-inf) = NaN, and tan(NaN) = NaN
86 if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
87 if (xbits.is_signaling_nan()) {
88 fputil::raise_except_if_required(FE_INVALID);
89 return FPBits::quiet_nan().get_val();
90 }
91 // x = +/-inf
92 if (x_abs == 0x7c00) {
93 fputil::set_errno_if_required(EDOM);
94 fputil::raise_except_if_required(FE_INVALID);
95 }
96
97 return x + FPBits::quiet_nan().get_val();
98 }
99
100 // Range reduction:
101 // For |x| > pi/32, we perform range reduction as follows:
102 // Find k and y such that:
103 // x = (k + y) * pi/32;
104 // k is an integer, |y| < 0.5
105 //
106 // This is done by performing:
107 // k = round(x * 32/pi)
108 // y = x * 32/pi - k
109 //
110 // Once k and y are computed, we then deduce the answer by the formula:
111 // tan(x) = sin(x) / cos(x)
112 // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
113 float sin_k, cos_k, sin_y, cosm1_y;
114 sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
115
116 // Note that, cosm1_y = cos_y - 1:
117 using fputil::multiply_add;
118 return fputil::cast<float16>(
119 multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) /
120 multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k)));
121}
122
123} // namespace LIBC_NAMESPACE_DECL
124

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