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