1 /* 2 * Single-precision scalar cospi 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 cospi(x) - cospif. 24 Maximum error: 2.64 ULP: 25 cospif(0x1.37e844p-4) got 0x1.f16b3p-1 26 want 0x1.f16b2ap-1. */ 27 float 28 cospif (float x) 29 { 30 if (isinf (x)) 31 return __math_invalidf (x); 32 33 float ax = asfloat (asuint (x) & ~0x80000000); 34 35 /* Edge cases for when cospif should be exactly +/- 1. (Integers) 36 0x1p23 is the limit for single precision to store any decimal places. */ 37 if (ax >= 0x1p24f) 38 return 1; 39 40 uint32_t m = roundf (ax); 41 if (m == ax) 42 return (m & 1) ? -1 : 1; 43 44 /* Any non-integer values >= 0x1p22f will be int +0.5. 45 These values should return exactly 0. */ 46 if (ax >= 0x1p22f) 47 return 0; 48 49 /* For very small inputs, squaring r causes underflow. 50 Values below this threshold can be approximated via cospi(x) ~= 1 - 51 (pi*x). */ 52 if (ax < 0x1p-31f) 53 return 1 - (C0 * x); 54 55 /* n = rint(|x|). */ 56 float n = ax + Shift; 57 uint32_t sign = asuint (n) << 31; 58 n = n - Shift; 59 60 /* We know that cospi(x) = sinpi(0.5 - x) 61 range reduction and offset into sinpi range -1/2 .. 1/2 62 r = 0.5 - |x - rint(x)|. */ 63 float r = 0.5f - fabs (ax - n); 64 65 /* y = sin(pi * r). */ 66 float r2 = r * r; 67 float y = fmaf (C5, r2, C4); 68 y = fmaf (y, r2, C3); 69 y = fmaf (y, r2, C2); 70 y = fmaf (y, r2, C1); 71 y = fmaf (y, r2, C0); 72 73 /* As all values are reduced to -1/2 .. 1/2, the result of cos(x) always be 74 positive, therefore, the sign must be introduced based upon if x rounds to 75 odd or even. */ 76 return asfloat (asuint (y * r) ^ sign); 77 } 78 79 PL_SIG (S, F, 1, cospi, -0.9, 0.9) 80 PL_TEST_ULP (cospif, 2.15) 81 PL_TEST_SYM_INTERVAL (cospif, 0, 0x1p-31, 5000) 82 PL_TEST_SYM_INTERVAL (cospif, 0x1p-31, 0.5, 10000) 83 PL_TEST_SYM_INTERVAL (cospif, 0.5, 0x1p22f, 10000) 84 PL_TEST_SYM_INTERVAL (cospif, 0x1p22f, inf, 10000) 85