1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Core approximation for single-precision vector sincos
3*f3087befSAndrew Turner *
4*f3087befSAndrew Turner * Copyright (c) 2023-2024, Arm Limited.
5*f3087befSAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*f3087befSAndrew Turner */
7*f3087befSAndrew Turner
8*f3087befSAndrew Turner #include "sv_math.h"
9*f3087befSAndrew Turner
10*f3087befSAndrew Turner const static struct sv_sincosf_data
11*f3087befSAndrew Turner {
12*f3087befSAndrew Turner float poly_sin[3], poly_cos[3], pio2[3], inv_pio2, shift, range_val;
13*f3087befSAndrew Turner } sv_sincosf_data = {
14*f3087befSAndrew Turner .poly_sin = { /* Generated using Remez, odd coeffs only, in [-pi/4, pi/4]. */
15*f3087befSAndrew Turner -0x1.555546p-3, 0x1.11076p-7, -0x1.994eb4p-13 },
16*f3087befSAndrew Turner .poly_cos = { /* Generated using Remez, even coeffs only, in [-pi/4, pi/4]. */
17*f3087befSAndrew Turner 0x1.55554ap-5, -0x1.6c0c1ap-10, 0x1.99e0eep-16 },
18*f3087befSAndrew Turner .pio2 = { 0x1.921fb6p+0f, -0x1.777a5cp-25f, -0x1.ee59dap-50f },
19*f3087befSAndrew Turner .inv_pio2 = 0x1.45f306p-1f,
20*f3087befSAndrew Turner .shift = 0x1.8p23,
21*f3087befSAndrew Turner .range_val = 0x1p20
22*f3087befSAndrew Turner };
23*f3087befSAndrew Turner
24*f3087befSAndrew Turner static inline svbool_t
check_ge_rangeval(svbool_t pg,svfloat32_t x,const struct sv_sincosf_data * d)25*f3087befSAndrew Turner check_ge_rangeval (svbool_t pg, svfloat32_t x, const struct sv_sincosf_data *d)
26*f3087befSAndrew Turner {
27*f3087befSAndrew Turner svbool_t in_bounds = svaclt (pg, x, d->range_val);
28*f3087befSAndrew Turner return svnot_z (pg, in_bounds);
29*f3087befSAndrew Turner }
30*f3087befSAndrew Turner
31*f3087befSAndrew Turner /* Single-precision vector function allowing calculation of both sin and cos in
32*f3087befSAndrew Turner one function call, using shared argument reduction and separate low-order
33*f3087befSAndrew Turner polynomials.
34*f3087befSAndrew Turner Worst-case error for sin is 1.67 ULP:
35*f3087befSAndrew Turner sv_sincosf_sin(0x1.c704c4p+19) got 0x1.fff698p-5 want 0x1.fff69cp-5
36*f3087befSAndrew Turner Worst-case error for cos is 1.81 ULP:
37*f3087befSAndrew Turner sv_sincosf_cos(0x1.e506fp+19) got -0x1.ffec6ep-6 want -0x1.ffec72p-6. */
38*f3087befSAndrew Turner static inline svfloat32x2_t
sv_sincosf_inline(svbool_t pg,svfloat32_t x,const struct sv_sincosf_data * d)39*f3087befSAndrew Turner sv_sincosf_inline (svbool_t pg, svfloat32_t x, const struct sv_sincosf_data *d)
40*f3087befSAndrew Turner {
41*f3087befSAndrew Turner /* n = rint ( x / (pi/2) ). */
42*f3087befSAndrew Turner svfloat32_t q = svmla_x (pg, sv_f32 (d->shift), x, d->inv_pio2);
43*f3087befSAndrew Turner q = svsub_x (pg, q, d->shift);
44*f3087befSAndrew Turner svint32_t n = svcvt_s32_x (pg, q);
45*f3087befSAndrew Turner
46*f3087befSAndrew Turner /* Reduce x such that r is in [ -pi/4, pi/4 ]. */
47*f3087befSAndrew Turner svfloat32_t r = x;
48*f3087befSAndrew Turner r = svmls_x (pg, r, q, d->pio2[0]);
49*f3087befSAndrew Turner r = svmls_x (pg, r, q, d->pio2[1]);
50*f3087befSAndrew Turner r = svmls_x (pg, r, q, d->pio2[2]);
51*f3087befSAndrew Turner
52*f3087befSAndrew Turner /* Approximate sin(r) ~= r + r^3 * poly_sin(r^2). */
53*f3087befSAndrew Turner svfloat32_t r2 = svmul_x (pg, r, r), r3 = svmul_x (pg, r, r2);
54*f3087befSAndrew Turner svfloat32_t s = svmla_x (pg, sv_f32 (d->poly_sin[1]), r2, d->poly_sin[2]);
55*f3087befSAndrew Turner s = svmad_x (pg, r2, s, d->poly_sin[0]);
56*f3087befSAndrew Turner s = svmla_x (pg, r, r3, s);
57*f3087befSAndrew Turner
58*f3087befSAndrew Turner /* Approximate cos(r) ~= 1 - (r^2)/2 + r^4 * poly_cos(r^2). */
59*f3087befSAndrew Turner svfloat32_t r4 = svmul_x (pg, r2, r2);
60*f3087befSAndrew Turner svfloat32_t p = svmla_x (pg, sv_f32 (d->poly_cos[1]), r2, d->poly_cos[2]);
61*f3087befSAndrew Turner svfloat32_t c = svmad_x (pg, sv_f32 (d->poly_cos[0]), r2, -0.5);
62*f3087befSAndrew Turner c = svmla_x (pg, c, r4, p);
63*f3087befSAndrew Turner c = svmad_x (pg, r2, c, 1);
64*f3087befSAndrew Turner
65*f3087befSAndrew Turner svuint32_t un = svreinterpret_u32 (n);
66*f3087befSAndrew Turner /* If odd quadrant, swap cos and sin. */
67*f3087befSAndrew Turner svbool_t swap = svcmpeq (pg, svlsl_x (pg, un, 31), 0);
68*f3087befSAndrew Turner svfloat32_t ss = svsel (swap, s, c);
69*f3087befSAndrew Turner svfloat32_t cc = svsel (swap, c, s);
70*f3087befSAndrew Turner
71*f3087befSAndrew Turner /* Fix signs according to quadrant.
72*f3087befSAndrew Turner ss = asfloat(asuint(ss) ^ ((n & 2) << 30))
73*f3087befSAndrew Turner cc = asfloat(asuint(cc) & (((n + 1) & 2) << 30)). */
74*f3087befSAndrew Turner svuint32_t sin_sign = svlsl_x (pg, svand_x (pg, un, 2), 30);
75*f3087befSAndrew Turner svuint32_t cos_sign = svlsl_x (
76*f3087befSAndrew Turner pg, svand_x (pg, svreinterpret_u32 (svadd_x (pg, n, 1)), 2), 30);
77*f3087befSAndrew Turner ss = svreinterpret_f32 (sveor_x (pg, svreinterpret_u32 (ss), sin_sign));
78*f3087befSAndrew Turner cc = svreinterpret_f32 (sveor_x (pg, svreinterpret_u32 (cc), cos_sign));
79*f3087befSAndrew Turner
80*f3087befSAndrew Turner return svcreate2 (ss, cc);
81*f3087befSAndrew Turner }
82