1 /* 2 * Single-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 "math_config.h" 9 #include "test_sig.h" 10 #include "test_defs.h" 11 #include "poly_scalar_f32.h" 12 13 /* Taylor series coefficents for sin(pi * x). */ 14 const static struct sincospif_data 15 { 16 float poly[6]; 17 } sincospif_data = { 18 /* Taylor series coefficents for sin(pi * x). */ 19 .poly = { 0x1.921fb6p1f, -0x1.4abbcep2f, 0x1.466bc6p1f, -0x1.32d2ccp-1f, 20 0x1.50783p-4f, -0x1.e30750p-8f }, 21 }; 22 23 /* Top 12 bits of the float representation with the sign bit cleared. */ 24 static inline uint32_t 25 abstop12 (float x) 26 { 27 return (asuint (x) >> 20) & 0x7ff; 28 } 29 30 /* Triages special cases into 4 categories: 31 -1 or +1 if iy represents half an integer 32 -1 if round(y) is odd. 33 +1 if round(y) is even. 34 -2 or +2 if iy represents and integer. 35 -2 if iy is odd. 36 +2 if iy is even. 37 The argument is the bit representation of a positive non-zero 38 finite floating-point value which is either a half or an integer. */ 39 static inline int 40 checkint (uint32_t iy) 41 { 42 int e = iy >> 23; 43 if (e > 0x7f + 23) 44 return 2; 45 if (iy & ((1 << (0x7f + 23 - e)) - 1)) 46 { 47 if ((iy - 1) & 2) 48 return -1; 49 else 50 return 1; 51 } 52 if (iy & (1 << (0x7f + 23 - e))) 53 return -2; 54 return 2; 55 } 56 57 /* Approximation for scalar single-precision sincospif(x). 58 Maximum error for sin: 3.04 ULP: 59 sincospif_sin(0x1.c597ccp-2) got 0x1.f7cd56p-1 want 0x1.f7cd5p-1. 60 Maximum error for cos: 3.18 ULP: 61 sincospif_cos(0x1.d341a8p-5) got 0x1.f7cd56p-1 want 0x1.f7cd5p-1. */ 62 void 63 arm_math_sincospif (float x, float *out_sin, float *out_cos) 64 { 65 66 const struct sincospif_data *d = ptr_barrier (&sincospif_data); 67 uint32_t sign = asuint (x) & 0x80000000; 68 69 /* abs(x) in [0, 0x1p22]. */ 70 if (likely (abstop12 (x) < abstop12 (0x1p22))) 71 { 72 /* ar_s = x - n (range reduction into -1/2 .. 1/2). */ 73 float ar_s = x - rintf (x); 74 /* We know that cospi(x) = sinpi(0.5 - x) 75 range reduction and offset into sinpi range -1/2 .. 1/2 76 ar_c = 0.5 - |x - n|. */ 77 float ar_c = 0.5f - fabsf (ar_s); 78 79 float ar2_s = ar_s * ar_s; 80 float ar2_c = ar_c * ar_c; 81 float ar4_s = ar2_s * ar2_s; 82 float ar4_c = ar2_c * ar2_c; 83 84 uint32_t cc_sign = lrintf (x) << 31; 85 uint32_t ss_sign = cc_sign; 86 if (ar_s == 0) 87 ss_sign = sign; 88 89 /* As all values are reduced to -1/2 .. 1/2, the result of cos(x) 90 always be positive, therefore, the sign must be introduced 91 based upon if x rounds to odd or even. For sin(x) the sign is 92 copied from x. */ 93 *out_sin = pw_horner_5_f32 (ar2_s, ar4_s, d->poly) 94 * asfloat (asuint (ar_s) ^ ss_sign); 95 *out_cos = pw_horner_5_f32 (ar2_c, ar4_c, d->poly) 96 * asfloat (asuint (ar_c) ^ cc_sign); 97 return; 98 } 99 else 100 { 101 /* When abs(x) > 0x1p22, the x will be either 102 - Half integer (relevant if abs(x) in [0x1p22, 0x1p23]) 103 - Odd integer (relevant if abs(x) in [0x1p22, 0x1p24]) 104 - Even integer (relevant if abs(x) in [0x1p22, inf]) 105 - Inf or NaN. */ 106 if (abstop12 (x) >= 0x7f8) 107 { 108 float inv_result = __math_invalidf (x); 109 *out_sin = inv_result; 110 *out_cos = inv_result; 111 return; 112 } 113 else 114 { 115 uint32_t ax = asuint (x) & 0x7fffffff; 116 int m = checkint (ax); 117 if (m & 1) 118 { 119 *out_sin = sign ? -m : m; 120 *out_cos = 0; 121 return; 122 } 123 else 124 { 125 *out_sin = asfloat (sign); 126 *out_cos = m >> 1; 127 return; 128 } 129 } 130 } 131 } 132 133 #if WANT_TRIGPI_TESTS 134 TEST_DISABLE_FENV (arm_math_sincospif_sin) 135 TEST_DISABLE_FENV (arm_math_sincospif_cos) 136 TEST_ULP (arm_math_sincospif_sin, 2.54) 137 TEST_ULP (arm_math_sincospif_cos, 2.68) 138 # define SINCOSPIF_INTERVAL(lo, hi, n) \ 139 TEST_SYM_INTERVAL (arm_math_sincospif_sin, lo, hi, n) \ 140 TEST_SYM_INTERVAL (arm_math_sincospif_cos, lo, hi, n) 141 SINCOSPIF_INTERVAL (0, 0x1p-31, 10000) 142 SINCOSPIF_INTERVAL (0x1p-31, 1, 50000) 143 SINCOSPIF_INTERVAL (1, 0x1p22f, 50000) 144 SINCOSPIF_INTERVAL (0x1p22f, inf, 10000) 145 #endif 146