1 /* 2 * Double-precision vector cos 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 HalfPi v_f64 (0x1.921fb54442d18p+0) 34 #define Pi1 v_f64 (0x1.921fb54442d18p+1) 35 #define Pi2 v_f64 (0x1.1a62633145c06p-53) 36 #define Pi3 v_f64 (0x1.c1cd129024e09p-106) 37 #define Shift v_f64 (0x1.8p52) 38 #define RangeVal v_f64 (0x1p23) 39 #define AbsMask v_u64 (0x7fffffffffffffff) 40 41 VPCS_ATTR 42 __attribute__ ((noinline)) static v_f64_t 43 specialcase (v_f64_t x, v_f64_t y, v_u64_t cmp) 44 { 45 return v_call_f64 (cos, x, y, cmp); 46 } 47 48 VPCS_ATTR 49 v_f64_t 50 V_NAME(cos) (v_f64_t x) 51 { 52 v_f64_t n, r, r2, y; 53 v_u64_t odd, cmp; 54 55 r = v_as_f64_u64 (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/2)/pi) - 0.5. */ 59 n = v_fma_f64 (InvPi, r + HalfPi, Shift); 60 odd = v_as_u64_f64 (n) << 63; 61 n -= Shift; 62 n -= v_f64 (0.5); 63 64 /* r = |x| - n*pi (range reduction into -pi/2 .. pi/2). */ 65 r = v_fma_f64 (-Pi1, n, r); 66 r = v_fma_f64 (-Pi2, n, r); 67 r = v_fma_f64 (-Pi3, n, r); 68 69 /* sin(r) poly approx. */ 70 r2 = r * r; 71 y = v_fma_f64 (C7, r2, C6); 72 y = v_fma_f64 (y, r2, C5); 73 y = v_fma_f64 (y, r2, C4); 74 y = v_fma_f64 (y, r2, C3); 75 y = v_fma_f64 (y, r2, C2); 76 y = v_fma_f64 (y, r2, C1); 77 y = v_fma_f64 (y * r2, r, r); 78 79 /* sign. */ 80 y = v_as_f64_u64 (v_as_u64_f64 (y) ^ odd); 81 82 if (unlikely (v_any_u64 (cmp))) 83 return specialcase (x, y, cmp); 84 return y; 85 } 86 VPCS_ALIAS 87 #endif 88