1 /* 2 * Double-precision vector sin function. 3 * 4 * Copyright (c) 2019, Arm Limited. 5 * SPDX-License-Identifier: MIT 6 */ 7 8 #include "mathlib.h" 9 #include "v_math.h" 10 #if V_SUPPORTED 11 12 static const double Poly[] = { 13 /* worst-case error is 3.5 ulp. 14 abs error: 0x1.be222a58p-53 in [-pi/2, pi/2]. */ 15 -0x1.9f4a9c8b21dc9p-41, 16 0x1.60e88a10163f2p-33, 17 -0x1.ae6361b7254e7p-26, 18 0x1.71de382e8d62bp-19, 19 -0x1.a01a019aeb4ffp-13, 20 0x1.111111110b25ep-7, 21 -0x1.55555555554c3p-3, 22 }; 23 24 #define C7 v_f64 (Poly[0]) 25 #define C6 v_f64 (Poly[1]) 26 #define C5 v_f64 (Poly[2]) 27 #define C4 v_f64 (Poly[3]) 28 #define C3 v_f64 (Poly[4]) 29 #define C2 v_f64 (Poly[5]) 30 #define C1 v_f64 (Poly[6]) 31 32 #define InvPi v_f64 (0x1.45f306dc9c883p-2) 33 #define Pi1 v_f64 (0x1.921fb54442d18p+1) 34 #define Pi2 v_f64 (0x1.1a62633145c06p-53) 35 #define Pi3 v_f64 (0x1.c1cd129024e09p-106) 36 #define Shift v_f64 (0x1.8p52) 37 #define RangeVal v_f64 (0x1p23) 38 #define AbsMask v_u64 (0x7fffffffffffffff) 39 40 VPCS_ATTR 41 __attribute__ ((noinline)) static v_f64_t 42 specialcase (v_f64_t x, v_f64_t y, v_u64_t cmp) 43 { 44 return v_call_f64 (sin, x, y, cmp); 45 } 46 47 VPCS_ATTR 48 v_f64_t 49 V_NAME(sin) (v_f64_t x) 50 { 51 v_f64_t n, r, r2, y; 52 v_u64_t sign, odd, cmp; 53 54 r = v_as_f64_u64 (v_as_u64_f64 (x) & AbsMask); 55 sign = v_as_u64_f64 (x) & ~AbsMask; 56 cmp = v_cond_u64 (v_as_u64_f64 (r) >= v_as_u64_f64 (RangeVal)); 57 58 /* n = rint(|x|/pi). */ 59 n = v_fma_f64 (InvPi, r, Shift); 60 odd = v_as_u64_f64 (n) << 63; 61 n -= Shift; 62 63 /* r = |x| - n*pi (range reduction into -pi/2 .. pi/2). */ 64 r = v_fma_f64 (-Pi1, n, r); 65 r = v_fma_f64 (-Pi2, n, r); 66 r = v_fma_f64 (-Pi3, n, r); 67 68 /* sin(r) poly approx. */ 69 r2 = r * r; 70 y = v_fma_f64 (C7, r2, C6); 71 y = v_fma_f64 (y, r2, C5); 72 y = v_fma_f64 (y, r2, C4); 73 y = v_fma_f64 (y, r2, C3); 74 y = v_fma_f64 (y, r2, C2); 75 y = v_fma_f64 (y, r2, C1); 76 y = v_fma_f64 (y * r2, r, r); 77 78 /* sign. */ 79 y = v_as_f64_u64 (v_as_u64_f64 (y) ^ sign ^ odd); 80 81 if (unlikely (v_any_u64 (cmp))) 82 return specialcase (x, y, cmp); 83 return y; 84 } 85 VPCS_ALIAS 86 #endif 87