1072a4ba8SAndrew Turner /*
2072a4ba8SAndrew Turner * Single-precision SVE sin(x) function.
3072a4ba8SAndrew Turner *
4072a4ba8SAndrew Turner * Copyright (c) 2019-2023, Arm Limited.
5072a4ba8SAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6072a4ba8SAndrew Turner */
7072a4ba8SAndrew Turner
8072a4ba8SAndrew Turner #include "sv_math.h"
9072a4ba8SAndrew Turner #include "pl_sig.h"
10072a4ba8SAndrew Turner #include "pl_test.h"
11072a4ba8SAndrew Turner
12*5a02ffc3SAndrew Turner static const struct data
13*5a02ffc3SAndrew Turner {
14*5a02ffc3SAndrew Turner float poly[4];
15*5a02ffc3SAndrew Turner /* Pi-related values to be loaded as one quad-word and used with
16*5a02ffc3SAndrew Turner svmla_lane. */
17*5a02ffc3SAndrew Turner float negpi1, negpi2, negpi3, invpi;
18*5a02ffc3SAndrew Turner float shift;
19*5a02ffc3SAndrew Turner } data = {
20*5a02ffc3SAndrew Turner .poly = {
21*5a02ffc3SAndrew Turner /* Non-zero coefficients from the degree 9 Taylor series expansion of
22*5a02ffc3SAndrew Turner sin. */
23*5a02ffc3SAndrew Turner -0x1.555548p-3f, 0x1.110df4p-7f, -0x1.9f42eap-13f, 0x1.5b2e76p-19f
24*5a02ffc3SAndrew Turner },
25*5a02ffc3SAndrew Turner .negpi1 = -0x1.921fb6p+1f,
26*5a02ffc3SAndrew Turner .negpi2 = 0x1.777a5cp-24f,
27*5a02ffc3SAndrew Turner .negpi3 = 0x1.ee59dap-49f,
28*5a02ffc3SAndrew Turner .invpi = 0x1.45f306p-2f,
29*5a02ffc3SAndrew Turner .shift = 0x1.8p+23f
30*5a02ffc3SAndrew Turner };
31072a4ba8SAndrew Turner
32*5a02ffc3SAndrew Turner #define RangeVal 0x49800000 /* asuint32 (0x1p20f). */
33*5a02ffc3SAndrew Turner #define C(i) sv_f32 (d->poly[i])
34072a4ba8SAndrew Turner
35*5a02ffc3SAndrew Turner static svfloat32_t NOINLINE
special_case(svfloat32_t x,svfloat32_t y,svbool_t cmp)36*5a02ffc3SAndrew Turner special_case (svfloat32_t x, svfloat32_t y, svbool_t cmp)
37072a4ba8SAndrew Turner {
38072a4ba8SAndrew Turner return sv_call_f32 (sinf, x, y, cmp);
39072a4ba8SAndrew Turner }
40072a4ba8SAndrew Turner
41072a4ba8SAndrew Turner /* A fast SVE implementation of sinf.
42072a4ba8SAndrew Turner Maximum error: 1.89 ULPs.
43072a4ba8SAndrew Turner This maximum error is achieved at multiple values in [-2^18, 2^18]
44072a4ba8SAndrew Turner but one example is:
45*5a02ffc3SAndrew Turner SV_NAME_F1 (sin)(0x1.9247a4p+0) got 0x1.fffff6p-1 want 0x1.fffffap-1. */
SV_NAME_F1(sin)46*5a02ffc3SAndrew Turner svfloat32_t SV_NAME_F1 (sin) (svfloat32_t x, const svbool_t pg)
47072a4ba8SAndrew Turner {
48*5a02ffc3SAndrew Turner const struct data *d = ptr_barrier (&data);
49072a4ba8SAndrew Turner
50*5a02ffc3SAndrew Turner svfloat32_t ax = svabs_x (pg, x);
51*5a02ffc3SAndrew Turner svuint32_t sign
52*5a02ffc3SAndrew Turner = sveor_x (pg, svreinterpret_u32 (x), svreinterpret_u32 (ax));
53*5a02ffc3SAndrew Turner svbool_t cmp = svcmpge (pg, svreinterpret_u32 (ax), RangeVal);
54*5a02ffc3SAndrew Turner
55*5a02ffc3SAndrew Turner /* pi_vals are a quad-word of helper values - the first 3 elements contain
56*5a02ffc3SAndrew Turner -pi in extended precision, the last contains 1 / pi. */
57*5a02ffc3SAndrew Turner svfloat32_t pi_vals = svld1rq (svptrue_b32 (), &d->negpi1);
58072a4ba8SAndrew Turner
59072a4ba8SAndrew Turner /* n = rint(|x|/pi). */
60*5a02ffc3SAndrew Turner svfloat32_t n = svmla_lane (sv_f32 (d->shift), ax, pi_vals, 3);
61*5a02ffc3SAndrew Turner svuint32_t odd = svlsl_x (pg, svreinterpret_u32 (n), 31);
62*5a02ffc3SAndrew Turner n = svsub_x (pg, n, d->shift);
63072a4ba8SAndrew Turner
64072a4ba8SAndrew Turner /* r = |x| - n*pi (range reduction into -pi/2 .. pi/2). */
65*5a02ffc3SAndrew Turner svfloat32_t r;
66*5a02ffc3SAndrew Turner r = svmla_lane (ax, n, pi_vals, 0);
67*5a02ffc3SAndrew Turner r = svmla_lane (r, n, pi_vals, 1);
68*5a02ffc3SAndrew Turner r = svmla_lane (r, n, pi_vals, 2);
69072a4ba8SAndrew Turner
70072a4ba8SAndrew Turner /* sin(r) approx using a degree 9 polynomial from the Taylor series
71072a4ba8SAndrew Turner expansion. Note that only the odd terms of this are non-zero. */
72*5a02ffc3SAndrew Turner svfloat32_t r2 = svmul_x (pg, r, r);
73*5a02ffc3SAndrew Turner svfloat32_t y;
74*5a02ffc3SAndrew Turner y = svmla_x (pg, C (2), r2, C (3));
75*5a02ffc3SAndrew Turner y = svmla_x (pg, C (1), r2, y);
76*5a02ffc3SAndrew Turner y = svmla_x (pg, C (0), r2, y);
77*5a02ffc3SAndrew Turner y = svmla_x (pg, r, r, svmul_x (pg, y, r2));
78072a4ba8SAndrew Turner
79072a4ba8SAndrew Turner /* sign = y^sign^odd. */
80*5a02ffc3SAndrew Turner sign = sveor_x (pg, sign, odd);
81072a4ba8SAndrew Turner
82072a4ba8SAndrew Turner if (unlikely (svptest_any (pg, cmp)))
83*5a02ffc3SAndrew Turner return special_case (x,
84*5a02ffc3SAndrew Turner svreinterpret_f32 (sveor_x (
85*5a02ffc3SAndrew Turner svnot_z (pg, cmp), svreinterpret_u32 (y), sign)),
86*5a02ffc3SAndrew Turner cmp);
87*5a02ffc3SAndrew Turner return svreinterpret_f32 (sveor_x (pg, svreinterpret_u32 (y), sign));
88072a4ba8SAndrew Turner }
89072a4ba8SAndrew Turner
90072a4ba8SAndrew Turner PL_SIG (SV, F, 1, sin, -3.1, 3.1)
91*5a02ffc3SAndrew Turner PL_TEST_ULP (SV_NAME_F1 (sin), 1.40)
92*5a02ffc3SAndrew Turner PL_TEST_SYM_INTERVAL (SV_NAME_F1 (sin), 0, 0x1p23, 1000000)
93*5a02ffc3SAndrew Turner PL_TEST_SYM_INTERVAL (SV_NAME_F1 (sin), 0x1p23, inf, 10000)
94