1//===----------------------------------------------------------------------===//
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/// \file
10/// Double-precision implementation of tanf.
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANF_DOUBLE_EVAL_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_TANF_DOUBLE_EVAL_H
16
17#include "src/__support/FPUtil/FEnvImpl.h"
18#include "src/__support/FPUtil/FPBits.h"
19#include "src/__support/FPUtil/PolyEval.h"
20#include "src/__support/FPUtil/except_value_utils.h"
21#include "src/__support/FPUtil/multiply_add.h"
22#include "src/__support/FPUtil/nearest_integer.h"
23#include "src/__support/macros/config.h"
24#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
25#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
26#include "src/__support/math/sincosf_utils.h"
27
28namespace LIBC_NAMESPACE_DECL {
29
30namespace math {
31namespace double_eval {
32
33namespace tanf_internal {
34
35#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
36// Exceptional cases for tanf.
37LIBC_INLINE_VAR constexpr size_t N_EXCEPTS = 6;
38
39LIBC_INLINE_VAR constexpr fputil::ExceptValues<float, N_EXCEPTS> TANF_EXCEPTS{.values: {
40 // (inputs, RZ output, RU offset, RD offset, RN offset)
41 // x = 0x1.ada6aap27, tan(x) = 0x1.e80304p-3 (RZ)
42 {.input: 0x4d56d355, .rnd_towardzero_result: 0x3e740182, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
43 // x = 0x1.862064p33, tan(x) = -0x1.8dee56p-3 (RZ)
44 {.input: 0x50431032, .rnd_towardzero_result: 0xbe46f72b, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 1},
45 // x = 0x1.af61dap48, tan(x) = 0x1.60d1c6p-2 (RZ)
46 {.input: 0x57d7b0ed, .rnd_towardzero_result: 0x3eb068e3, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
47 // x = 0x1.0088bcp52, tan(x) = 0x1.ca1edp0 (RZ)
48 {.input: 0x5980445e, .rnd_towardzero_result: 0x3fe50f68, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
49 // x = 0x1.f90dfcp72, tan(x) = 0x1.597f9cp-1 (RZ)
50 {.input: 0x63fc86fe, .rnd_towardzero_result: 0x3f2cbfce, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
51 // x = 0x1.a6ce12p86, tan(x) = -0x1.c5612ep-1 (RZ)
52 {.input: 0x6ad36709, .rnd_towardzero_result: 0xbf62b097, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
53}};
54#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
55
56} // namespace tanf_internal
57
58LIBC_INLINE float tanf(float x) {
59 using namespace sincosf_utils_internal;
60 using namespace tanf_internal;
61 using FPBits = typename fputil::FPBits<float>;
62 FPBits xbits(x);
63 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
64
65 // |x| < pi/32
66 if (LIBC_UNLIKELY(x_abs <= 0x3dc9'0fdbU)) {
67 double xd = static_cast<double>(x);
68
69 // |x| < 0x1.0p-12f
70 if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) {
71 if (LIBC_UNLIKELY(x_abs == 0U)) {
72 // For signed zeros.
73 return x;
74 }
75 // When |x| < 2^-12, the relative error of the approximation tan(x) ~ x
76 // is:
77 // |tan(x) - x| / |tan(x)| < |x^3| / (3|x|)
78 // = x^2 / 3
79 // < 2^-25
80 // < epsilon(1)/2.
81 // So the correctly rounded values of tan(x) are:
82 // = x + sign(x)*eps(x) if rounding mode = FE_UPWARD and x is positive,
83 // or (rounding mode = FE_DOWNWARD and x is
84 // negative),
85 // = x otherwise.
86 // To simplify the rounding decision and make it more efficient, we use
87 // fma(x, 2^-25, x) instead.
88 // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
89 // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
90 // |x| < 2^-125. For targets without FMA instructions, we simply use
91 // double for intermediate results as it is more efficient than using an
92 // emulated version of FMA.
93#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
94 return fputil::multiply_add(x, 0x1.0p-25f, x);
95#else
96 return static_cast<float>(fputil::multiply_add(x: xd, y: 0x1.0p-25, z: xd));
97#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
98 }
99
100 // |x| < pi/32
101 double xsq = xd * xd;
102
103 // Degree-9 minimax odd polynomial of tan(x) generated by Sollya with:
104 // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/32]);
105 double result =
106 fputil::polyeval(x: xsq, a0: 1.0, a: 0x1.555555553d022p-2, a: 0x1.111111ce442c1p-3,
107 a: 0x1.ba180a6bbdecdp-5, a: 0x1.69c0a88a0b71fp-6);
108 return static_cast<float>(xd * result);
109 }
110
111#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
112 bool x_sign = xbits.uintval() >> 31;
113 // Check for exceptional values
114 if (LIBC_UNLIKELY(x_abs == 0x3f8a1f62U)) {
115 // |x| = 0x1.143ec4p0
116 float sign = x_sign ? -1.0f : 1.0f;
117
118 // volatile is used to prevent compiler (gcc) from optimizing the
119 // computation, making the results incorrect in different rounding modes.
120 volatile float tmp = 0x1.ddf9f4p0f;
121 tmp = fputil::multiply_add(x: sign, y: tmp, z: sign * 0x1.1p-24f);
122
123 return tmp;
124 }
125#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
126
127 // |x| > 0x1.ada6a8p+27f
128 if (LIBC_UNLIKELY(x_abs > 0x4d56'd354U)) {
129 // Inf or NaN
130 if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
131 if (xbits.is_signaling_nan()) {
132 fputil::raise_except_if_required(FE_INVALID);
133 return FPBits::quiet_nan().get_val();
134 }
135
136 if (x_abs == 0x7f80'0000U) {
137 fputil::set_errno_if_required(EDOM);
138 fputil::raise_except_if_required(FE_INVALID);
139 }
140 return x + FPBits::quiet_nan().get_val();
141 }
142#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
143 // Other large exceptional values
144 if (auto r = TANF_EXCEPTS.lookup_odd(x_abs, sign: x_sign);
145 LIBC_UNLIKELY(r.has_value()))
146 return r.value();
147#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
148 }
149
150 // For |x| >= pi/32, we use the definition of tan(x) function:
151 // tan(x) = sin(x) / cos(x)
152 // The we follow the same computations of sin(x) and cos(x) as sinf, cosf,
153 // and sincosf.
154
155 double xd = static_cast<double>(x);
156 double sin_k, cos_k, sin_y, cosm1_y;
157
158 sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
159 // tan(x) = sin(x) / cos(x)
160 // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
161 using fputil::multiply_add;
162 return static_cast<float>(
163 multiply_add(x: sin_y, y: cos_k, z: multiply_add(x: cosm1_y, y: sin_k, z: sin_k)) /
164 multiply_add(x: sin_y, y: -sin_k, z: multiply_add(x: cosm1_y, y: cos_k, z: cos_k)));
165}
166
167} // namespace double_eval
168} // namespace math
169
170} // namespace LIBC_NAMESPACE_DECL
171
172#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANF_DOUBLE_EVAL_H
173