1 /* 2 * Double-precision SVE sinpi(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 "mathlib.h" 9 #include "sv_math.h" 10 #include "pl_sig.h" 11 #include "pl_test.h" 12 #include "poly_sve_f64.h" 13 14 static const struct data 15 { 16 double poly[10]; 17 } data = { 18 /* Polynomial coefficients generated using Remez algorithm, 19 see sinpi.sollya for details. */ 20 .poly = { 0x1.921fb54442d184p1, -0x1.4abbce625be53p2, 0x1.466bc6775ab16p1, 21 -0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4, -0x1.e30750a28c88ep-8, 22 0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16, 23 0x1.af86ae521260bp-21, -0x1.012a9870eeb7dp-25 }, 24 }; 25 26 /* A fast SVE implementation of sinpi. 27 Maximum error 3.10 ULP: 28 _ZGVsMxv_sinpi(0x1.df1a14f1b235p-2) got 0x1.fd64f541606cp-1 29 want 0x1.fd64f541606c3p-1. */ 30 svfloat64_t SV_NAME_D1 (sinpi) (svfloat64_t x, const svbool_t pg) 31 { 32 const struct data *d = ptr_barrier (&data); 33 34 /* range reduction into -1/2 .. 1/2) 35 with n = rint(x) and r = r - n. */ 36 svfloat64_t n = svrinta_x (pg, x); 37 svfloat64_t r = svsub_x (pg, x, n); 38 39 /* Result should be negated based on if n is odd or not. */ 40 svuint64_t intn = svreinterpret_u64 (svcvt_s64_x (pg, n)); 41 svuint64_t sign = svlsl_z (pg, intn, 63); 42 43 /* y = sin(r). */ 44 svfloat64_t r2 = svmul_x (pg, r, r); 45 svfloat64_t r4 = svmul_x (pg, r2, r2); 46 svfloat64_t y = sv_pw_horner_9_f64_x (pg, r2, r4, d->poly); 47 y = svmul_x (pg, y, r); 48 49 return svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (y), sign)); 50 } 51 52 PL_SIG (SV, D, 1, sinpi, -0.9, 0.9) 53 PL_TEST_ULP (SV_NAME_D1 (sinpi), 2.61) 54 PL_TEST_SYM_INTERVAL (SV_NAME_D1 (sinpi), 0, 0x1p-63, 5000) 55 PL_TEST_SYM_INTERVAL (SV_NAME_D1 (sinpi), 0x1p-63, 0.5, 10000) 56 PL_TEST_SYM_INTERVAL (SV_NAME_D1 (sinpi), 0.5, 0x1p51, 10000) 57 PL_TEST_SYM_INTERVAL (SV_NAME_D1 (sinpi), 0x1p51, inf, 10000) 58