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