1 //===-- Collection of utils for sinf16/cosf16 -------------------*- 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_MATH_GENERIC_SINCOSF16_UTILS_H
10 #define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H
11
12 #include "src/__support/FPUtil/FPBits.h"
13 #include "src/__support/FPUtil/PolyEval.h"
14 #include "src/__support/FPUtil/nearest_integer.h"
15 #include "src/__support/common.h"
16 #include "src/__support/macros/config.h"
17
18 namespace LIBC_NAMESPACE_DECL {
19
20 // Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
21 // Table is generated with Sollya as follows:
22 // > display = hexadecimmal;
23 // > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); };
24 constexpr float SIN_K_PI_OVER_32[64] = {
25 0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2,
26 0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1,
27 0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1,
28 0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1,
29 0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1,
30 0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1,
31 0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2,
32 0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4,
33 0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2,
34 -0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1,
35 -0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1,
36 -0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1,
37 -0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1,
38 -0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1,
39 -0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2,
40 -0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4};
41
range_reduction_sincospif16(float x,float & y)42 LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
43 float kf = fputil::nearest_integer(x * 32);
44 y = fputil::multiply_add<float>(x, 32.0, -kf);
45
46 return static_cast<int32_t>(kf);
47 }
48
sincospif16_eval(float xf,float & sin_k,float & cos_k,float & sin_y,float & cosm1_y)49 LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k,
50 float &sin_y, float &cosm1_y) {
51 float y;
52 int32_t k = range_reduction_sincospif16(xf, y);
53
54 sin_k = SIN_K_PI_OVER_32[k & 63];
55 cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
56
57 // Recall, after range reduction, -0.5 <= y <= 0.5. For very small values of
58 // y, calculating sin(y * p/32) can be inaccurate. Generating a polynomial for
59 // sin(y * p/32)/y instead significantly reduces the relative errors.
60 float ysq = y * y;
61
62 // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
63 // with:
64 // > Q = fpminimax(sin(y * pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]);
65 sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f,
66 0x1.a03354p-21f, -0x1.ad02d2p-20f);
67
68 // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya
69 // with:
70 // > P = fpminimax(cos(y * pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]);
71 cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f,
72 0x1.a6f7a2p-29f);
73 }
74
75 } // namespace LIBC_NAMESPACE_DECL
76
77 #endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H
78