1 /* 2 * Double-precision SVE 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 "sv_math.h" 9 #include "pl_sig.h" 10 #include "pl_test.h" 11 12 static const struct data 13 { 14 uint64_t tiny_bound, thres; 15 } data = { 16 .tiny_bound = 0x0c80000000000000, /* asuint (0x1p-102). */ 17 .thres = 0x7300000000000000, /* asuint (inf) - tiny_bound. */ 18 }; 19 20 static svfloat64_t NOINLINE 21 special_case (svfloat64_t sqsum, svfloat64_t x, svfloat64_t y, svbool_t pg, 22 svbool_t special) 23 { 24 return sv_call2_f64 (hypot, x, y, svsqrt_x (pg, sqsum), special); 25 } 26 27 /* SVE implementation of double-precision hypot. 28 Maximum error observed is 1.21 ULP: 29 _ZGVsMxvv_hypot (-0x1.6a22d0412cdd3p+352, 0x1.d3d89bd66fb1ap+330) 30 got 0x1.6a22d0412cfp+352 31 want 0x1.6a22d0412cf01p+352. */ 32 svfloat64_t SV_NAME_D2 (hypot) (svfloat64_t x, svfloat64_t y, svbool_t pg) 33 { 34 const struct data *d = ptr_barrier (&data); 35 36 svfloat64_t sqsum = svmla_x (pg, svmul_x (pg, x, x), y, y); 37 38 svbool_t special = svcmpge ( 39 pg, svsub_x (pg, svreinterpret_u64 (sqsum), d->tiny_bound), d->thres); 40 41 if (unlikely (svptest_any (pg, special))) 42 return special_case (sqsum, x, y, pg, special); 43 return svsqrt_x (pg, sqsum); 44 } 45 46 PL_SIG (SV, D, 2, hypot, -10.0, 10.0) 47 PL_TEST_ULP (SV_NAME_D2 (hypot), 0.71) 48 PL_TEST_INTERVAL2 (SV_NAME_D2 (hypot), 0, inf, 0, inf, 10000) 49 PL_TEST_INTERVAL2 (SV_NAME_D2 (hypot), 0, inf, -0, -inf, 10000) 50 PL_TEST_INTERVAL2 (SV_NAME_D2 (hypot), -0, -inf, 0, inf, 10000) 51 PL_TEST_INTERVAL2 (SV_NAME_D2 (hypot), -0, -inf, -0, -inf, 10000) 52