1*f3087befSAndrew Turner /* 2*f3087befSAndrew Turner * Single-precision scalar sinpi function. 3*f3087befSAndrew Turner * 4*f3087befSAndrew Turner * Copyright (c) 2023-2024, Arm Limited. 5*f3087befSAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6*f3087befSAndrew Turner */ 7*f3087befSAndrew Turner 8*f3087befSAndrew Turner #include "mathlib.h" 9*f3087befSAndrew Turner #include "math_config.h" 10*f3087befSAndrew Turner #include "test_sig.h" 11*f3087befSAndrew Turner #include "test_defs.h" 12*f3087befSAndrew Turner 13*f3087befSAndrew Turner /* Taylor series coefficents for sin(pi * x). */ 14*f3087befSAndrew Turner #define C0 0x1.921fb6p1f 15*f3087befSAndrew Turner #define C1 -0x1.4abbcep2f 16*f3087befSAndrew Turner #define C2 0x1.466bc6p1f 17*f3087befSAndrew Turner #define C3 -0x1.32d2ccp-1f 18*f3087befSAndrew Turner #define C4 0x1.50783p-4f 19*f3087befSAndrew Turner #define C5 -0x1.e30750p-8f 20*f3087befSAndrew Turner 21*f3087befSAndrew Turner #define Shift 0x1.0p+23f 22*f3087befSAndrew Turner 23*f3087befSAndrew Turner /* Approximation for scalar single-precision sinpi(x) - sinpif. 24*f3087befSAndrew Turner Maximum error: 2.48 ULP: 25*f3087befSAndrew Turner sinpif(0x1.d062b6p-2) got 0x1.fa8c06p-1 26*f3087befSAndrew Turner want 0x1.fa8c02p-1. */ 27*f3087befSAndrew Turner float 28*f3087befSAndrew Turner arm_math_sinpif (float x) 29*f3087befSAndrew Turner { 30*f3087befSAndrew Turner if (isinf (x) || isnan (x)) 31*f3087befSAndrew Turner return __math_invalidf (x); 32*f3087befSAndrew Turner 33*f3087befSAndrew Turner float r = asfloat (asuint (x) & ~0x80000000); 34*f3087befSAndrew Turner uint32_t sign = asuint (x) & 0x80000000; 35*f3087befSAndrew Turner 36*f3087befSAndrew Turner /* Edge cases for when sinpif should be exactly 0. (Integers) 37*f3087befSAndrew Turner 0x1p23 is the limit for single precision to store any decimal places. */ 38*f3087befSAndrew Turner if (r >= 0x1p23f) 39*f3087befSAndrew Turner return asfloat (sign); 40*f3087befSAndrew Turner 41*f3087befSAndrew Turner int32_t m = roundf (r); 42*f3087befSAndrew Turner if (m == r) 43*f3087befSAndrew Turner return asfloat (sign); 44*f3087befSAndrew Turner 45*f3087befSAndrew Turner /* For very small inputs, squaring r causes underflow. 46*f3087befSAndrew Turner Values below this threshold can be approximated via sinpi(x) ~= pi*x. */ 47*f3087befSAndrew Turner if (r < 0x1p-31f) 48*f3087befSAndrew Turner return C0 * x; 49*f3087befSAndrew Turner 50*f3087befSAndrew Turner /* Any non-integer values >= 0x1p22f will be int + 0.5. 51*f3087befSAndrew Turner These values should return exactly 1 or -1. */ 52*f3087befSAndrew Turner if (r >= 0x1p22f) 53*f3087befSAndrew Turner { 54*f3087befSAndrew Turner uint32_t iy = ((m & 1) << 31) ^ asuint (-1.0f); 55*f3087befSAndrew Turner return asfloat (sign ^ iy); 56*f3087befSAndrew Turner } 57*f3087befSAndrew Turner 58*f3087befSAndrew Turner /* n = rint(|x|). */ 59*f3087befSAndrew Turner float n = r + Shift; 60*f3087befSAndrew Turner sign ^= (asuint (n) << 31); 61*f3087befSAndrew Turner n = n - Shift; 62*f3087befSAndrew Turner 63*f3087befSAndrew Turner /* r = |x| - n (range reduction into -1/2 .. 1/2). */ 64*f3087befSAndrew Turner r = r - n; 65*f3087befSAndrew Turner 66*f3087befSAndrew Turner /* y = sin(pi * r). */ 67*f3087befSAndrew Turner float r2 = r * r; 68*f3087befSAndrew Turner float y = fmaf (C5, r2, C4); 69*f3087befSAndrew Turner y = fmaf (y, r2, C3); 70*f3087befSAndrew Turner y = fmaf (y, r2, C2); 71*f3087befSAndrew Turner y = fmaf (y, r2, C1); 72*f3087befSAndrew Turner y = fmaf (y, r2, C0); 73*f3087befSAndrew Turner 74*f3087befSAndrew Turner /* Copy sign of x to sin(|x|). */ 75*f3087befSAndrew Turner return asfloat (asuint (y * r) ^ sign); 76*f3087befSAndrew Turner } 77*f3087befSAndrew Turner 78*f3087befSAndrew Turner #if WANT_EXPERIMENTAL_MATH 79*f3087befSAndrew Turner float 80*f3087befSAndrew Turner sinpif (float x) 81*f3087befSAndrew Turner { 82*f3087befSAndrew Turner return arm_math_sinpif (x); 83*f3087befSAndrew Turner } 84*f3087befSAndrew Turner #endif 85*f3087befSAndrew Turner 86*f3087befSAndrew Turner #if WANT_TRIGPI_TESTS 87*f3087befSAndrew Turner TEST_ULP (arm_math_sinpif, 1.99) 88*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_sinpif, 0, 0x1p-31, 5000) 89*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_sinpif, 0x1p-31, 0.5, 10000) 90*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_sinpif, 0.5, 0x1p22f, 10000) 91*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_sinpif, 0x1p22f, inf, 10000) 92*f3087befSAndrew Turner #endif 93