1 /* 2 * Double-precision scalar sincospi function. 3 * 4 * Copyright (c) 2024, 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 "test_sig.h" 11 #include "test_defs.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 const static struct sincospi_data 20 { 21 double poly[10]; 22 } sincospi_data = { 23 /* Taylor series coefficents for sin(pi * x). */ 24 .poly = { 0x1.921fb54442d184p1, -0x1.2aef39896f94bp0, 0x1.466bc6775ab16p1, 25 -0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4, -0x1.e30750a28c88ep-8, 26 0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16, 27 0x1.af86ae521260bp-21, -0x1.012a9870eeb7dp-25 }, 28 }; 29 30 /* Top 12 bits of a double (sign and exponent bits). */ 31 static inline uint64_t 32 abstop12 (double x) 33 { 34 return (asuint64 (x) >> 52) & 0x7ff; 35 } 36 37 /* Triages special cases into 4 categories: 38 -1 or +1 if iy represents half an integer 39 -1 if round(y) is odd. 40 +1 if round(y) is even. 41 -2 or +2 if iy represents and integer. 42 -2 if iy is odd. 43 +2 if iy is even. 44 The argument is the bit representation of a positive non-zero 45 finite floating-point value which is either a half or an integer. */ 46 static inline int 47 checkint (uint64_t iy) 48 { 49 int e = iy >> 52; 50 if (e > 0x3ff + 52) 51 return 2; 52 if (iy & ((1ULL << (0x3ff + 52 - e)) - 1)) 53 { 54 if ((iy - 1) & 2) 55 return -1; 56 else 57 return 1; 58 } 59 if (iy & (1 << (0x3ff + 52 - e))) 60 return -2; 61 return 2; 62 } 63 64 /* Approximation for scalar double-precision sincospi(x). 65 Maximum error for sin: 3.46 ULP: 66 sincospif_sin(0x1.3d8a067cd8961p+14) got 0x1.ffe609a279008p-1 want 67 0x1.ffe609a27900cp-1. 68 Maximum error for cos: 3.66 ULP: 69 sincospif_cos(0x1.a0ec6997557eep-24) got 0x1.ffffffffffe59p-1 want 70 0x1.ffffffffffe5dp-1. */ 71 void 72 arm_math_sincospi (double x, double *out_sin, double *out_cos) 73 { 74 const struct sincospi_data *d = ptr_barrier (&sincospi_data); 75 uint64_t sign = asuint64 (x) & 0x8000000000000000; 76 77 if (likely (abstop12 (x) < abstop12 (0x1p51))) 78 { 79 /* ax = |x| - n (range reduction into -1/2 .. 1/2). */ 80 double ar_s = x - rint (x); 81 82 /* We know that cospi(x) = sinpi(0.5 - x) 83 range reduction and offset into sinpi range -1/2 .. 1/2 84 ax = 0.5 - |x - rint(x)|. */ 85 double ar_c = 0.5 - fabs (ar_s); 86 87 /* ss = sin(pi * ax). */ 88 double ar2_s = ar_s * ar_s; 89 double ar2_c = ar_c * ar_c; 90 double ar4_s = ar2_s * ar2_s; 91 double ar4_c = ar2_c * ar2_c; 92 93 uint64_t cc_sign = ((uint64_t) llrint (x)) << 63; 94 uint64_t ss_sign = cc_sign; 95 if (ar_s == 0) 96 ss_sign = sign; 97 98 double ss = pw_horner_9_f64 (ar2_s, ar4_s, d->poly); 99 double cc = pw_horner_9_f64 (ar2_c, ar4_c, d->poly); 100 101 /* As all values are reduced to -1/2 .. 1/2, the result of cos(x) 102 always be positive, therefore, the sign must be introduced 103 based upon if x rounds to odd or even. For sin(x) the sign is 104 copied from x. */ 105 *out_sin 106 = asdouble (asuint64 (fma (-4 * ar2_s, ar_s, ss * ar_s)) ^ ss_sign); 107 *out_cos 108 = asdouble (asuint64 (fma (-4 * ar2_c, ar_c, cc * ar_c)) ^ cc_sign); 109 } 110 else 111 { 112 /* When abs(x) > 0x1p51, the x will be either 113 - Half integer (relevant if abs(x) in [0x1p51, 0x1p52]) 114 - Odd integer (relevant if abs(x) in [0x1p52, 0x1p53]) 115 - Even integer (relevant if abs(x) in [0x1p53, inf]) 116 - Inf or NaN. */ 117 if (abstop12 (x) >= 0x7ff) 118 { 119 double inv_result = __math_invalid (x); 120 *out_sin = inv_result; 121 *out_cos = inv_result; 122 return; 123 } 124 else 125 { 126 uint64_t ax = asuint64 (x) & 0x7fffffffffffffff; 127 int m = checkint (ax); 128 /* The case where ax is half integer. */ 129 if (m & 1) 130 { 131 *out_sin = sign ? -m : m; 132 *out_cos = 0; 133 return; 134 } 135 /* The case where ax is integer. */ 136 else 137 { 138 *out_sin = asdouble (sign); 139 *out_cos = m >> 1; 140 return; 141 } 142 } 143 } 144 } 145 146 #if WANT_TRIGPI_TESTS 147 TEST_DISABLE_FENV (arm_math_sincospi_sin) 148 TEST_DISABLE_FENV (arm_math_sincospi_cos) 149 TEST_ULP (arm_math_sincospi_sin, 2.96) 150 TEST_ULP (arm_math_sincospi_cos, 3.16) 151 # define SINCOS_INTERVAL(lo, hi, n) \ 152 TEST_SYM_INTERVAL (arm_math_sincospi_sin, lo, hi, n) \ 153 TEST_SYM_INTERVAL (arm_math_sincospi_cos, lo, hi, n) 154 SINCOS_INTERVAL (0, 0x1p-63, 10000) 155 SINCOS_INTERVAL (0x1p-63, 0.5, 50000) 156 SINCOS_INTERVAL (0.5, 0x1p51, 50000) 157 SINCOS_INTERVAL (0x1p51, inf, 10000) 158 #endif 159