1 /* Single-precision vector (Advanced SIMD) tan function
3 Copyright (C) 2023-2024 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_advsimd_f32.h"
23 static const struct data
26 float32x4_t pi_consts
;
29 float32x4_t range_val
;
32 /* Coefficients generated using FPMinimax. */
33 .poly
= { V4 (0x1.55555p
-2f
), V4 (0x1.11166p
-3f
), V4 (0x1.b88a78p
-5f
),
34 V4 (0x1.7b5756p
-6f
), V4 (0x1.4ef4cep
-8f
), V4 (0x1.0e1e74p
-7f
) },
35 /* Stores constants: (-pi/2)_high, (-pi/2)_mid, (-pi/2)_low, and 2/pi. */
37 = { -0x1.921fb6p
+0f
, 0x1.777a5cp
-25f
, 0x1.ee59dap
-50f
, 0x1.45f306p
-1f
},
38 .shift
= V4 (0x1.8p
+23f
),
40 .range_val
= V4 (0x1p
15f
),
44 #define RangeVal v_u32 (0x47000000) /* asuint32(0x1p15f). */
45 #define TinyBound v_u32 (0x30000000) /* asuint32 (0x1p-31f). */
46 #define Thresh v_u32 (0x16000000) /* asuint32(RangeVal) - TinyBound. */
48 /* Special cases (fall back to scalar calls). */
49 static float32x4_t VPCS_ATTR NOINLINE
50 special_case (float32x4_t x
, float32x4_t y
, uint32x4_t cmp
)
52 return v_call_f32 (tanf
, x
, y
, cmp
);
55 /* Use a full Estrin scheme to evaluate polynomial. */
56 static inline float32x4_t
57 eval_poly (float32x4_t z
, const struct data
*d
)
59 float32x4_t z2
= vmulq_f32 (z
, z
);
61 /* Tiny z (<= 0x1p-31) will underflow when calculating z^4.
62 If fp exceptions are to be triggered correctly,
63 sidestep this by fixing such lanes to 0. */
65 = vcleq_u32 (vreinterpretq_u32_f32 (vabsq_f32 (z
)), TinyBound
);
66 if (__glibc_unlikely (v_any_u32 (will_uflow
)))
67 z2
= vbslq_f32 (will_uflow
, v_f32 (0), z2
);
69 float32x4_t z4
= vmulq_f32 (z2
, z2
);
70 return v_estrin_5_f32 (z
, z2
, z4
, d
->poly
);
73 /* Fast implementation of AdvSIMD tanf.
74 Maximum error is 3.45 ULP:
75 __v_tanf(-0x1.e5f0cap+13) got 0x1.ff9856p-1
76 want 0x1.ff9850p-1. */
77 float32x4_t VPCS_ATTR NOINLINE
V_NAME_F1 (tan
) (float32x4_t x
)
79 const struct data
*d
= ptr_barrier (&data
);
80 float32x4_t special_arg
= x
;
82 /* iax >= RangeVal means x, if not inf or NaN, is too large to perform fast
85 uint32x4_t iax
= vreinterpretq_u32_f32 (vabsq_f32 (x
));
86 /* If fp exceptions are to be triggered correctly, also special-case tiny
87 input, as this will load to overflow later. Fix any special lanes to 1 to
88 prevent any exceptions being triggered. */
89 uint32x4_t special
= vcgeq_u32 (vsubq_u32 (iax
, TinyBound
), Thresh
);
90 if (__glibc_unlikely (v_any_u32 (special
)))
91 x
= vbslq_f32 (special
, v_f32 (1.0f
), x
);
93 /* Otherwise, special-case large and special values. */
94 uint32x4_t special
= vcageq_f32 (x
, d
->range_val
);
97 /* n = rint(x/(pi/2)). */
98 float32x4_t q
= vfmaq_laneq_f32 (d
->shift
, x
, d
->pi_consts
, 3);
99 float32x4_t n
= vsubq_f32 (q
, d
->shift
);
100 /* Determine if x lives in an interval, where |tan(x)| grows to infinity. */
101 uint32x4_t pred_alt
= vtstq_u32 (vreinterpretq_u32_f32 (q
), v_u32 (1));
103 /* r = x - n * (pi/2) (range reduction into -pi./4 .. pi/4). */
105 r
= vfmaq_laneq_f32 (x
, n
, d
->pi_consts
, 0);
106 r
= vfmaq_laneq_f32 (r
, n
, d
->pi_consts
, 1);
107 r
= vfmaq_laneq_f32 (r
, n
, d
->pi_consts
, 2);
109 /* If x lives in an interval, where |tan(x)|
110 - is finite, then use a polynomial approximation of the form
111 tan(r) ~ r + r^3 * P(r^2) = r + r * r^2 * P(r^2).
112 - grows to infinity then use symmetries of tangent and the identity
113 tan(r) = cotan(pi/2 - r) to express tan(x) as 1/tan(-r). Finally, use
114 the same polynomial approximation of tan as above. */
116 /* Invert sign of r if odd quadrant. */
117 float32x4_t z
= vmulq_f32 (r
, vbslq_f32 (pred_alt
, v_f32 (-1), v_f32 (1)));
119 /* Evaluate polynomial approximation of tangent on [-pi/4, pi/4]. */
120 float32x4_t z2
= vmulq_f32 (r
, r
);
121 float32x4_t p
= eval_poly (z2
, d
);
122 float32x4_t y
= vfmaq_f32 (z
, vmulq_f32 (z
, z2
), p
);
124 /* Compute reciprocal and apply if required. */
125 float32x4_t inv_y
= vdivq_f32 (v_f32 (1.0f
), y
);
127 if (__glibc_unlikely (v_any_u32 (special
)))
128 return special_case (special_arg
, vbslq_f32 (pred_alt
, inv_y
, y
), special
);
129 return vbslq_f32 (pred_alt
, inv_y
, y
);
131 libmvec_hidden_def (V_NAME_F1 (tan
))
132 HALF_WIDTH_ALIAS_F1 (tan
)