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_MATH_GENERIC_SINCOSF_UTILS_H 10 #define LLVM_LIBC_SRC_MATH_GENERIC_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) 18 #include "range_reduction_fma.h" 19 // using namespace LIBC_NAMESPACE::fma; 20 using LIBC_NAMESPACE::fma::FAST_PASS_BOUND; 21 using LIBC_NAMESPACE::fma::large_range_reduction; 22 using LIBC_NAMESPACE::fma::small_range_reduction; 23 24 #else 25 #include "range_reduction.h" 26 // using namespace LIBC_NAMESPACE::generic; 27 using LIBC_NAMESPACE::generic::FAST_PASS_BOUND; 28 using LIBC_NAMESPACE::generic::large_range_reduction; 29 using LIBC_NAMESPACE::generic::small_range_reduction; 30 #endif // LIBC_TARGET_CPU_HAS_FMA 31 32 namespace 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)); }; 38 const 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 63 static 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 Sollya 79 // 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 91 LIBC_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, where 107 // k = round(x * 32) and y = (x * 32) - k. 108 // => pi * x = (k + y) * pi / 32 109 static LIBC_INLINE int64_t range_reduction_sincospi(double x, double &y) { 110 double kd = fputil::nearest_integer(x * 32); 111 y = fputil::multiply_add<double>(x, 32.0, -kd); 112 113 return static_cast<int64_t>(kd); 114 } 115 116 LIBC_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_DECL 124 125 #endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF_UTILS_H 126