126 lines · c
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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H10#define LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H11 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_FMA16 17#if defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)18#include "range_reduction_fma.h"19// using namespace LIBC_NAMESPACE::fma;20using LIBC_NAMESPACE::fma::FAST_PASS_BOUND;21using LIBC_NAMESPACE::fma::large_range_reduction;22using LIBC_NAMESPACE::fma::small_range_reduction;23 24#else25#include "range_reduction.h"26// using namespace LIBC_NAMESPACE::generic;27using LIBC_NAMESPACE::generic::FAST_PASS_BOUND;28using LIBC_NAMESPACE::generic::large_range_reduction;29using LIBC_NAMESPACE::generic::small_range_reduction;30#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE31 32namespace LIBC_NAMESPACE_DECL {33 34// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.35// Table is generated with Sollya as follow:36// > display = hexadecimal;37// > for k from 0 to 63 do { D(sin(k * pi/32)); };38const double SIN_K_PI_OVER_32[64] = {39 0x0.0000000000000p+0, 0x1.917a6bc29b42cp-4, 0x1.8f8b83c69a60bp-3,40 0x1.294062ed59f06p-2, 0x1.87de2a6aea963p-2, 0x1.e2b5d3806f63bp-2,41 0x1.1c73b39ae68c8p-1, 0x1.44cf325091dd6p-1, 0x1.6a09e667f3bcdp-1,42 0x1.8bc806b151741p-1, 0x1.a9b66290ea1a3p-1, 0x1.c38b2f180bdb1p-1,43 0x1.d906bcf328d46p-1, 0x1.e9f4156c62ddap-1, 0x1.f6297cff75cbp-1,44 0x1.fd88da3d12526p-1, 0x1.0000000000000p+0, 0x1.fd88da3d12526p-1,45 0x1.f6297cff75cbp-1, 0x1.e9f4156c62ddap-1, 0x1.d906bcf328d46p-1,46 0x1.c38b2f180bdb1p-1, 0x1.a9b66290ea1a3p-1, 0x1.8bc806b151741p-1,47 0x1.6a09e667f3bcdp-1, 0x1.44cf325091dd6p-1, 0x1.1c73b39ae68c8p-1,48 0x1.e2b5d3806f63bp-2, 0x1.87de2a6aea963p-2, 0x1.294062ed59f06p-2,49 0x1.8f8b83c69a60bp-3, 0x1.917a6bc29b42cp-4, 0x0.0000000000000p+0,50 -0x1.917a6bc29b42cp-4, -0x1.8f8b83c69a60bp-3, -0x1.294062ed59f06p-2,51 -0x1.87de2a6aea963p-2, -0x1.e2b5d3806f63bp-2, -0x1.1c73b39ae68c8p-1,52 -0x1.44cf325091dd6p-1, -0x1.6a09e667f3bcdp-1, -0x1.8bc806b151741p-1,53 -0x1.a9b66290ea1a3p-1, -0x1.c38b2f180bdb1p-1, -0x1.d906bcf328d46p-1,54 -0x1.e9f4156c62ddap-1, -0x1.f6297cff75cbp-1, -0x1.fd88da3d12526p-1,55 -0x1.0000000000000p+0, -0x1.fd88da3d12526p-1, -0x1.f6297cff75cbp-1,56 -0x1.e9f4156c62ddap-1, -0x1.d906bcf328d46p-1, -0x1.c38b2f180bdb1p-1,57 -0x1.a9b66290ea1a3p-1, -0x1.8bc806b151741p-1, -0x1.6a09e667f3bcdp-1,58 -0x1.44cf325091dd6p-1, -0x1.1c73b39ae68c8p-1, -0x1.e2b5d3806f63bp-2,59 -0x1.87de2a6aea963p-2, -0x1.294062ed59f06p-2, -0x1.8f8b83c69a60bp-3,60 -0x1.917a6bc29b42cp-4,61};62 63static LIBC_INLINE void sincosf_poly_eval(int64_t k, double y, double &sin_k,64 double &cos_k, double &sin_y,65 double &cosm1_y) {66 // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.67 // So k is an integer and -0.5 <= y <= 0.5.68 // Then sin(x) = sin((k + y)*pi/32)69 // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)70 71 sin_k = SIN_K_PI_OVER_32[k & 63];72 // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).73 // cos_k = cos(k * pi/32)74 cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];75 76 double ysq = y * y;77 78 // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya79 // with:80 // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);81 sin_y =82 y * fputil::polyeval(ysq, 0x1.921fb54442d18p-4, -0x1.4abbce625abb1p-13,83 0x1.466bc624f2776p-24, -0x1.32c3a619d4a7ep-36);84 // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:85 // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);86 // Note that cosm1_y = cos(y*pi/32) - 1.87 cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3cc9be430bp-8,88 0x1.03c1f070c2e27p-18, -0x1.55cc84bd942p-30);89}90 91LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,92 double &cos_k, double &sin_y, double &cosm1_y) {93 int64_t k;94 double y;95 96 if (LIBC_LIKELY(x_abs < FAST_PASS_BOUND)) {97 k = small_range_reduction(xd, y);98 } else {99 fputil::FPBits<float> x_bits(x_abs);100 k = large_range_reduction(xd, x_bits.get_exponent(), y);101 }102 103 sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);104}105 106// Return k and y, where107// k = round(x * 32) and y = (x * 32) - k.108// => pi * x = (k + y) * pi / 32109static LIBC_INLINE int64_t range_reduction_sincospi(double x, double &y) {110 double kd = fputil::nearest_integer(x * 32);111 y = fputil::multiply_add(x, 32.0, -kd);112 113 return static_cast<int64_t>(kd);114}115 116LIBC_INLINE void sincospif_eval(double xd, double &sin_k, double &cos_k,117 double &sin_y, double &cosm1_y) {118 double y;119 int64_t k = range_reduction_sincospi(xd, y);120 sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);121}122 123} // namespace LIBC_NAMESPACE_DECL124 125#endif // LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H126