]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/aarch64/fpu/tan_sve.c
1 /* Double-precision vector (SVE) tan function
3 Copyright (C) 2023 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
11 The GNU C Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, see
18 <https://www.gnu.org/licenses/>. */
21 #include "poly_sve_f64.h"
23 static const struct data
26 double half_pi_hi
, half_pi_lo
, inv_half_pi
, range_val
, shift
;
28 /* Polynomial generated with FPMinimax. */
29 .poly
= { 0x1.5555555555556p
-2, 0x1.1111111110a63p
-3, 0x1.ba1ba1bb46414p
-5,
30 0x1.664f47e5b5445p
-6, 0x1.226e5e5ecdfa3p
-7, 0x1.d6c7ddbf87047p
-9,
31 0x1.7ea75d05b583ep
-10, 0x1.289f22964a03cp
-11,
32 0x1.4e4fd14147622p
-12, },
33 .half_pi_hi
= 0x1.921fb54442d18p0
,
34 .half_pi_lo
= 0x1.1a62633145c07p
-54,
35 .inv_half_pi
= 0x1.45f306dc9c883p
-1,
40 static svfloat64_t NOINLINE
41 special_case (svfloat64_t x
, svfloat64_t y
, svbool_t special
)
43 return sv_call_f64 (tan
, x
, y
, special
);
46 /* Vector approximation for double-precision tan.
47 Maximum measured error is 3.48 ULP:
48 _ZGVsMxv_tan(0x1.4457047ef78d8p+20) got -0x1.f6ccd8ecf7dedp+37
49 want -0x1.f6ccd8ecf7deap+37. */
50 svfloat64_t
SV_NAME_D1 (tan
) (svfloat64_t x
, svbool_t pg
)
52 const struct data
*dat
= ptr_barrier (&data
);
54 /* Invert condition to catch NaNs and Infs as well as large values. */
55 svbool_t special
= svnot_z (pg
, svaclt (pg
, x
, dat
->range_val
));
57 /* q = nearest integer to 2 * x / pi. */
58 svfloat64_t shift
= sv_f64 (dat
->shift
);
59 svfloat64_t q
= svmla_x (pg
, shift
, x
, dat
->inv_half_pi
);
60 q
= svsub_x (pg
, q
, shift
);
61 svint64_t qi
= svcvt_s64_x (pg
, q
);
63 /* Use q to reduce x to r in [-pi/4, pi/4], by:
64 r = x - q * pi/2, in extended precision. */
66 svfloat64_t half_pi
= svld1rq (svptrue_b64 (), &dat
->half_pi_hi
);
67 r
= svmls_lane (r
, q
, half_pi
, 0);
68 r
= svmls_lane (r
, q
, half_pi
, 1);
69 /* Further reduce r to [-pi/8, pi/8], to be reconstructed using double angle
71 r
= svmul_x (pg
, r
, 0.5);
73 /* Approximate tan(r) using order 8 polynomial.
74 tan(x) is odd, so polynomial has the form:
75 tan(x) ~= x + C0 * x^3 + C1 * x^5 + C3 * x^7 + ...
76 Hence we first approximate P(r) = C1 + C2 * r^2 + C3 * r^4 + ...
77 Then compute the approximation by:
78 tan(r) ~= r + r^3 * (C0 + r^2 * P(r)). */
79 svfloat64_t r2
= svmul_x (pg
, r
, r
);
80 svfloat64_t r4
= svmul_x (pg
, r2
, r2
);
81 svfloat64_t r8
= svmul_x (pg
, r4
, r4
);
82 /* Use offset version coeff array by 1 to evaluate from C1 onwards. */
83 svfloat64_t p
= sv_estrin_7_f64_x (pg
, r2
, r4
, r8
, dat
->poly
+ 1);
84 p
= svmad_x (pg
, p
, r2
, dat
->poly
[0]);
85 p
= svmla_x (pg
, r
, r2
, svmul_x (pg
, p
, r
));
87 /* Recombination uses double-angle formula:
88 tan(2x) = 2 * tan(x) / (1 - (tan(x))^2)
89 and reciprocity around pi/2:
90 tan(x) = 1 / (tan(pi/2 - x))
91 to assemble result using change-of-sign and conditional selection of
92 numerator/denominator dependent on odd/even-ness of q (hence quadrant). */
94 = svcmpeq (pg
, svand_x (pg
, svreinterpret_u64 (qi
), 1), 0);
96 svfloat64_t n
= svmad_x (pg
, p
, p
, -1);
97 svfloat64_t d
= svmul_x (pg
, p
, 2);
99 n
= svneg_m (n
, use_recip
, d
);
100 d
= svsel (use_recip
, swap
, d
);
101 if (__glibc_unlikely (svptest_any (pg
, special
)))
102 return special_case (x
, svdiv_x (svnot_z (pg
, special
), n
, d
), special
);
103 return svdiv_x (pg
, n
, d
);