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___SUPPORT_MATH_SINCOSF16_UTILS_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF16_UTILS_H
11
12#include "src/__support/FPUtil/PolyEval.h"
13#include "src/__support/FPUtil/nearest_integer.h"
14#include "src/__support/common.h"
15#include "src/__support/macros/config.h"
16
17namespace LIBC_NAMESPACE_DECL {
18
19namespace math {
20
21namespace sincosf16_internal {
22
23// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
24// Table is generated with Sollya as follows:
25// > display = hexadecimmal;
26// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); };
27LIBC_INLINE_VAR constexpr float SIN_K_PI_OVER_32[64] = {
28 0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2,
29 0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1,
30 0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1,
31 0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1,
32 0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1,
33 0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1,
34 0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2,
35 0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4,
36 0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2,
37 -0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1,
38 -0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1,
39 -0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1,
40 -0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1,
41 -0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1,
42 -0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2,
43 -0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4};
44
45LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
46 float kf = fputil::nearest_integer(x: x * 32);
47 y = fputil::multiply_add(x, y: 32.0f, z: -kf);
48
49 return static_cast<int32_t>(kf);
50}
51
52// Recall, range reduction:
53// k = round(x * 32/pi)
54//
55// The precision choice of 'double' in the following function is to minimize
56// rounding errors in this initial scaling step,
57// preserving enough bits so errors accumulated while computing the subtraction:
58// y = x * 32/pi - round(x * 32/pi)
59// are beyond the least-significant bit of single-precision used during
60// further intermediate computation.
61LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) {
62 // Generated by Sollya with:
63 // > D(32/pi);
64 constexpr double THIRTYTWO_OVER_PI = 0x1.45f306dc9c883p3;
65
66 double prod = x * THIRTYTWO_OVER_PI;
67 double kd = fputil::nearest_integer(x: prod);
68 y = static_cast<float>(prod - kd);
69
70 return static_cast<int32_t>(kd);
71}
72
73LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k,
74 float &cos_k, float &sin_y,
75 float &cosm1_y) {
76
77 sin_k = SIN_K_PI_OVER_32[k & 63];
78 cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
79
80 // Recall, after range reduction, -0.5 <= y <= 0.5. For very small values of
81 // y, calculating sin(y * p/32) can be inaccurate. Generating a polynomial for
82 // sin(y * p/32)/y instead significantly reduces the relative errors.
83 float ysq = y * y;
84
85 // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
86 // with:
87 // > Q = fpminimax(sin(y * pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]);
88 sin_y = y * fputil::polyeval(x: ysq, a0: 0x1.921fb6p-4f, a: -0x1.4aeabcp-13f,
89 a: 0x1.a03354p-21f, a: -0x1.ad02d2p-20f);
90
91 // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya
92 // with:
93 // > P = fpminimax(cos(y * pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]);
94 cosm1_y = ysq * fputil::polyeval(x: ysq, a0: -0x1.3bd3ccp-8f, a: 0x1.03a61ap-18f,
95 a: 0x1.a6f7a2p-29f);
96}
97
98LIBC_INLINE void sincosf16_eval(float xf, float &sin_k, float &cos_k,
99 float &sin_y, float &cosm1_y) {
100 float y;
101 int32_t k = range_reduction_sincosf16(x: xf, y);
102
103 sincosf16_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
104}
105
106LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k,
107 float &sin_y, float &cosm1_y) {
108 float y;
109 int32_t k = range_reduction_sincospif16(x: xf, y);
110
111 sincosf16_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
112}
113
114} // namespace sincosf16_internal
115
116} // namespace math
117
118} // namespace LIBC_NAMESPACE_DECL
119
120#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF16_UTILS_H
121