1 //===-- Half-precision sinpif function ------------------------------------===// 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 #include "src/math/sinpif16.h" 10 #include "hdr/errno_macros.h" 11 #include "hdr/fenv_macros.h" 12 #include "sincosf16_utils.h" 13 #include "src/__support/FPUtil/FEnvImpl.h" 14 #include "src/__support/FPUtil/FPBits.h" 15 #include "src/__support/FPUtil/cast.h" 16 #include "src/__support/FPUtil/multiply_add.h" 17 18 namespace LIBC_NAMESPACE_DECL { 19 20 LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { 21 using FPBits = typename fputil::FPBits<float16>; 22 FPBits xbits(x); 23 24 uint16_t x_u = xbits.uintval(); 25 uint16_t x_abs = x_u & 0x7fff; 26 float xf = x; 27 28 // Range reduction: 29 // For |x| > 1/32, we perform range reduction as follows: 30 // Find k and y such that: 31 // x = (k + y) * 1/32 32 // k is an integer 33 // |y| < 0.5 34 // 35 // This is done by performing: 36 // k = round(x * 32) 37 // y = x * 32 - k 38 // 39 // Once k and y are computed, we then deduce the answer by the sine of sum 40 // formula: 41 // sin(x * pi) = sin((k + y) * pi/32) 42 // = sin(k * pi/32) * cos(y * pi/32) + 43 // sin(y * pi/32) * cos(k * pi/32) 44 45 // For signed zeros 46 if (LIBC_UNLIKELY(x_abs == 0U)) 47 return x; 48 49 // Numbers greater or equal to 2^10 are integers, or infinity, or NaN 50 if (LIBC_UNLIKELY(x_abs >= 0x6400)) { 51 // Check for NaN or infinity values 52 if (LIBC_UNLIKELY(x_abs >= 0x7c00)) { 53 // If value is equal to infinity 54 if (x_abs == 0x7c00) { 55 fputil::set_errno_if_required(EDOM); 56 fputil::raise_except_if_required(FE_INVALID); 57 } 58 59 return x + FPBits::quiet_nan().get_val(); 60 } 61 return FPBits::zero(xbits.sign()).get_val(); 62 } 63 64 float sin_k, cos_k, sin_y, cosm1_y; 65 sincospif16_eval(xf, sin_k, cos_k, sin_y, cosm1_y); 66 67 if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0)) 68 return FPBits::zero(xbits.sign()).get_val(); 69 70 // Since, cosm1_y = cos_y - 1, therefore: 71 // sin(x * pi) = cos_k * sin_y + sin_k + (cosm1_y * sin_k) 72 return fputil::cast<float16>(fputil::multiply_add( 73 sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k))); 74 } 75 76 } // namespace LIBC_NAMESPACE_DECL 77