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