1/*
2 * Single-precision vector hypot(x) function.
3 *
4 * Copyright (c) 2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8#include "v_math.h"
9#include "pl_sig.h"
10#include "pl_test.h"
11
12#if WANT_SIMD_EXCEPT
13static 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
21static 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
31static float32x4_t VPCS_ATTR NOINLINE
32special_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
44float32x4_t VPCS_ATTR 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
71float32x4_t VPCS_ATTR 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 = vcge_u16 (
78      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
88PL_SIG (V, F, 2, hypot, -10.0, 10.0)
89PL_TEST_ULP (V_NAME_F2 (hypot), 1.21)
90PL_TEST_EXPECT_FENV (V_NAME_F2 (hypot), WANT_SIMD_EXCEPT)
91PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), 0, inf, 0, inf, 10000)
92PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), 0, inf, -0, -inf, 10000)
93PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), -0, -inf, 0, inf, 10000)
94PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), -0, -inf, -0, -inf, 10000)
95