1 /* 2 * Double-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 #include "poly_scalar_f64.h" 13 14 /* Taylor series coefficents for sin(pi * x). 15 C2 coefficient (orginally ~=5.16771278) has been split into two parts: 16 C2_hi = 4, C2_lo = C2 - C2_hi (~=1.16771278) 17 This change in magnitude reduces floating point rounding errors. 18 C2_hi is then reintroduced after the polynomial approxmation. */ 19 static const double poly[] 20 = { 0x1.921fb54442d184p1, -0x1.2aef39896f94bp0, 0x1.466bc6775ab16p1, 21 -0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4, -0x1.e30750a28c88ep-8, 22 0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16, 0x1.af86ae521260bp-21, 23 -0x1.012a9870eeb7dp-25 }; 24 25 #define Shift 0x1.8p+52 26 27 /* Approximation for scalar double-precision cospi(x). 28 Maximum error: 3.13 ULP: 29 cospi(0x1.160b129300112p-21) got 0x1.fffffffffd16bp-1 30 want 0x1.fffffffffd16ep-1. */ 31 double 32 cospi (double x) 33 { 34 if (isinf (x)) 35 return __math_invalid (x); 36 37 double ax = asdouble (asuint64 (x) & ~0x8000000000000000); 38 39 /* Edge cases for when cospif should be exactly 1. (Integers) 40 0x1p53 is the limit for single precision to store any decimal places. */ 41 if (ax >= 0x1p53) 42 return 1; 43 44 /* If x is an integer, return +- 1, based upon if x is odd. */ 45 uint64_t m = (uint64_t) ax; 46 if (m == ax) 47 return (m & 1) ? -1 : 1; 48 49 /* For very small inputs, squaring r causes underflow. 50 Values below this threshold can be approximated via 51 cospi(x) ~= 1. */ 52 if (ax < 0x1p-63) 53 return 1; 54 55 /* Any non-integer values >= 0x1x51 will be int +0.5. 56 These values should return exactly 0. */ 57 if (ax >= 0x1p51) 58 return 0; 59 60 /* n = rint(|x|). */ 61 double n = ax + Shift; 62 uint64_t sign = asuint64 (n) << 63; 63 n = n - Shift; 64 65 /* We know that cospi(x) = sinpi(0.5 - x) 66 range reduction and offset into sinpi range -1/2 .. 1/2 67 r = 0.5 - |x - rint(x)|. */ 68 double r = 0.5 - fabs (ax - n); 69 70 /* y = sin(r). */ 71 double r2 = r * r; 72 double y = horner_9_f64 (r2, poly); 73 y = y * r; 74 75 /* Reintroduce C2_hi. */ 76 y = fma (-4 * r2, r, y); 77 78 /* As all values are reduced to -1/2 .. 1/2, the result of cos(x) always be 79 positive, therefore, the sign must be introduced based upon if x rounds to 80 odd or even. */ 81 return asdouble (asuint64 (y) ^ sign); 82 } 83 84 PL_SIG (S, D, 1, cospi, -0.9, 0.9) 85 PL_TEST_ULP (cospi, 2.63) 86 PL_TEST_SYM_INTERVAL (cospi, 0, 0x1p-63, 5000) 87 PL_TEST_SYM_INTERVAL (cospi, 0x1p-63, 0.5, 10000) 88 PL_TEST_SYM_INTERVAL (cospi, 0.5, 0x1p51f, 10000) 89 PL_TEST_SYM_INTERVAL (cospi, 0x1p51f, inf, 10000) 90