1 /* Single-precision vector (AdvSIMD) exp2 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 uint16x8_t special_bound
;
27 uint32x4_t off
, mantissa_mask
;
30 /* Coefficients generated using Remez algorithm approximate
31 log2(1+r)/r for r in [ -1/3, 1/3 ].
32 rel error: 0x1.c4c4b0cp-26. */
33 .poly
= { V4 (0x1.715476p0f
), /* (float)(1 / ln(2)). */
34 V4 (-0x1.715458p
-1f
), V4 (0x1.ec701cp
-2f
), V4 (-0x1.7171a4p
-2f
),
35 V4 (0x1.27a0b8p
-2f
), V4 (-0x1.e5143ep
-3f
), V4 (0x1.9d8ecap
-3f
),
36 V4 (-0x1.c675bp
-3f
), V4 (0x1.9e495p
-3f
) },
37 .min_norm
= V4 (0x00800000),
38 .special_bound
= V8 (0x7f00), /* asuint32(inf) - min_norm. */
39 .off
= V4 (0x3f2aaaab), /* 0.666667. */
40 .mantissa_mask
= V4 (0x007fffff),
43 static float32x4_t VPCS_ATTR NOINLINE
44 special_case (float32x4_t x
, float32x4_t n
, float32x4_t p
, float32x4_t r
,
47 /* Fall back to scalar code. */
48 return v_call_f32 (log2f
, x
, vfmaq_f32 (n
, p
, r
), vmovl_u16 (cmp
));
51 /* Fast implementation for single precision AdvSIMD log2,
52 relies on same argument reduction as AdvSIMD logf.
53 Maximum error: 2.48 ULPs
54 _ZGVnN4v_log2f(0x1.558174p+0) got 0x1.a9be84p-2
56 float32x4_t VPCS_ATTR NOINLINE
V_NAME_F1 (log2
) (float32x4_t x
)
58 const struct data
*d
= ptr_barrier (&data
);
59 uint32x4_t u
= vreinterpretq_u32_f32 (x
);
60 uint16x4_t special
= vcge_u16 (vsubhn_u32 (u
, d
->min_norm
),
61 vget_low_u16 (d
->special_bound
));
63 /* x = 2^n * (1+r), where 2/3 < 1+r < 4/3. */
64 u
= vsubq_u32 (u
, d
->off
);
65 float32x4_t n
= vcvtq_f32_s32 (
66 vshrq_n_s32 (vreinterpretq_s32_u32 (u
), 23)); /* signextend. */
67 u
= vaddq_u32 (vandq_u32 (u
, d
->mantissa_mask
), d
->off
);
68 float32x4_t r
= vsubq_f32 (vreinterpretq_f32_u32 (u
), v_f32 (1.0f
));
70 /* y = log2(1+r) + n. */
71 float32x4_t r2
= vmulq_f32 (r
, r
);
72 float32x4_t p
= v_pw_horner_8_f32 (r
, r2
, d
->poly
);
74 if (__glibc_unlikely (v_any_u16h (special
)))
75 return special_case (x
, n
, p
, r
, special
);
76 return vfmaq_f32 (n
, p
, r
);
78 libmvec_hidden_def (V_NAME_F1 (log2
))
79 HALF_WIDTH_ALIAS_F1 (log2
)