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/// Implementation header for cos.
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_COS_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_COS_H
16
17#include "range_reduction_double_common.h"
18#include "sincos_eval.h"
19#include "src/__support/FPUtil/FEnvImpl.h"
20#include "src/__support/FPUtil/FPBits.h"
21#include "src/__support/FPUtil/double_double.h"
22#include "src/__support/FPUtil/dyadic_float.h"
23#include "src/__support/FPUtil/except_value_utils.h"
24#include "src/__support/macros/config.h"
25#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
26#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
27
28#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
29#include "range_reduction_double_fma.h"
30#else
31#include "range_reduction_double_nofma.h"
32#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
33
34namespace LIBC_NAMESPACE_DECL {
35
36namespace math {
37
38#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
39LIBC_INLINE double
40cos_accurate(double x, uint16_t x_e, unsigned k,
41 const range_reduction_double_internal::LargeRangeReduction
42 &range_reduction_large) {
43 using namespace math::range_reduction_double_internal;
44 using FPBits = typename fputil::FPBits<double>;
45
46 DFloat128 u_f128, sin_u, cos_u;
47 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT))
48 u_f128 = range_reduction_small_f128(x);
49 else
50 u_f128 = range_reduction_large.accurate();
51
52 math::sincos_eval_internal::sincos_eval(u: u_f128, sin_u, cos_u);
53
54 auto get_sin_k = [](unsigned kk) -> DFloat128 {
55 unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
56 DFloat128 ans = SIN_K_PI_OVER_128_F128[idx];
57 if (kk & 128)
58 ans.sign = Sign::NEG;
59 return ans;
60 };
61
62 // -sin(k * pi/128) = sin((k + 128) * pi/128)
63 // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
64 DFloat128 msin_k_f128 = get_sin_k(k + 128);
65 DFloat128 cos_k_f128 = get_sin_k(k + 64);
66
67 // cos(x) = cos((k * pi/128 + u)
68 // = cos(u) * cos(k*pi/128) - sin(u) * sin(k*pi/128)
69 DFloat128 r = fputil::quick_add(a: fputil::quick_mul(a: cos_k_f128, b: cos_u),
70 b: fputil::quick_mul(a: msin_k_f128, b: sin_u));
71
72 // TODO: Add assertion if Ziv's accuracy tests fail in debug mode.
73 // https://github.com/llvm/llvm-project/issues/96452.
74
75 return static_cast<double>(r);
76}
77#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
78
79LIBC_INLINE double cos(double x) {
80 using namespace range_reduction_double_internal;
81 using FPBits = typename fputil::FPBits<double>;
82 FPBits xbits(x);
83
84 uint16_t x_e = xbits.get_biased_exponent();
85
86 DoubleDouble y;
87 unsigned k = 0;
88 LargeRangeReduction range_reduction_large;
89
90 // |x| < 2^16.
91 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) {
92 // |x| < 2^-4
93 if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 4)) {
94 // |x| < 2^-27
95 if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) {
96 // Signed zeros.
97 if (LIBC_UNLIKELY(x == 0.0))
98 return 1.0;
99
100 // For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2.
101 return fputil::round_result_slightly_down(value_rn: 1.0);
102 }
103 // No range reduction needed.
104
105 // Use degree-8 polynomial approximation:
106 // cos(x) ~ 1 + a1 * x^2 + a2 * x^4 + a3 * x^6 + a4 * x^8
107 // ~ 1 + x^2 * Q(x^2).
108 // > P = fpminimax(cos(x), [|0, 2, 4, 6, 8|], [|1, D...|], [0, 2^-4]);
109 // > dirtyinfnorm(cos(x) - P, [-2^-4, 2^-4]);
110 // 0x1.3cfe...p-70
111 // > P;
112 constexpr double COEFFS[] = {-0x1p-1, 0x1.5555555555262p-5,
113 -0x1.6c16c1508bff1p-10,
114 0x1.a00ffd769159ap-16};
115 double x_sq = x * x;
116 double c0 = fputil::multiply_add(x: x_sq, y: COEFFS[1], z: COEFFS[0]);
117 double c1 = fputil::multiply_add(x: x_sq, y: COEFFS[3], z: COEFFS[2]);
118 double x4 = x_sq * x_sq;
119 double r_lo = fputil::multiply_add(x: x4, y: c1, z: c0) * x_sq;
120
121#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
122 return 1.0 + r_lo;
123#else
124 // Overall errors <= 1.0 * ulp(x^2) + 2^-69 for default rounding mode
125 // <= 2.0 * ulp(x^2) + 2^-69 for directed rounding mode.
126 double err = fputil::multiply_add(x: x_sq, y: 0x1.0p-51, z: 0x1.0p-69);
127 double r_lo_u = r_lo + err;
128 double r_lo_l = r_lo - err;
129 double r_upper = 1.0 + r_lo_u;
130 double r_lower = 1.0 + r_lo_l;
131
132 if (LIBC_LIKELY(r_upper == r_lower))
133 return r_upper;
134
135 k = range_reduction_small(x, u&: y);
136 return cos_accurate(x, x_e, k, range_reduction_large);
137#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
138 } else {
139 // Small range reduction.
140 k = range_reduction_small(x, u&: y);
141 }
142 } else {
143 // Inf or NaN
144 if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
145 if (xbits.is_signaling_nan()) {
146 fputil::raise_except_if_required(FE_INVALID);
147 return FPBits::quiet_nan().get_val();
148 }
149 // cos(+-Inf) = NaN
150 if (xbits.get_mantissa() == 0) {
151 fputil::set_errno_if_required(EDOM);
152 fputil::raise_except_if_required(FE_INVALID);
153 }
154 return x + FPBits::quiet_nan().get_val();
155 }
156
157 // Large range reduction.
158 k = range_reduction_large.fast(x, u&: y);
159 }
160
161 DoubleDouble sin_y, cos_y;
162
163 [[maybe_unused]] double err =
164 math::sincos_eval_internal::sincos_eval(u: y, sin_u&: sin_y, cos_u&: cos_y);
165
166 // Look up sin(k * pi/128) and cos(k * pi/128)
167#ifdef LIBC_MATH_HAS_SMALL_TABLES
168 // Memory saving versions. Use 65-entry table.
169 auto get_idx_dd = [](unsigned kk) -> DoubleDouble {
170 unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
171 DoubleDouble ans = SIN_K_PI_OVER_128[idx];
172 if (kk & 128) {
173 ans.hi = -ans.hi;
174 ans.lo = -ans.lo;
175 }
176 return ans;
177 };
178 DoubleDouble msin_k = get_idx_dd(k + 128);
179 DoubleDouble cos_k = get_idx_dd(k + 64);
180#else
181 // Fast look up version, but needs 256-entry table.
182 // -sin(k * pi/128) = sin((k + 128) * pi/128)
183 // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
184 DoubleDouble msin_k = SIN_K_PI_OVER_128[(k + 128) & 255];
185 DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255];
186#endif // LIBC_MATH_HAS_SMALL_TABLES
187
188 // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128).
189 // So k is an integer and -pi / 256 <= y <= pi / 256.
190 // Then cos(x) = cos((k * pi/128 + y)
191 // = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128)
192 DoubleDouble cos_k_cos_y = fputil::quick_mult(a: cos_y, b: cos_k);
193 DoubleDouble msin_k_sin_y = fputil::quick_mult(a: sin_y, b: msin_k);
194 // When k != 64 mod 128,
195 // |cos( k * pi/128 )| > pi/128 - epsilon > |y| >= |sin(y)|,
196 // and cos(y) > 1 - pi/128. So we can use Fast2Sum for the subtraction:
197 // cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128).
198 DoubleDouble rr = fputil::exact_add(a: cos_k_cos_y.hi, b: msin_k_sin_y.hi);
199 rr.lo += msin_k_sin_y.lo + cos_k_cos_y.lo;
200
201#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
202 return rr.hi + rr.lo;
203#else
204 double rlp = rr.lo + err;
205 double rlm = rr.lo - err;
206
207 double r_upper = rr.hi + rlp; // (rr.lo + ERR);
208 double r_lower = rr.hi + rlm; // (rr.lo - ERR);
209
210 // Ziv's rounding test.
211 if (LIBC_LIKELY(r_upper == r_lower))
212 return r_upper;
213
214 return cos_accurate(x, x_e, k, range_reduction_large);
215#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
216}
217
218} // namespace math
219
220} // namespace LIBC_NAMESPACE_DECL
221
222#endif // LLVM_LIBC_SRC___SUPPORT_MATH_COS_H
223