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 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 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 44 float32x4_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 71 float32x4_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 88 PL_SIG (V, F, 2, hypot, -10.0, 10.0) 89 PL_TEST_ULP (V_NAME_F2 (hypot), 1.21) 90 PL_TEST_EXPECT_FENV (V_NAME_F2 (hypot), WANT_SIMD_EXCEPT) 91 PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), 0, inf, 0, inf, 10000) 92 PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), 0, inf, -0, -inf, 10000) 93 PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), -0, -inf, 0, inf, 10000) 94 PL_TEST_INTERVAL2 (V_NAME_F2 (hypot), -0, -inf, -0, -inf, 10000) 95