1//===-- Collection of utils for sinf/cosf/sincosf ---------------*- C++ -*-===//
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_SINCOSF_UTILS_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
11
12#include "src/__support/FPUtil/FPBits.h"
13#include "src/__support/FPUtil/PolyEval.h"
14#include "src/__support/macros/config.h"
15#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
16
17#if defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
18#include "range_reduction_fma.h"
19#else
20#include "range_reduction.h"
21#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
22
23namespace LIBC_NAMESPACE_DECL {
24
25namespace math {
26
27namespace sincosf_utils_internal {
28
29#if defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
30
31// using namespace LIBC_NAMESPACE::fma;
32using math::trigonometric_fma_utils_internal::FAST_PASS_BOUND;
33using math::trigonometric_fma_utils_internal::large_range_reduction;
34using math::trigonometric_fma_utils_internal::small_range_reduction;
35
36#else
37
38// using namespace LIBC_NAMESPACE::generic;
39using math::trigonometric_func_utils_internal::FAST_PASS_BOUND;
40using math::trigonometric_func_utils_internal::large_range_reduction;
41using math::trigonometric_func_utils_internal::small_range_reduction;
42
43#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
44
45// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
46// Table is generated with Sollya as follow:
47// > display = hexadecimal;
48// > for k from 0 to 63 do { D(sin(k * pi/32)); };
49LIBC_INLINE_VAR const double SIN_K_PI_OVER_32[64] = {
50 0x0.0000000000000p+0, 0x1.917a6bc29b42cp-4, 0x1.8f8b83c69a60bp-3,
51 0x1.294062ed59f06p-2, 0x1.87de2a6aea963p-2, 0x1.e2b5d3806f63bp-2,
52 0x1.1c73b39ae68c8p-1, 0x1.44cf325091dd6p-1, 0x1.6a09e667f3bcdp-1,
53 0x1.8bc806b151741p-1, 0x1.a9b66290ea1a3p-1, 0x1.c38b2f180bdb1p-1,
54 0x1.d906bcf328d46p-1, 0x1.e9f4156c62ddap-1, 0x1.f6297cff75cbp-1,
55 0x1.fd88da3d12526p-1, 0x1.0000000000000p+0, 0x1.fd88da3d12526p-1,
56 0x1.f6297cff75cbp-1, 0x1.e9f4156c62ddap-1, 0x1.d906bcf328d46p-1,
57 0x1.c38b2f180bdb1p-1, 0x1.a9b66290ea1a3p-1, 0x1.8bc806b151741p-1,
58 0x1.6a09e667f3bcdp-1, 0x1.44cf325091dd6p-1, 0x1.1c73b39ae68c8p-1,
59 0x1.e2b5d3806f63bp-2, 0x1.87de2a6aea963p-2, 0x1.294062ed59f06p-2,
60 0x1.8f8b83c69a60bp-3, 0x1.917a6bc29b42cp-4, 0x0.0000000000000p+0,
61 -0x1.917a6bc29b42cp-4, -0x1.8f8b83c69a60bp-3, -0x1.294062ed59f06p-2,
62 -0x1.87de2a6aea963p-2, -0x1.e2b5d3806f63bp-2, -0x1.1c73b39ae68c8p-1,
63 -0x1.44cf325091dd6p-1, -0x1.6a09e667f3bcdp-1, -0x1.8bc806b151741p-1,
64 -0x1.a9b66290ea1a3p-1, -0x1.c38b2f180bdb1p-1, -0x1.d906bcf328d46p-1,
65 -0x1.e9f4156c62ddap-1, -0x1.f6297cff75cbp-1, -0x1.fd88da3d12526p-1,
66 -0x1.0000000000000p+0, -0x1.fd88da3d12526p-1, -0x1.f6297cff75cbp-1,
67 -0x1.e9f4156c62ddap-1, -0x1.d906bcf328d46p-1, -0x1.c38b2f180bdb1p-1,
68 -0x1.a9b66290ea1a3p-1, -0x1.8bc806b151741p-1, -0x1.6a09e667f3bcdp-1,
69 -0x1.44cf325091dd6p-1, -0x1.1c73b39ae68c8p-1, -0x1.e2b5d3806f63bp-2,
70 -0x1.87de2a6aea963p-2, -0x1.294062ed59f06p-2, -0x1.8f8b83c69a60bp-3,
71 -0x1.917a6bc29b42cp-4,
72};
73
74LIBC_INLINE void sincosf_poly_eval(int64_t k, double y, double &sin_k,
75 double &cos_k, double &sin_y,
76 double &cosm1_y) {
77 // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
78 // So k is an integer and -0.5 <= y <= 0.5.
79 // Then sin(x) = sin((k + y)*pi/32)
80 // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
81
82 sin_k = SIN_K_PI_OVER_32[k & 63];
83 // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
84 // cos_k = cos(k * pi/32)
85 cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
86
87 double ysq = y * y;
88
89 // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
90 // with:
91 // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
92 sin_y =
93 y * fputil::polyeval(x: ysq, a0: 0x1.921fb54442d18p-4, a: -0x1.4abbce625abb1p-13,
94 a: 0x1.466bc624f2776p-24, a: -0x1.32c3a619d4a7ep-36);
95 // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
96 // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
97 // Note that cosm1_y = cos(y*pi/32) - 1.
98 cosm1_y = ysq * fputil::polyeval(x: ysq, a0: -0x1.3bd3cc9be430bp-8,
99 a: 0x1.03c1f070c2e27p-18, a: -0x1.55cc84bd942p-30);
100}
101
102LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
103 double &cos_k, double &sin_y, double &cosm1_y) {
104 int64_t k;
105 double y;
106
107 if (LIBC_LIKELY(x_abs < FAST_PASS_BOUND)) {
108 k = small_range_reduction(x: xd, y);
109 } else {
110 fputil::FPBits<float> x_bits(x_abs);
111 k = large_range_reduction(x: xd, x_exp: x_bits.get_exponent(), y);
112 }
113
114 sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
115}
116
117// Return k and y, where
118// k = round(x * 32) and y = (x * 32) - k.
119// => pi * x = (k + y) * pi / 32
120LIBC_INLINE int64_t range_reduction_sincospi(double x, double &y) {
121 double kd = fputil::nearest_integer(x: x * 32);
122 y = fputil::multiply_add(x, y: 32.0, z: -kd);
123
124 return static_cast<int64_t>(kd);
125}
126
127LIBC_INLINE void sincospif_eval(double xd, double &sin_k, double &cos_k,
128 double &sin_y, double &cosm1_y) {
129 double y;
130 int64_t k = range_reduction_sincospi(x: xd, y);
131 sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
132}
133
134} // namespace sincosf_utils_internal
135
136} // namespace math
137
138} // namespace LIBC_NAMESPACE_DECL
139
140#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
141