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