1 /*
2 * Single-precision vector hypot(x) function.
3 *
4 * Copyright (c) 2023-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "v_math.h"
9 #include "test_sig.h"
10 #include "test_defs.h"
11
12 #if WANT_SIMD_EXCEPT
13 static const struct data
14 {
15 uint32x4_t tiny_bound, thres;
16 } data = {
17 .tiny_bound = V4 (0x20000000), /* asuint (0x1p-63). */
18 .thres = V4 (0x3f000000), /* asuint (0x1p63) - tiny_bound. */
19 };
20 #else
21 static const struct data
22 {
23 uint32x4_t tiny_bound;
24 uint16x8_t thres;
25 } data = {
26 .tiny_bound = V4 (0x0C800000), /* asuint (0x1p-102). */
27 .thres = V8 (0x7300), /* asuint (inf) - tiny_bound. */
28 };
29 #endif
30
31 static float32x4_t VPCS_ATTR NOINLINE
special_case(float32x4_t x,float32x4_t y,float32x4_t sqsum,uint16x4_t special)32 special_case (float32x4_t x, float32x4_t y, float32x4_t sqsum,
33 uint16x4_t special)
34 {
35 return v_call2_f32 (hypotf, x, y, vsqrtq_f32 (sqsum), vmovl_u16 (special));
36 }
37
38 /* Vector implementation of single-precision hypot.
39 Maximum error observed is 1.21 ULP:
40 _ZGVnN4vv_hypotf (0x1.6a419cp-13, 0x1.82a852p-22) got 0x1.6a41d2p-13
41 want 0x1.6a41dp-13. */
42 #if WANT_SIMD_EXCEPT
43
V_NAME_F2(hypot)44 float32x4_t VPCS_ATTR NOINLINE V_NAME_F2 (hypot) (float32x4_t x, float32x4_t y)
45 {
46 const struct data *d = ptr_barrier (&data);
47
48 float32x4_t ax = vabsq_f32 (x);
49 float32x4_t ay = vabsq_f32 (y);
50
51 uint32x4_t ix = vreinterpretq_u32_f32 (ax);
52 uint32x4_t iy = vreinterpretq_u32_f32 (ay);
53
54 /* Extreme values, NaNs, and infinities should be handled by the scalar
55 fallback for correct flag handling. */
56 uint32x4_t specialx = vcgeq_u32 (vsubq_u32 (ix, d->tiny_bound), d->thres);
57 uint32x4_t specialy = vcgeq_u32 (vsubq_u32 (iy, d->tiny_bound), d->thres);
58 ax = v_zerofy_f32 (ax, specialx);
59 ay = v_zerofy_f32 (ay, specialy);
60 uint16x4_t special = vaddhn_u32 (specialx, specialy);
61
62 float32x4_t sqsum = vfmaq_f32 (vmulq_f32 (ax, ax), ay, ay);
63
64 if (unlikely (v_any_u16h (special)))
65 return special_case (x, y, sqsum, special);
66
67 return vsqrtq_f32 (sqsum);
68 }
69 #else
70
V_NAME_F2(hypot)71 float32x4_t VPCS_ATTR NOINLINE V_NAME_F2 (hypot) (float32x4_t x, float32x4_t y)
72 {
73 const struct data *d = ptr_barrier (&data);
74
75 float32x4_t sqsum = vfmaq_f32 (vmulq_f32 (x, x), y, y);
76
77 uint16x4_t special
78 = vcge_u16 (vsubhn_u32 (vreinterpretq_u32_f32 (sqsum), d->tiny_bound),
79 vget_low_u16 (d->thres));
80
81 if (unlikely (v_any_u16h (special)))
82 return special_case (x, y, sqsum, special);
83
84 return vsqrtq_f32 (sqsum);
85 }
86 #endif
87
88 HALF_WIDTH_ALIAS_F2 (hypot)
89
90 TEST_SIG (V, F, 2, hypot, -10.0, 10.0)
91 TEST_ULP (V_NAME_F2 (hypot), 1.21)
92 TEST_DISABLE_FENV_IF_NOT (V_NAME_F2 (hypot), WANT_SIMD_EXCEPT)
93 TEST_INTERVAL2 (V_NAME_F2 (hypot), 0, inf, 0, inf, 10000)
94 TEST_INTERVAL2 (V_NAME_F2 (hypot), 0, inf, -0, -inf, 10000)
95 TEST_INTERVAL2 (V_NAME_F2 (hypot), -0, -inf, 0, inf, 10000)
96 TEST_INTERVAL2 (V_NAME_F2 (hypot), -0, -inf, -0, -inf, 10000)
97