1*5a02ffc3SAndrew Turner /*
2*5a02ffc3SAndrew Turner * Core approximation for double-precision vector sincos
3*5a02ffc3SAndrew Turner *
4*5a02ffc3SAndrew Turner * Copyright (c) 2023, Arm Limited.
5*5a02ffc3SAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*5a02ffc3SAndrew Turner */
7*5a02ffc3SAndrew Turner
8*5a02ffc3SAndrew Turner #include "sv_math.h"
9*5a02ffc3SAndrew Turner #include "poly_sve_f64.h"
10*5a02ffc3SAndrew Turner
11*5a02ffc3SAndrew Turner static const struct sv_sincos_data
12*5a02ffc3SAndrew Turner {
13*5a02ffc3SAndrew Turner double sin_poly[7], cos_poly[6], pio2[3];
14*5a02ffc3SAndrew Turner double inv_pio2, shift, range_val;
15*5a02ffc3SAndrew Turner } sv_sincos_data = {
16*5a02ffc3SAndrew Turner .inv_pio2 = 0x1.45f306dc9c882p-1,
17*5a02ffc3SAndrew Turner .pio2 = { 0x1.921fb50000000p+0, 0x1.110b460000000p-26,
18*5a02ffc3SAndrew Turner 0x1.1a62633145c07p-54 },
19*5a02ffc3SAndrew Turner .shift = 0x1.8p52,
20*5a02ffc3SAndrew Turner .sin_poly = { /* Computed using Remez in [-pi/2, pi/2]. */
21*5a02ffc3SAndrew Turner -0x1.555555555547bp-3, 0x1.1111111108a4dp-7,
22*5a02ffc3SAndrew Turner -0x1.a01a019936f27p-13, 0x1.71de37a97d93ep-19,
23*5a02ffc3SAndrew Turner -0x1.ae633919987c6p-26, 0x1.60e277ae07cecp-33,
24*5a02ffc3SAndrew Turner -0x1.9e9540300a1p-41 },
25*5a02ffc3SAndrew Turner .cos_poly = { /* Computed using Remez in [-pi/4, pi/4]. */
26*5a02ffc3SAndrew Turner 0x1.555555555554cp-5, -0x1.6c16c16c1521fp-10,
27*5a02ffc3SAndrew Turner 0x1.a01a019cbf62ap-16, -0x1.27e4f812b681ep-22,
28*5a02ffc3SAndrew Turner 0x1.1ee9f152a57cdp-29, -0x1.8fb131098404bp-37 },
29*5a02ffc3SAndrew Turner .range_val = 0x1p23, };
30*5a02ffc3SAndrew Turner
31*5a02ffc3SAndrew Turner static inline svbool_t
check_ge_rangeval(svbool_t pg,svfloat64_t x,const struct sv_sincos_data * d)32*5a02ffc3SAndrew Turner check_ge_rangeval (svbool_t pg, svfloat64_t x, const struct sv_sincos_data *d)
33*5a02ffc3SAndrew Turner {
34*5a02ffc3SAndrew Turner svbool_t in_bounds = svaclt (pg, x, d->range_val);
35*5a02ffc3SAndrew Turner return svnot_z (pg, in_bounds);
36*5a02ffc3SAndrew Turner }
37*5a02ffc3SAndrew Turner
38*5a02ffc3SAndrew Turner /* Double-precision vector function allowing calculation of both sin and cos in
39*5a02ffc3SAndrew Turner one function call, using shared argument reduction and separate polynomials.
40*5a02ffc3SAndrew Turner Largest observed error is for sin, 3.22 ULP:
41*5a02ffc3SAndrew Turner v_sincos_sin (0x1.d70eef40f39b1p+12) got -0x1.ffe9537d5dbb7p-3
42*5a02ffc3SAndrew Turner want -0x1.ffe9537d5dbb4p-3. */
43*5a02ffc3SAndrew Turner static inline svfloat64x2_t
sv_sincos_inline(svbool_t pg,svfloat64_t x,const struct sv_sincos_data * d)44*5a02ffc3SAndrew Turner sv_sincos_inline (svbool_t pg, svfloat64_t x, const struct sv_sincos_data *d)
45*5a02ffc3SAndrew Turner {
46*5a02ffc3SAndrew Turner /* q = nearest integer to 2 * x / pi. */
47*5a02ffc3SAndrew Turner svfloat64_t q = svsub_x (pg, svmla_x (pg, sv_f64 (d->shift), x, d->inv_pio2),
48*5a02ffc3SAndrew Turner d->shift);
49*5a02ffc3SAndrew Turner svint64_t n = svcvt_s64_x (pg, q);
50*5a02ffc3SAndrew Turner
51*5a02ffc3SAndrew Turner /* Reduce x such that r is in [ -pi/4, pi/4 ]. */
52*5a02ffc3SAndrew Turner svfloat64_t r = x;
53*5a02ffc3SAndrew Turner r = svmls_x (pg, r, q, d->pio2[0]);
54*5a02ffc3SAndrew Turner r = svmls_x (pg, r, q, d->pio2[1]);
55*5a02ffc3SAndrew Turner r = svmls_x (pg, r, q, d->pio2[2]);
56*5a02ffc3SAndrew Turner
57*5a02ffc3SAndrew Turner svfloat64_t r2 = svmul_x (pg, r, r), r3 = svmul_x (pg, r2, r),
58*5a02ffc3SAndrew Turner r4 = svmul_x (pg, r2, r2);
59*5a02ffc3SAndrew Turner
60*5a02ffc3SAndrew Turner /* Approximate sin(r) ~= r + r^3 * poly_sin(r^2). */
61*5a02ffc3SAndrew Turner svfloat64_t s = sv_pw_horner_6_f64_x (pg, r2, r4, d->sin_poly);
62*5a02ffc3SAndrew Turner s = svmla_x (pg, r, r3, s);
63*5a02ffc3SAndrew Turner
64*5a02ffc3SAndrew Turner /* Approximate cos(r) ~= 1 - (r^2)/2 + r^4 * poly_cos(r^2). */
65*5a02ffc3SAndrew Turner svfloat64_t c = sv_pw_horner_5_f64_x (pg, r2, r4, d->cos_poly);
66*5a02ffc3SAndrew Turner c = svmad_x (pg, c, r2, -0.5);
67*5a02ffc3SAndrew Turner c = svmad_x (pg, c, r2, 1);
68*5a02ffc3SAndrew Turner
69*5a02ffc3SAndrew Turner svuint64_t un = svreinterpret_u64 (n);
70*5a02ffc3SAndrew Turner /* If odd quadrant, swap cos and sin. */
71*5a02ffc3SAndrew Turner svbool_t swap = svcmpeq (pg, svlsl_x (pg, un, 63), 0);
72*5a02ffc3SAndrew Turner svfloat64_t ss = svsel (swap, s, c);
73*5a02ffc3SAndrew Turner svfloat64_t cc = svsel (swap, c, s);
74*5a02ffc3SAndrew Turner
75*5a02ffc3SAndrew Turner /* Fix signs according to quadrant.
76*5a02ffc3SAndrew Turner ss = asdouble(asuint64(ss) ^ ((n & 2) << 62))
77*5a02ffc3SAndrew Turner cc = asdouble(asuint64(cc) & (((n + 1) & 2) << 62)). */
78*5a02ffc3SAndrew Turner svuint64_t sin_sign = svlsl_x (pg, svand_x (pg, un, 2), 62);
79*5a02ffc3SAndrew Turner svuint64_t cos_sign = svlsl_x (
80*5a02ffc3SAndrew Turner pg, svand_x (pg, svreinterpret_u64 (svadd_x (pg, n, 1)), 2), 62);
81*5a02ffc3SAndrew Turner ss = svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (ss), sin_sign));
82*5a02ffc3SAndrew Turner cc = svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (cc), cos_sign));
83*5a02ffc3SAndrew Turner
84*5a02ffc3SAndrew Turner return svcreate2 (ss, cc);
85*5a02ffc3SAndrew Turner }
86