| //===-- Collection of utils for sinf/cosf/sincosf ---------------*- C++ -*-===// |
| // |
| // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
| // See https://llvm.org/LICENSE.txt for license information. |
| // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
| // |
| //===----------------------------------------------------------------------===// |
| |
| #ifndef LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H |
| #define LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H |
| |
| #include "src/__support/FPUtil/FPBits.h" |
| #include "src/__support/FPUtil/PolyEval.h" |
| #include "src/__support/common.h" |
| |
| #if defined(LIBC_TARGET_HAS_FMA) |
| #include "range_reduction_fma.h" |
| // using namespace __llvm_libc::fma; |
| using __llvm_libc::fma::FAST_PASS_BOUND; |
| using __llvm_libc::fma::large_range_reduction; |
| using __llvm_libc::fma::small_range_reduction; |
| #else |
| #include "range_reduction.h" |
| // using namespace __llvm_libc::generic; |
| using __llvm_libc::generic::FAST_PASS_BOUND; |
| using __llvm_libc::generic::large_range_reduction; |
| using __llvm_libc::generic::small_range_reduction; |
| #endif // LIBC_TARGET_HAS_FMA |
| |
| namespace __llvm_libc { |
| |
| // Lookup table for sin(k * pi / 16) with k = 0, ..., 31. |
| // Table is generated with Sollya as follow: |
| // > display = hexadecimal; |
| // > for k from 0 to 31 do { D(sin(k * pi/16)); }; |
| const double SIN_K_PI_OVER_16[32] = { |
| 0x0.0000000000000p+0, 0x1.8f8b83c69a60bp-3, 0x1.87de2a6aea963p-2, |
| 0x1.1c73b39ae68c8p-1, 0x1.6a09e667f3bcdp-1, 0x1.a9b66290ea1a3p-1, |
| 0x1.d906bcf328d46p-1, 0x1.f6297cff75cb0p-1, 0x1.0000000000000p+0, |
| 0x1.f6297cff75cb0p-1, 0x1.d906bcf328d46p-1, 0x1.a9b66290ea1a3p-1, |
| 0x1.6a09e667f3bcdp-1, 0x1.1c73b39ae68c8p-1, 0x1.87de2a6aea963p-2, |
| 0x1.8f8b83c69a60bp-3, 0x0.0000000000000p+0, -0x1.8f8b83c69a60bp-3, |
| -0x1.87de2a6aea963p-2, -0x1.1c73b39ae68c8p-1, -0x1.6a09e667f3bcdp-1, |
| -0x1.a9b66290ea1a3p-1, -0x1.d906bcf328d46p-1, -0x1.f6297cff75cb0p-1, |
| -0x1.0000000000000p+0, -0x1.f6297cff75cb0p-1, -0x1.d906bcf328d46p-1, |
| -0x1.a9b66290ea1a3p-1, -0x1.6a09e667f3bcdp-1, -0x1.1c73b39ae68c8p-1, |
| -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60bp-3}; |
| |
| static inline void sincosf_eval(double xd, uint32_t x_abs, double &sin_k, |
| double &cos_k, double &sin_y, double &cosm1_y) { |
| int64_t k; |
| double y; |
| |
| if (likely(x_abs < FAST_PASS_BOUND)) { |
| k = small_range_reduction(xd, y); |
| } else { |
| fputil::FPBits<float> x_bits(x_abs); |
| k = large_range_reduction(xd, x_bits.get_exponent(), y); |
| } |
| |
| // After range reduction, k = round(x * 16 / pi) and y = (x * 16 / pi) - k. |
| // So k is an integer and -0.5 <= y <= 0.5. |
| // Then sin(x) = sin((k + y)*pi/16) |
| // = sin(y*pi/16) * cos(k*pi/16) + cos(y*pi/16) * sin(k*pi/16) |
| |
| sin_k = SIN_K_PI_OVER_16[k & 31]; |
| // cos(k * pi/16) = sin(k * pi/16 + pi/2) = sin((k + 8) * pi/16). |
| // cos_k = y * cos(k * pi/16) |
| cos_k = SIN_K_PI_OVER_16[(k + 8) & 31]; |
| |
| double ysq = y * y; |
| |
| // Degree-6 minimax even polynomial for sin(y*pi/16)/y generated by Sollya |
| // with: |
| // > Q = fpminimax(sin(y*pi/16)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]); |
| sin_y = y * fputil::polyeval(ysq, 0x1.921fb54442d17p-3, -0x1.4abbce6256adp-10, |
| 0x1.466bc5a5ac6b3p-19, -0x1.32bdcb4207562p-29); |
| // Degree-8 minimax even polynomial for cos(y*pi/16) generated by Sollya with: |
| // > P = fpminimax(cos(x*pi/16), [|0, 2, 4, 6, 8|], [|1, D...|], [0, 0.5]); |
| // Note that cosm1_y = cos(y*pi/16) - 1. |
| cosm1_y = |
| ysq * fputil::polyeval(ysq, -0x1.3bd3cc9be45dcp-6, 0x1.03c1f081b08ap-14, |
| -0x1.55d3c6fb0fb6ep-24, 0x1.e1d3d60f58873p-35); |
| } |
| |
| } // namespace __llvm_libc |
| |
| #endif // LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H |