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